text:
#------------------------------------
# HMC sampler
#--------------------------------------------
function hmc!(phik_curr, tildetx_hr_curr, parlr, parhr, scale_hmc = 1.0e-3)
    ePs = scale_hmc * rand()
    ulim =  30
    h0 = smooth_heavy(parlr.grd.r, 0.5, 1, 1500, 1/200) .* parlr.cPP ./ (parlr.grd.deltk^2)
    mk = 1.0e-2 ./ h0
    mk[parlr.pMaskBool] = 0.0

    phik_test = copy(phik_curr)
    rk   = white(parlr) .* sqrt(mk); # note that the variance of real(pk_init) and imag(pk_init) is mk/2
    grad, loglike   = ttk_grad_wlog(tildetx_hr_curr, phik_test, parlr, parhr)
    h_at_zero = 0.5 * sum( abs2(rk[!parlr.pMaskBool])./(2*mk[!parlr.pMaskBool]/2)) - loglike # the 0.5 is out front since only half the sum is unique

    for HMCCounter = 1:ulim
        loglike = lfrog!(phik_test, rk, tildetx_hr_curr, parlr, parhr, ePs, mk)
    end

    h_at_end = 0.5 * sum( abs2(rk[!parlr.pMaskBool])./(2*mk[!parlr.pMaskBool]/2) ) - loglike # the 0.5 is out front since only half the sum is unique
    prob_accept = minimum([1, exp(h_at_zero - h_at_end)])
    if rand() < prob_accept
        phik_curr[:] = phik_test
        println("Accept: prob_accept = $(round(prob_accept,4)), h_at_end = $(round(h_at_end)), h_at_zero = $(round(h_at_zero)), loglike = $(round(loglike))")
        return 1
    else
        println("Reject: prob_accept = $(round(prob_accept,4)), h_at_end = $(round(h_at_end)), h_at_zero = $(round(h_at_zero)), loglike = $(round(loglike))")
        return 0
    end
end

function lfrog!(phik_curr, rk, tildetx_hr_curr, parlr, parhr, ePs, mk)
    grad, loglike   = ttk_grad_wlog(tildetx_hr_curr, phik_curr, parlr, parhr)
    rk_halfstep =  rk +  ePs .* grad ./ 2.0
    inv_mk = 1./ (mk ./ 2.0)
    inv_mk[parlr.pMaskBool] = 0.0
    phik_curr[:] = phik_curr + ePs .* inv_mk .* rk_halfstep
    grad, loglike   = ttk_grad_wlog(tildetx_hr_curr, phik_curr, parlr, parhr)
    rk[:] = rk_halfstep + ePs .* grad ./ 2.0
    loglike
end

#-----------------------------
#  messenger algorithm
#------------------------------------
rwhite(ln, par::SpectrumGrids) = (par.grd.deltk/par.grd.deltx)*rft(randn(ln, ln), par)
function gibbspass_t!(sx, tx, phik_curr, ytx, maskvarx, parlr, parhr, coolingVec = [Inf for k=1:100])
    phidx1_hr_curr, phidx2_hr_curr, phidx1_lr_curr, phidx2_lr_curr = phi_easy_approx(phik_curr, parlr, parhr)
    dx, Nx = embedd(ytx, phidx1_lr_curr, phidx2_lr_curr, maskvarx, parlr, parhr)
    tk    = rft(tx, parhr)
    sk    = similar(tk)
    gpass!(sx, sk, tx, tk, dx, Nx, parhr, coolingVec)
    phidx1_hr_curr, phidx2_hr_curr
end
function gpass!(sx, sk, tx, tk, dx, Nx, par, cool)
  ln    = size(tx, 1)
  d2k   = par.grd.deltk * par.grd.deltk
  d2x   = par.grd.deltx * par.grd.deltx
  delt0 = 1 / d2k
  barNx = 0.99 * minimum(Nx)
  barNk = barNx * d2x * delt0 # var in fourier
  tldNx = Nx .- barNx
  Sk    = scale(par.cTT[1:size(sk, 1),1:size(sk, 2)], delt0)
  for uplim in cool
    λbarNk = (uplim > 8000.0) ? barNk : max(barNk, delt0 * par.CTell2d[int(uplim)])  # var in fourier
    λbarNx = λbarNk / (d2x * delt0)
    # --- update s
    tmpx   = 1 ./ (1 / λbarNk .+ 1 ./ Sk)
    sk[:]  = tmpx .* scale(tk, 1 / λbarNk)
    sk[:] += vec(rwhite(ln, par) .* √(tmpx))
    sx[:]  = irft(sk, par)
    # --- update t
    tmpx   = 1 ./ (1 ./ tldNx .+ 1 / λbarNx)
    tx[:]  = tmpx .* (dx ./ tldNx + sx ./ λbarNx)
    tx[:] += vec(randn(ln, ln) .* √(tmpx))
    tk[:]  = rft(tx, par)
  end
  nothing
end

##--------------------------------
# gradient of phi given  tildetx
#------------------------------------
function gradupdate!(phik_curr, tildetx_hr_curr, parlr, parhr, scale_grad=1.0e-3)
    for cntr = 1:30
        grad, loglike = ttk_grad_wlog(tildetx_hr_curr, phik_curr, parlr, parhr)
        println(loglike)
        phik_curr[:]  =  phik_curr + grad .* scale_grad .* smooth_heavy(parlr.grd.r, 1, 2, 1000, 1/200) .* parlr.cPP ./ (parlr.grd.deltk^2)
    end
end

function ttk_grad_wlog(tildetx_hr_sim, phik_curr, parlr, parhr)
    parlrgrddeltk2 = parlr.grd.deltk * parlr.grd.deltk
    phidx1_lr =  ifft2r(complex(0.0,1.0) .* parlr.grd.k1 .* phik_curr, parlr)
    phidx2_lr =  ifft2r(complex(0.0,1.0) .* parlr.grd.k2 .* phik_curr, parlr)
    #----- the following is slightly more efficient when working with real transforms
    tildetk_hr_sim   = rft(tildetx_hr_sim, parhr)
    rd1, rd2 = size(tildetk_hr_sim)
    tildetxd1   = irft(complex(0.0,1.0) .* parhr.grd.k1[1:rd1, 1:rd2] .* tildetk_hr_sim, parhr)
    tildetxd2   = irft(complex(0.0,1.0) .* parhr.grd.k2[1:rd1, 1:rd2] .* tildetk_hr_sim, parhr)

    tildetx_unlensed     = spline_interp2(parhr.grd.x, parhr.grd.y, tildetx_hr_sim   , parlr.grd.x - phidx1_lr, parlr.grd.y - phidx2_lr, 0.01)
    tildetxd1_unlensed   = spline_interp2(parhr.grd.x, parhr.grd.y, tildetxd1  , parlr.grd.x - phidx1_lr, parlr.grd.y - phidx2_lr, 0.01)
    tildetxd2_unlensed   = spline_interp2(parhr.grd.x, parhr.grd.y, tildetxd2  , parlr.grd.x - phidx1_lr, parlr.grd.y - phidx2_lr, 0.01)
    tildetk_unlensed     =  fft2(tildetx_unlensed  , parlr);  tildetk_unlensed[parlr.cMaskBool]   = 0.0
    tildetkd1_unlensed   =  fft2(tildetxd1_unlensed, parlr);  tildetkd1_unlensed[parlr.cMaskBool] = 0.0
    tildetkd2_unlensed   =  fft2(tildetxd2_unlensed, parlr);  tildetkd2_unlensed[parlr.cMaskBool] = 0.0
    #--------------- log likelihood
    loglike = 0.0
    maskforpsi = one(phik_curr)
    maskforpsi[parlr.pMaskBool] = 0.0
    for k=2:length(phik_curr)
        loglike += - 0.25 * abs2(maskforpsi[k] * phik_curr[k]) / (parlr.cPP[k] / (2.0 * parlrgrddeltk2))
        loglike += - 0.25 * abs2(tildetk_unlensed[k]) / (parlr.cTT[k] / (2.0 * parlrgrddeltk2))
    end
    #--------------------------------------
    Bk = tildetk_unlensed ./ (parlr.cTT)
    Bk[parlr.cMaskBool] = 0.0
    Bx = ifft2r(Bk, parlr)
    constantterm =  complex(0.0,-2.0) * parlrgrddeltk2 # the minus is here to convert psi gradient to phi gradient
    term  = constantterm * parlr.grd.k1 .* fft2(ifft2r(tildetkd1_unlensed, parlr) .* Bx, parlr)
    term += constantterm * parlr.grd.k2 .* fft2(ifft2r(tildetkd2_unlensed, parlr) .* Bx, parlr)
    term -= 2.0 * parlrgrddeltk2 *  phik_curr ./ parlr.cPP
    term[parlr.pMaskBool] = 0.0
    term, loglike
end

#----------------------------------------------
# quadratic estimate
#----------------------------------------------
function ttk_est(ytk, par)
    cTT = par.cTT
    cTTl_n_m = par.cTTLen + par.cNT .* par.cMaskInf
    cTTl_n = par.cTTLen + par.cNT
    tildetxd1_flt   = ifft2r(im * par.grd.k1 .* ytk .* cTT ./ cTTl_n_m, par)
    tildetxd2_flt   = ifft2r(im * par.grd.k2 .* ytk .* cTT ./ cTTl_n_m, par)
    Bx_flt     = ifft2r(ytk ./ cTTl_n_m, par)
    ttk  = complex(0.0,-2.0) * par.grd.k1 .* fft2(tildetxd1_flt .* Bx_flt, par)
    ttk += complex(0.0,-2.0) * par.grd.k2 .* fft2(tildetxd2_flt .* Bx_flt, par)
    ttk = ttk .* ttk_al(par)
    ttk[par.pMaskBool] = 0.0
    ttk
end
function ttk_al(parin)
    cTT = parin.cTT
    cTTl_n_m = parin.cTTLen + parin.cNT .* parin.cMaskInf
    cTTl_n = parin.cTTLen + parin.cNT

    gttsqk = zero(parin.grd.k1)

    tmp  = 2*ifft2r(im .* parin.grd.k1 .* cTTl_n .* cTT ./(cTTl_n_m.*cTTl_n_m),parin) .*  ifft2r(im.*parin.grd.k1.*cTTl_n.*(cTT)./(cTTl_n_m.*cTTl_n_m),parin)
    tmp += ifft2r(cTTl_n./(cTTl_n_m.*cTTl_n_m),parin) .*  ifft2r(im.*parin.grd.k1.*im.*parin.grd.k1.*cTTl_n.*(cTT.*cTT)./(cTTl_n_m.*cTTl_n_m),parin)
    tmp += ifft2r(cTTl_n./(cTTl_n_m.*cTTl_n_m),parin) .*  ifft2r(im.*parin.grd.k1.*im.*parin.grd.k1.*cTTl_n.*(cTT.*cTT)./(cTTl_n_m.*cTTl_n_m),parin)
    gttsqk += - 1.0 * parin.grd.k1 .* parin.grd.k1 .* (1./(2*pi)) .* fft2(tmp,parin)

    tmp  = 2*ifft2r(im .* parin.grd.k2 .* cTTl_n .* cTT ./(cTTl_n_m.*cTTl_n_m),parin) .*  ifft2r(im.*parin.grd.k1.*cTTl_n.*(cTT)./(cTTl_n_m.*cTTl_n_m),parin)
    tmp += ifft2r(cTTl_n./(cTTl_n_m.*cTTl_n_m),parin) .*  ifft2r(im.*parin.grd.k2.*im.*parin.grd.k1.*cTTl_n.*(cTT.*cTT)./(cTTl_n_m.*cTTl_n_m),parin)
    tmp += ifft2r(cTTl_n./(cTTl_n_m.*cTTl_n_m),parin) .*  ifft2r(im.*parin.grd.k2.*im.*parin.grd.k1.*cTTl_n.*(cTT.*cTT)./(cTTl_n_m.*cTTl_n_m),parin)
    gttsqk += - 2.0 * parin.grd.k2 .* parin.grd.k1 .* (1./(2*pi)) .* fft2(tmp,parin)

    tmp  = 2*ifft2r(im .* parin.grd.k2 .* cTTl_n .* cTT ./(cTTl_n_m.*cTTl_n_m),parin) .*  ifft2r(im.*parin.grd.k2.*cTTl_n.*(cTT)./(cTTl_n_m.*cTTl_n_m),parin)
    tmp += ifft2r(cTTl_n./(cTTl_n_m.*cTTl_n_m),parin) .*  ifft2r(im.*parin.grd.k2.*im.*parin.grd.k2.*cTTl_n.*(cTT.*cTT)./(cTTl_n_m.*cTTl_n_m),parin)
    tmp += ifft2r(cTTl_n./(cTTl_n_m.*cTTl_n_m),parin) .*  ifft2r(im.*parin.grd.k2.*im.*parin.grd.k2.*cTTl_n.*(cTT.*cTT)./(cTTl_n_m.*cTTl_n_m),parin)
    gttsqk += - 1.0 * parin.grd.k2 .* parin.grd.k2 .* (1./(2*pi)) .* fft2(tmp,parin)

    1.0 ./ abs(real(gttsqk))
end
# this printed out the code that follows
# for ka = ["parin.grd.k1", "parin.grd.k2"], kbin = ["parin.grd.k1", "parin.grd.k2"]
    # println("""
        # tmp  = 2*ifft2r(im .* $ka .* cTTl_n .* cTT ./(cTTl_n_m.*cTTl_n_m),par) .*  ifft2r(im.*$kb.*cTTl_n.*(cTT)./(cTTl_n_m.*cTTl_n_m),par)
        # tmp += ifft2r(cTTl_n./(cTTl_n_m.*cTTl_n_m),parin) .*  ifft2r(im.*$ka.*im.*$kb.*cTTl_n.*(cTT.*cTT)./(cTTl_n_m.*cTTl_n_m),parin)
        # tmp += ifft2r(cTTl_n./(cTTl_n_m.*cTTl_n_m),parin) .*  ifft2r(im.*$ka.*im.*$kb.*cTTl_n.*(cTT.*cTT)./(cTTl_n_m.*cTTl_n_m),parin)
        # gttsqk += - 1.0 * $ka .* $kb .* (1./(2*pi)) .* fft2(tmp,parin)
        # """)
# end
function ttk_n0(parin)
    tmp = 2.0 * ttk_al(parin)
    tmp[parin.cMaskBool] = Inf
    tmp
end

#----------------------------------
# simluation algorithm for generating synthetic data
#----------------------------------
function  simulate_start(par::SpectrumGrids; pad_proportion = 0.01)
    d = 2.0
    #    simulate noise for Q and U  on low res
    znt =  ((par.grd.deltk/par.grd.deltx)^(d/2.0))*fft2(randn(size(par.grd.x)),par.grd.deltx)
    ntk = znt.*sqrt(par.cNT)./(par.grd.deltk^(d/2.0))
    #   simulate stationary Q,U on high reslution
    hr_factor = 3
    parhr= setpar(par.grd.pixel_size_arcmin/hr_factor, par.grd.n*hr_factor, par.beamFWHM, par.nugget_at_each_pixel, 3000.0, 2.0);
    zt   = ((parhr.grd.deltk/parhr.grd.deltx)^(d/2.0)) * fft2(randn(size(parhr.grd.x)),parhr.grd.deltx)
    tk   = zt.*sqrt(parhr.cTT)./(parhr.grd.deltk^(d/2.0));
    tx   = ifft2r(tk, parhr.grd.deltk)
    zP   = ((parhr.grd.deltk/parhr.grd.deltx)^(d/2.0)) * fft2(randn(size(parhr.grd.x)),parhr.grd.deltx) #the real and imag parts should have var 1/2.
    phik = zP.*sqrt(parhr.cPP)./(parhr.grd.deltk^(d/2.0))
    phix = ifft2r(phik, parhr.grd.deltk)
    phidx1 = ifft2r(complex(0.0,1.0) .* parhr.grd.k1 .* phik, parhr)
    phidx2 = ifft2r(complex(0.0,1.0) .* parhr.grd.k2 .* phik, parhr)
    #   Lense Q,U
    tildetx = spline_interp2(parhr.grd.x, parhr.grd.y, tx, parhr.grd.x + phidx1, parhr.grd.y + phidx2, pad_proportion)
    #tildetx = ifft2r(parhr.beam.*fft2(tildetx,parhr.grd.deltx), parhr.grd.deltk) #h the beam is already a part of the noise...
    tildetk = fft2(tildetx[1:hr_factor:end, 1:hr_factor:end], par)
    ytk = ntk + tildetk
    ytk, tildetk, phix[1:hr_factor:end, 1:hr_factor:end], tx
end

#----------------------------------
# miscilanious function
#----------------------------------
white(par::SpectrumGrids) = (par.grd.deltk/par.grd.deltx)*fft2(randn(size(par.grd.x)),par.grd.deltx)
function upsample(psik_lr::Array{Complex{Float64},2} ,parlr, parhr)
    index =  (-parlr.grd.nyq-parlr.grd.deltk./10) .<= parhr.grd.kside .< (parlr.grd.nyq-parlr.grd.deltk./10)
    psik_hr = zeros(Complex{Float64},size(parhr.grd.k1))
    psik_hr[index,index]=psik_lr
    return psik_hr
end
function upsample(txlr::Array{Float64,2} ,parlr, parhr)
    tklr = fft2(txlr, parlr)
    index =  (-parlr.grd.nyq-parlr.grd.deltk./10) .<= parhr.grd.kside .< (parlr.grd.nyq-parlr.grd.deltk./10)
    tkhr = zeros(Complex{Float64},size(parhr.grd.k1))
    tkhr[index,index]=tklr
    ifft2r(tkhr, parhr)
end
function downsample(txhr::Array{Float64,2} ,parlr, parhr)
    freq = round(parlr.grd.deltx / parhr.grd.deltx)
    txhr[1:freq:end,1:freq:end]
end
function  phi_easy_approx(phik, parlr, parhr)
    phik_hr = upsample(phik, parlr, parhr)
    phidx1_hr = ifft2r(complex(0.0,1.0) .* parhr.grd.k1 .* phik_hr, parhr)
    phidx2_hr = ifft2r(complex(0.0,1.0) .* parhr.grd.k2 .* phik_hr, parhr)
    phidx1_lr = ifft2r(complex(0.0,1.0) .* parlr.grd.k1 .* phik, parlr)
    phidx2_lr = ifft2r(complex(0.0,1.0) .* parlr.grd.k2 .* phik, parlr)
    phidx1_hr, phidx2_hr, phidx1_lr, phidx2_lr
end
function embedd(ytx, phidx1_lr, phidx2_lr, Nx_lowres, parlr, parhr; pad_proportion = 0.1)
    index_matrix1 = reshape(1:length(parhr.grd.x),size(parhr.grd.x)); #this index matrix only needs to be made once
    I_obs_into_xyfull = int(nearest_interp2(parhr.grd.x ,parhr.grd.y, index_matrix1, parlr.grd.x + phidx1_lr, parlr.grd.y + phidx2_lr; pad_proportion = pad_proportion))
    dx = zeros(Float64,size(parhr.grd.x))
    dx[I_obs_into_xyfull]  = vec(ytx)
    NN = fill(Inf, size(dx))
    NN[I_obs_into_xyfull]= vec(Nx_lowres)
    dx, NN
end
function smooth_heavy(x,height1,height2,location, transition_sharpness)
    zeroto1 = 0.5 .+ 0.5 .* tanh(transition_sharpness .* (x.-location))
    (height2 .- height1) .* zeroto1 .+ height1
end
function radialAve(Xk,par)
    block_counter = 0.0
    block_delta = par.grd.deltk
    radXk = zero(Xk)
    while true
      k_index   = ((block_counter)*block_delta .< par.grd.r .<= (block_counter + 3)*block_delta) & (1 .< par.grd.r)
      all(~k_index) && break
      radXk[k_index] =  mean( Xk[k_index] )
      block_counter += 1.0
    end
    radXk
end
#  plt.semilogy(parglob.grd.r[:,1], radialAve(abs2(ttk1),parglob)[:,1]); plt.show()
#  plt.semilogy(parglob.grd.r[:,1], radialAve(ttk_al(parglob),parglob)[:,1]); plt.show()

function repstretch(start, stop, rep::Int, total::Int)
    tdr = int(total / rep)
    return vec(linspace(start, stop, tdr)' .* ones(rep, tdr))
end

function pos(x)
        y = real(x)
        y[y.<0.0] = 0.0
        y
end

max_stars_repo_name: EthanAnderes/BayesianCmbLensing

max_issues_repo_name: EthanAnderes/BayesianCmbLensing

max_forks_repo_name: EthanAnderes/BayesianCmbLensing


text:
using FlameGraphs, AbstractTrees, Colors
using Base.StackTraces: StackFrame
using Test, Profile

# useful for testing
stackframe(func, file, line; C=false) = StackFrame(Symbol(func), Symbol(file), line, nothing, C, false, 0)

@testset "flamegraph" begin
    backtraces = UInt64[0, 4, 3, 2, 1,   # order: calles then caller
                        0, 6, 5, 1,
                        0, 8, 7,
                        0, 4, 3, 2, 1,
                        0]
    lidict = Dict{UInt64,StackFrame}(1=>stackframe(:f1, :file1, 1),
                                     2=>stackframe(:f2, :file1, 5),
                                     3=>stackframe(:f3, :file2, 1),
                                     4=>stackframe(:f2, :file1, 15),
                                     5=>stackframe(:f4, :file1, 20),
                                     6=>stackframe(:f5, :file3, 1),
                                     7=>stackframe(:f1, :file1, 2),
                                     8=>stackframe(:f6, :file3, 10))
    g = flamegraph(backtraces; lidict=lidict)
    @test all(node->node.data.status == 0, PostOrderDFS(g))
    level1 = collect(g)
    @test length(level1) == 2
    n1, n2 = level1
    @test n1.data.sf.func === :f1
    @test n1.data.sf.line == 1
    @test n1.data.span == 1:3
    @test n2.data.sf.func === :f1
    @test n2.data.sf.line == 2
    @test n2.data.span == 4:4
    level2a = collect(n1)
    @test length(level2a) == 2
    n3, n4 = level2a
    @test n3.data.sf.func === :f2
    @test n3.data.sf.line == 5
    @test n3.data.span == 1:2
    @test n4.data.sf.func === :f4
    @test n4.data.sf.line == 20
    @test n4.data.span == 3:3
    level2b = collect(n2)
    @test length(level2b) == 1
    n5 = level2b[1]
    @test n5.data.sf.func === :f6
    @test n5.data.sf.line == 10
    @test n5.data.span == 4:4
    level3a = collect(n3)
    @test length(level3a) == 1
    n6 = level3a[1]
    @test n6.data.sf.func === :f3
    @test n6.data.sf.line == 1
    @test n6.data.span == 1:2
    level3b = collect(n4)
    @test length(level3b) == 1
    n7 = level3b[1]
    @test n7.data.sf.func === :f5
    @test n7.data.sf.line == 1
    @test n7.data.span == 3:3
    @test isempty(n5)
    level4a = collect(n6)
    @test length(level4a) == 1
    n8 = level4a[1]
    @test n8.data.sf.func === :f2
    @test n8.data.sf.line == 15
    @test n8.data.span == 1:2
    @test isempty(n7)
    @test isempty(n8)

    io = IOBuffer()
    print_tree(io, g)
    str = String(take!(io))
    @test str == """
FlameGraphs.NodeData(ip:0x0, 0x00, 1:4)
├─ FlameGraphs.NodeData(f1 at file1:1, 0x00, 1:3)
│  ├─ FlameGraphs.NodeData(f2 at file1:5, 0x00, 1:2)
│  │  └─ FlameGraphs.NodeData(f3 at file2:1, 0x00, 1:2)
│  │     └─ FlameGraphs.NodeData(f2 at file1:15, 0x00, 1:2)
│  └─ FlameGraphs.NodeData(f4 at file1:20, 0x00, 3:3)
│     └─ FlameGraphs.NodeData(f5 at file3:1, 0x00, 3:3)
└─ FlameGraphs.NodeData(f1 at file1:2, 0x00, 4:4)
   └─ FlameGraphs.NodeData(f6 at file3:10, 0x00, 4:4)
"""

    # pruning
    c = g.child.child
    @test c.child != c
    g = flamegraph(backtraces; lidict=lidict, pruned=[(:f3, "file2")])
    c = g.child.child
    @test c.child == c

    # Now make some of them C calls
    lidict = Dict{UInt64,StackFrame}(1=>stackframe(:f1, :file1, 1),
                                     2=>stackframe(:jl_f, :filec, 55; C=true),
                                     3=>stackframe(:jl_invoke, :file2, 1; C=true),
                                     4=>stackframe(:f2, :file1, 15),
                                     5=>stackframe(:f4, :file1, 20),
                                     6=>stackframe(:f5, :file3, 1),
                                     7=>stackframe(:f1, :file1, 2),
                                     8=>stackframe(:f6, :file3, 10))
    g = flamegraph(backtraces; lidict=lidict)
    level1 = collect(g)
    @test length(level1) == 2
    n1, n2 = level1
    @test n1.data.sf.func === :f1
    @test n1.data.sf.line == 1
    @test n1.data.span == 1:3
    # Note: only span 2:3 contributes to the rt dispatch, but since we mark the caller
    # rather than the callee, and 1:3 all are the same caller, the whole thing gets marked.
    @test n1.data.status == FlameGraphs.runtime_dispatch
    @test n2.data.sf.func === :f1
    @test n2.data.sf.line == 2
    @test n2.data.span == 4:4
    @test n2.data.status == 0
    level2a = collect(n1)
    @test length(level2a) == 2
    n3, n4 = level2a
    @test n3.data.sf.func === :f4
    @test n3.data.sf.line == 20
    @test n3.data.span == 1:1
    @test n4.data.sf.func === :f2
    @test n4.data.sf.line == 15
    @test n4.data.span == 2:3
    level2b = collect(n2)
    @test length(level2b) == 1
    n5 = level2b[1]
    @test n5.data.sf.func === :f6
    @test n5.data.sf.line == 10
    @test n5.data.span == 4:4
    level3a = collect(n3)
    @test length(level3a) == 1
    n6 = level3a[1]
    @test n6.data.sf.func === :f5
    @test n6.data.sf.line == 1
    @test n6.data.span == 1:1
    @test isempty(n4)
    @test isempty(n5)
    @test isempty(n6)

    # Pruning REPL code
    lidict = Dict{UInt64,StackFrame}(1=>stackframe(:f1, :file1, 1),
                                     2=>stackframe(:eval_user_input, "REPL.jl", 5),
                                     3=>stackframe(:f3, :file2, 1),
                                     4=>stackframe(:f2, :file1, 15),
                                     5=>stackframe(:f4, :file1, 20),
                                     6=>stackframe(:f5, :file3, 1),
                                     7=>stackframe(:f1, :file1, 2),
                                     8=>stackframe(:f6, :file3, 10))
    g = flamegraph(backtraces; lidict=lidict)
    sfc = [c.data.sf for c in g]
    @test lidict[3] ∈ sfc
    @test lidict[1] ∉ sfc
end

@testset "flamepixels" begin
    backtraces = UInt64[0, 4, 3, 2, 1,   # order: calles then caller
                        0, 6, 5, 1,
                        0, 8, 7,
                        0, 4, 3, 2, 1,
                        0]
    lidict = Dict{UInt64,StackFrame}(1=>stackframe(:f1, :file1, 1),
                                     2=>stackframe(:f2, :file1, 5),
                                     3=>stackframe(:f3, :file2, 1),
                                     4=>stackframe(:f2, :file1, 15),
                                     5=>stackframe(:f4, :file1, 20),
                                     6=>stackframe(:f5, :file3, 1),
                                     7=>stackframe(:f1, :file1, 2),
                                     8=>stackframe(:f6, :file3, 10))
    g = flamegraph(backtraces; lidict=lidict)
    img = flamepixels(g)
    fc = FlameGraphs.FlameColors()
    @test all(img[:,1] .== fc.colors[1])
    @test all(img[1:3,2] .== fc.colors[3])
    @test img[4,2] == fc.colors[4]
    @test all(img[1:2,3] .== fc.colors[1])
    @test img[3,3] == fc.colors[2]
    @test img[4,3] == fc.colors[1]
    @test all(img[1:2,4] .== fc.colors[3])
    @test img[3,4] == fc.colors[4]
    @test img[4,4] == fc.colorbg
    @test all(img[1:2,5] .== fc.colors[1])
    @test all(img[3:4,5] .== fc.colorbg)
    imgtags = flametags(g, img)
    @test axes(imgtags) == axes(img)
    @test imgtags[1,2] == lidict[1]
    @test imgtags[4,2] == lidict[7]
    @test imgtags[4,end] == StackTraces.UNKNOWN

    lidict = Dict{UInt64,StackFrame}(1=>stackframe(:f1, :file1, 1),
                                     2=>stackframe(:jl_f, :filec, 55; C=true),
                                     3=>stackframe(:jl_invoke, :file2, 1; C=true),
                                     4=>stackframe(:f2, :file1, 15),
                                     5=>stackframe(:f4, :file1, 20),
                                     6=>stackframe(:f5, :file3, 1),
                                     7=>stackframe(:f1, :file1, 2),
                                     8=>stackframe(:f6, :file3, 10))
    g = flamegraph(backtraces; lidict=lidict)
    img = flamepixels(g)
    @test all(img[:,1] .== fc.colors[1])
    # Note: only span 2:3 contributes to the rt dispatch, but since we mark the caller
    # rather than the callee, and 1:3 all are the same caller, the whole thing gets marked.
    @test all(img[1:3,2] .== fc.colorsrt[3])
    @test img[4,2] == fc.colors[4]
    @test img[1,3] == fc.colors[1]
    @test all(img[2:3,3] .== fc.colors[2])
    @test img[4,3] == fc.colors[1]
    @test img[1,4] == fc.colors[3]
    @test all(img[2:4,4] .== fc.colorbg)

    # Customizing FlameColors
    # the classic colors which were used in FlameGraphs v0.1 or ProfileView v0.5
    fc2 = FlameColors(
        parse.(RGB, ["#E870DD", "#32B44E", "#1AA2FF", "#00DEE6", "#FFA49C",
                     "#9E9E9E", "#A8A200", "#CDB9FF", "#00E5B2", "#FF5F82"]),
        colorant"white", colorant"black", [colorant"red"], [colorant"orange"])
    img = flamepixels(fc2, g)
    @test all(img[:,1] .== fc2.colors[1])
    @test all(img[1:3,2] .== fc2.colorsrt[1])
    @test img[4,2] == fc2.colors[7]
    @test img[1,3] == fc2.colors[1]
    @test all(img[2:3,3] .== fc2.colors[2])
    @test img[4,3] == fc2.colors[3]
    @test img[1,4] == fc2.colors[6]
    @test all(img[2:4,4] .== fc2.colorbg)
end

@testset "Profiling" begin
     A = randn(100, 100, 200)
     Profile.clear()
     @profile mapslices(sum, A; dims=2)
     g = flamegraph()
     @test FlameGraphs.depth(g) > 10
     img = flamepixels(StackFrameCategory(), flamegraph(C=true))
     @test any(img .== colorant"orange")
     A = [1,2,3]
     sum(A)
     Profile.clear()
     @profile sum(A)
     Sys.islinux() && @test_logs (:warn, r"There were no samples collected.") flamegraph() === nothing
end

ext: jl


text:
using Pkg;
Pkg.activate("././.")
using Revise
using ITensors

#
# DMRG calculation of the extended Hubbard model
# ground state wavefunction, and spin densities
#

let
  ## PARAMS ##

  # Model
  N = 8
  t = 1.0
  U = 8

  # DMRG
  sweeps = Sweeps(20)
  setmaxdim!(sweeps, 50, 100, 200, 400, 800, 800)
  setcutoff!(sweeps, 1E-12)
  @show sweeps

  # TEBD
  T = 10
  τ = 0.01
  TEBD_maxdim = 800
  TEBD_cutoff = 1E-14

  ## SET UP THE MODEL ##

  # Make N electron sites with conserved number and spin
  sites = siteinds("Electron", N; conserve_qns=true)

  # Make the Hamiltonian
  ampo = OpSum()
  for i in 1:(N - 1)
    ampo += -t, "Cdagup", i, "Cup", i + 1
    ampo += -t, "Cdagup", i + 1, "Cup", i
    ampo += -t, "Cdagdn", i, "Cdn", i + 1
    ampo += -t, "Cdagdn", i + 1, "Cdn", i

    ampo += U, "Nupdn", i
  end
  # Edge site
  ampo += U, "Nupdn", N
  H = MPO(ampo, sites)

  ## FIND GROUND STATE ##

  # Initialize wavefunction
  psi0 = productMPS(sites, n -> isodd(n) ? "Up" : "Dn")
  @show flux(psi0)

  # Start DMRG calculation
  energy, psi = dmrg(H, psi0, sweeps)

  ## TIME-EVOLVE ##

  # Make the trotter gates
  gates = ITensor[]
  for j=1:N-1
    s1 = sites[j] # site j
    s2 = sites[j+1] # site j+1

    hj = -t * (op("Cdagup*F",s1) * op("Cup",s2)  # t * (c^†_jσ c_{j+1}σ + h.c.)
                 -op("Cup*F",s1) * op("Cdagup",s2)
                 +op("Cdagdn*F",s1) * op("Cdn",s2)
                 -op("Cdn*F",s1) * op("Cdagdn",s2))

    Gj = exp(-1.0im * τ/2 * hj)
    push!(gates,Gj)
    hj2 = U * op("Nupdn",s1)
    Gj2 = exp(-1.0im * τ/2 * hj2)
    push!(gates,Gj2)
  end
  # End site
  hn = U*op("Nupdn",sites[N])
  Gn = exp(-1.0im * τ/2 * hn)
  push!(gates,Gn)
  # Append gates in reverse to complete Trotter formula
  append!(gates,reverse(gates))

  function compute_overlap(ψ1::MPS, ψ2::MPS)
      LinearAlgebra.norm(inner(ψ1, ψ2))
  end

  # Evolve
  nsteps = floor(T/τ) # Number of time steps for time evolution
  t = 0.0
  ϕ = copy(psi)
  for step in 1:nsteps
    print(floor(Int,step),"-")
    ϕ = apply(gates, ϕ; maxdim=TEBD_maxdim, cutoff=TEBD_cutoff)
    t += τ

    @show compute_overlap(ϕ, psi)

  end

end


text:
# using FNVHash

include("..\\generalized-parser\\source\\LogProbs.jl")
include("..\\generalized-parser\\source\\ProductSemiRings.jl")
include("..\\chinese-restaurants\\ChineseRestaurants.jl")
# include("JazzGrammars3.jl")
include("JazzGrammars5.jl")

# include("..\\generalized-parser\\source\\ChartParser2.jl")
# include("..\\generalized-parser\\source\\ChartParser5.jl")
include("..\\generalized-parser\\source\\ChartParser6.jl")

include("read_goldstandard.jl")
include("..\\tree-annotation-tool\\treebank_functions.jl")

###################
### Simple Test ###
###################

grammar = make_jazz_grammar()
# grammar = make_plain_PCFG_jazz_grammar()
@time forest = run_chartparser(map(Terminal, ["0maj","1maj","7maj","0maj"]), grammar)
# @time forest = run_chartparser(map(Terminal, ["0maj","0maj","0maj"]), grammar)
println(score(forest))
println(best_tree(forest))

####################
### Catalan Test ###
####################

# chords(n::Int) = [Terminal("0maj") for i in 1:n]
# number_parses(n::Int) = score(run_chartparser(chords(n), grammar))
# catalan_number(n::Int) = div(factorial(BigInt(2*n)), factorial(BigInt(n+1)) * factorial(BigInt(n)))
#
# for n in 1:20
#     println(n, "\t", number_parses(n), "\t", catalan_number(n-1))
# end

######################
### Profiling Test ###
######################

# corpus_dir = "C:\\Users\\Daniel\\Documents\\Bitbucket\\musical-grammars\\data\\Jazz-selection\\"
# testfile = "'Round Midnight_id_05848_midicons.txt"
# testfile2 = "After You've Gone_id_06924_allanah.txt"
#
# function translate_corpus_notation(chord::AbstractString)
#     m = match(r"([A-G]b*#*)(.*)(_)", chord)
#     string(pitch_class(m[1]), m[2])
# end
#
# function read_chords(file::AbstractString)
#     chords = map(x->x[2], map(split, map(strip, readlines(file)[2:end-1])))
#     map(translate_corpus_notation, chords)
# end

# chords = read_chords(corpus_dir*testfile2)
# @time forest = run_chartparser(chords, grammar)
# println(score(forest))
#
# using ProfileView
# Profile.clear()
# @profile run_chartparser(chords, grammar)
# ProfileView.view()

# workspace()


text:
using Redis
import DataStructures: OrderedSet
using Random
using Dates
using Test
using Base

include(joinpath(dirname(@__FILE__),"client_tests.jl"))
include(joinpath(dirname(@__FILE__),"redis_tests.jl"))


text:
module C
const NAMES = [
    "kCp0"
    "kCpB0"
    "kCpM0"
    "kCpBM0"
    "kCpS"
    "kCpBS"
    "kCpMS"
    "kCpBMS"
    "kmCpS"
    "kmCpBS"
    "kmCpMS"
    "kmCpBMS"
    "kCpu"
    "kCpBu"
    "kCpMu"
    "kCpBMu"
    "kmCpu"
    "kmCpBu"
    "kmCpMu"
    "kmCpBMu"

    "kBaM"
    "kCaB"
    "kCaM"
    "kCBaM"
    "kCMaB"
    "kCaBM"
    "kCpaBM"
    "kCpaB"
    "kCpaM"
    "kCpBaM"
    "kCpMaB"
    "kBdM"
    "kCdB"
    "kCdM"
    "kCBdM"
    "kCMdB"
    "kCpdB"
    "kCpdM"
    "kCpBdM"
    "kCpMdB"
    "kCdBM"
    "kCpdBM"

    "kTp0"
    "kCTpS"
    "kCBTpS"
    "kCMTpS"
    "kCBMTpS"
    "kmCTpS"
    "kmCBTpS"
    "kmCMTpS"
    "kmCBMTpS"
    "kTpIKK1"
    "kmTpIKK1"
    "kTpu"
    "kmTpu"
    "kTpC0"
    "kCTpCS"
    "kCBTpCS"
    "kCMTpCS"
    "kCBMTpCS"
    "kmCTpCS"
    "kmCBTpCS"
    "kmCMTpCS"
    "kmCBMTpCS"
    "kTpCIKK1"
    "kmTpCIKK1"
    "kTpCIKK2"
    "kmTpCIKK2"
    "kTpCIKK3"
    "kmTpCIKK3"
    "kTpCu"
    "kmTpCu"
    "kTaB"
    "kTaM"
    "kTaBM"
    "kTpaB"
    "kTpaM"
    "kTpaBM"
    "kTCd"
    "kTpCd"

    "kIp0"
    "kIpTAKp"
    "kmIpTAKp"
    "kIpTAKpC"
    "kmIpTAKpC"
    "kIpu"
    "kmIpu"
    "kICp0"
    "kICpTAKp"
    "kmICpTAKp"
    "kICpTAKpC"
    "kmICpTAKpC"
    "kIpCu"
    "kmIpCu"
    "kIaB"
    "kIaM"
    "kIaBM"
    "kIpaB"
    "kIpaM"
    "kIpaBM"
    "kICd"
    "kIpCd"
    "kIpCfaIKKpC"
    "kmIpCfaIKKpC"
    "kIpCfaIKKppC"
    "kmIpCfaIKKppC"
    "kIppCu"
    "kmIppCu"
    "kIppCd"
    "kIppChf"
    "kmIppChf"
    "kIpphf"
    "kmIpphf"
    "kIir"
    "kmIir"
    "kIpA20"
    "kmIpA20"
    "kIpCA20"
    "kmIpCA20"

    "kassaikk"
    "kassaikknfkb"
    "kassanfkbikk"
    "kdisaikk"
    "kdisaikknfkb"
    "kdisanfkbikk"
    "kassa"
    "kdisa"
    "kdegbounda"
    "kdegboundaIKK"
    "kdegfreea"
    "kdegfreeaIKK"
    "kshutboundikbain"
    "kshutboundikbaout"
    "kshutfreeikbain"
    "kshutfreeikbaout"
    "k0mrnaikba"
    "kdegmrnaikba"
    "khillprodmrnaikba"
    "kpikba"
    "kprodmrnaikba"

    "kassbikk"
    "kassbikknfkb"
    "kassbnfkbikk"
    "kdisbikk"
    "kdisbikknfkb"
    "kdisbnfkbikk"
    "kassb"
    "kdegboundb"
    "kdegboundbIKK"
    "kdegfreeb"
    "kdegfreebIKK"
    "kdisb"
    "kshutboundikbbin"
    "kshutboundikbbout"
    "kshutfreeikbbin"
    "kshutfreeikbbout"
    "k0mrnaikbb"
    "kdegmrnaikbb"
    "khillprodmrnaikbb"
    "kpikbb"
    "kprodmrnaikbb"
    "kasseikk"
    "kasseikknfkb"
    "kassenfkbikk"
    "kdiseikk"
    "kdiseikknfkb"
    "kdisenfkbikk"

    "kasse"
    "kdegbounde"
    "kdegboundeIKK"
    "kdegfreee"
    "kdegfreeeIKK"
    "kdise"
    "kshutboundikbein"
    "kshutboundikbeout"
    "kshutfreeikbein"
    "kshutfreeikbeout"
    "k0mrnaikbe"
    "kdegmrnaikbe"
    "khillprodmrnaikbe"
    "kpikbe"
    "kprodmrnaikbe"

    "kshutnfkbin"
    "kshutnfkbout"

    "k0mrnaa20"
    "kdegmrnaa20"
    "kprodmrnaa20"
    "kpa20"
    "kdega20"
    "khillprodmrnaa20"
    "delayrnae"
    "X"
    "term"
]

for (idx,name) in enumerate(NAMES)
    eval(Meta.parse("const $name = $idx"))
end

const NUM = length(NAMES)

end  # module

max_stars_repo_name: himoto/paramEstim

max_issues_repo_name: himoto/paramEstim


text:
macro except_str(expr, err_type)
    return quote
        let err = nothing
            try
                $(esc(expr))
            catch err
            end
            err === nothing && error("expected failure, but no exception thrown")
            @test typeof(err) === $(esc(err_type))
            buf = IOBuffer()
            showerror(buf, err)
            String(take!(buf))
        end
    end
end

@testset "error hints" begin
    @testset "ones" begin
        for T in (RGB, RGB{N0f8})
            err_str = @except_str one(T) MethodError
            @test occursin(r"MethodError: no method matching one\(::Type\{RGB.*\}", err_str)
            @test occursin("You may need to `using ColorVectorSpace`.", err_str)

            err_str = @except_str one(T(1, 1, 1)) MethodError
            @test occursin(r"MethodError: no method matching one\(::RGB\{.*\}", err_str)
            @test occursin("You may need to `using ColorVectorSpace`.", err_str)

            err_str = @except_str ones(T) MethodError
            @test occursin(r"MethodError: no method matching one\(::Type\{RGB.*\}", err_str)
            @test occursin("You may need to `using ColorVectorSpace`.", err_str)
        end
    end
    @testset "Math" begin
        gray = Gray(0.8)
        rgb = RGB{Float32}(1, 0, 0)
        err_str = @except_str gray + rgb MethodError
        @test occursin("no method matching +(::Gray{Float64}, ::RGB{Float32})", err_str)
        @test occursin("Math on colors is deliberately undefined in ColorTypes, but see the ColorVectorSpace package", err_str)

        err_str = @except_str gray * rgb MethodError
        @test occursin("no method matching *(::Gray{Float64}, ::RGB{Float32})", err_str)
        @test occursin("Math on colors is deliberately undefined in ColorTypes, but see the ColorVectorSpace package", err_str)
        @test occursin("You may also need `⋅`, `⊙`, or `⊗`.", err_str)

        err_str = @except_str rgb * 0.5 MethodError
        @test !occursin("You may also need `⋅`, `⊙`, or `⊗`.", err_str)
    end
end


text:
module BenchFillArrays

using SIMDIndexing

include("../examples/vdot.jl")
# using LinearAlgebra
using FillArrays
using BenchmarkTools

suite = BenchmarkGroup()

arr1 = randn(2^20)
arr2 = Ones{Float64}(length(arr1))

# Bump SIMD width to 16.  It looks like `@simd` macro (or LLVM?) likes
# to unroll loops a lot.  So I have to increase the vector size to
# make the performance comparable.
N = 16

# There is no way `vdot` is faster; but it must not be _slower_.
suite["vdot"] = @benchmarkable vdot($arr1, $arr2, $(Vec{N, Float64}))
suite["sum"] = @benchmarkable sum($arr1)
# suite["dot"] = @benchmarkable dot($arr1, $arr2)

function main()
    global results = run(suite; verbose = true)
    show(IOContext(stdout, :compact => false), results)

    println()
    println()
    println("Benchmark result is stored in $(@__MODULE__).results")

    println("vdot vs sum")
    display(judge(minimum.((results["vdot"], results["sum"]))...))
end

end  # module
BenchFillArrays.suite


text:
using Random
using Statistics
using StatsBase
using Plots
import Dates
import Base.+
import Base.-

include("Point_Struct.jl")
include("fun.jl")
include("constants.jl")
include("integrate.jl")

#################### Integration #####################

calls = Int64(1E3)
it    = 5
dimensions = 6.

Δₐ = 0.4
ss = 0.04 #0.02 rad  / sec

fₛ = ss / Δₐ

Φₛ = π/2         # start point
vₛ = 2*π / 60.   # speed rad/sec
dt = 0.3         # Time interval

Start_Time = 0   #parse(Float64, ARGS[1])
Stop_Time  = 35 #parse(Float64, ARGS[2])
file_name  = "ciao.txt" #ARGS[3]

Step_start = Start_Time / dt
Step_stop  = Stop_Time / dt

time_s = Array{Float64, 1}(undef, 0)
Coo  = Array{Float64, 1}(undef, 0)
std_err = Array{Float64,1}(undef, 0)
el_arr = Array{Float64,1}(undef, 0)
az_arr = Array{Float64,1}(undef, 0)

el = π/(4.5)

# Raster scan
# xu = Point{Float64}(4000.0,  π/(4.5) + θᵦ, Φₛ + θᵦ, undef, undef, undef)
# xl = Point{Float64}(0.,      π/(4.5) - θᵦ,  Φₛ - θᵦ, undef, undef, undef )

xu = Point{Float64}(2000.0,  el + θᵦ, Φₛ + θᵦ, undef, undef, undef)
xl = Point{Float64}(00.,     el - θᵦ, Φₛ - θᵦ, undef, undef, undef)

start = Dates.value(Dates.now())

#hit=Array{Float64}(undef,0)

for t = Step_start:Step_stop
    tempo = Float64(t*dt)
    # xu = Point{Float64}(2000.0,  el + θᵦ, Φₛ + θᵦ, undef, undef, undef)
    # xl = Point{Float64}(00.,     el - θᵦ, Φₛ - θᵦ, undef, undef, undef)
    # Rivedere gli estremi di integrazione.
    # xu_p = Point{Float64}(5035.,   2.8902652413026098,  (Φₛ + vₛ*t*dt) + θᵦ, undef, undef, undef)
    # xl_p = Point{Float64}(0.,     0.12566370614359174,  (Φₛ + vₛ*t*dt) - θᵦ , undef, undef ,undef)

    # Raster scan
    # xu_p = Point{Float64}(4000.,   π/(4.5) + θᵦ,  (Φₛ + (Δₐ/2)*sin(2*π*fₛ*t*dt)) + θᵦ, undef, undef, undef)
    # xl_p = Point{Float64}(0.,      π/(4.5) - θᵦ,  (Φₛ + (Δₐ/2)*sin(2*π*fₛ*t*dt)) - θᵦ , undef, undef ,undef)

    xu_p = Point{Float64}(2000.,   el + θᵦ,  (Φₛ + (Δₐ/2)*sin(2*π*fₛ*t*dt)) + θᵦ, undef, undef, undef)
    xl_p = Point{Float64}(0.,      el - θᵦ,  (Φₛ + (Δₐ/2)*sin(2*π*fₛ*t*dt)) - θᵦ, undef, undef ,undef)

    result, chi_square = Integrate(corr, xl, xu, xl_p, xu_p, calls, it, dimensions, tempo )
    error = (chi_square/(it+1))^0.5
    println(t*dt," ", result, " ± ", error)

    append!(time_s, t*dt)
    append!(Coo, result)
    append!(std_err,error)

    if t*dt % 10 == 0
        global el -= 2*θᵦ
    end
    println(el, " ", Φₛ + (Δₐ/2)*sin(2*π*fₛ*t*dt))
    append!(el_arr, el)
    append!(az_arr, Φₛ + (Δₐ/2)*sin(2*π*fₛ*t*dt))
end

ord = sortperm(time_s)
time_s = time_s[ord]
Coo  = Coo[ord]

#
# # Print to file
# f = open(file_name, "w")
# for j = 1:min(length(Coo), length(time_s))
#     println(f, time_s[j], " ", Coo[j])
# end
# close(f)
x = Array{Float64}(undef, 0)
y = Array{Float64}(undef, 0)

x_w = Array{Float64}(undef, 0)
y_w = Array{Float64}(undef, 0)

time_arr = Array{Float64}(undef, 0)
Coo_arr  = Array{Float64}(undef, 0)

phi_cir = ((0:1000) ./ 1000) .* π
punti_cir = [Point{Float64}(2500, π/4., i, undef, undef, undef)() for i in phi_cir]

x_c = [i.x/2000 for i in punti_cir ]
y_c = [i.y/2000 for i in punti_cir ]

@gif for i = 1:length(az_arr)

    pun = Point{Float64}(2500., π/4., az_arr[i], undef, undef, undef )()
    append!(x, pun.x/2000)
    append!(y, pun.y/2000)

    append!(x_w, -1-(45/2000)*i*cos((π/2)+(π/7)) )
    append!(y_w, 0+(45/2000)*i*sin((π/2)+(π/7)) )

    p1=plot(x, y, xlims=(-1000/2000,1000/2000), ylims=(-1,1), arrow=arrow(), legend=false, ylabel="Y", xlabel="X")
    p1=plot!(x_c,y_c, xlims=(-1000/2000,1000/2000), ylims=(-1,1))

    for j = -5:5
        for k = 0:8
            p1=plot!(x_w .+ k*(500/2000) , y_w .- j*(500/2000), arrow=arrow(), legend = false)
        end
    end
    append!(time_arr, time_s[i])
    append!(Coo_arr, Coo[i]/findmax(Coo)[1])

    p2 = plot(time_arr, Coo_arr, xlims=(0,35), ylims=(0,1.2), arrow=arrow(),  ylabel="Autocorrelation C_00^0t [normalized unit]", xlabel="time [s]")

    plot(p1,p2, layout =(1,2))

end every 1

#plot(time_s, Coo/findmax(Coo)[1], seriestype=:scatter, ylims=(0,1.1), title="Test", xlabel = "Time [seconds]", ylabel="Autocorrelation Level [Normalized Units]", m=(1, :circle, 1), bg=RGB(1,1,1))#, yerr=(std_err ./ findmax(Coo)[1]))
#plot(hit)
#plot!(hit, seriestype=:scatter)
#
# using DSP
# period = periodogram(Coo, fs=1)
# Freq   = freq(period)
# A      = power(period)

#stop = Dates.value(Dates.now())
#println("Performance: ", ( stop - start ) / (1000 * (Step_stop - Step_start)), " sec/step" )

max_stars_repo_name: gruffoc/AFO

max_issues_repo_name: gruffoc/AFO

max_forks_repo_name: gruffoc/AFO


text:
### FULL MATRIX LOGISTIC REGRESSION ###
using LinearAlgebra
using Random
# using ReverseDiff
using Optim
using Plots
using SparseArrays
using ZigZagBoomerang
const ZZB = ZigZagBoomerang
#generate data

# partial derivative
include("./fulldesign.jl")
Random.seed!(1)
p = 100
N = 1000
w = 0.95
ξtrue = 2*(rand(p).> w)
(X, y) = generate_data(ξtrue, p, N)

N, d = size(X)

### Finding the mode ###
# Energy function logistic
# smooth component of the prior measure
const prec = 0.01
function pri(x)
    sqrt(prec)*exp(-prec*dot(x,x)*0.5)/sqrt(2π)
end

# gradient of the smooth component of the prior
function ∇pri(x)
    -prec.*x
end

function ∇pri(x, i)
    -prec*x[i]
end

# Ψ such that μ(dx) = \exp(-ψ)∏(dx_i + k^(-1) δ_0(dx_i))
function ϕ(x, A, y)
    sum(log.(1 .+ exp.(A*x)) .- y.*A*x) #- pri(x)
end

# jth term of the kth partial derivative
function ∇ϕs(x, k, j, A, y)
    A[j,k]*exp(dot(A[j,:], x))/(1 + exp(dot(A[j,:], x))) - y[j]*A[j,k] #- ∇pri(x, j)
end

include("preprocessing_logistic.jl")
∇ϕξref = [∇ϕ(ξref, k, X, y) for k in 1:length(ξref)]

### STANDARD STICKY ZIG-ZAG WITH CONTROL VARIATES
function ∂ϕ!(ξ, k, ∇ϕξref, ξref, ∇pri, y, X, N)
    j = rand(1:length(y)) # subsampling over sample size
    E_tilde = ∇ϕs(ξ, k, j, X, y)
    E_tilde′ = ∇ϕs(ξref, k, j, X, y) # could be pre-computed
    return ∇ϕξref[k] + N*(E_tilde - E_tilde′) - ∇pri(ξ, k)
end

struct MyBoundLog
    c::Float64
end

c = [MyBoundLog(0.0) for i in 1:d]
for k in 1:d
    C_k = 0.0
    for j in 1:N
        C_k = max(C_k, N*0.25*abs.(X[j,k]).*norm(X[j,:]))
    end
    c[k] = MyBoundLog(C_k)
end

# a and b
function ZZB.ab(G, i, x, θ, c::Vector{MyBoundLog}, F::ZigZag, ∇ϕξ′, ξ′, precision, args...)
    #max(0, θ[i]*∇ϕξ′[i]) + abs(θ[i])*c[i]*norm(ξ - ξ′), abs(θ[i])*c[i]*norm(θ)
    # with prior
    # max(0, θ[i]*(∇ϕξ′[i] + precision*x[i])) + abs(θ[i])*c[i].c*norm(x - ξ′), abs(θ[i])*c[i].c*norm(θ) + θ[i]^2*precision
    max(0, θ[i]*(∇ϕξ′[i] - ∇pri(x, i))) + abs(θ[i])*c[i].c*norm(x - ξ′), abs(θ[i])*c[i].c*norm(θ) - ∇pri(θ, i)^2
    # NB: if ξ and θ are sparse vectors, this implemention is not efficient
    # To Check θ[i]*(∇ϕξ′[i] + precision*x[i]))
end

F = ZigZag(ones(d,d), zeros(d))
adapt = false
su = true #strong upperbounds which continue to be upperbounds after
κ = getkappa(w, pri)
t0, x0, θ0, T = 0.0, randn(d), ones(d), 100.0
Ξ, (t, x, θ), (acc, num), c = ZZB.sspdmp(∂ϕ!, t0, x0, θ0, T, c, F, κ,
                                    ∇ϕξref, ξref, ∇pri, y, X, N;
                                    strong_upperbounds = su ,
                                     adapt = adapt)

error("")

first2((t,x)) = t => x[1:2]
sep(x) = first.(x), last.(x)
ts1, xs1 = sep(discretize(Ξ, 0.05))
p1 = plot(ts1, getindex.(xs1,1), leg = false)
savefig(p1, "figure1" )
p2 = plot(ts1, getindex.(xs1,7), leg = false)
savefig(p2, "figure2")
save
#
#
# #############################################################################
# ####### SUBSAMPLING OVER DIMENSIONS WITH CONTROL VARIATES ###################
# #############################################################################
# function gradϕ!(ξ, ∇ϕξref, ξref, precision, y, X, N, d)
#     k = rand(1:d) # subsampling over dimension
#     j = rand(1:N) # subsampling over sample size
#     E_tilde = ∇ϕs(ξ, k, j, X, y)
#     E_tilde′ = ∇ϕs(ξref, k, j, X, y)
#     return (∇ϕξref[k] + d*N*(E_tilde - E_tilde′) + precision*(ξ[k] - ξref[k]), k)
# end
# struct MyBoundSub
#     c::Float64
# end
# c = MyBoundSub(0.0)
# C_k = 0.0
# for k in 1:d
#     for j in 1:N
#         global C_k = max(C_k, N*d*0.25*abs.(X[j,k]).*norm(X[j,:]))
#     end
#     global c = MyBoundSub(C_k)
# end
#
# function ZZB.ab(x, θ, c::MyBoundSub, F::ZigZag, ∇ϕξ′, ξ′, precision, y, X, N, d, args...)
#     #max(0, θ[i]*∇ϕξ′[i]) + abs(θ[i])*c[i]*norm(ξ - ξ′), abs(θ[i])*c[i]*norm(θ)
#     # with prior
#     ref = max([θ[i]*(∇ϕξ′[i] + precision*x[i]) for i in 1:d]..., 0)
#     θstar = max([abs(θ[i]) for i in 1:d]...)
#     ref + θstar*c.c*norm(x - ξ′), θstar*c.c*norm(θ) + θstar^2*precision # NB: if ξ and θ are sparse vectors, this implemention is not efficient
# end
#
# Ξ, (t, x, θ), (acc, num), c = pdmp_sub(gradϕ!, t0, x0, θ0, T, c, F, ∇ϕξref, ξref, precision, y, X, N, d; adapt=false, factor=2.0)


text:
using Temporal, Dates, Test

# constant variables
N = 252*3
K = 4
data = cumsum(randn(N,K), dims=1)
dates = today()-Day(N-1):Day(1):today()
times = DateTime(dates[1]):Hour(1):DateTime(dates[1])+Hour(N-1)
X = TS(data, dates)
Y = TS(data, times)

test_scripts = [
    "basic.jl",
    "subsetting.jl",
    "mutating.jl",
    "combinations.jl",
    "operations.jl",
    "filter.jl",
    "ohlc.jl",
    "collapse.jl",
    "conversion.jl",
    "data.jl",
]

for test_script in test_scripts
    include(test_script)
end


text:
module Pmm
using ...Types, ...Utils
using DifferentialEquations, LinearAlgebra
#include("aerodynamic_coeff_fragment.jl")
export timeOfFlight!, trajectoryPMM!

    function pmm(du,u,p,t)
        w_bar, A, R, g₀,ω_bar, m,target,coeffAero,proj  = p

        vx= u[4]
        vy=u[5]
        vz=u[6]
        X₁ = u[1]
        X₂ = u[2]
        X₃ = u[3]

       u_bar =[u[4], u[5], u[6]]
       v_bar= u_bar - w_bar
        v = norm(v_bar)
        du_bar = [du[4], du[5],du[6]]
        #println(X₂)
        ρ =  Utils.density(X₂)
        #println("ρ", " ", ρ)
        α = atan(vz/vx)
        β = asin(vy/v)

        mach = Utils.machNumber(v,X₂)

        CD₀ = coeffAero(proj,mach)

        CD = CD₀

        FDx = - ρ*A/2*(CD)*v*vx #drag
        FDy = - ρ*A/2*(CD)*v*vy #drag
        FDz = - ρ*A/2*(CD)*v*vz #drag

        FGx = - g₀*X₁/R*m
        FGy = -g₀*(1-2*X₂/R)*m
        FGz = - g₀*X₃/R*m

        FCx = -2*cross(ω_bar,u_bar)[1]*m
        FCy = -2*cross(ω_bar,u_bar)[2]*m
        FCz = -2*cross(ω_bar,u_bar)[3]*m

        Fx =  FDx  +FGx +FCx
        Fy =  FDy  +FGy +FCy
        Fz =  FDz  +FGz +FCz

        du[1] = u[4] #x
        du[2] = u[5] #y
        du[3] = u[6] #z
        du[4] =Fx/m #ax
        du[5] = Fy/m#ay
        du[6] = Fz/m#az

    end

    function timeOfFlight!(p::Projectile2D, target::Target2D) #PmmZEroDrag
        #t = 0.0
        range = target.position[1]
        flight_time = 0.0

    x0 = p.position[1]
    y0 = p.position[2]
    v0x = p.velocity_m[1]
    v0y = p.velocity_m[2]
    u0 = [x0,y0,v0x,v0y]

    tspan = (0.0,tmax)
    prob = ODEProblem(pmm,u0,tspan)
    #sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8)
    integrator = init(prob, Tsit5(), reltol=1e-8, abstol=1e-8)
    oldx = integrator.u[1]
    oldy = integrator.u[2]
    oldu = integrator.u[3]
    oldv = integrator.u[4]
    while range > p.position[1]
        oldx = integrator.u[1]
        oldy = integrator.u[2]
        oldu = integrator.u[3]
        oldv = integrator.u[4]
    step!(integrator,dt,true)
    p.position[1]= integrator.u[1]
    #println(integrator.t, " ",p.position[1], " ", range)
    end
    p.position[1] = oldx
    p.position[2] = oldy
    p.velocity_m[1] = oldu
    p.velocity_m[2] = oldv
    flight_time = integrator.t-dt
    p.yaw = atan(p.velocity_m[2]/p.velocity_m[1])

    return flight_time
end

function condition(u,t,integrator) # Event when event_f(u,t) == 0
  #u[2]
  #target.x-u[1]
  #dist-u[1]
  integrator.p[7][1]-u[1]
end

affect!(integrator) = terminate!(integrator)

function trajectoryPMM!(u0, tspan, p, proj::AbstractPenetrator)
    #global dist = target.position[1]
    prob = ODEProblem(pmm,u0,tspan, p)
#cb = DiscreteCallback(condition,affect!)
#println(target.x)
    cb = ContinuousCallback(condition,affect!)
    sol = solve(prob,Tsit5(),callback=cb)
#sol = solve(prob,Tsit5(),callback=cb)
#sol = solve(prob,BS5(),callback=cb)
#sol = solve(prob,AutoTsit5(Rosenbrock23()),callback=cb)
 proj.position[1] = sol.u[end][1]
 proj.position[2] = sol.u[end][2]
 proj.position[3] = sol.u[end][3]

 proj.velocity[1] = sol.u[end][4]
 proj.velocity[2] = sol.u[end][5]
 proj.velocity[3] = sol.u[end][6]

 proj.tof = sol.t[end]
end

end


text:
## -----------  Logistic

"""
	$(SIGNATURES)

Admission probability = pMin + (pMax - pMin) / (1 + Q exp(-B (x - M))).

`M` is essentially a left-right shifter.
`B` is a slope parameter. Needs to be high b/c logistic maps [-Inf, Inf] -> [0, 1].
`Q` is mostly redundant and can be fixed.

Parameters may vary by college (as indicated by switches). But some colleges are open admission. They have no parameters.

All students get into the lowest college. Parameters only need to be set for colleges 2+.
"""
mutable struct AdmProbFctLogisticSwitches{F1} <: AbstractAdmProbFctSwitches{F1}
    pvec :: ParamVector
    nColleges :: Int
    nOpenColleges :: Int
end

mutable struct AdmProbFctLogistic{F1} <: AbstractAdmProbFct{F1}
    objId :: ObjectId
    switches :: AdmProbFctLogisticSwitches{F1}
    pMinV :: Vector{F1}
    pMaxV :: Vector{F1}
    qV :: Vector{F1}
    bV :: Vector{F1}
    mV :: Vector{F1}
end

# The order matters here!
param_names(af :: AdmProbFctLogisticSwitches{F1}) where F1 =
    param_names(typeof(af));
param_names(::Type{AdmProbFctLogisticSwitches{F1}}) where F1 =
    [:pMinV, :pMaxV, :qV, :bV, :mV];

logistic(x, pMin, pMax, q, b, m) =
    pMin .+ (pMax .- pMin) ./ (1.0 .+ q .* exp.(-b .* (x .- m)));

# Lazy.@forward AdmProbFctLogistic.switches (
#     by_college, by_college!, college_index, get_param, get_params, param_names
#     );

## -------  Constructors

# nc: total no of colleges
function make_test_admprob_fct_logistic_switches(nc)
    init_admprob_fct_logistic_switches(ObjectId(:test), nc);
end

function make_test_admprob_fct_logistic(nc)
    switches = make_test_admprob_fct_logistic_switches(nc);
    af = init_admprob_fct(switches);
    @assert validate_admprob_fct(af);
    return af
end

"""
	$(SIGNATURES)

Set up switches for logistic admission probability functions.
"""
function init_admprob_fct_logistic_switches(
    objId :: ObjectId,  nc :: Integer;
    byCollegeV :: Vector{Symbol} = [:bV, :mV],
    calibratedV :: Vector{Symbol} = [:bV, :mV]
    )
    nOpen = 1;
    make_admprob_fct_switches(objId, AdmProbFctLogisticSwitches{Float64},
        nc, nOpen, byCollegeV, calibratedV);
end

# notation from symbol table +++++
# Min admission prob for each college.
# Input is no of colleges that need parameters.
function init_pmin(nc :: Integer, isCalibrated :: Bool)
    sz = (nc, );
    return Param(:pMinV, "Min prob", "pMinV",
        fill(0.05, sz), fill(0.01, sz),
        fill(0.005, sz), fill(0.45, sz), isCalibrated);
end

function init_pmax(nc :: Integer, isCalibrated :: Bool)
    sz = (nc, );
    return Param(:pMaxV, "Max prob", "pMaxV",
        fill(0.95, sz), fill(0.99, sz),
        fill(0.5, sz), fill(0.995, sz), isCalibrated);
end

# Q is equivalent to M. So can be set to 1.
function init_q(nc :: Integer, isCalibrated :: Bool)
    sz = (nc, );
    pName = :qV;
    v = fill(1.0, sz);
    return Param(pName, "Logistic $pName", string(pName), v, v,
        0.1 .* v, 10.0 .* v, isCalibrated);
end

# Slope
# Logistic maps [-Inf, Inf] -> [0, 1]. So B must be large.
function init_b(nc :: Integer, isCalibrated :: Bool)
    sz = (nc, );
    pName = :bV;
    v = fill(5.0, sz);
    return Param(pName, "Logistic $pName", string(pName), v, v,
        fill(0.1, sz), fill(20.0, sz), isCalibrated);
end

# The input is a percentile. So M should also be on the order of [0, 1].
# This can vary across colleges (shifts curve left/right).
function init_m(nc :: Integer, isCalibrated :: Bool)
    sz = (nc, );
    pName = :mV;
    if nc > 1
        v = collect(LinRange(-0.4, 0.4, nc));
    else
        v = fill(0.4, sz);
    end
    sz = size(v);
    Param(pName, "Logistic $pName", string(pName), v, v,
        fill(-0.8, sz), fill(0.8, sz), isCalibrated);
end

function init_param(::Type{AdmProbFctLogisticSwitches{F1}}, pName, nc, isCalibrated) where F1
    if pName == :pMinV
        return init_pmin(nc, isCalibrated);
    elseif pName == :pMaxV
        return init_pmax(nc, isCalibrated);
    elseif pName == :mV
        return init_m(nc, isCalibrated);
    elseif pName == :bV
        return init_b(nc, isCalibrated);
    elseif pName == :qV
        return init_q(nc, isCalibrated);
    else
        error("Invalid: $pName");
    end
end

"""
	$(SIGNATURES)

Make the object that holds the admission probability functions for a set of colleges.
"""
function init_admprob_fct(switches :: AdmProbFctLogisticSwitches{F1}) where F1
    pV = [param_default_value(switches, pName)  for pName in param_names(switches)];
    af = AdmProbFctLogistic(get_object_id(switches), switches, pV...);
        # param_default_value(switches, :pMinV), param_default_value(switches, :pMaxV),
        # param_default_value(switches, :qV), param_default_value(switches, :bV),
        # param_default_value(switches, :mV));
    @assert validate_admprob_fct(af);
    return af
end

# stub ++++++
function validate_admprob_fct_switches(switches :: AdmProbFctLogisticSwitches{F1}) where F1
    isValid = true;
    return isValid;
end

function validate_admprob_fct(switches :: AdmProbFctLogistic{F1}) where F1
    isValid = true;
    return isValid;
end

## ----------  Access

# Maps a parameter name such as `:pMinV` into "pMin"
# function sym_from_name(pName)
#     pStr = string(pName);
#     @assert last(pStr) == 'V'  "Should end in V: $pName";
#     pStr2 = pStr[1 : (end-1)];
#     return pStr2
# end

## ---------  Admission prob for one college

"""
	$(SIGNATURES)

Make admission probability function for one college. Maps endowment percentile into [0, 1].
"""
function make_admprob_function(af :: AdmProbFctLogistic{F1},
    iCollege :: Integer) where F1

    if iCollege <= n_open_colleges(af)
        fct = x -> 1.0;
    else
        paramV = get_params(af, iCollege);
        fct = x -> logistic(x, paramV...);
    end
    return fct
end

"""
	$(SIGNATURES)

Probability of being admitted into a specific college (just based on its admissions prob function).
"""
function prob_admit(af :: AdmProbFctLogistic{F1}, iCollege :: Integer, hsGpa) where F1
    if iCollege <= n_open_colleges(af)
        probV = prob_admit_open(hsGpa);
    else
        probV = logistic(hsGpa, get_params(af, iCollege)...);
    end
    return probV
end

# --------------------


text:
const LAYOUT_TESTS = [
"""@cstruct T {
}""",

"""@cstruct T {
	c::Cchar
}""",

"""@cstruct T {
	s::Cshort
}""",

"""@cstruct T {
	i::Cint
}""",

"""@cstruct T {
	l::Clong
}""",

"""@cstruct T {
	ll::Clonglong
}""",

"""@cstruct T {
	x::Cchar
	c::Cchar
}""",

"""@cstruct T {
	x::Cchar
	s::Cshort
}""",

"""@cstruct T {
	x::Cchar
	i::Cint
}""",

"""@cstruct T {
	x::Cchar
	l::Clong
}""",

"""@cstruct T {
	x::Cchar
	ll::Clonglong
}""",

"""@cstruct T {
	x::Cchar
	f::Cfloat
}""",

"""@cstruct T {
	x::Cchar
	d::Cdouble
}""",

"""@cstruct T {
	x::Cchar
	ld::Clongdouble
}""",

"""@cstruct T {
	c::Cchar:3
}""",

"""@cstruct T {
	s::Cshort:3
}""",

"""@cstruct T {
	i::Cint:3
}""",

"""@cstruct T {
	l::Clong:3
}""",

"""@cstruct T {
	ll::Clonglong:3
}""",

"""@cstruct T {
	x::Cchar:3
	c::Cchar
}""",

"""@cstruct T {
	x::Cchar:3
	s::Cshort
}""",

"""@cstruct T {
	x::Cchar:3
	i::Cint
}""",

"""@cstruct T {
	x::Cchar:3
	l::Clong
}""",

"""@cstruct T {
	x::Cchar:3
	l::Clonglong
}""",

"""@cstruct T {
	x::Cchar:3
	f::Cfloat
}""",

"""@cstruct T {
	x::Cchar:3
	d::Cdouble
}""",

"""@cstruct T {
	x::Cchar:3
	ld::Clongdouble
}""",

"""@cstruct T {
	c::Cchar:3
	x::Cchar
}""",

"""@cstruct T {
	s::Cshort:3
	x::Cchar
}""",

"""@cstruct T {
	i::Cint:3
	x::Cchar
}""",

"""@cstruct T {
	l::Clong:3
	x::Cchar
}""",

"""@cstruct T {
	l::Clonglong:3
	x::Cchar
}""",

"""@cstruct T {
	x::Cchar:3
	c::Cchar:3
}""",

"""@cstruct T {
	x::Cchar:3
	s::Cshort:3
}""",

"""@cstruct T {
	x::Cchar:3
	i::Cint:3
}""",

"""@cstruct T {
	x::Cchar:3
	l::Clong:3
}""",

"""@cstruct T {
	x::Cchar:3
	ll::Clonglong:3
}""",

"""@cstruct T {
	x::Cchar:3
	c::Cchar:5
}""",

"""@cstruct T {
	x::Cchar:3
	s::Cshort:13
}""",

"""@cstruct T {
	x::Cchar:3
	i::Cint:29
}""",

"""@cstruct T {
	x::Cchar:3
	l::Clong:29
}""",

"""@cstruct T {
	x::Cchar:3
	ll::Clonglong:61
}""",

"""@cstruct T {
	x::Cchar:3
	@cstruct {
		c::Cchar:3
	}
}""",

"""@cstruct T {
	x::Cchar:3
	@cstruct {
		s::Cshort:3
	}
}""",

"""@cstruct T {
	x::Cchar:3
	@cstruct {
		i::Cint:3
	}
}""",

"""@cstruct T {
	x::Cchar:3
	@cstruct {
		l::Clong:3
	}
}""",

"""@cstruct T {
	x::Cchar:3
	@cstruct {
		ll::Clonglong:3
	}
}""",

"""@cstruct T {
	x::Cchar:3
	@cstruct {
		c::Cchar:3
	}
}""",

"""@cstruct T {
	x::Cshort:3
	@cstruct {
		s::Cshort:3
	}
}""",

"""@cstruct T {
	x::Cint:3
	@cstruct {
		i::Cint:3
	}
}""",

"""@cstruct T {
	x::Clong:3
	@cstruct {
		l::Clong:3
	}
}""",

"""@cstruct T {
	x::Clonglong:3
	@cstruct {
		ll::Clonglong:3
	}
}""",

"""@cstruct T {
	x::Cchar:3
	@cunion {
		c::Cchar:3
	}
}""",

"""@cstruct T {
	x::Cchar:3
	@cunion {
		s::Cshort:3
	}
}""",

"""@cstruct T {
	x::Cchar:3
	@cunion {
		i::Cint:3
	}
}""",

"""@cstruct T {
	x::Cchar:3
	@cunion {
		l::Clong:3
	}
}""",

"""@cstruct T {
	x::Cchar:3
	@cunion {
		ll::Clonglong:3
	}
}""",

"""@cstruct T {
	x::Cchar:3
	@cunion {
		c::Cchar:3
	}
}""",

"""@cstruct T {
	x::Cshort:3
	@cunion {
		s::Cshort:3
	}
}""",

"""@cstruct T {
	x::Cint:3
	@cunion {
		i::Cint:3
	}
}""",

"""@cstruct T {
	x::Clong:3
	@cunion {
		l::Clong:3
	}
}""",

"""@cstruct T {
	x::Clonglong:3
	@cunion {
		ll::Clonglong:3
	}
}""",

"""@cstruct T {
	x::Cchar:3
	@calign sizeof(Cchar)*2
	@cstruct {
		c::Cchar:3
	}
}""",

"""@cstruct T {
	x::Cshort:3
	@calign sizeof(Cshort)*2
	@cstruct {
		s::Cshort:3
	}
}""",

"""@cstruct T {
	x::Cint:3
	@calign sizeof(Cint)*2
	@cstruct {
		i::Cint:3
	}
}""",

"""@cstruct T {
	x::Clong:3
	@calign sizeof(Clong)*2
	@cstruct {
		l::Clong:3
	}
}""",

"""@cstruct T {
	x::Clonglong:3
	@calign sizeof(Clonglong)*2
	@cstruct {
		ll::Clonglong:3
	}
}""",

"""@cstruct T {
	x::Cchar:3
	@cstruct {
		@calign sizeof(Cint)*2
		c::Cchar
	}
}""",

"""@cstruct T {
	x::Cshort:3
	@cstruct {
		@calign sizeof(Cint)*2
		c::Cchar
	}
}""",

"""@cstruct T {
	x::Cint:3
	@cstruct {
		@calign sizeof(Cint)*2
		c::Cchar
	}
}""",

"""@cstruct T {
	x::Clong:3
	@cstruct {
		@calign sizeof(Cint)*2
		c::Cchar
	}
}""",

"""@cstruct T {
	x::Clonglong:3
	@cstruct {
		@calign sizeof(Cint)*2
		c::Cchar
	}
}""",

"""@cstruct T {
	x::Cchar:3
	@cstruct {
		c::Cchar:3
	} __packed__
}""",

"""@cstruct T {
	x::Cchar:3
	@cstruct {
		s::Cshort:3
	} __packed__
}""",

"""@cstruct T {
	x::Cchar:3
	@cstruct {
		i::Cint:3
	} __packed__
}""",

"""@cstruct T {
	x::Cchar:3
	@cstruct {
		l::Clong:3
	} __packed__
}""",

"""@cstruct T {
	x::Cchar:3
	@cstruct {
		ll::Clonglong:3
	} __packed__
}""",

"""@cstruct T {
	x::Cchar:3
	@cstruct {
		c::Cchar:3
	} __packed__
}""",

"""@cstruct T {
	x::Cshort:3
	@cstruct {
		s::Cshort:3
	} __packed__
}""",

"""@cstruct T {
	x::Cint:3
	@cstruct {
		i::Cint:3
	} __packed__
}""",

"""@cstruct T {
	x::Clong:3
	@cstruct {
		l::Clong:3
	} __packed__
}""",

"""@cstruct T {
	x::Clonglong:3
	@cstruct {
		ll::Clonglong:3
	} __packed__
}""",

"""@cstruct T {
	x::Cchar:3
	@cunion {
		c::Cchar:3
	} __packed__
}""",

"""@cstruct T {
	x::Cchar:3
	@cunion {
		s::Cshort:3
	} __packed__
}""",

"""@cstruct T {
	x::Cchar:3
	@cunion {
		i::Cint:3
	} __packed__
}""",

"""@cstruct T {
	x::Cchar:3
	@cunion {
		l::Clong:3
	} __packed__
}""",

"""@cstruct T {
	x::Cchar:3
	@cunion {
		ll::Clonglong:3
	} __packed__
}""",

"""@cstruct T {
	x::Cchar:3
	@cunion {
		c::Cchar:3
	} __packed__
}""",

"""@cstruct T {
	x::Cshort:3
	@cunion {
		s::Cshort:3
	} __packed__
}""",

"""@cstruct T {
	x::Cint:3
	@cunion {
		i::Cint:3
	} __packed__
}""",

"""@cstruct T {
	x::Clong:3
	@cunion {
		l::Clong:3
	} __packed__
}""",

"""@cstruct T {
	x::Clonglong:3
	@cunion {
		ll::Clonglong:3
	} __packed__
}""",

"""@cstruct T {
	x::Cchar:3
	@calign sizeof(Cchar)*2
	@cstruct {
		c::Cchar:3
	} __packed__
}""",

"""@cstruct T {
	x::Cshort:3
	@calign sizeof(Cshort)*2
	@cstruct {
		s::Cshort:3
	} __packed__
}""",

"""@cstruct T {
	x::Cint:3
	@calign sizeof(Cint)*2
	@cstruct {
		i::Cint:3
	} __packed__
}""",

"""@cstruct T {
	x::Clong:3
	@calign sizeof(Clong)*2
	@cstruct {
		l::Clong:3
	} __packed__
}""",

"""@cstruct T {
	x::Clonglong:3
	@calign sizeof(Clonglong)*2
	@cstruct {
		ll::Clonglong:3
	} __packed__
}""",

"""@cstruct T {
	x::Cchar:3
	@cstruct {
		@calign sizeof(Cint)*2
		c::Cchar
	} __packed__
}""",

"""@cstruct T {
	x::Cshort:3
	@cstruct {
		@calign sizeof(Cint)*2
		c::Cchar
	} __packed__
}""",

"""@cstruct T {
	x::Cint:3
	@cstruct {
		@calign sizeof(Cint)*2
		c::Cchar
	} __packed__
}""",

"""@cstruct T {
	x::Clong:3
	@cstruct {
		@calign sizeof(Cint)*2
		c::Cchar
	} __packed__
}""",

"""@cstruct T {
	x::Clonglong:3
	@cstruct {
		@calign sizeof(Cint)*2
		c::Cchar
	} __packed__
}""",

"""@cstruct T {
	e::@cenum {
		VAL1,
		VAL2,
		VAL3,
	}
}""",

"""@cstruct T {
	e::@cenum {
		VAL1,
		VAL2 = VAL1,
		VAL3 = 1 << 10,
	}
}""",

"""@cstruct T {
	e::@cenum {
		VAL1 = -1,
		VAL2 = 0,
		VAL3 = 1,
	}
}""",

"""@cstruct T {
	e::@cenum {
		VAL1 = -1,
		VAL2 = 0,
		VAL3 = 0xffffffff,
	}
}""",

"""@cstruct T {
	x::Cchar
	e::@cenum {
		VAL1,
		VAL2,
		VAL3,
	}
}""",

"""@cstruct T {
	x::Cchar
	e::@cenum {
		VAL1,
		VAL2 = VAL1,
		VAL3 = 1 << 10,
	}
}""",

"""@cstruct T {
	x::Cchar
	e::@cenum {
		VAL1 = -1,
		VAL2 = 0,
		VAL3 = 1,
	}
}""",

"""@cstruct T {
	x::Cchar
	e::@cenum {
		VAL1 = -1,
		VAL2 = 0,
		VAL3 = 0xffffffff,
	}
}""",

"""@cstruct T {
	x::Cchar:3
	e::@cenum {
		VAL1,
		VAL2,
		VAL3,
	}
}""",

"""@cstruct T {
	x::Cchar:3
	e::@cenum {
		VAL1,
		VAL2 = VAL1,
		VAL3 = 1 << 10,
	}
}""",

"""@cstruct T {
	x::Cchar:3
	e::@cenum {
		VAL1 = -1,
		VAL2 = 0,
		VAL3 = 1,
	}
}""",

"""@cstruct T {
	x::Cchar:3
	e::@cenum {
		VAL1 = -1,
		VAL2 = 0,
		VAL3 = 0xffffffff,
	}
}""",

"""@cstruct T {
	e::@cenum {
		VAL1,
		VAL2,
		VAL3,
	} __packed__
}""",

"""@cstruct T {
	e::@cenum {
		VAL1,
		VAL2 = VAL1,
		VAL3 = 1 << 10,
	} __packed__
}""",

"""@cstruct T {
	e::@cenum {
		VAL1 = -1,
		VAL2 = 0,
		VAL3 = 1,
	} __packed__
}""",

"""@cstruct T {
	e::@cenum {
		VAL1 = -1,
		VAL2 = 0,
		VAL3 = 0xffffffff,
	} __packed__
}""",

"""@cstruct T {
	x::Cchar
	e::@cenum {
		VAL1,
		VAL2,
		VAL3,
	} __packed__
}""",

"""@cstruct T {
	x::Cchar
	e::@cenum {
		VAL1,
		VAL2 = VAL1,
		VAL3 = 1 << 10,
	} __packed__
}""",

"""@cstruct T {
	x::Cchar
	e::@cenum {
		VAL1 = -1,
		VAL2 = 0,
		VAL3 = 1,
	} __packed__
}""",

"""@cstruct T {
	x::Cchar
	e::@cenum {
		VAL1 = -1,
		VAL2 = 0,
		VAL3 = 0xffffffff,
	} __packed__
}""",

"""@cstruct T {
	x::Cchar:3
	e::@cenum {
		VAL1,
		VAL2,
		VAL3,
	} __packed__
}""",

"""@cstruct T {
	x::Cchar:3
	e::@cenum {
		VAL1,
		VAL2 = VAL1,
		VAL3 = 1 << 10,
	} __packed__
}""",

"""@cstruct T {
	x::Cchar:3
	e::@cenum {
		VAL1 = -1,
		VAL2 = 0,
		VAL3 = 1,
	} __packed__
}""",

"""@cstruct T {
	x::Cchar:3
	e::@cenum {
		VAL1 = -1,
		VAL2 = 0,
		VAL3 = 0xffffffff,
	} __packed__
}""",
]


text:
@testset "KLIEP -- $optlib" for optlib in [OptimLib, ConvexLib]
  for (i, (pair, rtol)) in enumerate([(pair₁, 2e-1), (pair₂, 4e-1)])
    d_nu, d_de = pair
    Random.seed!(123)
    x_nu, x_de = rand(d_nu, 1_000), rand(d_de, 500)

    # estimated density ratio
    σ, b = 1.0, 100
    r̂ = densratio(x_nu, x_de, KLIEP(σ=σ, b=b), optlib=optlib)

    # simplex constraints
    @test abs(mean(r̂) - 1) ≤ 1e-2
    @test all(r̂ .≤ Inf)

    r = pdf.(d_nu, x_de) ./ pdf.(d_de, x_de)
    @test r ≈ r̂ rtol=rtol

    if visualtests
      gr(size=(800, 800))
      @plottest plot_d_nu(pair, x_de, r̂) joinpath(datadir, "KLIEP-$optlib-$i.png") !istravis
    end
  end
end


text:
module Data
export List, Fun
include("TypedFn.jl")
@info TypedFn
using MLStyle.Data.TypedFn

include("List.jl")
end


text:
###############################################################################
#  0.14.0
###############################################################################

@deprecate(
    TimeArray(d::AbstractVector{D}, v::AbstractArray{T,N},
              c::Vector{S}, m::Any = nothing;
              args...) where {T,N,D<:TimeType,S<:AbstractString},
    TimeArray(d, v, Symbol.(c), m; args...)
)

@deprecate(
    TimeArray(d::D, v::AbstractArray{T,N},
              c::Vector{S}, m::Any = nothing;
              args...) where {T,N,D<:TimeType,S<:AbstractString},
    TimeArray(d, v, Symbol.(c), m; args...)
)

@deprecate getindex(ta::TimeArray, s::AbstractString)     getindex(ta, Symbol(s))
@deprecate getindex(ta::TimeArray, ss::AbstractString...) getindex(ta, Symbol.(ss)...)

@deprecate rename(ta::TimeArray, col::String)         rename(ta::TimeArray, Symbol(col))
@deprecate rename(ta::TimeArray, col::Vector{String}) rename(ta, Symbol.(col))

@deprecate(
    merge(ta1::TimeArray{T}, ta2::TimeArray, method::Symbol; kw...) where {T},
    merge(ta1, ta2; method = method, kw...))


text:
"""
    IMERGHalfHourly{ST<:AbstractString, DT<:TimeType} <: IMERGDataset

Object containing information on Half-Hourly IMERG datasets to be downloaded
"""
struct IMERGHalfHourly{ST<:AbstractString, DT<:TimeType} <: IMERGDataset
	npdID :: ST
	lname :: ST
	doi   :: ST
    dtbeg :: DT
    dtend :: DT
    sroot :: ST
    hroot :: ST
    fpref :: ST
    fsuff :: ST
end

"""
    IMERGEarlyHH(
        ST = String,
        DT = Date;
        dtbeg :: TimeType,
        dtend :: TimeType,
        sroot :: AbstractString = homedir(),
    ) -> npd :: IMERGHalfHourly{ST,DT}

Creates a `IMERGHalfHourly` dataset `npd` to retrieve datasets from the Near Real-Time Early processing runs for Half-Hourly output

Keyword Arguments
=================
- `dtbeg` : Date at which download / analysis of the dataset begins
- `dtend` : Date at which download / analysis of the dataset ends
- `sroot` : The directory in which the folder `imergearlyhh` will be created for data downloads, storage and analysis, default is the home directoy called by `homedir()`

The following fields in `npd` will be fixed as below:
- `npdID` : imergearlyhh
- `lname` : Early IMERG Half-Hourly
- `doi`   : 10.5067/GPM/IMERG/3B-HH-E/06
- `hroot` : https://gpm1.gesdisc.eosdis.nasa.gov/opendap/GPM_L3/GPM_3IMERGHHE.06
- `fpref` : 3B-HHR-E.MS.MRG.3IMERG
- `fsuff` : V06B.HDF5
"""
function IMERGEarlyHH(
    ST = String,
    DT = Date;
    dtbeg :: TimeType,
    dtend :: TimeType,
    sroot :: AbstractString = homedir(),
)

	@info "$(modulelog()) - Setting up data structure containing information on Early IMERG Half-Hourly data to be downloaded"

    fol = joinpath(sroot,"imergearlyhh"); if !isdir(fol); mkpath(fol) end

    return IMERGHalfHourly{ST,DT}(
		"imergearlyhh", "Early IMERG Half-Hourly", "10.5067/GPM/IMERG/3B-HH-E/06",
        dtbeg, dtend,
		joinpath(sroot,"imergearlyhh"),
        "https://gpm1.gesdisc.eosdis.nasa.gov/opendap/GPM_L3/GPM_3IMERGHHE.06",
        "3B-HHR-E.MS.MRG.3IMERG", "V06B.HDF5",
    )

end

"""
    IMERGLateHH(
        ST = String,
        DT = Date;
        dtbeg :: TimeType,
        dtend :: TimeType,
        sroot :: AbstractString = homedir(),
    ) -> npd :: IMERGHalfHourly{ST,DT}

Creates a `IMERGHalfHourly` dataset `npd` to retrieve datasets from the Near Real-Time Late processing runs for Half-Hourly output

Keyword Arguments
=================
- `dtbeg` : Date at which download / analysis of the dataset begins
- `dtend` : Date at which download / analysis of the dataset ends
- `sroot` : The directory in which the folder `imerglatehh` will be created for data downloads, storage and analysis, default is the home directoy called by `homedir()`

The following fields in `npd` will be fixed as below:
- `npdID` : imerglatehh
- `lname` : Late IMERG Half-Hourly
- `doi`   : 10.5067/GPM/IMERG/3B-HH-L/06
- `hroot` : https://gpm1.gesdisc.eosdis.nasa.gov/opendap/GPM_L3/GPM_3IMERGHHL.06
- `fpref` : 3B-HHR-L.MS.MRG.3IMERG
- `fsuff` : V06B.HDF5
"""
function IMERGLateHH(
    ST = String,
    DT = Date;
    dtbeg :: TimeType,
    dtend :: TimeType,
    sroot :: AbstractString = homedir(),
)

	@info "$(modulelog()) - Setting up data structure containing information on Late IMERG Half-Hourly data to be downloaded"

    fol = joinpath(sroot,"imerglatehh"); if !isdir(fol); mkpath(fol) end

    return IMERGHalfHourly{ST,DT}(
		"imerglatehh", "Late IMERG Half-Hourly", "10.5067/GPM/IMERG/3B-HH-L/06",
        dtbeg, dtend,
		joinpath(sroot,"imerglatehh"),
        "https://gpm1.gesdisc.eosdis.nasa.gov/opendap/GPM_L3/GPM_3IMERGHHL.06",
        "3B-HHR-L.MS.MRG.3IMERG", "V06B.HDF5",
    )

end

"""
    IMERGFinalHH(
        ST = String,
        DT = Date;
        dtbeg :: TimeType,
        dtend :: TimeType,
        sroot :: AbstractString = homedir(),
    ) -> npd :: IMERGHalfHourly{ST,DT}

Creates a `IMERGHalfHourly` dataset `npd` to retrieve datasets from the final post-processing runs for Half-Hourly output

Keyword Arguments
=================
- `dtbeg` : Date at which download / analysis of the dataset begins
- `dtend` : Date at which download / analysis of the dataset ends
- `sroot` : The directory in which the folder `imergfinalhh` will be created for data downloads, storage and analysis, default is the home directoy called by `homedir()`

The following fields in `npd` will be fixed as below:
- `npdID` : imergfinalhh
- `lname` : Final IMERG Half-Hourly
- `doi`   : 10.5067/GPM/IMERG/3B-HH/06
- `hroot` : https://gpm1.gesdisc.eosdis.nasa.gov/opendap/GPM_L3/GPM_3IMERGHH.06
- `fpref` : 3B-HHR.MS.MRG.3IMERG
- `fsuff` : V06B.HDF5
"""
function IMERGFinalHH(
    ST = String,
    DT = Date;
    dtbeg :: TimeType,
    dtend :: TimeType,
    sroot :: AbstractString = homedir(),
)

	@info "$(modulelog()) - Setting up data structure containing information on Final IMERG Half-Hourly data to be downloaded"

    fol = joinpath(sroot,"imergfinalhh"); if !isdir(fol); mkpath(fol) end

    return IMERGHalfHourly{ST,DT}(
		"imergfinalhh", "Final IMERG Half-Hourly", "10.5067/GPM/IMERG/3B-HH/06",
        dtbeg, dtend,
		joinpath(sroot,"imergfinalhh"),
        "https://gpm1.gesdisc.eosdis.nasa.gov/opendap/GPM_L3/GPM_3IMERGHH.06",
        "3B-HHR.MS.MRG.3IMERG", "V06B.HDF5",
    )

end

function show(io::IO, npd::IMERGHalfHourly{ST,DT}) where {ST<:AbstractString, DT<:TimeType}
    print(
		io,
		"The NASA Precipitation Dataset {$ST,$DT} has the following properties:\n",
		"    Dataset ID      (npdID) : ", npd.npdID, '\n',
		"    Logging Name    (lname) : ", npd.lname, '\n',
		"    DOI URL          (doi)  : ", npd.doi,   '\n',
		"    Data Directory  (sroot) : ", npd.sroot, '\n',
		"    Date Begin      (dtbeg) : ", npd.dtbeg, '\n',
		"    Date End        (dtend) : ", npd.dtend, '\n',
		"    Timestep                : 30 minutes\n",
        "    Data Resolution         : 0.1º\n",
        "    Data Server     (hroot) : ", npd.hroot, '\n',
	)
end

max_stars_repo_name: natgeo-wong/NASAPrecipitation

max_issues_repo_name: natgeo-wong/NASAPrecipitation

max_forks_repo_name: natgeo-wong/NASAPrecipitation


text:
using DifferentialEquations
using Roots
using PyCall
pygui(:qt)
using PyPlot
pygui(true)

include("YfindNP.jl")

Energy=.25
# Energy=0.1251;
# Energy=0.1249;
H=(2*Energy)^2
h=sqrt(H)
t_end=10000
barrier=20
max_hit=500

Hamil(XX,YY,QQ,PP)=( (QQ-XX)^2+(PP-YY)^2 )*( (QQ+XX)^2+(PP+YY)^2 )/((PP^4+2*PP^2*(XX^2-1)+(1+XX^2)^2 )*(QQ^4+2*QQ^2*(YY^2+1)+(YY^2-1)^2 ))
H_test(u)=( (u[3]-u[1])^2+(u[2]-u[4])^2 )*( (u[3]+u[1])^2+(u[2]+u[4])^2 )/ ((u[2]^4+2*u[2]^2*(u[1]^2-1)+(1+u[1]^2)^2 )*(u[3]^4+2*u[3]^2*(u[4]^2+1)+(u[4]^2-1)^2 ))

condition(u,t,integrator)= u[5]>p[1] || maximum([abs(u[1]),abs(u[2]),abs(u[3]),abs(u[4])])>p[2]
function condition2(u,t,integrator) # Event when event_f(u,t) == 0
 u[1]
end
affect!(integrator) = terminate!(integrator)
function affect2!(integrator)
integrator.u[5]=integrator.u[5]+1
 end
cb2 = ContinuousCallback(condition2,affect2!,nothing)
cb1 = DiscreteCallback(condition,affect!)

cb=CallbackSet(cb2,cb1)

ODE1(z,w)=conj(  im * w.*( 1 ./(z.^2-w.^2)+1 ./(1+w.^2) ))
ODE2(z,w)=conj(  im * z.*( 1 ./(w.^2-z.^2)+1 ./(1+z.^2) ))
function Eq_of_M(du,u,p,t)
    du[1]=real(ODE1(u[1]+im*u[2],u[3]+im*u[4]))
    du[2]=imag(ODE1(u[1]+im*u[2],u[3]+im*u[4]))
    du[3]=real(ODE2(u[1]+im*u[2],u[3]+im*u[4]))
    du[4]=imag(ODE2(u[1]+im*u[2],u[3]+im*u[4]))
    du[5]=0
    return
end
# refernce
# z=X+iP W=Q+iY . Plane Q=0, P=0
# z=u[1]+i u[2] ; w=u[3]+i u[4]

#Q=0, have P find Y

Q=.1
P=.1
Y=YfindNP(Q,P,H)
u0=zeros(5)
u0[1]=0 #X
u0[2]=P #P
u0[3]=Q #Q
u0[4]=Y[1] #Y
u0[5]=0
TEST1_0=H_test(u0)/H
p=zeros(2)
p[1]=max_hit
p[2]=barrier
prob = ODEProblem(Eq_of_M,u0,(0., t_end),p)
# t,A=solve(prob, RK4(),reltol=1e-15,abstol=1e-16,maxiters=1e15)

sol=solve(prob, Vern9(),reltol=1e-13,abstol=1e-15,maxiters=1e20,callback=cb)
figure()
# =1e-16,
# figure()
A=sol[:,:]
plot3D(A[3,:],A[2,:],A[1,:],linewidth=.5)

#testing if energy is conserved at end
uf=zeros(4)
uf[1]=A[1,end]
uf[2]=A[2,end]
uf[3]=A[3,end]
uf[4]=A[4,end]

#finding maxes for exit criteria
u1max=maximum(abs.(A[1,:]))
u2max=maximum(abs.(A[2,:]))
u3max=maximum(abs.(A[3,:]))
u4max=maximum(abs.(A[4,:]))
N=length(A[1,:])
H_v=zeros(N)
for k=1:N
    um=zeros(4)
    um[1]=A[1,k]
    um[2]=A[2,k]
    um[3]=A[3,k]
    um[4]=A[4,k]
    H_v[k]=H_test(um)/H
end

TEST1_f=H_test(uf)/H
# figure()
# plot(1:N, log.(H_v))


text:
using Documenter
using MultipleTesting

makedocs(
    modules = [MultipleTesting],
    format = :html,
    sitename = "MultipleTesting.jl",
    doctest = true,
    checkdocs = :exports,
    linkcheck = true,
    html_disable_git = true, # disable source and edit links to github
    html_canonical = "https://juliangehring.github.io/MultipleTesting.jl/stable/",
    pages = [
        "Home" => "index.md",
        "Library" => Any[
            "adjustment.md",
            "combination.md",
            "pi0.md",
            "higher-criticism.md",
            "models.md",
        ]
    ]
)

deploydocs(
    repo = "github.com/juliangehring/MultipleTesting.jl.git",
    julia = "0.7",
    target = "build",
    deps   = nothing,
    make   = nothing
)


text:
module VirtualTestDrive

using LibExpat
using Compat
using DataFrames

#import Base: TcpSocket, ==
import Base: TCPSocket, == #Xiaobai

macro lintpragma(s)
end

export
    ViresConnection,
    SCPMessage,
    DriverBehavior,

    BRISK_DRIVER,
    COMFORTABLE_DRIVER,
    DEFAULT_DRIVER,
    HASTY_DRIVER,
    INSECURE_DRIVER,

    find_vtd_directory,
    start_vires_vtd_tasks,
    stop_vires_components,
    rdbsniffer,

    load_and_init,
    record_scenario_run,
    set_driver_behavior,
    convert_recording_to_csv,
    register_to_authorize_framesteps!,
    unregister_to_authorize_framesteps!,

    create_camera_directly_above_ego,
    create_camera_relative_to_ego,

    step_with_udp_checking_frameid,
    step_with_udp_authorization,
    step_with_udp_explicit_stepwidth,

    write_and_wait_for_mirrored_message,

    get_xml_traffic_action_autonomous,

    idle_and_print_messages

    include("viresconnection.jl")
    include("constants.jl")
    include("rdb_definitions.jl")
    include("rdb_package_gen.jl")
    include("rdb.jl")
    include("system.jl")
    include("driver_definition.jl")
    include("io.jl")
    include("xml_utils.jl")
    include("scp_xml.jl")
    include("outward_facing.jl")
    include("virescsv.jl")

    include("rdb_process.jl")
end # module

max_stars_repo_name: sisl/VirtualTestDrive

max_issues_repo_name: sisl/VirtualTestDrive

max_forks_repo_name: sisl/VirtualTestDrive


text:
using Optim

struct HW
	m::Int64
    α::Float64
    β::Float64
    γ::Float64
    l0::Float64
    b0::Float64
    s0::Array{Float64, 1}

    HW(m::Int64) = new(m, 0, 0, 0, 0, 0, repeat([0], m))
	HW(m::Int64, α::Number, β::Number, γ::Number, l0::Number, b0::Number, s0::Array{Float64, 1}) =
		m == length(s0) ? new(m, Float64(α), Float64(β), Float64(γ), Float64(l0), Float64(b0),
		Float64.(s0)) : error("length(s0) must be equal to m")
end

function loss(model::HW, time_series)
	α = model.α
	β = model.β
	γ = model.γ
	l0 = model.l0
	b0 = model.b0
	m = model.m
	s = model.s0

	N = length(time_series)
	l_t = 0
	b_t = 0
	l_t_ = 0 #Variable to save l(t-1)
	b_t_ = 0 #Variable to save b(t-1)
	s_ = 0

	loss = 0

	for i in 0:(N - 1)
		if i == 0
			l_t = l0
			b_t = b0
		else
			l_t = (time_series[i] - s_) * α + (l_t_ + b_t_) * (1 - α)
			b_t = β * (l_t - l_t_) + (1 - β) * b_t_
		end
		l_t_ = l_t
		b_t_ = b_t
		s_ = s[i%m + 1]

		y_pred = l_t + b_t + s[i%m + 1]

		s[i%m + 1] = γ * (time_series[i + 1] - l_t_ - b_t_) + (1 - γ) * s[i%m + 1]

		loss += (time_series[i + 1] - y_pred)^2
	end

	return loss
end

function forecast(model::HW, time_serie, n_pred)
	α = model.α
	β = model.β
	γ = model.γ
	l0 = model.l0
	b0 = model.b0
	m = model.m
	s = model.s0

	N = length(time_serie)
	l_t = 0
	b_t = 0
	l_t_ = 0 #Variable to save l(t-1)
	b_t_ = 0 #Variable to save b(t-1)
	s_ = 0

	pred = []

	for i in 0:(N - 1)
		if i == 0
			l_t = l0
			b_t = b0
		else
			l_t = (time_serie[i] - s_) * α + (l_t_ + b_t_) * (1 - α)
			b_t = β * (l_t - l_t_) + (1 - β) * b_t_
		end
		l_t_ = l_t
		b_t_ = b_t
		s_ = s[i%m + 1]

		s[i%m + 1] = γ * (time_serie[i + 1] - l_t_ - b_t_) + (1 - γ) * s[i%m + 1]
	end

	l_t = (time_serie[end] - s_) * α + (l_t + b_t) * (1 - α)
	b_t = β * (l_t - l_t_) + (1 - β) * b_t_

	for i in N:(N+n_pred - 1)
		y_pred = l_t + b_t*(i-N+1) + s[i%m + 1]
		#The trend has to be added as many times as periods we want to forecast.
		push!(pred, y_pred)
	end

	return pred
end


text:
module SubdirTests
import ..Pkg # ensure we are using the correct Pkg

using Pkg, UUIDs, Test
using Pkg.REPLMode: pkgstr
using Pkg.Types: PackageSpec

using ..Utils

# Derived from RegistryTools' gitcmd.
function gitcmd(path::AbstractString)
    Cmd(["git", "-C", path, "-c", "user.name=RegistratorTests",
         "-c", "user.email=ci@juliacomputing.com"])
end

# Create a repository containing two packages in different
# subdirectories, `Package` and `Dep`, where the former depends on the
# latter. Return the tree hashes for the two packages.
function setup_packages_repository(dir)
    package_dir = joinpath(dir, "julia")
    mkpath(joinpath(package_dir, "src"))
    write(joinpath(package_dir, "Project.toml"), """
        name = "Package"
        uuid = "408b23ff-74ea-48c4-abc7-a671b41e2073"
        version = "1.0.0"

        [deps]
        Dep = "d43cb7ef-9818-40d3-bb27-28fb4aa46cc5"
        """)
    write(joinpath(package_dir, "src", "Package.jl"), """
        module Package end
        """)

    dep_dir = joinpath(dir, "dependencies", "Dep")
    mkpath(joinpath(dep_dir, "src"))
    write(joinpath(dep_dir, "Project.toml"), """
        name = "Dep"
        uuid = "d43cb7ef-9818-40d3-bb27-28fb4aa46cc5"
        version = "1.0.0"
        """)
    write(joinpath(dep_dir, "src", "Dep.jl"), """
        module Dep end
        """)

    git = gitcmd(dir)
    run(`$git init -q`)
    run(`$git add .`)
    run(`$git commit -qm 'Create repository.'`)
    package_tree_hash = readchomp(`$git rev-parse HEAD:julia`)
    dep_tree_hash = readchomp(`$git rev-parse HEAD:dependencies/Dep`)
    return package_tree_hash, dep_tree_hash
end

# Convert a path into a file URL.
function make_file_url(path)
    # Turn the slashes on Windows. In case the path starts with a
    # drive letter, an extra slash will be needed in the file URL.
    path = replace(path, "\\" => "/")
    if !startswith(path, "/")
        path = "/" * path
    end
    return "file://$(path)"
end

# Create a registry with the two packages `Package` and `Dep`.
function setup_registry(dir, packages_dir_url, package_tree_hash, dep_tree_hash)
    package_path = joinpath(dir, "P", "Package")
    dep_path = joinpath(dir, "D", "Dep")
    mkpath(package_path)
    mkpath(dep_path)
    write(joinpath(dir, "Registry.toml"), """
        name = "Registry"
        uuid = "cade28e2-3b52-4f58-aeba-0b1386f9894b"
        repo = "https://github.com"
        [packages]
        408b23ff-74ea-48c4-abc7-a671b41e2073 = { name = "Package", path = "P/Package" }
        d43cb7ef-9818-40d3-bb27-28fb4aa46cc5 = { name = "Dep", path = "D/Dep" }
        """)
    write(joinpath(package_path, "Package.toml"), """
        name = "Package"
        uuid = "408b23ff-74ea-48c4-abc7-a671b41e2073"
        repo = "$(packages_dir_url)"
        subdir = "julia"
        """)
    write(joinpath(package_path, "Versions.toml"), """
        ["1.0.0"]
        git-tree-sha1 = "$(package_tree_hash)"
        """)
    write(joinpath(package_path, "Deps.toml"), """
        [1]
        Dep = "d43cb7ef-9818-40d3-bb27-28fb4aa46cc5"
        """)

    write(joinpath(dep_path, "Package.toml"), """
        name = "Dep"
        uuid = "d43cb7ef-9818-40d3-bb27-28fb4aa46cc5"
        repo = "$(packages_dir_url)"
        subdir = "dependencies/Dep"
        """)
    write(joinpath(dep_path, "Versions.toml"), """
        ["1.0.0"]
        git-tree-sha1 = "$(dep_tree_hash)"
        """)

    git = gitcmd(dir)
    run(`$git init -q`)
    run(`$git add .`)
    run(`$git commit -qm 'Create repository.'`)
end

@testset "subdir" begin
    temp_pkg_dir() do depot
        # Apparently the working directory can turn out to be a
        # removed directory when getting here, which doesn't go well
        # with the `pkg"add ..."` calls. Just set it to something that
        # exists.
        cd(@__DIR__)
        # Setup a repository with two packages and a registry where
        # these packages are registered.
        packages_dir = mktempdir()
        registry_dir = mktempdir()
        packages_dir_url = make_file_url(packages_dir)
        tree_hashes = setup_packages_repository(packages_dir)
        setup_registry(registry_dir, packages_dir_url, tree_hashes...)
        pkgstr("registry add $(registry_dir)")
        dep = (name="Dep", uuid=UUID("d43cb7ef-9818-40d3-bb27-28fb4aa46cc5"))

        # Ordinary add from registry.
        pkg"add Package"
        @test isinstalled("Package")
        @test !isinstalled("Dep")
        @test isinstalled(dep)
        pkg"rm Package"

        pkg"add Dep"
        @test !isinstalled("Package")
        @test isinstalled("Dep")
        pkg"rm Dep"

        # Add version from registry.
        pkg"add Package@1.0.0"
        @test isinstalled("Package")
        @test !isinstalled("Dep")
        @test isinstalled(dep)
        pkg"rm Package"

        pkg"add Dep@1.0.0"
        @test !isinstalled("Package")
        @test isinstalled("Dep")
        pkg"rm Dep"

        # Add branch from registry.
        pkg"add Package#master"
        @test isinstalled("Package")
        @test !isinstalled("Dep")
        @test isinstalled(dep)
        pkg"rm Package"

        pkg"add Dep#master"
        @test !isinstalled("Package")
        @test isinstalled("Dep")
        pkg"rm Dep"

        # Develop from registry.
        pkg"develop Package"
        @test isinstalled("Package")
        @test !isinstalled("Dep")
        @test isinstalled(dep)
        pkg"rm Package"

        pkg"develop Dep"
        @test !isinstalled("Package")
        @test isinstalled("Dep")
        pkg"rm Dep"

        # Add from path.
        Pkg.add(Pkg.PackageSpec(path=packages_dir, subdir="julia"))
        @test isinstalled("Package")
        @test !isinstalled("Dep")
        @test isinstalled(dep)
        pkg"rm Package"

        Pkg.add(Pkg.PackageSpec(path=packages_dir, subdir="dependencies/Dep"))
        @test !isinstalled("Package")
        @test isinstalled("Dep")
        pkg"rm Dep"

        # Add from path, REPL subdir syntax.
        pkgstr("add $(packages_dir):julia")
        @test isinstalled("Package")
        @test !isinstalled("Dep")
        @test isinstalled(dep)
        pkg"rm Package"

        pkgstr("add $(packages_dir):dependencies/Dep")
        @test !isinstalled("Package")
        @test isinstalled("Dep")
        pkg"rm Dep"

        # Add from path at branch.
        Pkg.add(Pkg.PackageSpec(path=packages_dir, subdir="julia", rev="master"))
        @test isinstalled("Package")
        @test !isinstalled("Dep")
        @test isinstalled(dep)
        pkg"rm Package"

        Pkg.add(Pkg.PackageSpec(path=packages_dir, subdir="dependencies/Dep", rev="master"))
        @test !isinstalled("Package")
        @test isinstalled("Dep")
        pkg"rm Dep"

        # Add from path at branch, REPL subdir syntax
        pkgstr("add $(packages_dir):julia#master")
        @test isinstalled("Package")
        @test !isinstalled("Dep")
        @test isinstalled(dep)
        pkg"rm Package"

        pkgstr("add $(packages_dir):dependencies/Dep#master")
        @test !isinstalled("Package")
        @test isinstalled("Dep")
        pkg"rm Dep"

        # Develop from path.
        Pkg.develop(Pkg.PackageSpec(path=packages_dir, subdir="julia"))
        @test isinstalled("Package")
        @test !isinstalled("Dep")
        @test isinstalled(dep)
        pkg"rm Package"

        Pkg.develop(Pkg.PackageSpec(path=packages_dir, subdir="dependencies/Dep"))
        @test !isinstalled("Package")
        @test isinstalled("Dep")
        pkg"rm Dep"

        # Develop from path, REPL subdir syntax.
        pkgstr("develop $(packages_dir):julia")
        @test isinstalled("Package")
        @test !isinstalled("Dep")
        @test isinstalled(dep)
        pkg"rm Package"

        pkgstr("develop $(packages_dir):dependencies/Dep")
        @test !isinstalled("Package")
        @test isinstalled("Dep")
        pkg"rm Dep"

        # Add from url.
        Pkg.add(Pkg.PackageSpec(url=packages_dir_url, subdir="julia"))
        @test isinstalled("Package")
        @test !isinstalled("Dep")
        @test isinstalled(dep)
        pkg"rm Package"

        Pkg.add(Pkg.PackageSpec(url=packages_dir_url, subdir="dependencies/Dep"))
        @test !isinstalled("Package")
        @test isinstalled("Dep")
        pkg"rm Dep"

        # Add from url, REPL subdir syntax.
        pkgstr("add $(packages_dir_url):julia")
        @test isinstalled("Package")
        @test !isinstalled("Dep")
        @test isinstalled(dep)
        pkg"rm Package"

        pkgstr("add $(packages_dir_url):dependencies/Dep")
        @test !isinstalled("Package")
        @test isinstalled("Dep")
        pkg"rm Dep"

        # Add from url at branch.
        Pkg.add(Pkg.PackageSpec(url=packages_dir_url, subdir="julia",
                                rev="master"))
        @test isinstalled("Package")
        @test !isinstalled("Dep")
        @test isinstalled(dep)
        pkg"rm Package"

        Pkg.add(Pkg.PackageSpec(url=packages_dir_url, subdir="dependencies/Dep", rev="master"))
        @test !isinstalled("Package")
        @test isinstalled("Dep")
        pkg"rm Dep"

        # Add from url at branch, REPL subdir syntax.
        pkgstr("add $(packages_dir_url):julia#master")
        @test isinstalled("Package")
        @test !isinstalled("Dep")
        @test isinstalled(dep)
        pkg"rm Package"

        pkgstr("add $(packages_dir_url):dependencies/Dep#master")
        @test !isinstalled("Package")
        @test isinstalled("Dep")
        pkg"rm Dep"

        # Develop from url.
        Pkg.develop(Pkg.PackageSpec(url=packages_dir_url, subdir="julia"))
        @test isinstalled("Package")
        @test !isinstalled("Dep")
        @test isinstalled(dep)
        pkg"rm Package"

        Pkg.develop(Pkg.PackageSpec(url=packages_dir_url, subdir="dependencies/Dep"))
        @test !isinstalled("Package")
        @test isinstalled("Dep")
        pkg"rm Dep"

        # Develop from url, REPL subdir syntax.
        pkgstr("develop $(packages_dir_url):julia")
        @test isinstalled("Package")
        @test !isinstalled("Dep")
        @test isinstalled(dep)
        pkg"rm Package"

        pkgstr("develop $(packages_dir_url):dependencies/Dep")
        @test !isinstalled("Package")
        @test isinstalled("Dep")
        pkg"rm Dep"
    end
end

end # module


text:
using Base.Test
using SparseRF
using ScikitLearnBase

srand(2)
n,m = 10^3, 5 ;
features = rand(n,m);
weights = rand(-1:1,m);
labels = features * weights;

model = fit!(DecisionTreeRegressor(maxlabels=5, pruning_purity_threshold=0.1), features, labels)
@test R2(labels, predict(model, features)) > 0.8

model = fit!(RandomForestRegressor(ntrees=10, maxlabels=5, nsubfeatures=2), features, labels)
@test R2(labels, predict(model, features)) > 0.8

srand(2)
N = 3000
X = randn(N, 10)
y = randn(N)
maxdepth = 5
model = fit!(DecisionTreeRegressor(maxdepth=maxdepth), X, y)
@test depth(model) == maxdepth

## Test that the RNG arguments work as expected
srand(2)
X = randn(100, 10)
y = randn(100)
@test fit_predict!(RandomForestRegressor(; rng=10), X, y) ==
    fit_predict!(RandomForestRegressor(; rng=10), X, y)
@test fit_predict!(RandomForestRegressor(; rng=10), X, y) !=
    fit_predict!(RandomForestRegressor(; rng=22), X, y)


text:
module QuantumCircuitOpt

import JuMP
import LinearAlgebra
import Memento

const LA = LinearAlgebra
const QCO = QuantumCircuitOpt
const kron_symbol = 'x'
const qubit_separator = '_'

# Create our module level logger (this will get precompiled)
const _LOGGER = Memento.getlogger(@__MODULE__)

# Register the module level logger at runtime so that folks can access the logger via `getlogger(QuantumCircuitOpt)`
# NOTE: If this line is not included then the precompiled `QuantumCircuitOpt._LOGGER` won't be registered at runtime.
__init__() = Memento.register(_LOGGER)

"Suppresses information and warning messages output by QuantumCircuitOpt, for fine grained control use of the Memento package"
function silence()
    Memento.info(_LOGGER, "Suppressing information and warning messages for the rest of this session.  Use the Memento package for more fine-grained control of logging.")
    Memento.setlevel!(Memento.getlogger(QuantumCircuitOpt), "error")
end

"allows the user to set the logging level without the need to add Memento"
function logger_config!(level)
    Memento.config!(Memento.getlogger("QuantumCircuitOpt"), level)
end

include("data.jl")
include("gates.jl")
include("types.jl")
include("utility.jl")
include("chull.jl")
include("relaxations.jl")
include("qc_model.jl")
include("variables.jl")
include("constraints.jl")
include("objective.jl")
include("log.jl")
include("solution.jl")

end # module


text:
########################################################################
# Server for synchronous learning
########################################################################
mutable struct Server{T1<:Int64, T2<:Float64, T3<:Vector{T2}, T4<:Vector{T1}, T5<:Matrix{T2}, T6<:Array{T2, 3}, T7<:Vector{Client}}
    Ytrain::T4                       # training label
    Ytest::T4                        # test label
    num_classes::T1                  # number of classes
    num_clients::T1                  # number of clients
    num_epoches::T1                  # number of epoches
    batch_size::T1                   # batch size
    learning_rate::T2                # learning rate
    clients::T7                      # set of clients
    b::T3                            # server model
    embeddings::T6                   # embeddings for batches
    train_embeddings::T6             # embeddings for all training data (used for final evaluation)
    test_embeddings::T6              # embeddings for all test data (used for final evaluation)
    batch::T4                        # mini-batch
    grads::T5                        # gradient information
    function Server(Ytrain::Vector{Int64}, Ytest::Vector{Int64}, config::Dict{String, Union{Int64, Float64, String}})
        num_classes = config["num_classes"]
        num_clients = config["num_clients"]
        num_epoches = config["num_epoches"]
        batch_size = config["batch_size"]
        learning_rate = config["learning_rate"]
        clients = Vector{Client}(undef, num_clients)
        b = zeros(Float64, num_classes)
        embeddings = zeros(Float64, num_clients, num_classes, batch_size)
        train_embeddings = zeros(Float64, num_clients, num_classes, length(Ytrain))
        test_embeddings = zeros(Float64, num_clients, num_classes, length(Ytest))
        batch = zeros(Int64, batch_size)
        grads = zeros(Float64, num_classes, batch_size)
        new{Int64, Float64, Vector{Float64}, Vector{Int64}, Matrix{Float64}, Array{Float64, 3}, Vector{Client}}(Ytrain, Ytest, num_classes, num_clients, num_epoches, batch_size, learning_rate, clients, b, embeddings, train_embeddings, test_embeddings, batch, grads)
    end
end

# connect with client
function connect(s::Server, c::Client)
    s.clients[c.id] = c
end

# update batch
function update_batch(s::Server, batch::Vector{Int64})
    s.batch .= batch
end

# send embeddings to server
function send_embedding(c::Client, s::Server; tag = "batch")
    if tag == "batch"
        Xbatch = c.Xtrain[:, c.batch]
        if typeof(c.W) <: Matrix{Float64}
            embedding = c.W * Xbatch
        else
            embedding, c.back = Zygote.pullback(()->c.W(Xbatch), params(c.W))
        end
    elseif tag == "training"
        if typeof(c.W) <: Matrix{Float64}
            embedding = c.W * c.Xtrain
        else
            embedding = c.W(c.Xtrain)
        end
    else
        if typeof(c.W) <: Matrix{Float64}
            embedding = c.W * c.Xtest
        else
            embedding = c.W(c.Xtest)
        end
    end
    update_embedding(s, c.id, embedding, tag=tag)
end

# update embedding
function update_embedding(s::Server, id::Int64, embedding::Matrix{Float64}; tag = "batch")
    if tag == "batch"
        s.embeddings[id,:,:] .= embedding
    elseif tag == "training"
        s.train_embeddings[id,:,:] .= embedding
    else
        s.test_embeddings[id,:,:] .= embedding
    end
end

# compute mini-batch gradient
function compute_mini_batch_gradient(s::Server)
    batch_size = s.batch_size
    num_classes = s.num_classes
    sum_embeddings = reshape( sum( s.embeddings, dims=1), num_classes, batch_size )
    loss = 0.0
    grads = zeros( num_classes, batch_size )
    # compute mini-batch gradient
    for i = 1:batch_size
        y = s.Ytrain[ s.batch[i] ]
        emb = sum_embeddings[:, i] + s.b
        pred = softmax(emb)
        loss += neg_log_loss(pred, y)
        grads[:, i] .= pred
        grads[y, i] -= 1.0
    end
    # update local gradient information
    s.grads .= grads
    # send gradient information to clients
    for c in s.clients
        update_grads(c, grads)
    end
    # return mini-batch loss
    return loss / batch_size
end

# update server model b
function update_model(s::Server)
    bgrad = sum(s.grads, dims=2) ./ s.batch_size
    s.b .-= s.learning_rate * bgrad[:]
end

########################################################################
# Server for asynchronous learning
########################################################################
mutable struct AsynServer{T1<:Int64, T2<:Float64, T3<:Vector{T2}, T4<:Vector{T1}, T5<:SharedArray{T2, 3}, T6<:Vector{AsynClient}}
    Ytrain::T4                       # training label
    Ytest::T4                        # test label
    num_classes::T1                  # number of classes
    num_clients::T1                  # number of clients
    learning_rate::T2                # learning rate
    clients::T6                      # set of clients
    b::T3                            # server model
    embeddings::T5                   # latest embeddings
    train_embeddings::T5             # embeddings for all training data (used for final evaluation)
    test_embeddings::T5              # embeddings for all test data (used for final evaluation)
    function AsynServer(Ytrain::Vector{Int64}, Ytest::Vector{Int64}, config::Dict{String, Union{Int64, Float64, String}})
        num_classes = config["num_classes"]
        num_clients = config["num_clients"]
        learning_rate = config["learning_rate"]
        clients = Vector{AsynClient}(undef, num_clients)
        b = zeros(Float64, num_classes)
        embeddings = SharedArray{Float64}(num_clients, num_classes, length(Ytrain))
        train_embeddings = SharedArray{Float64}(num_clients, num_classes, length(Ytrain))
        test_embeddings = SharedArray{Float64}(num_clients, num_classes, length(Ytest))
        new{Int64, Float64, Vector{Float64}, Vector{Int64}, SharedArray{Float64, 3}, Vector{AsynClient}}(Ytrain, Ytest, num_classes, num_clients, learning_rate, clients, b, embeddings, train_embeddings, test_embeddings)
    end
end

# connect with client
function connect(s::AsynServer, c::AsynClient)
    s.clients[c.id] = c
end

# send embeddings to server
function send_embedding(c::AsynClient, s::AsynServer; tag = "batch")
    if tag == "batch"
        num_data = length(s.Ytrain)
        if c.num_commu == 0
            batch = collect(1:num_data)
        else
            batch = sample(collect(1:num_data), c.batch_size, replace=false)
        end
        Xbatch = c.Xtrain[:, batch]
        if typeof(c.W) <: Matrix{Float64}
            embedding = c.W * Xbatch
        else
            embedding, c.back = Zygote.pullback(()->c.W(Xbatch), params(c.W))
        end
        update_embedding(s, c.id, embedding, batch)
        @printf "Client %i finish uploading embedding \n" c.id
        return batch
    elseif tag == "training"
        if typeof(c.W) <: Matrix{Float64}
            embedding = c.W * c.Xtrain
        else
            embedding = c.W(c.Xtrain)
        end
        s.train_embeddings[c.id,:,:] .= embedding
    else
        if typeof(c.W) <: Matrix{Float64}
            embedding = c.W * c.Xtest
        else
            embedding = c.W(c.Xtest)
        end
        s.test_embeddings[c.id,:,:] .= embedding
    end
end

# update embedding
function update_embedding(s::AsynServer, id::Int64, embedding::Matrix{Float64}, batch::Vector{Int64})
    s.embeddings[id,:,batch] .= embedding
end

# compute gradient
function send_gradient(s::AsynServer, id::Int64, batch::Vector{Int64})
    batch_size = length(batch)
    num_classes = s.num_classes
    sum_embeddings = reshape( sum( s.embeddings[:,:,batch], dims=1), num_classes, batch_size )
    loss = 0.0
    grads = zeros( num_classes, batch_size )
    # compute mini-batch gradient
    for i = 1:batch_size
        y = s.Ytrain[ batch[i] ]
        emb = sum_embeddings[:, i] + s.b
        pred = softmax(emb)
        loss += neg_log_loss(pred, y)
        grads[:, i] .= pred
        grads[y, i] -= 1.0
    end

    # send gradient information to client
    update_grads(s.clients[id], grads)

    # return mini-batch loss
    return loss / batch_size
end

# Compute training and test accuracy
function eval(s::Union{Server, AsynServer})
    train_size = length(s.Ytrain)
    test_size = length(s.Ytest)
    num_classes = s.num_classes
    sum_train_embeddings = reshape( sum( s.train_embeddings, dims=1), num_classes, train_size )
    sum_test_embeddings = reshape( sum( s.test_embeddings, dims=1), num_classes, test_size )
    train_loss = 0.0
    train_acc = 0.0
    test_loss = 0.0
    test_acc = 0.0
    for i = 1:train_size
        y = s.Ytrain[i]
        emb = sum_train_embeddings[:, i] + s.b
        pred = softmax(emb)
        train_loss += neg_log_loss(pred, y)
        if argmax(pred) == y
            train_acc += 1.0
        end
    end
    for i = 1:test_size
        y = s.Ytest[i]
        emb = sum_test_embeddings[:, i] + s.b
        pred = softmax(emb)
        test_loss += neg_log_loss(pred, y)
        if argmax(pred) == y
            test_acc += 1.0
        end
    end
    return train_loss/train_size, train_acc/train_size, test_loss/test_size, test_acc/test_size
end