text:
using SeisAcoustic

# the path to segy file which include 3 sourec lines 4 component OBC data
work_dir = "/home/wgao1/BP_valhall/line1";
path_sgy = joinpath(work_dir, "line_666_668.segy");

# read text file header
text_header = read_text_header(path_sgy);
print_text_header(text_header)

# read file header
file_header = read_file_header(path_sgy);

# total number of traces
tmp = filesize(path_sgy);
num_traces = convert(Int64, (tmp - 3600) / (240+sizeof(Float32) * file_header.ns))

# read all the trace header
traces_header = extract_traces_header(path_sgy; print_interval=100000);

# create quick lookup table for shot gathers
path_tab = joinpath(work_dir, "lookup_table_obc.txt");
create_lookup_table(path_tab, path_sgy);

# creat acquisition geometry from shot gather
(sx, sy, gx, gy) = create_acquisition_geometry(traces_header; num_component=4);
i = 1; figure(figsize=(8,8));
scatter(gx[i], gy[i], s=1, c="g");
scatter(sx[i], sy[i], s=1, c="r");

# read one shot gather from segy file
idx = 317; (th, shot) = read_one_shot(path_sgy, path_tab, idx);

# dimensions of data cube
dt = 0.004
num_samples   = 2000
num_receivers = 2414
num_sources   = 315

# prepare one common source line shot 317 - 631
idx_l=317; idx_u=631;
(traces_header, data) = forward(path_sgy, path_tab, idx_l, idx_u, num_samples, num_receivers);

# write this source line as segy file
path_out = joinpath(work_dir, "shot_line1.segy");
write_segy_file(path_out, text_header, file_header, traces_header, data);

# prepare one common source line shot 946-632
idx_l=632; idx_u=946;
(traces_header, data) = backward(path_sgy, path_tab, idx_l, idx_u, num_samples, num_receivers);
path_out = joinpath(work_dir, "shot_line2.segy");
write_segy_file(path_out, text_header, file_header, traces_header, data);

# read the first generated segy file
path_out1 = joinpath(work_dir, "shot_line1.segy");
(text_header1, file_header1, traces_header1, data1) = read_segy_file(path_out1);
data1 = reshape(data1, num_samples, num_receivers, num_sources);

# read the second generated segy file
path_out2 = joinpath(work_dir, "shot_line2.segy");
(text_header2, file_header2, traces_header2, data2) = read_segy_file(path_out2);
data2 = reshape(data2, num_samples, num_receivers, num_sources);

# clipping the common receiver gathers simultaneously
dc1 = copy(data1);
dc2 = copy(data2);
for i = 1 : num_receivers
    tmp1 = dc1[:,i,:]
    tmp2 = dc2[:,i,:]
    val  = quantile(abs.(vcat(vec(tmp1), vec(tmp2))), 0.90)
    dc1[:,i,:] = clamp.(tmp1, -val, val)
    dc2[:,i,:] = clamp.(tmp2, -val, val)
    println("$i")
end

# save the clipped data
path_out1 = joinpath(work_dir, "clip_line1.segy");
write_segy_file(path_out1, text_header1, file_header1, traces_header1, dc1);

path_out2 = joinpath(work_dir, "clip_line2.segy");
write_segy_file(path_out2, text_header2, file_header2, traces_header2, dc2);

# plotting
idx = 2414
figure(); imshow(dc1[:,idx,:], cmap="gray", aspect=0.2)
figure(); imshow(dc2[:,idx,:], cmap="gray", aspect=0.2)

# begin for pseudo-blending
# random time shift the second shot line
tau1 = 4.0 * rand(num_sources);
# make the time delay is the integer times of time sample interval
for i = 1 : num_sources
    tmp = round(Int64, tau1[i]/dt)
    tau1[i] = tmp * dt
end
path_rts = joinpath(work_dir, "random_time_shift1.bin")
fp = open(path_rts, "w");
write(fp, num_sources);
write(fp, tau1); close(fp);

# read the random time shift
tau2 = 4.0 * rand(num_sources);
# make the time delay is the integer times of time sample interval
for i = 1 : num_sources
    tmp = round(Int64, tau2[i]/dt)
    tau2[i] = tmp * dt
end
path_rts = joinpath(work_dir, "random_time_shift2.bin")
fp = open(path_rts, "w");
write(fp, num_sources);
write(fp, tau2); close(fp);

# read the clipped seismic data
path_clip1 = joinpath(work_dir, "clip_line1.segy");
(text_header1, file_header1, traces_header1, d_clip1) = read_segy_file(path_clip1);
d_clip1 = reshape(d_clip1, num_samples, num_receivers, num_sources);

path_clip2 = joinpath(work_dir, "clip_line2.segy");
(text_header2, file_header2, traces_header2, d_clip2) = read_segy_file(path_clip2);
d_clip2 = reshape(d_clip2, num_samples, num_receivers, num_sources);

# show the clipped seismic data
idx = 50
figure(); imshow(d_clip1[:,:,idx], cmap="gray")
figure(); imshow(d_clip1[:,:,idx], cmap="gray")

# pseudo blending, align the shots from first source line
path_tau1 = joinpath(work_dir, "random_time_shift1.bin")
fp = open(path, "r")
num_tau = rand(fp, Int64)
tau1    = zeros(Float64, num_tau)
read!(fp, tau1); close(fp)
d_blend1 = pseudo_blending(d_clip1, d_clip2, tau1);

# pseudo blending, align the shots from second source line
path_tau2 = joinpath(work_dir, "random_time_shift2.bin")
fp = open(path_tau2, "r")
num_tau = read(fp, Int64)
tau2    = zeros(Float64, num_tau)
read!(fp, tau2); close(fp)
d_blend2 = pseudo_blending(d_clip2, d_clip1, tau2);

# write the blended shot gathers into segy file
path_out = joinpath(work_dir, "blend_line1.segy");
write_segy_file(path_out, text_header1, file_header1, traces_header1, d_blend1);

path_out = joinpath(work_dir, "blend_line2.segy");
write_segy_file(path_out, text_header2, file_header2, traces_header2, d_blend2);

# write the data into binary file can be accessed for training
path_in = joinpath(work_dir, "input_train1.bin")
fp = open(path_in, "w")
write(fp, num_samples, num_receivers, num_sources);
write(fp, vec(d_blend1)); close(fp);

path_in = joinpath(work_dir, "output_train1.bin")
fp = open(path_in, "w")
write(fp, num_samples, num_receivers, num_sources);
write(fp, vec(d_clip1)); close(fp);

# for the second part of training data
path_in = joinpath(work_dir, "input_train2.bin")
fp = open(path_in, "w")
write(fp, num_samples, num_receivers, num_sources);
write(fp, vec(d_blend2)); close(fp);

path_in = joinpath(work_dir, "output_train2.bin")
fp = open(path_in, "w")
write(fp, num_samples, num_receivers, num_sources);
write(fp, vec(d_clip2)); close(fp);

# examine the figure
idx = 1111
figure(); imshow(d_clip1[:,idx,:], cmap="gray", aspect=0.2)
figure(); imshow(d_clip2[:,idx,:], cmap="gray", aspect=0.2)

# show blended common receiver gathers
idx = 100
figure(); imshow(d_clip1[:,idx,:], cmap="gray", aspect=0.2)
figure(); imshow(d_blend1[:,idx,:], cmap="gray", aspect=0.2)

figure(); imshow(d_clip2[:,idx,:], cmap="gray", aspect=0.2)
figure(); imshow(d_blend2[:,idx,:], cmap="gray", aspect=0.2)

"""
   pseudo blending function
"""
function pseudo_blending(d1::Array{Tv,3}, d2::Array{Tv,3}, tau) where {Tv<:AbstractFloat}

    # size of data cube
    (num_samples, num_receivers, num_sources) = size(d1)

    # check the size of another cucbe
    if num_samples != size(d2,1) || num_receivers != size(d2,2) || num_sources != size(d2,3) || length(tau) != num_sources
       error("size mismatch")
    end

    # copy one cube
    db = copy(d1)

    # pseudo belending of two sources
    for i = 1 : num_sources
        istart = round(Int64, tau[i]/dt) + 1

        # loop over receivers
        for j = 1 : num_receivers

            #loop over samples
            for k = istart : num_samples
                idx = k - istart + 1
                db[k,j,i] = db[k,j,i] + d2[idx,j,i]
            end
        end

        println("finished $i sources")
    end

    return db
end

"""
   agc based on rms amplitude
"""
function rms_agc(data; hwl=0.25)

    n = size(data, 2)
    dout = copy(data)
    for i = 1 : n
        dout[:,i] = agc(dout[:,i]; half_window_length=hwl)
    end
    return dout
end

"""
   t-square amplitude gaining to compensate geometrical spreading
"""
function tsquare(d; dt=0.004)

    (nt, n1) = size(d)
    dout= copy(d)

    for i = 1 : nt
        t = (i-1) * dt

        for j = 1 : n1
            dout[i,j] = dout[i,j] * t^2
        end
    end
    return dout
end

"""
   automatic gaining control applied to one trace, hwl is the half window length in time,
"""
function agc(offset, trace::Vector{Tv}; half_window_length=0.5, vel=1500, dt=0.004, drms=1.0) where {Tv<:Real}

    nt   = length(trace)
    data = copy(trace)

    dtime = offset / vel
    istart= floor(Int64, dtime / dt)
    hwl   = round(Int64, half_window_length/dt)
    wl    = 2 * hwl + 1
    ratio = 1.0
    crms  = 0.0

    if istart - hwl < 1
       istart = hwl + 1
    end

    iend = nt - hwl
    if iend <= istart
       println("window size is too big")
       return data
    end

    # compute desired rms amplitude from seismic data
    if drms == 0.0
       idx  = istart
       il   = istart - hwl
       iu   = istart + hwl
       for i = il : iu
           drms = drms + trace[i]^2
       end
       drms = sqrt(drms / wl)
       crms = drms
    end

    # boost the middle part of the data
    for idx = istart+1 : iend
        crms= sqrt((crms^2 * wl - trace[idx-hwl-1]^2 + trace[idx+hwl]^2)/wl)
        ratio = drms / crms
        data[idx] = data[idx] * ratio
    end

    # boost the last hwl samples
    for idx = iend+1 : nt
        data[idx] = data[idx] * ratio
    end

    return data

end

# start from the first sample
function agc(trace::Vector{Tv}; half_window_length=0.5, dt=0.004, drms=1.0, delta=1e-16) where {Tv<:Real}

    nt   = length(trace)
    data = copy(trace)

    hwl = round(Int64, half_window_length/dt)
    wl  = 2 * hwl + 1

    if wl > nt
       println("window size is bigger than the length of input data")
       return data
    end
    # the rest of the code can handle the case when nt = wl

    # the start and end index of window center
    istart= hwl + 1
    iend  = nt  - hwl

    # set crms as the rms amplitude of the first window
    crms  = 0.0
    il = istart - hwl
    iu = istart + hwl
    for i = il : iu
        crms = crms + trace[i]^2
    end
    crms = sqrt(crms / wl)

    # the ratio between drms and crms
    ratio = drms / (crms+delta)

    # apply to the samples in the first half of window
    for i = 1 : istart
        data[i] = data[i] * ratio
    end

    # loop over from the second to the last window
    for idx = istart+1 : iend
        crms= sqrt((crms^2 * wl - trace[idx-hwl-1]^2 + trace[idx+hwl]^2)/wl)
        ratio = drms / (crms+delta)
        data[idx] = data[idx] * ratio
    end

    # boost the last hwl samples
    for idx = iend+1 : nt
        data[idx] = data[idx] * ratio
    end

    return data

end

function forward(path_sgy, path_tab, idx_l, idx_u, num_samples, num_receivers)

    # number of shot in this source line
    num_shots = idx_u - idx_l + 1

    # allocate memory for all the traces
    data = zeros(Float32, num_samples, num_receivers, num_shots);

    # trace header
    traces_header = Vector{TraceHeader}(undef, 0);

    ncount = 1
    for idx = idx_l : idx_u

        # read one shot gather
        (th, d) = read_one_shot(path_sgy, path_tab, idx)

        # take out the fourth component
        th = th[4:4:end]
        d  = d[:,4:4:end]

        # put to the main one
        traces_header = vcat(traces_header, th)
        data[:,:,ncount] = d

        println("$ncount")
        #jump to next shot
        ncount = ncount + 1
    end

    return traces_header, data
end

function backward(path_sgy, path_tab, idx_l, idx_u, num_samples, num_receivers)

    # number of shot in this source line
    num_shots = idx_u - idx_l + 1

    # allocate memory for all the traces
    data = zeros(Float32, num_samples, num_receivers, num_shots);

    # trace header
    traces_header = Vector{TraceHeader}(undef, 0);

    ncount = 1
    for idx = idx_u : -1 : idx_l

        # read one shot gather
        (th, d) = read_one_shot(path_sgy, path_tab, idx)

        # take out the fourth component
        th = th[4:4:end]
        d  = d[:,4:4:end]

        # put to the main one
        traces_header = vcat(traces_header, th)
        data[:,:,ncount] = d

        println("$ncount")
        #jump to next shot
        ncount = ncount + 1
    end

    return traces_header, data
end

# plot commands via seismic unix
tmp = "/home/wgao1/BP_valhall/line2/figure/p1.bin"; tmpfid = open(tmp, "w"); write(tmpfid, vec(data1[:,:,10])); close(tmpfid);
tmp = "/home/wgao1/BP_valhall/line2/figure/p2.bin"; tmpfid = open(tmp, "w"); write(tmpfid, vec(data1[:,:,110])); close(tmpfid);
tmp = "/home/wgao1/BP_valhall/line2/figure/p3.bin"; tmpfid = open(tmp, "w"); write(tmpfid, vec(data1[:,:,210])); close(tmpfid);
tmp = "/home/wgao1/BP_valhall/line2/figure/p4.bin"; tmpfid = open(tmp, "w"); write(tmpfid, vec(data1[:,:,310])); close(tmpfid);

psimage < /home/wgao1/BP_valhall/line2/figure/p1.bin xbox=0.5 ybox=0.5 height=8 width=12 labelsize=15 n1=2000 label1="Time" n2=2414 d2num=200 f2num=200 label2="Trace" perc=90 > /home/wgao1/BP_valhall/line2/figure/p1.eps
psimage < /home/wgao1/BP_valhall/line2/figure/p2.bin xbox=0.5 ybox=0.5 height=8 width=12 labelsize=15 n1=2000 label1="Time" n2=2414 d2num=200 f2num=200 label2="Trace" perc=90 > /home/wgao1/BP_valhall/line2/figure/p2.eps
psimage < /home/wgao1/BP_valhall/line2/figure/p3.bin xbox=0.5 ybox=0.5 height=8 width=12 labelsize=15 n1=2000 label1="Time" n2=2414 d2num=200 f2num=200 label2="Trace" perc=90 > /home/wgao1/BP_valhall/line2/figure/p3.eps
psimage < /home/wgao1/BP_valhall/line2/figure/p4.bin xbox=0.5 ybox=0.5 height=8 width=12 labelsize=15 n1=2000 label1="Time" n2=2414 d2num=200 f2num=200 label2="Trace" perc=90 > /home/wgao1/BP_valhall/line2/figure/p4.eps

# plot common receiver gather
tmp = "/home/wgao1/BP_valhall/line2/figure/r1.bin"; tmpfid = open(tmp, "w"); write(tmpfid, vec(data1[:,200,:])); close(tmpfid);
tmp = "/home/wgao1/BP_valhall/line2/figure/r2.bin"; tmpfid = open(tmp, "w"); write(tmpfid, vec(data1[:,400,:])); close(tmpfid);
tmp = "/home/wgao1/BP_valhall/line2/figure/r3.bin"; tmpfid = open(tmp, "w"); write(tmpfid, vec(data1[:,600,:])); close(tmpfid);
tmp = "/home/wgao1/BP_valhall/line2/figure/r4.bin"; tmpfid = open(tmp, "w"); write(tmpfid, vec(data1[:,800,:])); close(tmpfid);
tmp = "/home/wgao1/BP_valhall/line2/figure/r5.bin"; tmpfid = open(tmp, "w"); write(tmpfid, vec(data1[:,1000,:])); close(tmpfid);
tmp = "/home/wgao1/BP_valhall/line2/figure/r6.bin"; tmpfid = open(tmp, "w"); write(tmpfid, vec(data1[:,1200,:])); close(tmpfid);
tmp = "/home/wgao1/BP_valhall/line2/figure/r7.bin"; tmpfid = open(tmp, "w"); write(tmpfid, vec(data1[:,1400,:])); close(tmpfid);
tmp = "/home/wgao1/BP_valhall/line2/figure/r8.bin"; tmpfid = open(tmp, "w"); write(tmpfid, vec(data1[:,1600,:])); close(tmpfid);
tmp = "/home/wgao1/BP_valhall/line2/figure/r9.bin"; tmpfid = open(tmp, "w"); write(tmpfid, vec(data1[:,1800,:])); close(tmpfid);
tmp = "/home/wgao1/BP_valhall/line2/figure/r10.bin"; tmpfid = open(tmp, "w"); write(tmpfid, vec(data1[:,2000,:])); close(tmpfid);
tmp = "/home/wgao1/BP_valhall/line2/figure/r11.bin"; tmpfid = open(tmp, "w"); write(tmpfid, vec(data1[:,2200,:])); close(tmpfid);
tmp = "/home/wgao1/BP_valhall/line2/figure/r12.bin"; tmpfid = open(tmp, "w"); write(tmpfid, vec(data1[:,2400,:])); close(tmpfid);

psimage < /home/wgao1/BP_valhall/line2/figure/r1.bin xbox=0.5 ybox=0.5 height=8 width=12 labelsize=15 n1=2000 label1="Time" n2=315 d2num=50 f2num=50 label2="Trace" perc=90 > /home/wgao1/BP_valhall/line2/figure/r1.eps
psimage < /home/wgao1/BP_valhall/line2/figure/r2.bin xbox=0.5 ybox=0.5 height=8 width=12 labelsize=15 n1=2000 label1="Time" n2=315 d2num=50 f2num=50 label2="Trace" perc=90 > /home/wgao1/BP_valhall/line2/figure/r2.eps
psimage < /home/wgao1/BP_valhall/line2/figure/r3.bin xbox=0.5 ybox=0.5 height=8 width=12 labelsize=15 n1=2000 label1="Time" n2=315 d2num=50 f2num=50 label2="Trace" perc=90 > /home/wgao1/BP_valhall/line2/figure/r3.eps
psimage < /home/wgao1/BP_valhall/line2/figure/r4.bin xbox=0.5 ybox=0.5 height=8 width=12 labelsize=15 n1=2000 label1="Time" n2=315 d2num=50 f2num=50 label2="Trace" perc=90 > /home/wgao1/BP_valhall/line2/figure/r4.eps
psimage < /home/wgao1/BP_valhall/line2/figure/r5.bin xbox=0.5 ybox=0.5 height=8 width=12 labelsize=15 n1=2000 label1="Time" n2=315 d2num=50 f2num=50 label2="Trace" perc=90 > /home/wgao1/BP_valhall/line2/figure/r5.eps
psimage < /home/wgao1/BP_valhall/line2/figure/r6.bin xbox=0.5 ybox=0.5 height=8 width=12 labelsize=15 n1=2000 label1="Time" n2=315 d2num=50 f2num=50 label2="Trace" perc=90 > /home/wgao1/BP_valhall/line2/figure/r6.eps
psimage < /home/wgao1/BP_valhall/line2/figure/r7.bin xbox=0.5 ybox=0.5 height=8 width=12 labelsize=15 n1=2000 label1="Time" n2=315 d2num=50 f2num=50 label2="Trace" perc=90 > /home/wgao1/BP_valhall/line2/figure/r7.eps
psimage < /home/wgao1/BP_valhall/line2/figure/r8.bin xbox=0.5 ybox=0.5 height=8 width=12 labelsize=15 n1=2000 label1="Time" n2=315 d2num=50 f2num=50 label2="Trace" perc=90 > /home/wgao1/BP_valhall/line2/figure/r8.eps
psimage < /home/wgao1/BP_valhall/line2/figure/r9.bin xbox=0.5 ybox=0.5 height=8 width=12 labelsize=15 n1=2000 label1="Time" n2=315 d2num=50 f2num=50 label2="Trace" perc=90 > /home/wgao1/BP_valhall/line2/figure/r9.eps
psimage < /home/wgao1/BP_valhall/line2/figure/r10.bin xbox=0.5 ybox=0.5 height=8 width=12 labelsize=15 n1=2000 label1="Time" n2=315 d2num=50 f2num=50 label2="Trace" perc=90 > /home/wgao1/BP_valhall/line2/figure/r10.eps
psimage < /home/wgao1/BP_valhall/line2/figure/r11.bin xbox=0.5 ybox=0.5 height=8 width=12 labelsize=15 n1=2000 label1="Time" n2=315 d2num=50 f2num=50 label2="Trace" perc=90 > /home/wgao1/BP_valhall/line2/figure/r11.eps
psimage < /home/wgao1/BP_valhall/line2/figure/r12.bin xbox=0.5 ybox=0.5 height=8 width=12 labelsize=15 n1=2000 label1="Time" n2=315 d2num=50 f2num=50 label2="Trace" perc=90 > /home/wgao1/BP_valhall/line2/figure/r12.eps

# plot the second shot line
tmp = "/home/wgao1/BP_valhall/line2/figure/r21.bin"; tmpfid = open(tmp, "w"); write(tmpfid, vec(data2[:,200,:])); close(tmpfid);
tmp = "/home/wgao1/BP_valhall/line2/figure/r22.bin"; tmpfid = open(tmp, "w"); write(tmpfid, vec(data2[:,400,:])); close(tmpfid);
tmp = "/home/wgao1/BP_valhall/line2/figure/r23.bin"; tmpfid = open(tmp, "w"); write(tmpfid, vec(data2[:,600,:])); close(tmpfid);
tmp = "/home/wgao1/BP_valhall/line2/figure/r24.bin"; tmpfid = open(tmp, "w"); write(tmpfid, vec(data2[:,800,:])); close(tmpfid);
tmp = "/home/wgao1/BP_valhall/line2/figure/r25.bin"; tmpfid = open(tmp, "w"); write(tmpfid, vec(data2[:,1000,:])); close(tmpfid);
tmp = "/home/wgao1/BP_valhall/line2/figure/r26.bin"; tmpfid = open(tmp, "w"); write(tmpfid, vec(data2[:,1200,:])); close(tmpfid);
tmp = "/home/wgao1/BP_valhall/line2/figure/r27.bin"; tmpfid = open(tmp, "w"); write(tmpfid, vec(data2[:,1400,:])); close(tmpfid);
tmp = "/home/wgao1/BP_valhall/line2/figure/r28.bin"; tmpfid = open(tmp, "w"); write(tmpfid, vec(data2[:,1600,:])); close(tmpfid);
tmp = "/home/wgao1/BP_valhall/line2/figure/r29.bin"; tmpfid = open(tmp, "w"); write(tmpfid, vec(data2[:,1800,:])); close(tmpfid);
tmp = "/home/wgao1/BP_valhall/line2/figure/r210.bin"; tmpfid = open(tmp, "w"); write(tmpfid, vec(data2[:,2000,:])); close(tmpfid);
tmp = "/home/wgao1/BP_valhall/line2/figure/r211.bin"; tmpfid = open(tmp, "w"); write(tmpfid, vec(data2[:,2200,:])); close(tmpfid);
tmp = "/home/wgao1/BP_valhall/line2/figure/r212.bin"; tmpfid = open(tmp, "w"); write(tmpfid, vec(data2[:,2400,:])); close(tmpfid);

psimage < /home/wgao1/BP_valhall/line2/figure/r21.bin xbox=0.5 ybox=0.5 height=8 width=12 labelsize=15 n1=2000 label1="Time" n2=315 d2num=50 f2num=50 label2="Trace" perc=90 > /home/wgao1/BP_valhall/line2/figure/r21.eps
psimage < /home/wgao1/BP_valhall/line2/figure/r22.bin xbox=0.5 ybox=0.5 height=8 width=12 labelsize=15 n1=2000 label1="Time" n2=315 d2num=50 f2num=50 label2="Trace" perc=90 > /home/wgao1/BP_valhall/line2/figure/r22.eps
psimage < /home/wgao1/BP_valhall/line2/figure/r23.bin xbox=0.5 ybox=0.5 height=8 width=12 labelsize=15 n1=2000 label1="Time" n2=315 d2num=50 f2num=50 label2="Trace" perc=90 > /home/wgao1/BP_valhall/line2/figure/r23.eps
psimage < /home/wgao1/BP_valhall/line2/figure/r24.bin xbox=0.5 ybox=0.5 height=8 width=12 labelsize=15 n1=2000 label1="Time" n2=315 d2num=50 f2num=50 label2="Trace" perc=90 > /home/wgao1/BP_valhall/line2/figure/r24.eps
psimage < /home/wgao1/BP_valhall/line2/figure/r25.bin xbox=0.5 ybox=0.5 height=8 width=12 labelsize=15 n1=2000 label1="Time" n2=315 d2num=50 f2num=50 label2="Trace" perc=90 > /home/wgao1/BP_valhall/line2/figure/r25.eps
psimage < /home/wgao1/BP_valhall/line2/figure/r26.bin xbox=0.5 ybox=0.5 height=8 width=12 labelsize=15 n1=2000 label1="Time" n2=315 d2num=50 f2num=50 label2="Trace" perc=90 > /home/wgao1/BP_valhall/line2/figure/r26.eps
psimage < /home/wgao1/BP_valhall/line2/figure/r27.bin xbox=0.5 ybox=0.5 height=8 width=12 labelsize=15 n1=2000 label1="Time" n2=315 d2num=50 f2num=50 label2="Trace" perc=90 > /home/wgao1/BP_valhall/line2/figure/r27.eps
psimage < /home/wgao1/BP_valhall/line2/figure/r28.bin xbox=0.5 ybox=0.5 height=8 width=12 labelsize=15 n1=2000 label1="Time" n2=315 d2num=50 f2num=50 label2="Trace" perc=90 > /home/wgao1/BP_valhall/line2/figure/r28.eps
psimage < /home/wgao1/BP_valhall/line2/figure/r29.bin xbox=0.5 ybox=0.5 height=8 width=12 labelsize=15 n1=2000 label1="Time" n2=315 d2num=50 f2num=50 label2="Trace" perc=90 > /home/wgao1/BP_valhall/line2/figure/r29.eps
psimage < /home/wgao1/BP_valhall/line2/figure/r210.bin xbox=0.5 ybox=0.5 height=8 width=12 labelsize=15 n1=2000 label1="Time" n2=315 d2num=50 f2num=50 label2="Trace" perc=90 > /home/wgao1/BP_valhall/line2/figure/r210.eps
psimage < /home/wgao1/BP_valhall/line2/figure/r211.bin xbox=0.5 ybox=0.5 height=8 width=12 labelsize=15 n1=2000 label1="Time" n2=315 d2num=50 f2num=50 label2="Trace" perc=90 > /home/wgao1/BP_valhall/line2/figure/r211.eps
psimage < /home/wgao1/BP_valhall/line2/figure/r212.bin xbox=0.5 ybox=0.5 height=8 width=12 labelsize=15 n1=2000 label1="Time" n2=315 d2num=50 f2num=50 label2="Trace" perc=90 > /home/wgao1/BP_valhall/line2/figure/r212.eps

# download the figure to local machine
scp -r wgao1@saig-ml.physics.ualberta.ca:/home/wgao1/BP_valhall/line2/figure /Users/wenlei/Desktop

# 1. gaining control
# 2. time shift in frequency domain
"""
Shifting the shots according to the fire time, double the time samples padding at head
d: 3D Array, first dimension is time, second dimension is channel, third dimension is shot
"""
function DataShiftPadHead(d::Array{Tv,3}, shift::Vector{Tv}, dt::Tv) where Tv<:Float64
    (nt, nr, ns) = size(d)
    # double the time samples
    df = vcat(zeros(Float64, nt, nr, ns), d)
    df = fft(df, 1)
    wn = collect(0:nt)*2*pi/(2*nt*dt)
    ps = zeros(Complex{Float64}, nt+1)
    for i2 = 1 : ns   # loop over shots
        for iw = 1 : nt+1
            ps[iw] = exp(-im * wn[iw] * shift[i2])
        end
        for i1 = 1 : nr
            for iw = 1 : nt+1
                df[iw,i1,i2] = df[iw,i1,i2] * ps[iw]
            end
            for iw = 2 : nt
                idx = 2*nt - iw + 2
                df[idx,i1,i2] = conj(df[iw,i1,i2])
            end
        end
    end
    ds = real(ifft(df,1))[nt+1:end,:,:]
    return ds
end

"""
Shifting the shots according to the fire time, padding zeros at tail to make the number of
time samples equal to nextpow2.
d: 3D Array, first dimension is time, second dimension is channel, third dimension is shot
"""
function DataShiftPadTail(d::Array{Tv,3}, shift::Vector{Tv}, dt::Tv) where Tv<:Float64
    (nt, nr, ns) = size(d)
    nf = nextpow2(nt); hf = floor(Int64, nf/2);
    df = vcat(d, zeros(Float64, nf-nt, nr, ns)); df = fft(df, 1);
    wn = collect(0:hf)*2*pi/(nf*dt)
    ps = zeros(Complex{Float64}, nf+1)
    for i2 = 1 : ns   # loop over shots
        for iw = 1 : hf+1
            ps[iw] = exp(-im * wn[iw] * shift[i2])
        end
        for i1 = 1 : nr
            for iw = 1 : hf+1
                df[iw,i1,i2] = df[iw,i1,i2] * ps[iw]
            end
            for iw = 2 : hf
                idx = nf - iw + 2
                df[idx,i1,i2] = conj(df[iw,i1,i2])
            end
        end
    end
    ds = real(ifft(df,1))[1:nt,:,:]
    return ds
end

"""
Shifting the shots according to the fire time, no padding is applied
d: 3D Array, first dimension is time, second dimension is channel, third dimension is shot
"""
function DataShiftNoPad(d::Array{Tv,3}, shift::Vector{Tv}, dt::Tv) where Tv<:Float64
    (nt, nr, ns) = size(d)
    nf = floor(Int64, nt/2)  # assume the number of samples is even
    df = fft(d, 1)
    wn = collect(0:nf)*2*pi/(nt*dt)
    ps = zeros(Complex{Float64}, nf+1)
    for i2 = 1 : ns   # loop over shots
        for iw = 1 : nf+1
            ps[iw] = exp(-im * wn[iw] * shift[i2])
        end
        for i1 = 1 : nr
            for iw = 1 : nf+1
                df[iw,i1,i2] = df[iw,i1,i2] * ps[iw]
            end
            for iw = 2 : nf
                idx = nt - iw + 2
                df[idx,i1,i2] = conj(df[iw,i1,i2])
            end
        end
    end
    ds = real(ifft(df,1))
    return ds
end


text:
include("../../utility/sample.jl")

function increaseOrder(para, state, groups, rng = Random.GLOBAL_RNG)
    # random select a new group
    new, old = rand(rng, groups), state.group
    (new.order != old.order + 1) && return

    β, Kf = para.β, para.Kf
    δKnum = new.innerK[2] - old.innerK[2]
    δTnum = new.innerT[2] - old.innerT[2]
    propK, propT = 1.0, 1.0

    for ik in old.innerK[2]:new.innerK[2]
        propK = newK!(state.K[ik], Kf, Kf / 2.0, rng)
    end
    for it in old.innerT[2]:new.innerT[2]
        prop, state.T[new.innerT[2]] = newTau(β, rng)
        propT *= prop
    end

    state.propose[:increaseOrder] += 1
    newAbsWeight = abs(evaluate(para, state))
    R = propT * paraK * newAbsWeight * new.reweight / state.absWeight / old.reweight
    if rand(rng) < R
        state.accept[:increaseOrder] += 1
        state.group = new
        state.absWeight = newAbsWeight
    end
end

function decreaseOrder(para, state, )

ext: jl

max_stars_repo_name: kunyuan/DiagramCalculator

max_issues_repo_name: kunyuan/DiagramCalculator

max_forks_repo_name: kunyuan/DiagramCalculator


text:
function find_direction(B, ∇f, scheme::QuasiNewton{<:Direct})
   return -(B\∇f)
end
function find_direction(A, ∇f, scheme::QuasiNewton{<:Inverse})
   return -A*∇f
end

function find_direction!(d, B, ∇f, scheme::QuasiNewton{<:Direct})
   d .= -(B\∇f)
   d
end
function find_direction!(d, A, ∇f, scheme::QuasiNewton{<:Inverse})
   rmul!(mul!(d, A, ∇f), -1)
   d
end


text:
# This file was generated, do not modify it. # hide
confmat(mode.(ŷ), y[validation])

max_stars_repo_name: alan-turing-institute/MLJTutorials

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
push!(LOAD_PATH,"../src/")
using Documenter, AndersonMoore

makedocs(modules=[AndersonMoore],
        doctest=true)

deploydocs(deps   = Deps.pip("mkdocs", "python-markdown-math"),
    repo = "github.com/es335mathwiz/AndersonMoore.jl.git",
    julia  = "0.6.0",

    osname = "linux")


text:
"""
Shorthand for `isnothing(optional) ? fallback : optional`
"""
@inline ifnothing(optional, fallback) = isnothing(optional) ? fallback : optional

function round_to_IRect2D(r::Rect{2})
    newori = round.(Int, minimum(r))
    othercorner = round.(Int, maximum(r))
    newwidth = othercorner .- newori
    Rect{2, Int}(newori, newwidth)
end

function sceneareanode!(finalbbox, limits, aspect)

    scenearea = Node(IRect(0, 0, 100, 100))

    onany(finalbbox, limits, aspect) do bbox, limits, aspect

        w = width(bbox)
        h = height(bbox)
        # as = mw / mh
        as = w / h
        mw, mh = w, h

        if aspect isa AxisAspect
            aspect = aspect.aspect
        elseif aspect isa DataAspect
            aspect = limits.widths[1] / limits.widths[2]
        end

        if !isnothing(aspect)
            if as >= aspect
                # too wide
                mw *= aspect / as
            else
                # too high
                mh *= as / aspect
            end
        end

        restw = w - mw
        resth = h - mh

        # l = left(bbox) + alignment[1] * restw
        # b = bottom(bbox) + alignment[2] * resth
        l = left(bbox) + 0.5f0 * restw
        b = bottom(bbox) + 0.5f0 * resth

        newbbox = BBox(l, l + mw, b, b + mh)

        # only update scene if pixel positions change
        new_scenearea = round_to_IRect2D(newbbox)
        if new_scenearea != scenearea[]
            scenearea[] = new_scenearea
        end
    end

    scenearea
end

function roundedrectvertices(rect, cornerradius, cornersegments)
    cr = cornerradius
    csegs = cornersegments

    cr = min(width(rect) / 2, height(rect) / 2, cr)

    # inner corners
    ictl = topleft(rect) .+ Point2(cr, -cr)
    ictr = topright(rect) .+ Point2(-cr, -cr)
    icbl = bottomleft(rect) .+ Point2(cr, cr)
    icbr = bottomright(rect) .+ Point2(-cr, cr)

    # check if corners touch so we can remove one vertex that is doubled
    wtouching = width(rect) / 2 == cr
    htouching = height(rect) / 2 == cr

    cstr = if wtouching
        anglepoint.(Ref(ictr), LinRange(0, pi/2, csegs), cr)
    else
        anglepoint.(Ref(ictr), LinRange(0, pi/2, csegs)[1:end-1], cr)
    end
    cstl = if htouching
        anglepoint.(Ref(ictl), LinRange(pi/2, pi, csegs), cr)
    else
        anglepoint.(Ref(ictl), LinRange(pi/2, pi, csegs)[1:end-1], cr)
    end
    csbl = if wtouching
        anglepoint.(Ref(icbl), LinRange(pi, 3pi/2, csegs), cr)
    else
        anglepoint.(Ref(icbl), LinRange(pi, 3pi/2, csegs)[1:end-1], cr)
    end
    csbr = if htouching
        anglepoint.(Ref(icbr), LinRange(3pi/2, 2pi, csegs), cr)
    else
        anglepoint.(Ref(icbr), LinRange(3pi/2, 2pi, csegs)[1:end-1], cr)
    end
    arr = [cstr; cstl; csbl; csbr]
end

"""
    tightlimits!(la::Axis)

Sets the autolimit margins to zero on all sides.
"""
function tightlimits!(la::Axis)
    la.xautolimitmargin = (0, 0)
    la.yautolimitmargin = (0, 0)
    autolimits!(la)
end

"""
    tightlimits!(la::Axis, sides::Union{Left, Right, Bottom, Top}...)

Sets the autolimit margins to zero on all given sides.

Example:

```
tightlimits!(laxis, Bottom())
```
"""
function tightlimits!(la::Axis, sides::Union{Left, Right, Bottom, Top}...)
    for s in sides
        tightlimits!(la, s)
    end
end

function tightlimits!(la::Axis, ::Left)
    la.xautolimitmargin = Base.setindex(la.xautolimitmargin[], 0.0, 1)
    autolimits!(la)
end

function tightlimits!(la::Axis, ::Right)
    la.xautolimitmargin = Base.setindex(la.xautolimitmargin[], 0.0, 2)
    autolimits!(la)
end

function tightlimits!(la::Axis, ::Bottom)
    la.yautolimitmargin = Base.setindex(la.yautolimitmargin[], 0.0, 1)
    autolimits!(la)
end

function tightlimits!(la::Axis, ::Top)
    la.yautolimitmargin = Base.setindex(la.yautolimitmargin[], 0.0, 2)
    autolimits!(la)
end

"""
    layoutscene(padding = 30; kwargs...)

Create a `Scene` in `campixel!` mode and a `GridLayout` aligned to the scene's pixel area with `alignmode = Outside(padding)`.
"""
function layoutscene(padding = 30; kwargs...)
    scene = Scene(; camera = campixel!, kwargs...)
    gl = GridLayout(scene, alignmode = Outside(padding))
    scene, gl
end

"""
    layoutscene(nrows::Int, ncols::Int, padding = 30; kwargs...)

Create a `Scene` in `campixel!` mode and a `GridLayout` aligned to the scene's pixel area with size `nrows` x `ncols` and `alignmode = Outside(padding)`.
"""
function layoutscene(nrows::Int, ncols::Int, padding = 30; kwargs...)
    scene = Scene(; camera = campixel!, kwargs...)
    gl = GridLayout(scene, nrows, ncols, alignmode = Outside(padding))
    scene, gl
end

GridLayoutBase.GridLayout(scene::Scene, args...; kwargs...) = GridLayout(args...; bbox = lift(x -> FRect2D(x), pixelarea(scene)), kwargs...)

function axislines!(scene, rect, spinewidth, topspinevisible, rightspinevisible,
    leftspinevisible, bottomspinevisible, topspinecolor, leftspinecolor,
    rightspinecolor, bottomspinecolor)

    bottomline = lift(rect, spinewidth) do r, sw
        y = bottom(r)
        p1 = Point2(left(r) - 0.5sw, y)
        p2 = Point2(right(r) + 0.5sw, y)
        [p1, p2]
    end

    leftline = lift(rect, spinewidth) do r, sw
        x = left(r)
        p1 = Point2(x, bottom(r) - 0.5sw)
        p2 = Point2(x, top(r) + 0.5sw)
        [p1, p2]
    end

    topline = lift(rect, spinewidth) do r, sw
        y = top(r)
        p1 = Point2(left(r) - 0.5sw, y)
        p2 = Point2(right(r) + 0.5sw, y)
        [p1, p2]
    end

    rightline = lift(rect, spinewidth) do r, sw
        x = right(r)
        p1 = Point2(x, bottom(r) - 0.5sw)
        p2 = Point2(x, top(r) + 0.5sw)
        [p1, p2]
    end

    (lines!(scene, bottomline, linewidth = spinewidth, show_axis = false,
        visible = bottomspinevisible, color = bottomspinecolor),
    lines!(scene, leftline, linewidth = spinewidth, show_axis = false,
        visible = leftspinevisible, color = leftspinecolor),
    lines!(scene, rightline, linewidth = spinewidth, show_axis = false,
        visible = rightspinevisible, color = rightspinecolor),
    lines!(scene, topline, linewidth = spinewidth, show_axis = false,
        visible = topspinevisible, color = topspinecolor))
end

function interleave_vectors(vec1::Vector{T}, vec2::Vector{T}) where T
    n = length(vec1)
    @assert n == length(vec2)

    vec = Vector{T}(undef, 2 * n)
    @inbounds for i in 1:n
        k = 2(i - 1)
        vec[k + 1] = vec1[i]
        vec[k + 2] = vec2[i]
    end
    vec
end

"""
Take a sequence of variable definitions with docstrings above each and turn
them into Attributes, a Dict of varname => docstring pairs and a Dict of
varname => default_value pairs.

# Example

    attrs, docdict, defaultdict = @documented_attributes begin
        "The width."
        width = 10
        "The height."
        height = 20 + x
    end

    attrs == Attributes(
        width = 10,
        height = 20
    )

    docdict == Dict(
        width => "The width.",
        height => "The height."
    )

    defaultdict == Dict(
        width => "10",
        height => "20 + x"
    )
"""
macro documented_attributes(exp)
    if exp.head != :block
        error("Not a block")
    end

    expressions = filter(x -> !(x isa LineNumberNode), exp.args)

    vars_and_exps = map(expressions) do e
        if e.head == :macrocall && e.args[1] == GlobalRef(Core, Symbol("@doc"))
            varname = e.args[4].args[1]
            var_exp = e.args[4].args[2]
            str_exp = e.args[3]
        elseif e.head == Symbol("=")
            varname = e.args[1]
            var_exp = e.args[2]
            str_exp = "no description"
        else
            error("Neither docstringed variable nor normal variable: $e")
        end
        varname, var_exp, str_exp
    end

    # make a dictionary of :variable_name => docstring_expression
    exp_docdict = Expr(:call, :Dict,
        (Expr(:call, Symbol("=>"), QuoteNode(name), strexp)
            for (name, _, strexp) in vars_and_exps)...)

    # make a dictionary of :variable_name => docstring_expression
    defaults_dict = Expr(:call, :Dict,
        (Expr(:call, Symbol("=>"), QuoteNode(name), exp isa String ? "\"$exp\"" : string(exp))
            for (name, exp, _) in vars_and_exps)...)

    # make an Attributes instance with of variable_name = variable_expression
    exp_attrs = Expr(:call, :Attributes,
        (Expr(:kw, name, exp)
            for (name, exp, _) in vars_and_exps)...)

    esc(quote
        ($exp_attrs, $exp_docdict, $defaults_dict)
    end)
end

"""
Turn a combination of docdict and defaultdict from `@documented_attributes`
into a string to insert into a docstring.
"""
function docvarstring(docdict, defaultdict)
    buffer = IOBuffer()
    maxwidth = maximum(length ∘ string, keys(docdict))
    for (var, doc) in sort(collect(pairs(docdict)))
        print(buffer, "`$var`\\\nDefault: `$(defaultdict[var])`\\\n$doc\n\n")
    end
    String(take!(buffer))
end

function subtheme(scene, key::Symbol)
    sub = haskey(theme(scene), key) ? theme(scene, key) : Attributes()
    if !(sub isa Attributes)
        error("Subtheme is not of type Attributes but is $sub")
    end
    sub
end

"""
    xaxis_top!(la::Axis)

Move the x-axis to the top, while correctly aligning the tick labels at the bottom.
"""
function xaxis_top!(la::Axis)
    la.xaxisposition = :top
    la.xticklabelalign = (la.xticklabelalign[][1], :bottom)
    nothing
end

"""
    xaxis_bottom!(la::Axis)

Move the x-axis to the bottom, while correctly aligning the tick labels at the top.
"""
function xaxis_bottom!(la::Axis)
    la.xaxisposition = :bottom
    la.xticklabelalign = (la.xticklabelalign[][1], :top)
    nothing
end

"""
    yaxis_left!(la::Axis)

Move the y-axis to the left, while correctly aligning the tick labels at the right.
"""
function yaxis_left!(la::Axis)
    la.yaxisposition = :left
    la.yticklabelalign = (:right, la.yticklabelalign[][2])
    nothing
end

"""
    yaxis_right!(la::Axis)

Move the y-axis to the right, while correctly aligning the tick labels at the left.
"""
function yaxis_right!(la::Axis)
    la.yaxisposition = :right
    la.yticklabelalign = (:left, la.yticklabelalign[][2])
    nothing
end

"""
    labelslider!(scene, label, range; format = string, sliderkw = Dict(), labelkw = Dict(), valuekw = Dict(), layoutkw...)

Construct a horizontal GridLayout with a label, a slider and a value label in `scene`.

Returns a `NamedTuple`:

`(slider = slider, label = label, valuelabel = valuelabel, layout = layout)`

Specify a format function for the value label with the `format` keyword.
The slider is forwarded the keywords from `sliderkw`.
The label is forwarded the keywords from `labelkw`.
The value label is forwarded the keywords from `valuekw`.
All other keywords are forwarded to the `GridLayout`.

Example:

```
ls = labelslider!(scene, "Voltage:", 0:10; format = x -> "\$(x)V")
layout[1, 1] = ls.layout
```
"""
function labelslider!(scene, label, range; format = string,
        sliderkw = Dict(), labelkw = Dict(), valuekw = Dict(), layoutkw...)
    slider = Slider(scene; range = range, sliderkw...)
    label = Label(scene, label; labelkw...)
    valuelabel = Label(scene, lift(format, slider.value); valuekw...)
    layout = hbox!(label, slider, valuelabel; layoutkw...)
    (slider = slider, label = label, valuelabel = valuelabel, layout = layout)
end

"""
    labelslidergrid!(scene, labels, ranges; formats = [string],
        sliderkw = Dict(), labelkw = Dict(), valuekw = Dict(), layoutkw...)

Construct a GridLayout with a column of label, a column of sliders and a column of value labels in `scene`.
The argument values are broadcast, so you can use scalars if you want to keep labels, ranges or formats constant across rows.

Returns a `NamedTuple`:

`(sliders = sliders, labels = labels, valuelabels = valuelabels, layout = layout)`

Specify format functions for the value labels with the `formats` keyword.
The sliders are forwarded the keywords from `sliderkw`.
The labels are forwarded the keywords from `labelkw`.
The value labels are forwarded the keywords from `valuekw`.
All other keywords are forwarded to the `GridLayout`.

Example:

```
ls = labelslidergrid!(scene, ["Voltage", "Ampere"], Ref(0:0.1:100); format = x -> "\$(x)V")
layout[1, 1] = ls.layout
```
"""
function labelslidergrid!(scene, labels, ranges; formats = [string],
        sliderkw = Dict(), labelkw = Dict(), valuekw = Dict(), layoutkw...)

    elements = broadcast(labels, ranges, formats) do label, range, format
        slider = Slider(scene; range = range, sliderkw...)
        label = Label(scene, label; halign = :left, labelkw...)
        valuelabel = Label(scene, lift(format, slider.value); halign = :right, valuekw...)
        (; slider = slider, label = label, valuelabel = valuelabel)
    end

    sliders = map(x -> x.slider, elements)
    labels = map(x -> x.label, elements)
    valuelabels = map(x -> x.valuelabel, elements)

    layout = grid!(hcat(labels, sliders, valuelabels); layoutkw...)

    (sliders = sliders, labels = labels, valuelabels = valuelabels, layout = layout)
end

# helper function to create either h or vlines depending on `direction`
# this works only with LAxes because it needs to react to limit changes
function hvlines!(ax::Axis, direction::Int, datavals, axmins, axmaxs; attributes...)

    datavals, axmins, axmaxs = map(x -> x isa Observable ? x : Observable(x), (datavals, axmins, axmaxs))

    linesegs = lift(ax.limits, ax.scene.px_area, datavals, axmins, axmaxs) do lims, pxa,
            datavals, axmins, axmaxs

        xlims = (minimum(lims)[direction], maximum(lims)[direction])
        xfrac(f) = xlims[1] + f * (xlims[2] - xlims[1])
        segs = broadcast(datavals, axmins, axmaxs) do dataval, axmin, axmax
            if direction == 1
                (Point2f0(xfrac(axmin), dataval), Point2f0(xfrac(axmax), dataval))
            elseif direction == 2
                (Point2f0(dataval, xfrac(axmin)), Point2f0(dataval, xfrac(axmax)))
            else
                error("direction must be 1 or 2")
            end
        end
        # handle case that none of the inputs is an array, but we need an array for linesegments!
        if segs isa Tuple
            segs = [segs]
        end
        segs
    end

    linesegments!(ax, linesegs; xautolimits = direction == 2, yautolimits = direction == 1, attributes...)
end

"""
    hlines!(ax::Axis, ys; xmin = 0.0, xmax = 1.0, attrs...)

Create horizontal lines across `ax` at `ys` in data coordinates and `xmin` to `xmax`
in axis coordinates (0 to 1). All three of these can have single or multiple values because
they are broadcast to calculate the final line segments.
"""
hlines!(ax::Axis, ys; xmin = 0.0, xmax = 1.0, attrs...) =
    hvlines!(ax, 1, ys, xmin, xmax; attrs...)

"""
    vlines!(ax::Axis, xs; ymin = 0.0, ymax = 1.0, attrs...)

Create vertical lines across `ax` at `xs` in data coordinates and `ymin` to `ymax`
in axis coordinates (0 to 1). All three of these can have single or multiple values because
they are broadcast to calculate the final line segments.
"""
vlines!(ax::Axis, xs; ymin = 0.0, ymax = 1.0, attrs...) =
    hvlines!(ax, 2, xs, ymin, ymax; attrs...)


text:
### A Pluto.jl notebook ###
# v0.18.1

using Markdown
using InteractiveUtils

# ╔═╡ 49cb409b-e564-47aa-9dae-9bc5bffa991d
using UUIDs

# ╔═╡ b0350bd0-5dd2-4c73-b301-f076123144c2
using FileWatching

# ╔═╡ 729ae3bb-79c2-4fcd-8645-7e0071365537
md"""
# PlutoLinks.jl

Yayyyyyy
"""

# ╔═╡ 968f741b-e70f-4bc1-94b7-ebe8eff868ab
import PlutoHooks: @use_deps

# ╔═╡ c82c8aa9-46a9-4110-88af-8638625222e3
import PlutoHooks: @use_ref

# ╔═╡ 1df0a586-3692-11ec-0171-0b48a4a1c4bd
import PlutoHooks: @use_state

# ╔═╡ 89b3f807-2e24-4454-8f4c-b2a98aee571e
import PlutoHooks: @use_effect

# ╔═╡ bc0e4219-a40b-46f5-adb2-f164d8a9bbdb
import PlutoHooks: @use_memo

# ╔═╡ 9ec99592-955a-41bd-935a-b34f37bb5977
"""
Wraps a `Task` with the current cell. When the cell state is reset, sends an `InterruptException` to the underlying `Task`.
```julia
@use_task([]) do
	while true
		sleep(2.)
		@info "this is updating"
	end
end
```
It can be combined with `@use_state` for background updating of values.

I'm still wondering if it is best to have `deps=nothing` as a default, or have `deps=[]` or maybe even require deps explicitly so people are forced to know what they are doing.
"""
macro use_task(f, deps)
	quote
		try
			@use_deps($(esc(deps))) do
				_, refresh = @use_state(nothing)
				task_ref = @use_ref(Task($(esc(f))))

				@use_effect([]) do
					task = task_ref[]

					schedule(Task() do
						try
							fetch(task)
						finally
							refresh(nothing)
						end
					end)

					schedule(task)

					return function()
						if !istaskdone(task)
							try
								Base.schedule(task, InterruptException(), error=true)
							catch error
								nothing
							end
						end
					end
				end

				task_ref[]
			end
		catch e
			@warn "Got an error in use_task" e
		end
	end
end

# ╔═╡ 74b8a338-9eff-453c-9013-c11d5b57833b
"""
Watches a file and returns its content.

```julia
file_content = @use_file("my_dataset.csv") # file_content will contains the content of my_dataset.csv as a string
```

"""
macro use_file(filename)
	quote
		filename = $(esc(filename))
		update_time = @use_file_change(filename)
		@use_memo([update_time]) do
			read(filename, String)
		end
	end
end

# ╔═╡ 2ae53102-c6f8-4f84-8020-e3b28425240f
macro use_file_change(filename)
	quote
		filename = $(esc(filename))

		@use_deps([filename]) do
			last_update_time, set_last_update_time = @use_state(time())

			@use_task([]) do
				while true
					watch_file(filename)
					set_last_update_time(time())
				end
			end

			last_update_time
		end
	end
end

# ╔═╡ 480dd46c-cc31-46b5-bc2d-2e1680d5c682
function ingredients(path::String)
	# this is from the Julia source code (evalfile in base/loading.jl)
	# but with the modification that it returns the module instead of the last object
	name = Symbol(basename(path))
	m = Module(name)
	Core.eval(m,
        Expr(:toplevel,
             :(eval(x) = $(Expr(:core, :eval))($name, x)),
             :(include(x) = $(Expr(:top, :include))($name, x)),
             :(include(mapexpr::Function, x) = $(Expr(:top, :include))(mapexpr, $name, x)),
             :(include($path))))
	m
end

# ╔═╡ d84f47ba-7c18-4d6c-952c-c9a5748a51f8
"""
Watch a Julia source file and loads its content as a module. The loaded module is then reloaded automatically whenever the file changes.

```julia
MyModule = @ingredients("my_exported_function.jl")
```
"""
macro ingredients(filename)
	quote
		let
			filename = $(esc(filename))
			update_time = @use_file_change(filename)
			@use_memo([update_time]) do
				ingredients(filename)
			end
		end
	end
end

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
FileWatching = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee"
PlutoHooks = "0ff47ea0-7a50-410d-8455-4348d5de0774"
UUIDs = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[compat]
PlutoHooks = "~0.0.3"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[FileWatching]]
uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[PlutoHooks]]
deps = ["FileWatching", "InteractiveUtils", "Markdown", "UUIDs"]
git-tree-sha1 = "f297787f7d7507dada25f6769fe3f08f6b9b8b12"
uuid = "0ff47ea0-7a50-410d-8455-4348d5de0774"
version = "0.0.3"

[[Random]]
deps = ["SHA", "Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
"""

# ╔═╡ Cell order:
# ╟─729ae3bb-79c2-4fcd-8645-7e0071365537
# ╠═968f741b-e70f-4bc1-94b7-ebe8eff868ab
# ╠═c82c8aa9-46a9-4110-88af-8638625222e3
# ╠═1df0a586-3692-11ec-0171-0b48a4a1c4bd
# ╠═89b3f807-2e24-4454-8f4c-b2a98aee571e
# ╠═bc0e4219-a40b-46f5-adb2-f164d8a9bbdb
# ╠═49cb409b-e564-47aa-9dae-9bc5bffa991d
# ╠═b0350bd0-5dd2-4c73-b301-f076123144c2
# ╠═9ec99592-955a-41bd-935a-b34f37bb5977
# ╠═74b8a338-9eff-453c-9013-c11d5b57833b
# ╠═2ae53102-c6f8-4f84-8020-e3b28425240f
# ╠═480dd46c-cc31-46b5-bc2d-2e1680d5c682
# ╠═d84f47ba-7c18-4d6c-952c-c9a5748a51f8
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
using GRF
using Test

@testset "GRF.jl" begin
    # Write your tests here.
end


text:
@json struct Namespace
    id::Int
    name::String
    path::String
    kind::String
    full_path::String
end

@json struct Access
    access_level::Int
    notification_level::Int
end

@json struct Permissions
    project_access::Access
    group_access::Access
end

@json struct License
    key::String
    name::String
    nickname::String
    html_url::String
    source_url::String
end

@json struct ProjectGroup
    group_id::Int
    group_name::String
    group_full_path::String
    group_access_level::Int
end

@json struct Statistics
    commit_count::Int
    storage_size::Int
    repository_size::Int
    lfs_objects_size::Int
    job_artifacts_size::Int
end

@json struct Links
    self::String
    issues::String
    merge_requests::String
    repo_branches::String
    labels::String
    events::String
    members::String
end

@json struct Project
    id::Int
    description::String
    default_branch::String
    visibility::String
    ssh_url_to_repo::String
    http_url_to_repo::String
    web_url::String
    readme_url::String
    tag_list::Vector{String}
    owner::User
    name::String
    name_with_namespace::String
    path::String
    path_with_namespace::String
    issues_enabled::Bool
    open_issues_count::Int
    merge_requests_enabled::Bool
    jobs_enabled::Bool
    wiki_enabled::Bool
    snippets_enabled::Bool
    resolve_outdated_diff_discussions::Bool
    container_registry_enabled::Bool
    created_at::DateTime
    last_activity_at::DateTime
    creator_id::Int
    namespace::Namespace
    import_status::String
    import_error::String
    permissions::Permissions
    archived::Bool
    avatar_url::String
    license_url::String
    license::License
    shared_runners_enabled::Bool
    forks_count::Int
    star_count::Int
    runners_token::String
    public_jobs::Bool
    shared_with_groups::Vector{ProjectGroup}
    only_allow_merge_if_pipeline_succeeds::Bool
    only_allow_merge_if_all_discussions_are_resolved::Bool
    printing_merge_requests_link_enabled::Bool
    request_access_enabled::Bool
    merge_method::String
    statistics::Statistics
    _links => links::Links
end

@json struct FileContents
    file_name::String
    file_path::String
    size::Int
    encoding::String
    content::String
    content_sha256::String
    ref::String
    blob_id::String
    commit_id::String
    last_commit_id::String
end

endpoint(::GitLabAPI, ::typeof(get_user_repos)) =
    Endpoint(:GET, "/projects"; query=Dict("membership" => true))
endpoint(::GitLabAPI, ::typeof(get_user_repos), name::AStr) =
    Endpoint(:GET, "/users/$name/repos")
into(::GitLabAPI, ::typeof(get_user_repos)) = Vector{Project}

endpoint(::GitLabAPI, ::typeof(get_repo), owner_repo::AStr) =
    Endpoint(:GET, "/projects/$(encode(owner_repo))")
endpoint(::GitLabAPI, ::typeof(get_repo), owner::AStr, repo::AStr) =
    Endpoint(:GET, "/projects/$(encode(owner, repo))")
endpoint(::GitLabAPI, ::typeof(get_repo), owner::AStr, subgroup::AStr, repo::AStr) =
    Endpoint(:GET, "/projects/$(encode(owner, subgroup, repo))")
into(::GitLabAPI, ::typeof(get_repo)) = Project

endpoint(::GitLabAPI, ::typeof(create_repo)) =
    Endpoint(:POST, "/projects")
endpoint(::GitLabAPI, ::typeof(create_repo), id::Integer) =
    Endpoint(:POST, "/projects/user/$id")
into(::GitLabAPI, ::typeof(create_repo)) = Project

endpoint(::GitLabAPI, ::typeof(is_collaborator), owner::AStr, repo::AStr, id::Integer) =
    Endpoint(:GET, "/projects/$(encode(owner, repo))/members/$id"; allow_404=true)
postprocessor(::GitLabAPI, ::typeof(is_collaborator)) = DoSomething(ismember)
into(::GitLabAPI, ::typeof(is_collaborator)) = Bool

endpoint(::GitLabAPI, ::typeof(get_file_contents), id::Integer, path::AStr) =
    Endpoint(:GET, "/projects/$id/repository/files/$path"; query=Dict(:ref => "master"))
endpoint(::GitLabAPI, ::typeof(get_file_contents), owner::AStr, repo::AStr, path::AStr) =
    Endpoint(:GET, "/projects/$(encode(owner, repo))/repository/files/$path"; query=Dict(:ref => "master"))
into(::GitLabAPI, ::typeof(get_file_contents)) = FileContents


text:
using DataDeps

register(DataDep(
  "ailerons",
  "https://www.openml.org/d/296",
  "https://www.openml.org/data/get_csv/52060/ailerons.arff",
  "a866e46f5f5d34cd10583ef366bcc183f04cafda773fa55619d17fd9d70dd68d",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=41,
        feature_indices=1:40,
        header=true,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
"""
Finite difference schemes used with momentum descent:
* forward differences
xⁿ⁺¹ = (2-hb)xⁿ + (hb-1)xⁿ⁻¹ - h²∇E(xⁿ)

* backward differences
(1+hb)xⁿ⁺¹ = (2+hb)xⁿ - xⁿ⁻¹ - h²∇E(xⁿ)

* central differences
(2+hb)xⁿ⁺¹ = 4xⁿ + (hb-2)xⁿ⁻¹ - 2h²∇E(xⁿ)
"""

function finite_diff(scheme, X, Fend, b, h, canonical)
   H = canonical?h*h:h
   if scheme == :forward
      return X[end]*(2-h*b) + X[end-1]*(b*h-1) + H*Fend
   elseif scheme == :backward
      return X[end]*(2+h*b)/(1+h*b) + X[end-1]/(-1-h*b) + H*Fend/(b*h+1)
   elseif scheme == :central
       return (4*X[end] + X[end-1]*(b*h-2) + 2*H*Fend)/(2+b*h)
   else
      error("`momentum descent`: unknown finite differences scheme $(scheme)")
   end
end

function criterion(scheme, λh², hb)
   cond_p = nothing; cond_m = nothing;
   if scheme == :forward
      cond_p = abs(2-hb-λh² + sqrt(complex(hb*hb + (λh²)^2 - 4λh² + λh²*hb))) <2
      cond_m = abs(2-hb-λh² - sqrt(complex(hb*hb + (λh²)^2 - 4λh² + λh²*hb))) <2
   elseif scheme == :backward
      cond_p = abs((2+hb-λh² + sqrt(complex((2+hb-λh²)^2-4(1+hb)))) / (1+hb)) <2
      cond_m = abs((2+hb-λh² - sqrt(complex((2+hb-λh²)^2-4(1+hb)))) / (1+hb)) <2
   elseif scheme == :central
      cond_p = abs((2-λh² + sqrt(complex( (2-λh²)^2 - 4 + hb*hb ))) / (2+hb)) <1
      cond_m = abs((2-λh² - sqrt(complex( (2-λh²)^2 - 4 + hb*hb ))) / (2+hb)) <1
   else
      error("`momentum descent`: unknown finite differences scheme $(scheme)")
   end
   return cond_p && cond_m
end

function objective(b, λ)
    return real(b-sqrt(complex(b*b-4*λ)))
end

function b_choice(λ)
    η = 4*real(λ); ξ = 4*imag(λ)
    f = b -> objective(b, λ)
    b = 2*sqrt(real(λ))
    if ξ!=0
        a = 4η; b = -3ξ^2; c = 6*η*ξ^2; d = (η*ξ)^2-4ξ^4
        Δ₀ = b^2-3a*c; Δ₁ = 2b^3 - 9a*b*c + 27d*a^2
        C = ( complex(0.5(Δ₁+sign(Δ₁)*sqrt(complex(Δ₁^2-4Δ₀^3)))) )^(1/3)
        x₁ = -(1/(3a))*(b+C+Δ₀/C)
        k = -0.5-0.5*sqrt(3)*im
        x₂ = -(1/(3a))*(b + k*C + Δ₀/(k*C))
        x₃ = -(1/(3a))*(b + k*k*C + Δ₀/(k*k*C))
        x = [x₁; x₂; x₃]
        x_real = real(x[findmin(abs.(imag(x)))[2]])
        b = sqrt(x_real+η)
    end

#     f_b = maximum([f(bᵢ) for bᵢ in b_all])
#     b = b_all[indmax([f(bᵢ) for bᵢ in b_all])]
    return b
end

function stability(λ)
    η = 4*real(λ); ξ = 4*imag(λ)
    f = b -> objective(b, λ)
    b = b_choice(λ)
    f_b = f(b)

    return f_b, b
end

function forward_criterion(λh², hb)
    condition_p = abs(2-hb-λh² + sqrt(complex(hb*hb + (λh²)^2 - 4λh² + λh²*hb))) <2
    condition_m = abs(2-hb-λh² - sqrt(complex(hb*hb + (λh²)^2 - 4λh² + λh²*hb))) <2
    return condition_p && condition_m
end

function backward_criterion(λh², hb)
    condition_p = abs( (2+hb-λh² + sqrt(complex((2+hb-λh²)^2-4(1+hb)))) / (1+hb) ) <2
    condition_m = abs( (2+hb-λh² - sqrt(complex((2+hb-λh²)^2-4(1+hb)))) / (1+hb) ) <2
    return condition_p && condition_m
end

function central_criterion(λh², hb)
    condition_p = abs((2-λh² + sqrt(complex( (2-λh²)^2 - 4 + hb*hb )))/(2+hb)) <1
    condition_m = abs((2-λh² - sqrt(complex( (2-λh²)^2 - 4 + hb*hb )))/(2+hb)) <1
    return condition_p && condition_m
end

length_n(x, n, P) = in(n,2:length(x)-1) ? norm(P(n), x[n+1] - x[n-1]) :
             error("l_n not defined for n<2 and n>N-1")
t_n(x, n, P) = in(n,2:length(x)-1) ? (x[n+1] - x[n-1])/length_n(x, n, P) : zeros(x[n])

F(x, n, ∇E, P) = - ∇E(x[n]) + dot(t_n(x, n, P), ∇E(x[n])) * (P(n) * t_n(x, n, P))

∂Fⁿ⁺(x, n, ∇E, P) = (1/length_n(x, n, P)) * ( dot(t_n(x, n, P), ∇E(x[n])) * (I - kron(t_n(x, n, P), (P(n) * t_n(x, n, P))')) - kron(t_n(x, n, P), F(x, n, ∇E, P)') )
∂Fⁿ⁻(x, n, ∇E, P) = - ∂Fⁿ⁺(x, n, ∇E, P)

δFⁿ(x, n, ∇∇E, P, P⁻¹∇∇E) = - ( P⁻¹∇∇E(n) - kron(t_n(x, n, P), (∇∇E(n) * t_n(x, n, P))') )

∂Sⁿ(κ, x, n, P) = 2 * κ * kron(t_n(x, n, P), t_n(x, n, P)')

∂Sⁿ⁺(κ, x, n, P) = -κ * (kron(t_n(x, n, P), (P(n) * t_n(x, n, P))') +
(1/length_n(x, n, P)) * kron( P(n)*(x[n+1]-2*x[n]+x[n-1]) - dot(P(n)*(x[n+1]-2*x[n]+x[n-1]), t_n(x, n, P))*t_n(x, n, P), (P(n)*t_n(x, n,P))' ) +
(1/length_n(x, n, P)) * dot(P(n)*(x[n+1]-2*x[n]+x[n-1]), t_n(x, n, P)) * (I - kron(t_n(x, n, P), (P(n)*t_n(x, n, P))')))

∂Sⁿ⁻(κ, x, n, P) = -κ * (kron(t_n(x, n, P), t_n(x, n, P)') -
(1/length_n(x, n, P)) * kron( P(n)*(x[n+1]-2*x[n]+x[n-1]) - dot(P(n)*(x[n+1]-2*x[n]+x[n-1]), t_n(x, n, P))*t_n(x, n, P), (P(n)*t_n(x, n, P))' ) -
(1/length_n(x, n, P)) * dot(P(n)*(x[n+1]-2*x[n]+x[n-1]), t_n(x, n, P)) * (I - kron(t_n(x, n, P), (P(n)*t_n(x, n, P))')))

∂ₓFˣ(x, v, P0) = -P0 + 2 * kron(P0*v, v')

∂ᵥFˣ(x, v, ∇E0) = 2 * dot(∇E0, v) * (I - kron(v, v')) +
                        2 * (kron(∇E0, v') - dot(∇E0, v)*kron(v, v'))

∂ₓFᵛ(x, v, h, P0, Pv) = (- (Pv - P0) + kron((Pv-P0) * v, v') )/h

∂ᵥFᵛ(x, v, Hv, Pv) = -Pv + kron(v, (Pv*v)') +
                           dot(v, Hv)*(I-2*kron(v, v')) +
                           kron(Hv, v') +
                           kron(Pv*v, v') - (v'*Pv*v)*kron(v, v')

function ref{T}(A::Array{Array{T,2},2})
    return cat(1,[cat(2,A[n,:]...) for n=1:size(A,1)]...)
end


text:
using nGraph
using CxxWrap
using Test
using Flux
#using Random, Distributions
#using Flux, BenchmarkTools, Zygote

# Testing Strategy
# 1) Type Tests
# 2) Basic op tests without the Flux backend
#    These tests will be pretty basic because constructing graphps without the Flux backend
#    is pretty tedious and error prone.
#
# 3) Flux Backend

# Set OMP_NUM_THREADS to 1 because the sizes of everything used in the test suite is small
# enough that the overhead of parallelization on large systems hurts performance.
ENV["OMP_NUM_THREADS"] = 1

include("types.jl")
include("ops.jl")
# include("loss.jl")
# include("flux.jl")
# include("backprop.jl")
# include("train.jl")
# include("models/mnist.jl")
# include("models/inception.jl")

# include("models/resnet.jl")
# include("models/inception_v4.jl")


text:
# ---
# title: 650. 2 Keys Keyboard
# id: problem650
# author: Indigo
# date: 2021-06-28
# difficulty: Medium
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/2-keys-keyboard/description/>
# hidden: true
# ---
#
# Initially on a notepad only one character 'A' is present. You can perform two
# operations on this notepad for each step:
#
#   1. `Copy All`: You can copy all the characters present on the notepad (partial copy is not allowed).
#   2. `Paste`: You can paste the characters which are copied **last time**.
#
#
#
# Given a number `n`. You have to get **exactly** `n` 'A' on the notepad by
# performing the minimum number of steps permitted. Output the minimum number of
# steps to get `n` 'A'.
#
# **Example 1:**
#
#
#
#     Input: 3
#     Output: 3
#     Explanation:
#     Intitally, we have one character 'A'.
#     In step 1, we use **Copy All** operation.
#     In step 2, we use **Paste** operation to get 'AA'.
#     In step 3, we use **Paste** operation to get 'AAA'.
#
#
#
#
# **Note:**
#
#   1. The `n` will be in the range [1, 1000].
#
#
#
#
## @lc code=start
using LeetCode

function min_steps650(n::Int)
    res, k = 0, 2
    while n != 1
        while n % k == 0
            res += k
            n ÷= k
        end
        k += 1
    end
    return res
end
## @lc code=end


text:
import MathOptInterface
const MOI = MathOptInterface

import Linda:
    Column,
    Oracle, update!, optimize!, get_columns, get_dual_bound, get_objective_value
const House = DER.House

mutable struct HHOracle <: Oracle

    h::House
    model::MOI.ModelLike

    var2idx::Dict{Tuple{Symbol, Int, Symbol, Int}, MOI.VariableIndex}    # Variable indices in the model
    con2idx::Dict{Tuple{Symbol, Int, Symbol, Int}, MOI.ConstraintIndex}  # Constraint indices in the model

    # Shadow prices
    farkas::Bool
    π::Vector{Float64}
    σ::Float64
end

function update!(o::HHOracle, farkas::Bool, π::Vector{Float64}, σ::Float64)

    # Sanity checks
    T = o.h.num_timesteps
    length(π) == T || throw(DimensionMismatch("π has length $(length(π)) but should have length 2xT=$(2*T)."))

    # Update shadow prices
    o.farkas = farkas
    o.π .= π
    o.σ  = σ

    # Update objective coefficients
    pnet = [
        o.var2idx[(:house, o.h.index, :pnet, t)]
        for t in 1:T
    ]
    for t in 1:T
        MOI.modify(o.model,
            MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
            MOI.ScalarCoefficientChange(pnet[t], (!farkas) * o.h.price[t] - o.π[t])
        )
    end

    return nothing
end

function optimize!(o::HHOracle)
    MOI.optimize!(o.model)
    return nothing
end

function get_columns(o::HHOracle)

    T = o.h.num_timesteps

    # TODO: check solution status
    pnet = [
        MOI.get(o.model, MOI.VariablePrimal(),
            o.var2idx[(:house, o.h.index, :pnet, t)]
        )
        for t in 1:T
    ]

    # Compute real cost
    cost = dot(o.h.price, pnet)

    # Get reduced cost
    rc = MOI.get(o.model, MOI.ObjectiveValue()) - o.σ

    # Create column
    col = Column(
        cost,
        o.h.index,
        collect(1:2*T), pnet,
        true
    )

    return [(col, rc)]
end

function get_dual_bound(o::HHOracle)
    if o.farkas
        return -Inf
    else
        return MOI.get(o.model, MOI.ObjectiveValue()) - o.σ
    end
end

# Create an oracle from a house
function build_oracle(m::MOI.ModelLike, h::House)
    # Build pricing problem
    _, var2idx, con2idx = DR.DER.build_model!(h, m)

    # Setup oracle
    T = h.num_timesteps
    return HHOracle(h, m, var2idx, con2idx, false, zeros(T), 0.0)
end


text:
using DataConvenience
using DataFrames
df = DataFrame(col = rand(1_000_000), col1 = rand(1_000_000), col2 = rand(1_000_000))

fsort(df, :col) # sort by `:col`
fsort(df, [:col1, :col2]) # sort by `:col1` and `:col2`
fsort!(df, :col) # sort by `:col` # sort in-place by `:col`
fsort!(df, [:col1, :col2]) # sort in-place by `:col1` and `:col2`

df = DataFrame(col = rand(1_000_000), col1 = rand(1_000_000), col2 = rand(1_000_000))

using BenchmarkTools
fsort_1col = @belapsed fsort($df, :col) # sort by `:col`
fsort_2col = @belapsed fsort($df, [:col1, :col2]) # sort by `:col1` and `:col2`

sort_1col = @belapsed sort($df, :col) # sort by `:col`
sort_2col = @belapsed sort($df, [:col1, :col2]) # sort by `:col1` and `:col2`

using Plots
bar(["DataFrames.sort 1 col","DataFrames.sort 2 col2", "DataCon.sort 1 col","DataCon.sort 2 col2"],
    [sort_1col, sort_2col, fsort_1col, fsort_2col],
    title="DataFrames sort performance comparison",
    label = "seconds")

using DataFrames
using CSV

df = DataFrame(a = rand(1_000_000), b = rand(Int8, 1_000_000), c = rand(Int8, 1_000_000))

filepath = tempname()*".csv"
CSV.write(filepath, df)

for chunk in CsvChunkIterator(filepath)
  print(describe(chunk))
end

# read all column as String
for chunk in CsvChunkIterator(filepath, type=String)
    print(describe(chunk))
end

# read a three colunms csv where the column types are String, Int, Float32
for chunk in CsvChunkIterator(filepath, types=[String, Int, Float32])
  print(describe(chunk))
end

@replicate 10 8

x = Vector{Union{Missing, Int}}(undef, 10_000_000)

cmx = count_missing(x) # this is faster

cmx2 = countmissing(x) # this is faster

cimx = count(ismissing, x) # the way available at base

cmx == cimx # true

max_stars_repo_name: xiaodaigh/DataConvenience

max_issues_repo_name: xiaodaigh/DataConvenience

max_forks_repo_name: xiaodaigh/DataConvenience


text:
using MPI
using p4est_wrapper
using Test

# Initialize MPI if not initialized yet
if !MPI.Initialized()
    MPI.Init()
end

#############################################################################
# Definition of data structures and function callbacks
#############################################################################

# Refine callback
function my_refine( ::Ptr{p4est_t}, which_tree::p4est_topidx_t, quadrant::Ptr{p4est_quadrant_t})
    @assert which_tree == 0
    q = quadrant[]
    ((q.level == 0) || (q.level < 2 && p4est_quadrant_child_id(quadrant) == 1)) && return Cint(1)
    return Cint(0)
end

# C-callable refine callback
const my_refine_c = @cfunction(my_refine, Cint, (Ptr{p4est_t}, p4est_topidx_t, Ptr{p4est_quadrant_t}))

#############################################################################
# Main program
#############################################################################

mpicomm = p4est_wrapper.P4EST_ENABLE_MPI ? MPI.COMM_WORLD : Cint(0)

# Create a connectivity structure for the unit square.
unitsquare_connectivity = p4est_connectivity_new_unitsquare()
@assert Bool(p4est_connectivity_is_valid(unitsquare_connectivity))

# Create a new forest
unitsquare_forest = p4est_new(mpicomm, unitsquare_connectivity, sizeof(p4est_quadrant_t), C_NULL, C_NULL)
@test unitsquare_forest != C_NULL

# Register callback function to decide for refinement.
p4est_refine(unitsquare_forest, 1, my_refine_c, C_NULL)

# Equally partition the forest. The partition can be by element count or by a user-defined weight.
p4est_partition(unitsquare_forest, 0, C_NULL);

# Build the ghost layer.
ptr_to_p4est_ghost = p4est_ghost_new(unitsquare_forest, p4est_wrapper.P4EST_CONNECT_FULL)
p4est_ghost = ptr_to_p4est_ghost[]

# Obtain ghost quadrants
##p4est_quadrant_t * ptr_ghost_quadrants = (p4est_quadrant_t *) p4est_ghost->ghosts.array;
ptr_ghost_quadrants = Ptr{p4est_quadrant_t}(p4est_ghost.ghosts.array)
proc_offsets = unsafe_wrap(Array, p4est_ghost.proc_offsets, p4est_ghost.mpisize+1)

for i=1:p4est_ghost.mpisize
    for j=proc_offsets[i]:proc_offsets[i+1]-1
        quadrant = ptr_ghost_quadrants[j+1]
        piggy3 = quadrant.p.piggy3
        @test Bool(p4est_quadrant_is_valid(ptr_ghost_quadrants+(sizeof(p4est_quadrant_t)*j)))
        p4est_quadrant_print(p4est_wrapper.SC_LP_INFO, ptr_ghost_quadrants+(sizeof(p4est_quadrant_t)*j))
        print("(i, j, local_num, which_tree) ($(i),$(j),$(piggy3.local_num),$(piggy3.which_tree)) \n")
    end
end

# Destroy the ghost structure
p4est_ghost_destroy(ptr_to_p4est_ghost)
# Destroy the forest
p4est_destroy(unitsquare_forest)
# Destroy the connectivity
p4est_connectivity_destroy(unitsquare_connectivity)

# Finalize MPI if initialized and session is not interactive
if (MPI.Initialized() && !isinteractive())
    MPI.Finalize()
end


text:
# ----------------------------------------------------------------
"""
    Take an expression and spawn it
"""
function _gw_spawn(gw::GitWorker, ex::Expr, tname;
        follow = true, tout = 120.0,
        wt = 5.0, verb = false
    )

    !Meta.isexpr(ex, [:let, :block]) && error("A 'begin or let' block is expected. Ex: '@gw_spawn begin println(\"Hi\")' end")

    # task id
    tid = _task_id(tname)

    # readme
    readme = _readme(;tid,
        src = _expr_src(ex),
        desc = "Hi, no description provided!!!",
        lang = "julia"
    )

    ok_flag = _upload_jltask(gw,
        tid, ex; readme,
        vtime = _GW_DFLT_TASK_VTIME
    )

    if ok_flag
        @info("Task uploaded")
    else
        @error("Task uploading fail")
    end

    return :(nothing)
end

_gw_spawn(f::Function, taskid; kwargs...) = _gw_spawn(f(), taskid; kwargs...)

"""
    Take the return f()::Expr and spawn it
"""
gw_spawn(ex::Expr; kwargs...) = _gw_spawn(ex, _task_id("SPANW"); kwargs...)
gw_spawn(f::Function; kwargs...) = _gw_spawn(f(), _task_id("SPANW"); kwargs...)

# ----------------------------------------------------------------
"""
    Capture the expression and spawn it
"""
macro gw_spawn(ex)
    GitWorkers.gw_spawn(ex; follow = true, tout = 120.0, wt = 5.0)
    return :(nothing)
end

max_stars_repo_name: josePereiro/GitWorkers

max_issues_repo_name: josePereiro/GitWorkers

max_forks_repo_name: josePereiro/GitWorkers


text:
# loaders.jl
# utilities for loading in data / images

using FileIO: load

"""
    data = ir_load_brainweb_t1_256()

Load brainweb T1-weighted MRI slice of size `256 × 256`
"""
function ir_load_brainweb_t1_256()
	dirmirt = normpath(dirname(pathof(MIRT)), "..")
	file = joinpath(dirmirt, "data", "mri", "brainweb_t1_256.fld")
	return Float32.(load(file))' # transpose!
end


text:
function newton_solver(res, x0; iterations=500, tolerance=1e-6)
    x = x0
    for i = 1:iterations
        r = res(x)
        (norm(r, Inf) < tolerance) && break
        ∇r = FiniteDiff.finite_difference_derivative(x -> res(x), x)
        Δx = - r / ∇r
        α = 1.0
        r_cand = Inf
        while norm(r_cand, Inf) >= norm(r, Inf)
            r_cand = res(x + α * Δx)
            α /= 2
        end
        x += α * Δx
        (i == iterations) && (@show "solver failure")
    end
    return x
end

function catenary_parameters(x_start, x_goal, rope_length; iterations=500, tolerance=1e-6,
        a_guess=0.0, dx_guess=0.0)

    rope_length =  max(rope_length, (1+1e-5)*norm(x_goal - x_start))
    h = x_goal[1] - x_start[1]
    v = x_goal[2] - x_start[2]

    # find a i.e. the shape of the catenary
    res(b) = 1/sqrt(sqrt(rope_length^2 - v^2)/h - 1) - 1/sqrt(2b * sinh(1/(2b)) - 1)
    b = newton_solver(res, a_guess/h, iterations=iterations, tolerance=tolerance)
    a = b * h

    # find x offset
    function caten(x; a=a)
        return a * cosh(x/a)
    end
    res = dx -> caten(x_goal[1] + dx) - caten(x_start[1] + dx) - (x_goal[2] - x_start[2])
    dx = newton_solver(res, dx_guess, iterations=iterations, tolerance=tolerance)

    # find y offset
    dy = x_goal[2] - caten(x_goal[1] + dx)

    # final function
    return a, dx, dy
end

function catenary(x; a=0.0, dx=0.0, dy=0.0)
    return a * cosh((x+dx)/a) + dy
end

function link_transform(start, goal)
    # transforms a vertical line of length 1 into a line between start and goal
    v1 = [0.0, 0.0, 1.0]
    v2 = goal[1:3] - start[1:3]
    normalize!(v2)
    ax = cross(v1, v2)
    ang = acos(v1' * v2)
    q = axis_angle_to_quaternion(ang * normalize!(ax))

    rope_length = norm(goal - start)
    scaling = LinearMap(I * Diagonal([1.0, 1.0, rope_length]))
    transform = compose(Translation(start), LinearMap(q))
    return transform, scaling
end

function build_rope(vis::Visualizer; N::Int=10, color=Colors.RGBA(0,0,0,1),
        rope_type::Symbol=:cylinder, rope_radius=0.02, name=:rope)

    if rope_type == :line
        line = MeshCat.Line([xa, xb])
        material = LineBasicMaterial(color=color, wireframeLinewidth=10, linewidth=10)
    elseif rope_type == :cylinder
        line = MeshCat.Cylinder(Point(0,0,0.0), Point(0,0,1.0), rope_radius)
        material = MeshPhongMaterial(color=color)
    end

    for i = 1:N
        setobject!(vis[name][Symbol(i)][:scaled], line, material)
    end

    return vis
end

function set_straight_rope(vis::Visualizer, x_start, x_goal; N::Int=10, name::Symbol=:rope)
    Λ = range(0,1,length=N+1)
    x = [(1-λ)*x_start + λ * x_goal for λ in Λ]
    for i = 1:N
        transform, scaling = link_transform(x[i], x[i+1])
        settransform!(vis[name][Symbol(i)][:scaled], scaling)
        settransform!(vis[name][Symbol(i)], transform)
    end
    return nothing
end

function set_loose_rope(vis::Visualizer, x_start, x_goal; N::Int=10,
        rope_length=2norm(x_goal - x_start), min_altitude=-Inf, a_guess=1.0, dx_guess=0.0, name::Symbol=:rope)
    v = x_goal - x_start
    shadow_rope_length = norm(v[1:2])
    θ = atan(v[2], v[1])
    R = rotationmatrix(RotZ(-θ))
    v̄ = R * v # rotated into the xz plane

    a, dx, dy = catenary_parameters(zeros(2), v̄[[1,3]], rope_length, a_guess=a_guess, dx_guess=dx_guess)
    Λ = shadow_rope_length * range(0,1,length=N+1)
    x = []
    for i = 1:N+1
        xi = x_start + R' * [Λ[i], 0, catenary(Λ[i], a=a, dx=dx, dy=dy)]
        xi[3] = max(xi[3], min_altitude)
        push!(x, xi)
    end
    for i = 1:N
        transform, scaling = link_transform(x[i], x[i+1])
        settransform!(vis[name][Symbol(i)][:scaled], scaling)
        settransform!(vis[name][Symbol(i)], transform)
    end
    return a, dx
end

function animate_straight_rope(vis::Visualizer, start_traj::Vector, goal_traj::Vector;
        anim::Animation=MeshCat.Animation(100), N::Int=50, name=:rope)
    M = length(start_traj)

    for i = 1:M
        atframe(anim, i) do
            xa = start_traj[i]
            xb = goal_traj[i]
            set_straight_rope(vis, xa, xb, N=N, name=name)
        end
    end
    setanimation!(vis, anim)
    return vis, anim
end

function animate_loose_rope(vis::Visualizer, start_traj::Vector, goal_traj::Vector;
        anim::Animation=MeshCat.Animation(100), rope_length=30.0, N::Int=50, min_altitude=-Inf, name=:rope)
    M = length(start_traj)

    a_guess = 1.0
    dx_guess = 0.0

    for i = 1:M
        atframe(anim, i) do
            xa = start_traj[i]
            xb = goal_traj[i]
            a_guess, dx_guess = set_loose_rope(vis, xa, xb, rope_length=rope_length,
                N=N, min_altitude=min_altitude, a_guess=a_guess, dx_guess=dx_guess, name=name)
            # set_straight_rope(vis, xa, xb, N=N, name=name)
        end
    end
    setanimation!(vis, anim)
    return vis, anim
end


text:
# test to extract quartet from a "good" 5taxon network (Case G)
# Claudia November 2014
# also test to extract all quartets from a Data object
# Claudia January 2015

# Case G ---------

include("../examples/case_g_example.jl");
net.names

q1 = Quartet(1,["6","7","4","8"],[0.5,0.4,0.1]);
qnet = extractQuartet!(net,q1);
printEdges(qnet)
printNodes(qnet)

q2 = Quartet(2,["6","7","10","8"],[0.5,0.4,0.1]);
qnet = extractQuartet!(net,q2);
printEdges(qnet)
printNodes(qnet)

q3 = Quartet(3,["10","7","4","8"],[0.5,0.4,0.1]);
qnet = extractQuartet!(net,q3);
printEdges(qnet)
printNodes(qnet)

q4 = Quartet(4,["6","10","4","8"],[0.5,0.4,0.1]);
qnet = extractQuartet!(net,q4);
printEdges(qnet)
printNodes(qnet)

q5 = Quartet(5,["6","7","4","10"],[0.5,0.4,0.1]);
qnet = extractQuartet!(net,q5);
printEdges(qnet)
printNodes(qnet)

[n.number for n in qnet.leaf]
qnet.numTaxa

# Bad triangle ------------

## include("../bad_triangle_example.jl");
## printEdges(net)
## printNodes(net)
## net.names

## q1 = Quartet(1,["6","7","1","5"],[0.5,0.4,0.1]);
## qnet = extractQuartet!(net,q1)
## printEdges(qnet)
## printNodes(qnet)

## q2 = Quartet(2,["6","5","1","8"],[0.5,0.4,0.1]);
## qnet = extractQuartet!(net,q2);
## printEdges(qnet)
## printNodes(qnet)

## q3 = Quartet(3,["1","7","5","8"],[0.5,0.4,0.1]);
## qnet = extractQuartet!(net,q3);
## printEdges(qnet)
## printNodes(qnet)
## # correct, but does not update length of external edge 8 (should be 0.2, at is 0.1)

## q4 = Quartet(4,["6","1","7","8"],[0.5,0.4,0.1]);
## qnet = extractQuartet!(net,q4);
## printEdges(qnet)
## printNodes(qnet)
## # correct, but does not update length of external edge 9 (should be 0.2, at is 0.1)

## q5 = Quartet(5,["6","7","5","8"],[0.5,0.4,0.1]);
## qnet = extractQuartet!(net,q5);
## printEdges(qnet)
## printNodes(qnet)

## [n.number for n in qnet.leaf]
## qnet.numTaxa

# Bad diamond ------------

include("../examples/case_f_example.jl")
printEdges(net)
printNodes(net)
net.names
parameters!(net)
net.numht

parameters!(qnet,net)

q1 = Quartet(1,["6","7","4","8"],[0.5,0.4,0.1]);
qnet = extractQuartet!(net,q1);
printEdges(qnet)
printNodes(qnet)
qnet.indexht
qnet.hasEdge

q2 = Quartet(2,["6","7","10","8"],[0.5,0.4,0.1]);
qnet = extractQuartet!(net,q2);
printEdges(qnet)
printNodes(qnet)

q3 = Quartet(3,["10","7","4","8"],[0.5,0.4,0.1]);
qnet = extractQuartet!(net,q3);
printEdges(qnet)
printNodes(qnet)

q4 = Quartet(4,["6","10","4","8"],[0.5,0.4,0.1]);
qnet = extractQuartet!(net,q4);
printEdges(qnet)
printNodes(qnet)

q5 = Quartet(5,["6","7","4","10"],[0.5,0.4,0.1]);
qnet = extractQuartet!(net,q5);
printEdges(qnet)
printNodes(qnet)

[n.number for n in qnet.leaf]
qnet.numTaxa

# ------------------------------------
# extract all quartets
include("../examples/case_g_example.jl");

q1 = Quartet(1,["6","7","4","8"],[0.5,0.4,0.1]);
q2 = Quartet(2,["6","7","10","8"],[0.5,0.4,0.1]);
q3 = Quartet(3,["10","7","4","8"],[0.5,0.4,0.1]);
q4 = Quartet(4,["6","10","4","8"],[0.5,0.4,0.1]);
q5 = Quartet(5,["6","7","4","10"],[0.5,0.4,0.1]);

d = DataCF([q1,q2,q3,q4,q5]);
extractQuartet!(net,d);
calculateExpCFAll!(d)


text:
export Instruction, unsafe_delete!,
       metadata, metadata!,
       opcode,
       predicate_int, predicate_real

# forward definition of Instruction in src/core/value/constant.jl
identify(::Type{Value}, ::Val{API.LLVMInstructionValueKind}) = Instruction

identify(::Type{Instruction}, ref::API.LLVMValueRef) =
    identify(Instruction, Val{API.LLVMGetInstructionOpcode(ref)}())
identify(::Type{Instruction}, ::Val{K}) where {K} = error("Unknown instruction kind $K")

@inline function refcheck(::Type{T}, ref::API.LLVMValueRef) where T<:Instruction
    ref==C_NULL && throw(UndefRefError())
    T′ = identify(Instruction, ref)
    if T != T′
        error("invalid conversion of $T′ instruction reference to $T")
    end
end

# Construct a concretely typed instruction object from an abstract value ref
function Instruction(ref::API.LLVMValueRef)
    ref == C_NULL && throw(UndefRefError())
    T = identify(Instruction, ref)
    return T(ref)
end

Instruction(inst::Instruction) =
    Instruction(API.LLVMInstructionClone(inst))

unsafe_delete!(::BasicBlock, inst::Instruction) =
    API.LLVMInstructionEraseFromParent(inst)
Base.delete!(::BasicBlock, inst::Instruction) =
    API.LLVMInstructionRemoveFromParent(inst)

parent(inst::Instruction) =
    BasicBlock(API.LLVMGetInstructionParent(inst))

opcode(inst::Instruction) = API.LLVMGetInstructionOpcode(inst)

predicate_int(inst::Instruction) = API.LLVMGetICmpPredicate(inst)
predicate_real(inst::Instruction) = API.LLVMGetFCmpPredicate(inst)

## metadata iteration
# TODO: doesn't actually iterate, since we can't list the available keys

@enum(MD, MD_dbg = 0,
          MD_tbaa = 1,
          MD_prof = 2,
          MD_fpmath = 3,
          MD_range = 4,
          MD_tbaa_struct = 5,
          MD_invariant_load = 6,
          MD_alias_scope = 7,
          MD_noalias = 8,
          MD_nontemporal = 9,
          MD_mem_parallel_loop_access = 10,
          MD_nonnull = 11,
          MD_dereferenceable = 12,
          MD_dereferenceable_or_null = 13,
          MD_make_implicit = 14,
          MD_unpredictable = 15,
          MD_invariant_group = 16,
          MD_align = 17,
          MD_loop = 18,
          MD_type = 19,
          MD_section_prefix = 20,
          MD_absolute_symbol = 21,
          MD_associated = 22)

export InstructionMetadataDict

struct InstructionMetadataDict <: AbstractDict{MD,MetadataAsValue}
    inst::Instruction
end

metadata(inst::Instruction) = InstructionMetadataDict(inst)

Base.isempty(md::InstructionMetadataDict) =
  !convert(Core.Bool, API.LLVMHasMetadata(md.inst))

Base.haskey(md::InstructionMetadataDict, kind::MD) =
  API.LLVMGetMetadata(md.inst, kind) != C_NULL

function Base.getindex(md::InstructionMetadataDict, kind::MD)
    objref = API.LLVMGetMetadata(md.inst, kind)
    objref == C_NULL && throw(KeyError(name))
    return MetadataAsValue(objref)
  end

Base.setindex!(md::InstructionMetadataDict, node::MetadataAsValue, kind::MD) =
    API.LLVMSetMetadata(md.inst, kind, node)

Base.delete!(md::InstructionMetadataDict, kind::MD) =
    API.LLVMSetMetadata(md.inst, kind, C_NULL)

## instruction types

const opcodes = [:Ret, :Br, :Switch, :IndirectBr, :Invoke, :Unreachable, :Add, :FAdd, :Sub,
                 :FSub, :FNeg, :Mul, :FMul, :UDiv, :SDiv, :FDiv, :URem, :SRem, :FRem, :Shl, :LShr,
                 :AShr, :And, :Or, :Xor, :Alloca, :Load, :Store, :GetElementPtr, :Trunc,
                 :ZExt, :SExt, :FPToUI, :FPToSI, :UIToFP, :SIToFP, :FPTrunc, :FPExt,
                 :PtrToInt, :IntToPtr, :BitCast, :AddrSpaceCast, :ICmp, :FCmp, :PHI, :Call,
                 :Select, :UserOp1, :UserOp2, :VAArg, :ExtractElement, :InsertElement,
                 :ShuffleVector, :ExtractValue, :InsertValue, :Fence, :AtomicCmpXchg,
                 :AtomicRMW, :Resume, :LandingPad, :CleanupRet, :CatchRet, :CatchPad,
                 :CleanupPad, :CatchSwitch]

for op in opcodes
    typename = Symbol(op, :Inst)
    enum = Symbol(:LLVM, op)
    @eval begin
        @checked struct $typename <: Instruction
            ref::API.LLVMValueRef
        end
        identify(::Type{Instruction}, ::Val{API.$enum}) = $typename
    end
end

## call sites and invocations

export callconv, callconv!,
       istailcall, tailcall!,
       called_value

callconv(inst::Instruction) = API.LLVMGetInstructionCallConv(inst)
callconv!(inst::Instruction, cc) =
    API.LLVMSetInstructionCallConv(inst, cc)

istailcall(inst::Instruction) = convert(Core.Bool, API.LLVMIsTailCall(inst))
tailcall!(inst::Instruction, bool) = API.LLVMSetTailCall(inst, convert(Bool, bool))

called_value(inst::Instruction) = Value(API.LLVMGetCalledValue(inst))

## terminators

export isterminator, isconditional, condition, condition!, default_dest

isterminator(inst::Instruction) = API.LLVMIsATerminatorInst(inst) != C_NULL

isconditional(br::Instruction) = convert(Core.Bool, API.LLVMIsConditional(br))

condition(br::Instruction) =
    Value(API.LLVMGetCondition(br))
condition!(br::Instruction, cond::Value) =
    API.LLVMSetCondition(br, cond)

default_dest(switch::Instruction) =
    BasicBlock(API.LLVMGetSwitchDefaultDest(switch))

# successor iteration

export successors

struct TerminatorSuccessorSet
    term::Instruction
end

successors(term::Instruction) = TerminatorSuccessorSet(term)

Base.eltype(::TerminatorSuccessorSet) = BasicBlock

function Base.getindex(iter::TerminatorSuccessorSet, i)
    @boundscheck 1 <= i <= length(iter) || throw(BoundsError(iter, i))
    return BasicBlock(API.LLVMGetSuccessor(iter.term, i-1))
end

Base.setindex!(iter::TerminatorSuccessorSet, bb::BasicBlock, i) =
    API.LLVMSetSuccessor(iter.term, i-1, bb)

function Base.iterate(iter::TerminatorSuccessorSet, i=1)
    i >= length(iter) + 1 ? nothing : (iter[i], i+1)
end

Base.length(iter::TerminatorSuccessorSet) = API.LLVMGetNumSuccessors(iter.term)

Base.lastindex(iter::TerminatorSuccessorSet) = length(iter)

## phi nodes

# incoming iteration

export PhiIncomingSet

struct PhiIncomingSet
    phi::Instruction
end

incoming(phi::Instruction) = PhiIncomingSet(phi)

Base.eltype(::PhiIncomingSet) = Tuple{Value,BasicBlock}

function Base.getindex(iter::PhiIncomingSet, i)
    @boundscheck 1 <= i <= length(iter) || throw(BoundsError(iter, i))
    return tuple(Value(API.LLVMGetIncomingValue(iter.phi, i-1)),
                       BasicBlock(API.LLVMGetIncomingBlock(iter.phi, i-1)))
end

function Base.append!(iter::PhiIncomingSet, args::Vector{Tuple{V, BasicBlock}} where V <: Value)
    vals, blocks = zip(args...)
    API.LLVMAddIncoming(iter.phi, collect(vals),
                        collect(blocks), length(args))
end

Base.push!(iter::PhiIncomingSet, args::Tuple{<:Value, BasicBlock}) = append!(iter, [args])

Base.length(iter::PhiIncomingSet) = API.LLVMCountIncoming(iter.phi)


text:
using CUDAnative, CuArrays, BenchmarkTools, StaticArrays

function tiled_gpu_transpose(a_in, a_out, tile_res, dim1, dim2)
    i = (blockIdx().x-1) * blockDim().x + threadIdx().x
    @inbounds cartesian_in = CartesianIndices(a_in)[i]

    tile = @cuDynamicSharedMem(Float32, tile_res)

    @inbounds tile[threadIdx().x] = a_in[i]

    sync_threads()

    # Changing index order for transpose.
    # Note: the only way I can see to do this is to transfer it into an array
    #       for a swap, and then transfer back
    tmp_idx = MVector(cartesian_in.I)
    @inbounds tmp_idx[dim1], tmp_idx[dim2] = tmp_idx[dim2], tmp_idx[dim1]
    cartesian_out = CartesianIndex(tmp_idx.data)

    @inbounds a_out[cartesian_out] = tile[threadIdx().x]

    return nothing
end

function naive_gpu_transpose(a_in, a_out, dim1, dim2)
    i = (blockIdx().x-1) * blockDim().x + threadIdx().x
    @inbounds cartesian_in = CartesianIndices(a_in)[i]

    # Changing index order for transpose.
    # Note: the only way I can see to do this is to transfer it into an array
    #       for a swap, and then transfer back
    tmp_idx = MVector(cartesian_in.I)
    @inbounds tmp_idx[dim1], tmp_idx[dim2] = tmp_idx[dim2], tmp_idx[dim1]
    cartesian_out = CartesianIndex(tmp_idx.data)

    @inbounds a_out[cartesian_in] = a_in[cartesian_out]
    return nothing
end

function gpu_copy(a_in, a_out)
    i = (blockIdx().x-1) * blockDim().x + threadIdx().x

    @inbounds a_out[i] = a_in[i]
    return nothing
end

function tiled_gpu_copy(a_in, a_out, tile_res)
    i = (blockIdx().x-1) * blockDim().x + threadIdx().x

    tile = @cuDynamicSharedMem(Float32, tile_res)

    @inbounds tile[threadIdx().x] = a_in[i]

    sync_threads()

    @inbounds a_out[i] = tile[threadIdx().x]

    return nothing
end

function square_tiled_gpu_copy(a_in, a_out, tile_res)
    i = (blockIdx().x-1) * blockDim().x + threadIdx().x
    @inbounds cartesian_in = CartesianIndices(a_in)[i]

    tile = @cuDynamicSharedMem(Float32, (sqrt(tile_res), sqrt(tile_res)))

    @inbounds tile[threadIdx().x] = a_in[i]

    sync_threads()

    @inbounds a_out[i] = tile[threadIdx().x]

    return nothing
end

function main()

    res = 8192
    # blocks must be defined based on tile size that can fit into shared mem
    tile_res = 1024
    #res = 1024

    a = round.(rand(Float32, (res, res))*100)
    #a = ones(Float32, (res, res))
    d_a = CuArray(a)
    d_b = CuArray(zeros(Float32, res, res))

    println("Copy time is:")
    CuArrays.@time @cuda(threads = (tile_res, 1, 1),
                         blocks = (div(res*res,tile_res),1,1),
                         gpu_copy(d_a, d_b))

    println("Transpose time is:")
    CuArrays.@time @cuda(threads = (tile_res, 1, 1),
                         blocks = (div(res*res,tile_res),1,1),
                         naive_gpu_transpose(d_a, d_b, 1, 2))

    println("Tiled copy time is:")
    CuArrays.@time @cuda(threads = (tile_res, 1, 1),
                         blocks = (div(res*res,tile_res),1,1),
                         shmem = sizeof(Float32)*tile_res,
                         tiled_gpu_copy(d_a, d_b, tile_res))

    println("Tiled transpose time is:")
    CuArrays.@time @cuda(threads = (tile_res, 1, 1),
                         blocks = (div(res*res,tile_res),1,1),
                         shmem = sizeof(Float32)*tile_res,
                         tiled_gpu_transpose(d_a, d_b, tile_res, 1, 2))

    a = Array(d_a)
    b = Array(d_b)

    if (a == transpose(b))
    #if (a == b)
        println("Good job, man")
    else
        println("You failed. Sorry.")
        return a .- b
    end

    return nothing
end

main()