text:
@info("Load packages...")
windows = false
cd("/remote/ceph/group/gedet/data/sim/2020/2020-02-06_8380701d_st_ffischer")
include("gerda_data_processing/init.jl");
include("gerda_data_processing/fct.jl");

cal = JSON.parsefile("data_pulses/FEP_pulses/calibration/cal.json");

plots_base_path = "data_pulses/baselines_filtered/plots/RMS/"

pre_cut = 50        # keV
multiplicity = 1;
std_intersect = 2;
number_of_pulses = 1e4

# Find all detectors in baselines_raw
folders = glob("data_pulses/baselines_raw/*/*")
channels = []
for folder in folders
    ch = parse(Int64, split(split(folder, "/")[end], "-")[1])
    ch in channels ? "Channel is already in the list" : push!(channels, ch)
end
# Search in baselines_raw for all files for detector X
for eval_channel in channels
    Base.run(`clear`)
    ch_str = lpad(eval_channel, 2, "0");
    raw_files = glob("data_pulses/baselines_raw/*/" * ch_str * "-" * channel_to_name[eval_channel] * "/*")
    @info("Starting with Ch"*ch_str*" | "*channel_to_name[eval_channel]*" including "*string(length(raw_files))*" files:")
    data = Table( energy     = [],
                multiplicity = [],
                timestamp    = [],
                channel      = [],
                waveform     = [])

    # Read in all files -> data
    pro = Progress(length(raw_files), dt=0.5,
            barglyphs=BarGlyphs('|','█', ['▁' ,'▂' ,'▃' ,'▄' ,'▅' ,'▆', '▇'],' ','|',),
            barlen=10)
    for file in raw_files
        append!(data, HDF5.h5open(file, "r") do h5f
            LegendHDF5IO.readdata(h5f, "data")
        end)
        next!(pro)
    end
    @info(string(size(data, 1))*" raw baselines found")
    @info("Start X-Talk filtering")

    # Energy calibration
    waveforms = []
    for wf in data.waveform
        pulse = wf.value .- sum(wf.value)/length(wf.value)
        pulse .*= cal[string(eval_channel)]["value"]
        push!(waveforms, RDWaveform(wf.time, pulse))
    end
    # Filter for X-Talk
    rms = []
    for wf in waveforms
        pulse = wf.value
        if maximum(abs.(pulse)) < pre_cut
            temp = sqrt(sum(pulse.^2) / length(pulse))
            push!(rms, temp)
        end
    end

    f = StatsBase.fit(Histogram, rms, closed=:left, minimum(rms):0.03:maximum(rms));
    x_hist = float.(f.edges[1][1:end-1]);
    y_hist = float.(f.weights);
    fit = curve_fit(model, x_hist, y_hist, [mean(rms), 1, maximum(y_hist)], lower=[minimum(x_hist), 0.01, maximum(y_hist)/2])
    x_fit = minimum(rms):0.01:maximum(rms)
    y_fit = model(x_fit, fit.param)
    max_x = findall(x -> x == maximum(y_fit), y_fit)[1]
    @info("Creating histogram and fit")
    p = plot(x_hist, y_hist, label="Baseline RMS");
    p = plot!(x_fit, y_fit, label="Fit");
    p = scatter!([x_fit[max_x]], [y_fit[max_x]], label="Maximum at "*string(round(x_fit[max_x], digits=2)))
    p = vline!([x_fit[max_x] - 2 * fit.param[2], x_fit[max_x] + 2 * fit.param[2]], label="2 sigma", legend = :outertopright)
    p = plot!(xlim=[x_fit[max_x] - 5 * fit.param[2], x_fit[max_x] + 5 * fit.param[2]])
    if !isdir(plots_base_path)
        mkpath(plots_base_path)
    end
    ch_str = lpad(eval_channel, 2, "0");
    @suppress savefig(p, plots_base_path*ch_str*"-"*channel_to_name[eval_channel]*".png");
    @suppress savefig(p, plots_base_path*ch_str*"-"*channel_to_name[eval_channel]*".pdf");

    bl_index = []
    for bl in 1:1:length(waveforms)
        wf = waveforms[bl]
        temp = sqrt(sum(wf.value.^2) / length(wf.value))
        if temp >= fit.param[1] - 2*fit.param[2] && temp <= fit.param[1] + 2*fit.param[2]
            push!(bl_index, bl)
        end
    end

    ch_str = lpad(eval_channel, 2, "0");
    filename = "data_pulses/baselines_filtered/"*ch_str*"-"*channel_to_name[eval_channel]*"/"
    if !isdir(filename)
        mkpath(filename)
    end
    bl_parts = []
    if length(bl_index) > number_of_pulses
        for i in 1:1:Int(ceil(length(bl_index)/number_of_pulses))
            if i*number_of_pulses <= length(bl_index)
                part = ((i-1)*number_of_pulses + 1):1:i*number_of_pulses
            else
                part = ((i-1)*number_of_pulses + 1):1:length(bl_index)
            end
            push!(bl_parts, bl_index[Int.(part)])
        end
    else
        push!(bl_parts, bl_index)
    end
    @info("Storing data")
    pro = Progress(length(bl_parts), dt=0.5,
            barglyphs=BarGlyphs('|','█', ['▁' ,'▂' ,'▃' ,'▄' ,'▅' ,'▆', '▇'],' ','|',),
            barlen=10)
    for part in bl_parts
        part_num  = findall(x->x == part, bl_parts)[1]
        if part_num != length(bl_parts)
            init  = (part_num-1)*number_of_pulses +1
            final = part_num*number_of_pulses
        else
            init = (part_num-1)*number_of_pulses +1
            final = length(bl_index)
        end
        str_start = lpad( Int(init) , 7, "0")
        str_end   = lpad( Int(final), 7, "0")
        filename = "data_pulses/baselines_filtered/"*ch_str*"-"*channel_to_name[eval_channel]*"/"
        filename *= str_start*"-"*str_end*".h5"
        HDF5.h5open(filename, "w") do h5f
            LegendHDF5IO.writedata( h5f, "data", Table( energy       = float.(data.energy[part]),
                                                        multiplicity = float.(data.multiplicity[part]),
                                                        timestamp    = float.(data.timestamp[part]),
                                                        channel      = float.(data.channel[part]),
                                                        waveform     = StructArray{RDWaveform}(Array{RDWaveform,1}(waveforms[part]))
                                                        ))
        end
        next!(pro)
    end

end


text:
using garc

push!(ARGS, "arg")
garc.julia_main()

ext: jl


text:
import JuLIP.Potentials: @pot, @D, evaluate, evaluate_d, PairPotential

@pot type Morseold <: PairPotential
   e0::Float64
   A::Float64
   r0::Float64
end
evaluate(p::Morseold, r) = p.e0 * (exp(-2*p.A*(r/p.r0-1.0)) - 2.0*exp(-p.A*(r/p.r0-1.0)))
evaluate_d(p::Morseold, r) = -2.0*p.e0*p.A/p.r0*(exp(-2*p.A*(r/p.r0-1.0)) - exp(-p.A*(r/p.r0-1.0)))

const e0 = 0.99
const r0 = 1.05
const A = 4.1

lj = JuLIP.Potentials.AnalyticPotential(
                     "$e0 * (exp(-2*$A*(r/$r0-1.0)) - 2.0*exp(-$A*(r/$r0-1.0)))",
                     id="Morse(e0=$e0,r0=$r0,A=$A)")
@show typeof(lj)
ljold = Morseold(e0, A, r0)

rr = collect(linspace(0.9, 2.1, 100))
lj_r = [lj(r) for r in rr]
ljold_r = ljold(rr)
dlj_r = [(@D lj(r)) for r in rr]
dljold_r = @D ljold(rr)

@show vecnorm(lj_r - ljold_r, Inf)
@show vecnorm(dlj_r - dljold_r, Inf)

function test(N, p)
   r = 1.234
   s = 0.0
   for n = 1:N
      s += p(r)
   end
   return s
end

function testgrad(N, p)
   r = 1.234
   s = 0.0
   for n = 1:N
      s += @D p(r)
   end
   return s
end

test(10, lj);
test(10, ljold);
testgrad(10, lj);
testgrad(10, ljold);

println("1M evaluations of Analytic Potential")
@time test(1_000_000, lj)
println("1M evaluations of old Potential")
@time test(1_000_000, ljold)
println("1M ∇ evaluations of Analytic Potential")
@time testgrad(1_000_000, lj)
println("1M ∇ evaluations of old Potential")
@time testgrad(1_000_000, ljold)


text:
fish = parse.(Int, split(readline("day6/fish.txt"), ','))

function forecast_fish(init, days)
    curr = merge(
        Dict((i, 0) for i in 0:8),
        Dict((i, count(==(i), init)) for i in unique(init))
    )
    for _ in 1:days
        prev = copy(curr)
        for i in 8:-1:1
            curr[i - 1] = prev[i]
        end
        curr[6] += prev[0]
        curr[8] = prev[0]
    end
    curr
end

p1 = sum(values(forecast_fish(fish, 80)))

@show p1

p2 = sum(values(forecast_fish(fish, 256)))

@show p2


text:
import DataFrames
import SQLite

Base.eval(:(have_color = true))
include("validate_alone.jl")
#include("move.jl")

#main for standard and mini
function standard_mini(board, sourcex, sourcey, targetx, targety, option)
    #checks if is in valid line for promotion
    #println(option)
    if option == "!" || option == '!'
       # println("entered option")
        if !promotion_check(board, sourcex, sourcey, targetx, targety)
            return false
        end
    end

    #println("starts looking for name 2 : ", board[sourcey][sourcex].name)
    if board[sourcey][sourcex].name == "king"
           #println("enters 0")
        if king(board, sourcex, sourcey, targetx, targety)
            return true
        end
    elseif board[sourcey][sourcex].name ==  "rook"
        if rook(board, sourcex, sourcey, targetx, targety)
            return true
        end
    elseif board[sourcey][sourcex].name ==  "bishop"
        if bishop(board, sourcex, sourcey, targetx, targety)
            return true
        end
    elseif board[sourcey][sourcex].name ==  "dragon"
        if dragon_king(board, sourcex, sourcey, targetx, targety)
            return true
        end
    elseif board[sourcey][sourcex].name ==  "horse"
        if dragon_horse(board, sourcex, sourcey, targetx, targety)
            return true
        end
    elseif board[sourcey][sourcex].name == "gold"
        #if move number (aka i) is odd then black
        if board[sourcey][sourcex].side == 1
            if gold_black(board, sourcex, sourcey, targetx, targety)
                return true
            end
        else
            if gold_white(board, sourcex, sourcey, targetx, targety)
                return true
            end
        end
    elseif board[sourcey][sourcex].name ==  "silver"
        if board[sourcey][sourcex].side == 1
            if silver_black(board, sourcex, sourcey, targetx, targety)
                return true
            end
        else
            if silver_white(board, sourcex, sourcey, targetx, targety)
                return true
            end
        end
    elseif board[sourcey][sourcex].name ==  "knight"
        if board[sourcey][sourcex].side == 1
            if knight_black(board, sourcex, sourcey, targetx, targety)
                return true
            end
        else
            if knight_white(board, sourcex, sourcey, targetx, targety)
                return true
            end
        end
    elseif board[sourcey][sourcex].name ==  "lance"
        if board[sourcey][sourcex].side == 1
            if lance_black(board, sourcex, sourcey, targetx, targety)
                return true
            end
        else
            if lance_white(board, sourcex, sourcey, targetx, targety)
                return true
            end
        end
    elseif board[sourcey][sourcex].name ==  "pawn"
        #println("enters 1")
        if board[sourcey][sourcex].side == 1
                    #println("enters 2")
            if pawn_black(board, sourcex, sourcey, targetx, targety)
                return true
            end
        else
            if pawn_white(board, sourcex, sourcey, targetx, targety)
                return true
            end
        end
    else
        return false
    end
    return false
end

function dropping(board, targetx, targety, option, side)
    boardLength = length(board)

    if targetx > boardLength || targety > boardLength || targetx < 1|| targety < 1
        return false
    end

    #check that they have the piece that is being droppped (call before actually dropping)
    if !before_drop_check(option, side, capturedByOdd, capturedByEven)
        return false
    end

    if option != "" && option != "!"
        if drop_check(board, option, targetx, targety, side, capturedByOdd, capturedByEven)
            return true
        else
            return false
        end
    else
        return false
    end
end

#main for chu
function chu(board, sourcex, sourcey, targetx, targety, targetx2, targety2, option)
    #checks if is in valid line for promotion
    if option == "!" || option == '!'
        if targetx2 != -1 && targety2 != -1
            if !(promotion_check(board, sourcex, sourcey, targetx, targety) || promotion_check(board, targetx, targety, targetx2, targety2))
                return false
            end
        else
            if !promotion_check(board, sourcex, sourcey, targetx, targety)
                return false
            end
        end
    end

    if board[sourcey][sourcex].name == "cobra"
        if go_between(board, sourcex, sourcey, targetx, targety)
            return true
        end
    elseif board[sourcey][sourcex].name == "elephant"
        if board[sourcey][sourcex].side == 1
            if elephant_black(board, sourcex, sourcey, targetx, targety)
                return true
            end
        else
            if elephant_white(board, sourcex, sourcey, targetx, targety)
                return true
            end
        end
    elseif board[sourcey][sourcex].name ==  "pawn"
        if board[sourcey][sourcex].side == 1
            if pawn_black(board, sourcex, sourcey, targetx, targety)
                return true
            end
        else
            if pawn_white(board, sourcex, sourcey, targetx, targety)
                return true
            end
        end
    elseif board[sourcey][sourcex].name == "gold"
        if board[sourcey][sourcex].side == 1
            if gold_black(board, sourcex, sourcey, targetx, targety)
                return true
            end
        else
            if gold_white(board, sourcex, sourcey, targetx, targety)
                return true
            end
        end
    elseif board[sourcey][sourcex].name == "smover"
        if side_mover(board, sourcex, sourcey, targetx, targety)
            return true
        end
    elseif board[sourcey][sourcex].name == "boar"
        #promoted side mover is free boar
        if p_side_mover(board, sourcex, sourcey, targetx, targety)
            return true
        end
    elseif board[sourcey][sourcex].name == "vmover"
        if vertical_mover(board, sourcex, sourcey, targetx, targety)
            return true
        end
    elseif board[sourcey][sourcex].name == "ox"
        #promoted vertical mover is flying ox
        if p_vertical_mover(board, sourcex, sourcey, targetx, targety)
            return true
        end
    elseif board[sourcey][sourcex].name ==  "bishop"
        if bishop(board, sourcex, sourcey, targetx, targety)
            return true
        end
    elseif board[sourcey][sourcex].name ==  "horse"
        if dragon_horse(board, sourcex, sourcey, targetx, targety)
            return true
        end
    elseif board[sourcey][sourcex].name ==  "rook"
        if rook(board, sourcex, sourcey, targetx, targety)
            return true
        end
    elseif board[sourcey][sourcex].name ==  "dragon"
        if dragon_king(board, sourcex, sourcey, targetx, targety)
            return true
        end
    elseif board[sourcey][sourcex].name == "falcon"
        if board[sourcey][sourcex].side == 1
            if falcon_black(board, sourcex, sourcey, targetx, targety, targetx2, targety2)
                return true
            end
        else
            if falcon_white(board, sourcex, sourcey, targetx, targety, targetx2, targety2)
                return true
            end
        end
    elseif board[sourcey][sourcex].name == "soaring"
        if board[sourcey][sourcex].side  == 1
            if soaring_black(board, sourcex, sourcey, targetx, targety, targetx2, targety2)
                return true
            end
        else
            if soaring_white(board, sourcex, sourcey, targetx, targety, targetx2, targety2)
                return true
            end
        end
    elseif board[sourcey][sourcex].name ==  "lance"
        if board[sourcey][sourcex].side  == 1
            if lance_black(board, sourcex, sourcey, targetx, targety)
                return true
            end
        else
            if lance_white(board, sourcex, sourcey, targetx, targety)
                return true
            end
        end
    elseif board[sourcey][sourcex].name ==  "white"
        #white horse == promoted lance
        if board[sourcey][sourcex].side == 1
            if p_lance_black(board, sourcex, sourcey, targetx, targety)
                return true
            end
        else
            if p_lance_white(board, sourcex, sourcey, targetx, targety)
                return true
            end
        end
    elseif board[sourcey][sourcex].name ==  "chariot"
        if chariot(board, sourcex, sourcey, targetx, targety)
            return true
        end
    elseif board[sourcey][sourcex].name ==  "whale"
        #promoted reverse chariot == whale
        if board[sourcey][sourcex].side == 1
            if p_chariot_black(board, sourcex, sourcey, targetx, targety)
                return true
            end
        else
            if p_chariot_white(board, sourcex, sourcey, targetx, targety)
                return true
            end
        end
    elseif board[sourcey][sourcex].name ==  "tiger"
        if board[sourcey][sourcex].side == 1
            if tiger_black(board, sourcex, sourcey, targetx, targety)
                return true
            end
        else
            if tiger_white(board, sourcex, sourcey, targetx, targety)
                return true
            end
        end
    elseif board[sourcey][sourcex].name ==  "stag"
        #flying stag == promoted blind tiger
        if p_tiger(board, sourcex, sourcey, targetx, targety)
            return true
        end
    elseif board[sourcey][sourcex].name ==  "leopard"
        if leopard(board, sourcex, sourcey, targetx, targety)
            return true
        end
    elseif board[sourcey][sourcex].name ==  "copper"
        if board[sourcey][sourcex].side == 1
            if copper_black(board, sourcex, sourcey, targetx, targety)
                return true
            end
        else
            if copper_white(board, sourcex, sourcey, targetx, targety)
                return true
            end
        end
    elseif board[sourcey][sourcex].name ==  "silver"
        if board[sourcey][sourcex].side == 1
            if silver_black(board, sourcex, sourcey, targetx, targety)
                return true
            end
        else
            if silver_white(board, sourcex, sourcey, targetx, targety)
                return true
            end
        end
    elseif board[sourcey][sourcex].name ==  "king"
        #promoted drunk elephant == king == prince
        if king(board, sourcex, sourcey, targetx, targety)
            return true
        end
    elseif board[sourcey][sourcex].name ==  "kirin"
        #n in chu shogi is kirin
        if kirin(board, sourcex, sourcey, targetx, targety)
            return true
        end
    elseif board[sourcey][sourcex].name ==  "lion"
        #promoted kirin == lion
        if lion(board, sourcex, sourcey, targetx, targety, targetx2, targety2)
            return true
        end
    elseif board[sourcey][sourcex].name ==  "phoenix"
        if phoenix(board, sourcex, sourcey, targetx, targety)
            return true
        end
    elseif board[sourcey][sourcex].name ==  "queen"
        if queen(board, sourcex, sourcey, targetx, targety)
            return true
        end
    else
        return false
    end
    return false
end

function ten(board, sourcex, sourcey, targetx, targety, targetx2, targety2, targetx3, targety3, option)
       #checks if is in valid line for promotion
    if option == "!" || option == '!'
        if targetx2 != -1 && targety2 != -1  && targetx3 != -1 && targety3 != -1
            if !(promotion_check(board, sourcex, sourcey, targetx, targety) || promotion_check(board, targetx, targety, targetx2, targety2) || promotion_check(board, targetx2, targety2, targetx3, targety3))
                return false
            end
        elseif targetx2 != -1 && targety2 != -1
            if !(promotion_check(board, sourcex, sourcey, targetx, targety) || promotion_check(board, targetx, targety, targetx2, targety2))
                return false
            end
        else
            if !promotion_check(board, sourcex, sourcey, targetx, targety)
                return false
            end
        end
    end

    if board[sourcey][sourcex].name == "elephant"
        if board[sourcey][sourcex].side == 1
            if elephant_black(board, sourcex, sourcey, targetx, targety)
                return true
            end
        else
            if elephant_white(board, sourcex, sourcey, targetx, targety)
                return true
            end
        end
    elseif board[sourcey][sourcex].name ==  "pawn"
        if board[sourcey][sourcex].side == 1
            if pawn_black(board, sourcex, sourcey, targetx, targety)
                return true
            end
        else
            if pawn_white(board, sourcex, sourcey, targetx, targety)
                return true
            end
        end
    elseif board[sourcey][sourcex].name == "gold"
        if board[sourcey][sourcex].side == 1
            if gold_black(board, sourcex, sourcey, targetx, targety)
                return true
            end
        else
            if gold_white(board, sourcex, sourcey, targetx, targety)
                return true
            end
        end
    elseif board[sourcey][sourcex].name == "smover"
        if side_mover(board, sourcex, sourcey, targetx, targety)
            return true
        end
    elseif board[sourcey][sourcex].name == "boar"
        #promoted side mover is free boar
        if p_side_mover(board, sourcex, sourcey, targetx, targety)
            return true
        end
    elseif board[sourcey][sourcex].name == "vmover"
        if vertical_mover(board, sourcex, sourcey, targetx, targety)
            return true
        end
    elseif board[sourcey][sourcex].name == "ox"
        #promoted vertical mover is flying ox
        if p_vertical_mover(board, sourcex, sourcey, targetx, targety)
            return true
        end
    elseif board[sourcey][sourcex].name ==  "bishop"
        if bishop(board, sourcex, sourcey, targetx, targety)
            return true
        end
    elseif board[sourcey][sourcex].name ==  "horse"
        if dragon_horse(board, sourcex, sourcey, targetx, targety)
            return true
        end
    elseif board[sourcey][sourcex].name ==  "rook"
        if rook(board, sourcex, sourcey, targetx, targety)
            return true
        end
    elseif board[sourcey][sourcex].name ==  "dragon"
        if dragon_king(board, sourcex, sourcey, targetx, targety)
            return true
        end
    elseif board[sourcey][sourcex].name == "falcon"
        if board[sourcey][sourcex].side == 1
            if falcon_black(board, sourcex, sourcey, targetx, targety, targetx2, targety2)
                return true
            end
        else
            if falcon_white(board, sourcex, sourcey, targetx, targety, targetx2, targety2)
                return true
            end
        end
    elseif board[sourcey][sourcex].name == "soaring"
        if board[sourcey][sourcex].side  == 1
            if soaring_black(board, sourcex, sourcey, targetx, targety, targetx2, targety2)
                return true
            end
        else
            if soaring_white(board, sourcex, sourcey, targetx, targety, targetx2, targety2)
                return true
            end
        end
    elseif board[sourcey][sourcex].name ==  "lance"
        if board[sourcey][sourcex].side  == 1
            if lance_black(board, sourcex, sourcey, targetx, targety)
                return true
            end
        else
            if lance_white(board, sourcex, sourcey, targetx, targety)
                return true
            end
        end
    elseif board[sourcey][sourcex].name ==  "white"
        #white horse == promoted lance
        if board[sourcey][sourcex].side == 1
            if p_lance_black(board, sourcex, sourcey, targetx, targety)
                return true
            end
        else
            if p_lance_white(board, sourcex, sourcey, targetx, targety)
                return true
            end
        end
    elseif board[sourcey][sourcex].name ==  "chariot"
        if chariot(board, sourcex, sourcey, targetx, targety)
            return true
        end
    elseif board[sourcey][sourcex].name ==  "whale"
        #promoted reverse chariot == whale
        if board[sourcey][sourcex].side == 1
            if p_chariot_black(board, sourcex, sourcey, targetx, targety)
                return true
            end
        else
            if p_chariot_white(board, sourcex, sourcey, targetx, targety)
                return true
            end
        end
    elseif board[sourcey][sourcex].name ==  "tiger"
        if board[sourcey][sourcex].side == 1
            if tiger_black(board, sourcex, sourcey, targetx, targety)
                return true
            end
        else
            if tiger_white(board, sourcex, sourcey, targetx, targety)
                return true
            end
        end
    elseif board[sourcey][sourcex].name ==  "stag"
        #flying stag == promoted blind tiger
        if p_tiger(board, sourcex, sourcey, targetx, targety)
            return true
        end
    elseif board[sourcey][sourcex].name ==  "leopard"
        if leopard(board, sourcex, sourcey, targetx, targety)
            return true
        end
    elseif board[sourcey][sourcex].name ==  "copper"
        if board[sourcey][sourcex].side == 1
            if copper_black(board, sourcex, sourcey, targetx, targety)
                return true
            end
        else
            if copper_white(board, sourcex, sourcey, targetx, targety)
                return true
            end
        end
    elseif board[sourcey][sourcex].name ==  "silver"
        if board[sourcey][sourcex].side == 1
            if silver_black(board, sourcex, sourcey, targetx, targety)
                return true
            end
        else
            if silver_white(board, sourcex, sourcey, targetx, targety)
                return true
            end
        end
    elseif board[sourcey][sourcex].name ==  "king"
        #promoted drunk elephant == king == prince
        if king(board, sourcex, sourcey, targetx, targety)
            return true
        end
    elseif board[sourcey][sourcex].name ==  "kirin"
        #n in chu shogi is kirin
        if kirin(board, sourcex, sourcey, targetx, targety)
            return true
        end
    elseif board[sourcey][sourcex].name ==  "lion"
        #promoted kirin == lion
        if lion(board, sourcex, sourcey, targetx, targety, targetx2, targety2)
            return true
        end
    elseif board[sourcey][sourcex].name ==  "phoenix"
        if phoenix(board, sourcex, sourcey, targetx, targety)
            return true
        end
    elseif board[sourcey][sourcex].name ==  "queen"
        if queen(board, sourcex, sourcey, targetx, targety)
            return true
        end
    elseif board[sourcey][sourcex].name ==  "knight"
        if board[sourcey][sourcex].side == 1
            if knight_black(board, sourcex, sourcey, targetx, targety)
                return true
            end
        else
            if knight_white(board, sourcex, sourcey, targetx, targety)
                return true
            end
        end
    elseif board[sourcey][sourcex].name ==  "prince"
        if prince(board, sourcex, sourcey, targetx, targety)
            return true
        end
    elseif board[sourcey][sourcex].name ==  "vice"
        if vice(board, sourcex, sourcey, targetx, targety, targetx2, targety2, targetx3, targety3)
            return true
        end
    elseif board[sourcey][sourcex].name ==  "great"
        if great(board, sourcex, sourcey, targetx, targety)
            return true
        end
    elseif board[sourcey][sourcex].name ==  "bgeneral"
        if bgeneral(board, sourcex, sourcey, targetx, targety)
            return true
        end
    elseif board[sourcey][sourcex].name ==  "rgeneral"
        if rgeneral(board, sourcex, sourcey, targetx, targety)
            return true
        end
    elseif board[sourcey][sourcex].name ==  "demon"
        if demon(board, sourcex, sourcey, targetx, targety, targetx2, targety2, targetx3, targety3)
            return true
        end
    elseif board[sourcey][sourcex].name ==  "tetrarch"
        if tetrarch(board, sourcex, sourcey, targetx, targety, targetx2, targety2, targetx3, targety3)
            return true
        end
    elseif board[sourcey][sourcex].name ==  "buffalo"
        if buffalo(board, sourcex, sourcey, targetx, targety)
            return true
        end
    elseif board[sourcey][sourcex].name ==  "csoldier"
        if csoldier(board, sourcex, sourcey, targetx, targety)
            return true
        end
    elseif board[sourcey][sourcex].name ==  "ssoldier"
        if board[sourcey][sourcex].side == 1
            if ssoldier_black(board, sourcex, sourcey, targetx, targety)
                return true
            end
        else
            if ssoldier_white(board, sourcex, sourcey, targetx, targety)
                return true
            end
        end
    elseif board[sourcey][sourcex].name ==  "vsoldier"
        if board[sourcey][sourcex].side == 1
            if vsoldier_black(board, sourcex, sourcey, targetx, targety)
                return true
            end
        else
            if vsoldier_white(board, sourcex, sourcey, targetx, targety)
                return true
            end
        end
    elseif board[sourcey][sourcex].name ==  "iron"
        if board[sourcey][sourcex].side == 1
            if iron_black(board, sourcex, sourcey, targetx, targety)
                return true
            end
        else
            if iron_white(board, sourcex, sourcey, targetx, targety)
                return true
            end
        end
    elseif board[sourcey][sourcex].name ==  "eagle"
        if eagle(board, sourcex, sourcey, targetx, targety, targetx2, targety2)
           return true
        end
    elseif board[sourcey][sourcex].name ==  "hawk"
        if hawk(board, sourcex, sourcey, targetx, targety, targetx2, targety2)
            return true
        end
    elseif board[sourcey][sourcex].name ==  "multi"
        if board[sourcey][sourcex].side == 1
            if multi_black(board, sourcex, sourcey, targetx, targety)
                return true
            end
        else
            if multi_white(board, sourcex, sourcey, targetx, targety)
                return true
            end
        end
    elseif board[sourcey][sourcex].name ==  "dog"
        if board[sourcey][sourcex].side == 1
            if dog_black(board, sourcex, sourcey, targetx, targety)
                return true
            end
        else
            if dog_white(board, sourcex, sourcey, targetx, targety)
                return true
            end
        end
    else
        return false
    end
    #burning(board, board[sourcey][sourcex].side)
    return false
end

function validate(board, amove, capturedByOdd, capturedByEven)
    println("move number: $(amove.move_number)") # for debugging

    boardLength = length(board)

    if boardLength == 5
        gameType = "minishogi"
    elseif boardLength == 9
        gameType = "standard"
    elseif boardLength == 12
        gameType = "chu"
    else
        gameType = "ten"
    end

    i = amove.move_number
    moveType = amove.move_type
    sourcex = amove.sourcex
    sourcey = amove.sourcey
    targetx = amove.targetx
    targety = amove.targety
    option = amove.option
    targetx2 = amove.targetx2
    targety2 = amove.targety2
    targetx3 = amove.targetx3
    targety3 = amove.targety3

    side = i % 2

    if moveType == "resign"
        return true
    end

    if moveType == "move"
        #println("it is a move")
        if board[sourcey][sourcex].side == 1 && !isodd(i)
            return false
        end
        if board[sourcey][sourcex].side == 0 && !iseven(i)
            return false
        end
        #check if the target is valid or if it is actually moving
        if !target_valid(board, targetx, targety, board[sourcey][sourcex].side) || same_coordinates(sourcex, sourcey, targetx, targety) || !valid_boundaries(board, sourcex, sourcey, targetx, targety)
                            #println("generals")
            return false
        end

        if gameType == "chu"
            if chu(board, sourcex, sourcey, targetx, targety, targetx2, targety2, option)
                return true
            end
        elseif gameType == "ten"
            if ten(board, sourcex, sourcey, targetx, targety, targetx2, targety2, targetx3, targety3, option)
                return true
            end
        else
                #println("stand mini 1")
                #println("starts looking for name: ", board[sourcey][sourcex].name)
            if standard_mini(board, sourcex, sourcey, targetx, targety, option)
               # println("stand mini 2")
                return true
            end
        end

    elseif moveType == "drop"
        if gameType == "chu" || gameType == "ten"
            return false
        end
        #arrays?
        if dropping(board, targetx, targety, option, board[sourcey][sourcex].side, capturedByOdd, capturedByEven)
            return true
        end
    else
        return false
    end
    return false
end

##############
#### MAIN ####
##############

#=#sets initial board and boardLength

gameFile = SQLite.DB(ARGS[1])
metaTable, movesTable, gameType, seed, is_legal, is_timed = gameSetup(gameFile)
oldBoard, boardLength, numOfPromotionRanks = boardSetup(gameType)
board, capturedByEven, capturedByOdd, gameActive = updateBoard(oldBoard, movesTable)

amove = move(1,"move", 8, 1, 9, 1, "", false, -1, -1, -1, -1)

if validate(board, amove, capturedByOdd, capturedByEven)
    println("valid move") #debugging
else
    println("cheater") #debugging
end
=#

#=
println(elephant_black(board, sourcex, sourcey, targetx, targety))
println(elephant_white(board, sourcex, sourcey, targetx, targety))
println(pawn_black(board, sourcex, sourcey, targetx, targety))
println(pawn_white(board, sourcex, sourcey, targetx, targety))
println(gold_black(board, sourcex, sourcey, targetx, targety))
println(gold_white(board, sourcex, sourcey, targetx, targety))
println(side_mover(board, sourcex, sourcey, targetx, targety))
println(p_side_mover(board, sourcex, sourcey, targetx, targety))
println(vertical_mover(board, sourcex, sourcey, targetx, targety))
println(p_vertical_mover(board, sourcex, sourcey, targetx, targety))
println(bishop(board, sourcex, sourcey, targetx, targety))
println(dragon_horse(board, sourcex, sourcey, targetx, targety))
println(rook(board, sourcex, sourcey, targetx, targety))
println(dragon_king(board, sourcex, sourcey, targetx, targety))
println(falcon_black(board, sourcex, sourcey, targetx, targety, targetx2, targety2))
println(falcon_white(board, sourcex, sourcey, targetx, targety, targetx2, targety2))
println(soaring_black(board, sourcex, sourcey, targetx, targety, targetx2, targety2))
println(soaring_white(board, sourcex, sourcey, targetx, targety, targetx2, targety2))
println(lance_black(board, sourcex, sourcey, targetx, targety))
println(lance_white(board, sourcex, sourcey, targetx, targety))
println(p_lance_black(board, sourcex, sourcey, targetx, targety))
println(p_lance_white(board, sourcex, sourcey, targetx, targety))
println(chariot(board, sourcex, sourcey, targetx, targety))
println(p_chariot_black(board, sourcex, sourcey, targetx, targety))
println(p_chariot_white(board, sourcex, sourcey, targetx, targety))
println(tiger_black(board, sourcex, sourcey, targetx, targety))
println(tiger_white(board, sourcex, sourcey, targetx, targety))
println(p_tiger(board, sourcex, sourcey, targetx, targety))
println(leopard(board, sourcex, sourcey, targetx, targety))
println(copper_black(board, sourcex, sourcey, targetx, targety))
println(copper_white(board, sourcex, sourcey, targetx, targety))
println(silver_black(board, sourcex, sourcey, targetx, targety))
println(silver_white(board, sourcex, sourcey, targetx, targety))
println(king(board, sourcex, sourcey, targetx, targety))
println(kirin(board, sourcex, sourcey, targetx, targety))
println(lion(board, sourcex, sourcey, targetx, targety, targetx2, targety2))
println(phoenix(board, sourcex, sourcey, targetx, targety))
println(queen(board, sourcex, sourcey, targetx, targety))
println(knight_black(board, sourcex, sourcey, targetx, targety))
println(knight_white(board, sourcex, sourcey, targetx, targety))
println(prince(board, sourcex, sourcey, targetx, targety))
println(vice(board, sourcex, sourcey, targetx, targety, targetx2, targety2, targetx3, targety3))
println(great(board, sourcex, sourcey, targetx, targety))
println(bgeneral(board, sourcex, sourcey, targetx, targety))
println(rgeneral(board, sourcex, sourcey, targetx, targety))
println(demon(board, sourcex, sourcey, targetx, targety, targetx2, targety2, targetx3, targety3))
println(tetrarch(board, sourcex, sourcey, targetx, targety, targetx2, targety2, targetx3, targety3))
println(buffalo(board, sourcex, sourcey, targetx, targety))
println(csoldier(board, sourcex, sourcey, targetx, targety))
println(ssoldier_black(board, sourcex, sourcey, targetx, targety))
println(ssoldier_white(board, sourcex, sourcey, targetx, targety))
println(vsoldier_black(board, sourcex, sourcey, targetx, targety))
println(vsoldier_white(board, sourcex, sourcey, targetx, targety))
println(iron_black(board, sourcex, sourcey, targetx, targety))
println(iron_white(board, sourcex, sourcey, targetx, targety))
println(eagle(board, sourcex, sourcey, targetx, targety, targetx2, targety2))
println(hawk(board, sourcex, sourcey, targetx, targety, targetx2, targety2))
println(multi_black(board, sourcex, sourcey, targetx, targety))
println(multi_white(board, sourcex, sourcey, targetx, targety))
println(dog_black(board, sourcex, sourcey, targetx, targety))
println(dog_white(board, sourcex, sourcey, targetx, targety))

    return false=#

max_stars_repo_name: DarienFlamont/Shogi

max_issues_repo_name: DarienFlamont/Shogi

max_forks_repo_name: DarienFlamont/Shogi


text:
module SplitApplyCombine

using Base: @propagate_inbounds, @pure, promote_op
using Indexing

# Syntax
export @_

# collections -> scalar
export only

# collections -> collections
import Base: merge, merge!
export mapmany

# collections -> collections of collections
export group, groupinds, Groups, groupview, groupreduce, splitdims, splitdimsview, SplitDimsArray

# colletions of collections -> collections
export flatten #, flattenview

# collections of collections -> collections of collections
export innerjoin, ⨝, leftgroupjoin
# `pivot` or similar - like transpose, but nested.

include("only.jl")
include("merge.jl")
include("map.jl")
include("group.jl")
include("innerjoin.jl")
include("leftgroupjoin.jl")
include("splitdims.jl")

# Syntax
include("underscore.jl")

# Silly definitions missing from Base
# ===================================
# this should always work
Base.haskey(a, i) = i ∈ keys(a)

# mini-compat (more for my knowledge than anything)
if VERSION < v"0.7-"
    Base.keys(v::AbstractVector) = indices(v)[1]
    Base.keys(a::AbstractArray) = CartesianRange(indices(a)...)
    Base.keys(::NTuple{N,Any}) where {N} = Base.OneTo(N)
    Base.keys(::Number) = Base.OneTo(1)

    # A Nullable is a container with 0 or 1 values... so...
    Base.start(::Nullable) = false
    Base.done(n::Nullable, i::Bool) = isnull(n) | i
    Base.next(n::Nullable, i::Bool) = (n.valie, true)
    Base.first(n::Nullable) = get(n)
    Base.last(n::Nullable) = get(n)
    @propagate_inbounds function Base.getindex(n::Nullable)
    @boundscheck if !n.hasvalue
        return NullException()
    end
    return n.value
end

end

end # module

# Random thoughts:

# merge! many other mutating ops returns collection, so:
# why does setindex! return the set value(s), not the collection???


text:
# Input/Output support for fasta/fastq/bam/sam/bed and other formats
# IO functions are provided as core featuers of OpenGene, so it is not inside a sub-module
# ===========

export opengene_open,
	fastq_open,
	fastq_read,
	fastq_write,
	fasta_open,
	fasta_read,
	fasta_write,
	FastqPairStream,
	fastq_open_pair,
	fastq_flush_pair,
	fastq_close_pair,
	fastq_read_pair,
	fastq_write_pair,
	bed_read_intervals,
	bed_write_intervals

export vcf_read,
	vcf_write,
	gtf_read,
	gtf_read_row,
	gtf_read_row!,
	gtf_write

include("common.jl")

# for every file format, we create a folder
include("fastq/fastq.jl")
include("fasta/fasta.jl")
include("beds/bed.jl")
include("vcf/vcf.jl")
include("gtf/gtf.jl")


text:
# automatically generated by the FlatBuffers compiler, do not modify

include("namespace_test1_generated.jl")
if !isdefined(@__MODULE__(), :NamespaceA) @__MODULE__().eval(:(module NamespaceA __precompile__(false); import FlatBuffers end)) end
if !isdefined(NamespaceA, :NamespaceB) NamespaceA.eval(:(module NamespaceB __precompile__(false); import FlatBuffers end)) end
if !isdefined(@__MODULE__(), :NamespaceC) @__MODULE__().eval(:(module NamespaceC __precompile__(false); import FlatBuffers end)) end
include("NamespaceA/TableInFirstNS.jl")
include("NamespaceA/SecondTableInA.jl")

max_stars_repo_name: rjkat/flatbuffers-julia

max_issues_repo_name: rjkat/flatbuffers-julia

max_forks_repo_name: rjkat/flatbuffers-julia


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl
=#

insts = OrderedDict()
insts["minimal"] = [
    # rootdet
    ((6, MatNegGeom()),),
    ((6, MatNegGeomEFExp()),),
    ((6, MatNegGeomEFPow()),),
    # tr inv
    ((4, MatInv()),),
    ((4, MatInvEigOrd()),),
    ((4, MatInvDirect()),),
    # tr neglog
    ((4, MatNegLog()),),
    ((4, MatNegLogEigOrd()),),
    ((4, MatNegLogDirect()),),
    # tr negentropy
    ((4, MatNegEntropy()),),
    ((4, MatNegEntropyEigOrd()),),
    # tr power12
    ((5, MatPower12(1.5)),),
    ((5, MatPower12EigOrd(1.5)),),
    # tr neg2sqrt
    ((4, MatNeg2Sqrt()),),
    ((4, MatNeg2SqrtEigOrd()),),
    ]
insts["fast"] = [
    ((100, MatNegGeom()),),
    ((25, MatNegGeomEFExp()),),
    ((25, MatNegGeomEFPow()),),
    ((80, MatInv()),),
    ((12, MatInvEigOrd()),),
    ((50, MatInvDirect()),),
    ((150, MatNegLog()),),
    ((14, MatNegLogEigOrd()),),
    ((30, MatNegLogDirect()),),
    ((120, MatNegEntropy()),),
    ((15, MatNegEntropyEigOrd()),),
    ((100, MatPower12(1.5)),),
    ((10, MatPower12EigOrd(1.5)),),
    ]
insts["various"] = [
    ((800, MatNegGeom()),),
    ((140, MatNegGeomEFExp()),),
    ((80, MatNegGeomEFPow()), nothing, (default_tol_relax = 100,)),
    ((400, MatInv()),),
    ((20, MatInvEigOrd()),),
    ((250, MatInvDirect()),),
    ((300, MatNegLog()),),
    ((24, MatNegLogEigOrd()),),
    ((150, MatNegLogDirect()),),
    ((300, MatNegEntropy()),),
    ((20, MatNegEntropyEigOrd()),),
    ((30, MatPower12(1.5)),),
    ((18, MatPower12EigOrd(1.5)), nothing, (default_tol_relax = 100,)),
    ]
return (ExperimentDesignJuMP, insts)


text:
@nospecialize

"""
    scatter(x,y)
    scatter!(x,y)

Make a scatter plot of y vs x.

# Examples
```julia-repl
julia> scatter([1,2,3],[4,5,6],markersize=[3,4,5],markercolor=[:red,:green,:blue])
julia> scatter([(1,4),(2,5),(3,6)])
```
"""
@shorthands scatter

"""
    bar(x,y)
    bar!(x,y)

Make a bar plot of y vs x.

# Arguments

- $(_document_argument("bar_position"))
- $(_document_argument("bar_width"))
- $(_document_argument("bar_edges"))
- $(_document_argument("orientation"))

# Examples
```julia-repl
julia> bar([1,2,3],[4,5,6],fillcolor=[:red,:green,:blue],fillalpha=[0.2,0.4,0.6])
julia> bar([(1,4),(2,5),(3,6)])
```
"""
@shorthands bar

@shorthands barh

"""
    histogram(x)
    histogram!(x)

Plot a histogram.

# Arguments

- `x`: AbstractVector of values to be binned
- $(_document_argument("bins"))
- `weights`: Vector of weights for the values in `x`, for weighted bin counts
- $(_document_argument("normalize"))
- $(_document_argument("bar_position"))
- $(_document_argument("bar_width"))
- $(_document_argument("bar_edges"))
- $(_document_argument("orientation"))

# Example
```julia-repl
julia> histogram([1,2,1,1,4,3,8],bins=0:8)
julia> histogram([1,2,1,1,4,3,8],bins=0:8,weights=weights([4,7,3,9,12,2,6]))
```
"""
@shorthands histogram

"""
    barhist(x)
    barhist!(x)

Make a histogram bar plot. See `histogram`.
"""
@shorthands barhist

"""
    stephist(x)
    stephist!(x)

Make a histogram step plot (bin counts are represented using horizontal lines
instead of bars). See `histogram`.
"""
@shorthands stephist

"""
    scatterhist(x)
    scatterhist!(x)

Make a histogram scatter plot (bin counts are represented using points
instead of bars). See `histogram`.
"""
@shorthands scatterhist

"""
    histogram2d(x,y)
    histogram2d!(x,y)

Plot a two-dimensional histogram.

# Arguments

- `bins`: Number of bins (if an `Integer`) or bin edges (if an `AbtractVector`)
- `weights`: Vector of weights for the values in `x`. Each entry of x contributes
             its weight to the height of its bin.

# Example
```julia-repl
julia> histogram2d(randn(10_000),randn(10_000))
```
"""
@shorthands histogram2d

"""
    density(x)
    density!(x)

Make a line plot of a kernel density estimate of x.

# Arguments

- `x`: AbstractVector of samples for probability density estimation

# Example
```julia-repl
julia> using StatsPlots
julia> density(randn(100_000))
```
"""
@shorthands density

"""
    heatmap(x,y,z)
    heatmap!(x,y,z)

Plot a heatmap of the rectangular array `z`.

# Example
```julia-repl
julia> heatmap(randn(10,10))
```
"""
@shorthands heatmap
@shorthands plots_heatmap

"""
    hexbin(x,y)
    hexbin!(x,y)

Make a hexagonal binning plot (a histogram of the observations `(x[i],y[i])`
with hexagonal bins)

# Example
```julia-repl
julia> hexbin(randn(10_000), randn(10_000))
```
"""
@shorthands hexbin

"""
    sticks(x,y)
    sticks!(x,y)

Draw a stick plot of y vs x.

# Example
```julia-repl
julia> sticks(1:10)
```
"""
@shorthands sticks

"""
    hline(y)
    hline!(y)

Draw horizontal lines at positions specified by the values in
the AbstractVector `y`

# Example
```julia-repl
julia> hline([-1,0,2])
```
"""
@shorthands hline

"""
    vline(x)
    vline!(x)

Draw vertical lines at positions specified by the values in
the AbstractVector `x`

# Example
```julia-repl
julia> vline([-1,0,2])
```
"""
@shorthands vline

"""
    hspan(y)

Draw a rectangle between the horizontal line at position `y[1]`
and the horizontal line at position `y[2]`. If `length(y) ≥ 4`,
then further rectangles are drawn between `y[3]` and `y[4]`,
`y[5]` and `y[6]`, and so on. If `length(y)` is odd, then the
last entry of `y` is ignored.
# Example
```julia-repl
julia> hspan(1:6)
```
"""
@shorthands hspan

"""
    vspan(x)

Draw a rectangle between the vertical line at position `x[1]`
and the vertical line at position `x[2]`. If `length(x) ≥ 4`,
then further rectangles are drawn between `x[3]` and `x[4]`,
`x[5]` and `x[6]`, and so on. If `length(x)` is odd, then the
last entry of `x` is ignored.
# Example
```julia-repl
julia> vspan(1:6)
```
"""
@shorthands vspan

"""
    ohlc(x,y::Vector{OHLC})
    ohlc!(x,y::Vector{OHLC})

Make open-high-low-close plot. Each entry of y is represented by a vertical
segment extending from the low value to the high value, with short horizontal
segments on the left and right indicating the open and close values, respectively.

# Example
```julia-repl
julia> meanprices = cumsum(randn(100))
julia> y = OHLC[(p+rand(),p+1,p-1,p+rand()) for p in meanprices]
julia> ohlc(y)
```
"""
@shorthands ohlc

"""
    contour(x,y,z)
    contour!(x,y,z)

Draw contour lines of the `Surface` z.

# Arguments
- `levels`: Contour levels (if `AbstractVector`) or number of levels (if `Integer`)
- `fill`: Bool. Fill area between contours or draw contours only (false by default)

# Example
```julia-repl
julia> x = y = range(-20, stop = 20, length = 100)
julia> contour(x, y, (x, y) -> x^2 + y^2)
```
"""
@shorthands contour

"An alias for `contour` with fill = true."
@shorthands contourf

@shorthands contour3d

"""
    surface(x,y,z)
    surface!(x,y,z)

Draw a 3D surface plot.

# Example
```julia-repl
julia> using LinearAlgebra
julia> x = y = range(-3, stop = 3, length = 100)
julia> surface(x, y, (x, y) -> sinc(norm([x, y])))
```
"""
@shorthands surface

"""
    wireframe(x,y,z)
    wireframe!(x,y,z)

Draw a 3D wireframe plot.

# Example
```julia-repl
julia> wireframe(1:10,1:10,randn(10,10))
```
"""
@shorthands wireframe

"""
    path3d(x,y,z)
    path3d!(x,y,z)

Plot a 3D path from `(x[1],y[1],z[1])` to `(x[2],y[2],z[2])`,
..., to `(x[end],y[end],z[end])`.

# Example
```julia-repl
julia> path3d([0,1,2,3],[0,1,4,9],[0,1,8,27])
```
"""
@shorthands path3d

"""
    scatter3d(x,y,z)
    scatter3d!(x,y,z)

Make a 3D scatter plot.

# Example
```julia-repl
julia> scatter3d([0,1,2,3],[0,1,4,9],[0,1,8,27])
```
"""
@shorthands scatter3d

"""
    mesh3d(x,y,z)
    mesh3d(x,y,z; connections)

Plot a 3d mesh. On Plotly the triangles can be specified using the connections argument.

# Example
```Julia
x=[0, 1, 2, 0]
y=[0, 0, 1, 2]
z=[0, 2, 0, 1]

i=[0, 0, 0, 1]
j=[1, 2, 3, 2]
k=[2, 3, 1, 3]

plot(x,y,z,seriestype=:mesh3d;connections=(i,j,k))
```
"""
@shorthands mesh3d

"""
    boxplot(x, y)
    boxplot!(x, y)

Make a box and whisker plot.

# Keyword arguments
- `notch`: Bool. Notch the box plot? (false)
- `whisker_range`: Real. Whiskers extend `whisker_range`*IQR below the first quartile
           and above the third quartile. Values outside this range are shown as outliers (1.5)
- `outliers`: Bool. Show outliers? (true)
- `whisker_width`: Real or Symbol. Length of whiskers (:match)

# Example
```julia-repl
julia> using StatsPlots
julia> boxplot(repeat([1,2,3],outer=100),randn(300))
```
"""
@shorthands boxplot

"""
    violin(x,y,z)
    violin!(x,y,z)

Make a violin plot.

# Example
```julia-repl
julia> violin(repeat([1,2,3],outer=100),randn(300))
```
"""
@shorthands violin

"""
    quiver(x,y,quiver=(u,v))
    quiver!(x,y,quiver=(u,v))

Make a quiver (vector field) plot. The `i`th vector extends
from `(x[i],y[i])` to `(x[i] + u[i], y[i] + v[i])`.

# Example
```julia-repl
julia> quiver([1,2,3],[3,2,1],quiver=([1,1,1],[1,2,3]))
```
"""
@shorthands quiver

"""
    curves(x,y)
    curves!(x,y)

Draw a Bezier curve from `(x[1],y[1])` to `(x[end],y[end])`
with control points `(x[2],y[2]), ..., (x[end-1],y[end]-1)`

# Example
```julia-repl
julia> curves([1,2,3,4],[1,1,2,4])
```
"""
@shorthands curves

"Plot a pie diagram"
@shorthands pie

"Plot with seriestype :path3d"
plot3d(args...; kw...) = plot(args...; kw..., seriestype = :path3d)
plot3d!(args...; kw...) = plot!(args...; kw..., seriestype = :path3d)

"Add title to an existing plot"
title!(s::AbstractString; kw...) = plot!(; title = s, kw...)

"Add xlabel to an existing plot"
xlabel!(s::AbstractString; kw...) = plot!(; xlabel = s, kw...)

"Add ylabel to an existing plot"
ylabel!(s::AbstractString; kw...) = plot!(; ylabel = s, kw...)

"Set xlims for an existing plot"
xlims!(lims::Tuple; kw...) = plot!(; xlims = lims, kw...)

"Set ylims for an existing plot"
ylims!(lims::Tuple; kw...) = plot!(; ylims = lims, kw...)

"Set zlims for an existing plot"
zlims!(lims::Tuple; kw...) = plot!(; zlims = lims, kw...)

xlims!(xmin::Real, xmax::Real; kw...) = plot!(; xlims = (xmin, xmax), kw...)
ylims!(ymin::Real, ymax::Real; kw...) = plot!(; ylims = (ymin, ymax), kw...)
zlims!(zmin::Real, zmax::Real; kw...) = plot!(; zlims = (zmin, zmax), kw...)

"Set xticks for an existing plot"
xticks!(v::TicksArgs; kw...) = plot!(; xticks = v, kw...)

"Set yticks for an existing plot"
yticks!(v::TicksArgs; kw...) = plot!(; yticks = v, kw...)

xticks!(ticks::AVec{T}, labels::AVec{S}; kw...) where {T<:Real,S<:AbstractString} =
    plot!(; xticks = (ticks, labels), kw...)
yticks!(ticks::AVec{T}, labels::AVec{S}; kw...) where {T<:Real,S<:AbstractString} =
    plot!(; yticks = (ticks, labels), kw...)

"""
    annotate!(anns...)

Add annotations to an existing plot.

# Arguments

- `anns`: An `AbstractVector` of tuples of the form `(x,y,text)`. The `text` object
          can be a `String`, `PlotText` PlotText (created with `text(args...)`),
          or a tuple of arguments to `text` (e.g., `("Label", 8, :red, :top)`).

# Example
```julia-repl
julia> plot(1:10)
julia> annotate!([(7,3,"(7,3)"),(3,7,text("hey", 14, :left, :top, :green))])
julia> annotate!([(4, 4, ("More text", 8, 45.0, :bottom, :red))])
```
"""
annotate!(anns...; kw...) = plot!(; annotation = anns, kw...)
annotate!(anns::Tuple...; kw...)      = plot!(; annotation = collect(anns), kw...)
annotate!(anns::AVec{<:Tuple}; kw...) = plot!(; annotation = anns, kw...)

"Flip the current plots' x axis"
xflip!(flip::Bool = true; kw...) = plot!(; xflip = flip, kw...)

"Flip the current plots' y axis"
yflip!(flip::Bool = true; kw...) = plot!(; yflip = flip, kw...)

"Specify x axis attributes for an existing plot"
xaxis!(args...; kw...) = plot!(; xaxis = args, kw...)

"Specify y axis attributes for an existing plot"
yaxis!(args...; kw...) = plot!(; yaxis = args, kw...)
xgrid!(args...; kw...) = plot!(; xgrid = args, kw...)
ygrid!(args...; kw...) = plot!(; ygrid = args, kw...)

@specialize


text:
export LatLonSeparableDriver

type LatLonSeparableDriver <: DriverModel{LatLonAccel, IntegratedContinuous}
    C::IntegratedContinuous
    mlat::LateralDriverModel
    mlon::LongitudinalDriverModel
end

get_name(model::LatLonSeparableDriver) = @sprintf("%s + %s", get_name(model.mlat), get_name(model.mlon))
action_context(model::LatLonSeparableDriver) = model.C
function set_desired_speed!(model::LatLonSeparableDriver, v_des::Float64)
    set_desired_speed!(model.mlon, v_des)
    model
end
function observe!(model::LatLonSeparableDriver, scene::Scene, roadway::Roadway, egoid::Int)
    observe!(model.mlat, scene, roadway, egoid)
    observe!(model.mlon, scene, roadway, egoid)
    model
end
function Base.rand(model::LatLonSeparableDriver)
    alat = rand(model.mlat)
    alon = rand(model.mlon)
    LatLonAccel(alat, alon)
end
Distributions.pdf(model::LatLonSeparableDriver, a::LatLonAccel) = pdf(model.mlat, a.a_lat) * pdf(model.mlon, a.a_lon)
Distributions.logpdf(model::LatLonSeparableDriver, a::LatLonAccel) = logpdf(model.mlat, a.a_lat) + logpdf(model.mlon, a.a_lon)


text:
using ArgParse
using JSON
using Printf
using Dates
using Bridge
using StaticArrays
using BridgeSDEInference
using BridgeSDEInference: indices
using Random, LinearAlgebra, Distributions
import Distributions.logpdf
using NPZ
using PyCall
const State = SArray{Tuple{2},T,1,2} where {T}

settings = ArgParseSettings(
    description="Run FitzHugh-Nagumo model experiment (noisy observations, BridgeSDEInference)"
)

@add_arg_table! settings begin
    "--output-root-dir"
        default = "experiments"
        arg_type = String
        help = "Root directory to make experiment output subdirectory in"
    "--observation-noise-std"
        help = "Standard deviation of observation noise"
        arg_type = Float64
        default = 0.01
    "--num-steps-per-obs"
        help = "Number of time steps per interobservation interval to use in inference"
        arg_type = Int
        default = 40
    "--seed-chains"
        help = "Seed for random number generator used to simulate chains"
        arg_type = Int
        default = 20200710
    "--seed-data"
        help = "Seed for random number generator used to simulate observed data"
        arg_type = Int
        default = 4
    "--num-chain"
        help = "Number of independent chains to sample"
        arg_type = Int
        default = 4
    "--num-warm-up-iter"
        help = "Number of chain iterations in adaptive warm-up sampling stage"
        arg_type = Int
        default = 100000
    "--num-main-iter"
        help = "Number of chain iterations in main sampling stage"
        arg_type = Int
        default = 500000
    "--num-iter-per-readjust"
        help = "Number of chain iterations between each adaptation (readjustment) of sampler parameters during warm-up sampling stage"
        arg_type = Int
        default = 100
    "--pcn-imputation-update-initial-persistence"
        help = "Initial value of persistence parameter (in [0, 1]) for pre-conditioned Crank-Nicoloson update to paths"
        arg_type = Float64
        default = 0.975
    "--rwm-parameter-update-initial-step-size"
        help = "Initial value of step size parameter (positive) for random-walk Metropolis update to parameters"
        arg_type = Float64
        default = 0.5
    "--num-iter-per-print-summary"
        help = "Number of chain iterations between printing a summary of chain progress"
        arg_type = Int
        default = 500
    "--num-iter-per-save-path"
        help = "Number of chain iterations between each record of diffusion path"
        arg_type = Int
        default = 500
end

args = parse_args(ARGS, settings)

function logpdf(p::Normal, coord_idx, θ::Array{Float64})
    θi = θ[[indices(coord_idx)...]]
    @assert length(θi) == 1
    logpdf.(p, θi[1])
end

function logpdf(p::LogNormal, coord_idx, θ::Array{Float64})
    θi = θ[[indices(coord_idx)...]]
    @assert length(θi) == 1
    logpdf.(p, θi[1])
end

# Output directory

dir_name = @sprintf("σ_%.2g_%s", args["observation-noise-std"], string(DateTime(Dates.now())))
output_dir = joinpath(args["output-root-dir"], "fhn_noisy_bridge", dir_name)

if !isdir(output_dir)
    mkpath(output_dir)
end

open(joinpath(output_dir, "args.json"), "w") do f
  JSON.print(f, args, 2)
end

# Reference 'true' values

θ_ref = [
    0.1, # ϵ
    -0.8, # s
    1.5, # γ
    0.0, # β
    0.3, # σ
]
x0_ref = State(-0.5, -0.6)

# Prior distributions
θ_priors = Dict(
    "ϵ" => LogNormal(0, 1),
    "s" => Normal(0, 1),
    "γ" => LogNormal(0, 1),
    "β" => Normal(0, 0),
    "σ" => LogNormal(0, 1),
)
x0_prior = MultivariateNormal([0, 0], I)

# Diffusion

param = :regular
P = FitzhughDiffusion(param, θ_ref...)

# MCMC update information for θ coords

θ_coord_names = ["ϵ", "s", "γ", "β", "σ"]
θ_name_to_coord_dict = Dict(name => coord for (coord, name) in enumerate(θ_coord_names))
θ_update_coords = [1, 2, 3, 5]
θ_is_nonnegative = Dict(name => minimum(θ_priors[name]) == 0 for name in θ_coord_names)

# Observation scheme

num_obs = 100
L = @SMatrix [1. 0.]
σ_y = args["observation-noise-std"]
Σ = @SMatrix [σ_y^2]

# Time grid

dt = 1/50000
T = 20.0
tt = 0.0:dt:T

# Imputation grid

dt_inference = T / (num_obs * args["num-steps-per-obs"])

# Generate simulated data

Random.seed!(args["seed-data"])
X, _ = simulate_segment(0.0, x0_ref, P, tt);
num_obs = 100
num_steps_per_obs = div(length(tt), num_obs)
obs_noise = [rand(Normal(0, 1)) for i in 1:num_obs]

obs = (
    times = X.tt[num_steps_per_obs+1:num_steps_per_obs:end],
    values = [
        ((L*x)[1] + σ_y * n)
        for (x, n) in zip(X.yy[num_steps_per_obs+1:num_steps_per_obs:end], obs_noise)
    ]
)

# Generate chain initializations

Random.seed!(args["seed-chains"])
n_chain = args["num-chain"]
θ_inits = Array{Float64, 2}(undef, n_chain, length(θ_ref))
x0_inits = Array{Float64, 2}(undef, n_chain, length(x0_ref))
for c in 1:n_chain
    θ_inits[c, :] = [rand(θ_priors[name]) for name in θ_coord_names]
    x0_inits[c, :] = rand(x0_prior)
end

## Run chains

println("Starting chains for σ_y = $σ_y")
println('='^80)
mcmc_outputs = Array{Any, 1}(undef, n_chain)
mcmc_times = Array{Float64}(undef, n_chain)
for c in 1:n_chain
    println("Starting sampling chain $c of $n_chain...")
    println('-'^80)
    x0_init = State(x0_inits[c, :]...)
    θ_init = θ_inits[c, :]
    P_trgt = FitzhughDiffusion(param, θ_init...)
    P_aux = [
        FitzhughDiffusionAux(param, θ_init..., t₀, u, T, v)
        for (t₀, T, u, v) in zip(obs.times[1:end-1], obs.times[2:end],
                                 obs.values[1:end-1], obs.values[2:end])
    ]
    model_setup = DiffusionSetup(P_trgt, P_aux, PartObs())
    set_imputation_grid!(model_setup, dt_inference)
    set_observations!(
        model_setup, [L for _ in P_aux], [Σ for _ in P_aux],
        obs.values, obs.times)
    set_x0_prior!(model_setup, GsnStartingPt(mean(x0_prior), cov(x0_prior)), x0_init)
    initialise!(eltype(x0_init), model_setup, Vern7(), false, NoChangePt(100))
    set_auxiliary!(model_setup; skip_for_save=1, adaptive_prop=NoAdaptation())
    θ_updates = [
        ParamUpdate(
            MetropolisHastingsUpdt(), i, θ_init,
            UniformRandomWalk(args["rwm-parameter-update-initial-step-size"], θ_is_nonnegative[name]),
            θ_priors[name],
            UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P_aux, i))
        )
        for (i, name) in zip(θ_update_coords, θ_coord_names[θ_update_coords])
    ]
    mcmc_setup = MCMCSetup(
        Imputation(NoBlocking(), args["pcn-imputation-update-initial-persistence"],
        Vern7()), θ_updates...
    )
    schedule = MCMCSchedule(
        args["num-warm-up-iter"] + args["num-main-iter"], [[1, 2, 3, 4, 5]],
        (save=args["num-iter-per-save-path"], verbose=args["num-iter-per-print-summary"], warm_up=0,
        readjust=(x->((x%args["num-iter-per-readjust"]==0) && (x < args["num-warm-up-iter"]))), fuse=(x->false))
    )
    GC.gc()
    mcmc_times[c] = @elapsed mcmc_outputs[c] = mcmc(mcmc_setup, schedule, model_setup)
    println("...finshed sampling chain $c of $n_chain, time taken: $(mcmc_times[c]) seconds.")
    println('-'^80)
end

num_path = size(mcmc_outputs[1][1].paths, 1)
paths = permutedims(
    cat([
            cat([
                    hcat([Array(x) for x in mcmc_outputs[c][1].paths[p]]...)
                    for p in 1:num_path
                ]..., dims=3)
            for c in 1:n_chain
        ]..., dims=4
    ), (4, 3, 2, 1)
)
path_times = vcat(mcmc_outputs[1][1].time...)
traces = Dict(
    name=>hcat([[c[i] for c in output[2].θ_chain] for output in mcmc_outputs]...)'
    for (i, name) in zip(θ_update_coords, θ_coord_names[θ_update_coords]))
npzwrite(
    joinpath(output_dir, "traces.npz"),
    merge(traces, Dict("chain_times" => mcmc_times, "paths" => paths, "path_times" => path_times))
)

arviz = pyimport("arviz")
summary = arviz.summary(traces)
summary_dict = summary.to_dict()
summary_dict["total_sampling_time"] = sum(mcmc_times)

open(joinpath(output_dir, "summary.json"), "w") do f
  JSON.print(f, summary_dict, 2)
end

show(summary)

println("\nSaved chain traces and summary to " * output_dir)
println('='^80)

max_stars_repo_name: thiery-lab/manifold-mcmc-for-diffusions

max_issues_repo_name: thiery-lab/manifold-mcmc-for-diffusions

max_forks_repo_name: thiery-lab/manifold-mcmc-for-diffusions


text:
function MarkCellsBySteps(r)
    side = Ost
    moves_Sud_from_start = movementsandcountmoves(r, Sud)
    moves_West_from_start = movementsandcountmoves(r, West)
    counter = movementsandcountmoves(r, Ost)
    movements(r, West)

    markrowbyfixmoves(r, side, counter)
    putmarker!(r)
    movements(r, inverse(side))
    counter -= 1
    while isborder(r, Nord) == false
        move!(r, Nord)
        markrowbyfixmoves(r, side, counter)
        if counter >= 0
            putmarker!(r)
        end
        movements(r, inverse(side))
        counter -= 1
    end

    movements(r, Sud)
    movements(r, West)
    movementsbyfixmoves(r, Ost, moves_West_from_start)
    movementsbyfixmoves(r, Nord, moves_Sud_from_start)
end

function inverse(side)
    return HorizonSide((Int(side) + 2) % 4)
end

function markrowbyfixmoves(r, side, counter)
    while (isborder(r, side) == false) && (counter > 0)
        putmarker!(r)
        move!(r, side)
        counter -= 1
    end
end

function movements(r, side)
    while isborder(r, side) == false
        move!(r, side)
    end
end

function movementsandcountmoves(r, side)
    moves = 0
    while isborder(r, side) == false
        move!(r, side)
        moves += 1
    end
    return moves
end

function movementsbyfixmoves(r, side, moves)
    while moves != 0
        move!(r, side)
        moves -= 1
    end
end


text:
if(Pkg.installed("PyPlot") == nothing)
    Pkg.add("PyPlot")
end
using PyPlot

if !isfile("matiss/matiss_bc2016.oifits")
    if !isdir("matiss")
        mkdir("matiss")
    end
    download("http://www.opticalinterferometry.com/s/Object2-v3.oifits","matiss/matiss_bc2016.oifits")
end

function plotfunction(PDATA::PAINTER.PAINTER_Data,OIDATA::PAINTER.PAINTER_Input)
    nx = OIDATA.nx
    nw = OIDATA.nw
    wvl = OIDATA.wvl
    FOV = OIDATA.FOV
    x = PDATA.x
    w = PDATA.w
    if sum(w)==0
      w = 1.
    end
    indpix = linspace(-(FOV / 2), (FOV / 2), nx)
    pos = round(Int, [1, round(Int,nx / 4), round(Int,nx / 2), round(Int,nx * 3 / 4), nx])

    count_y = 0
    count_x = 0
    SubRow  = 2
    SubColumn = 7

    X3D = x .* max(0, w)
    VMIN = minimum(X3D)
    VMAX = maximum(X3D)

    figure("3D Image Reconstruction ")
    if nw == 250
      set = collect(1:10:250)
    elseif nw == 25
      set = collect(1:25)
    end
    n2 = 0
    for n in set
        n2+=1
        subplot(SubColumn, SubRow, n2)
        imshow( sqrt(sqrt( X3D[:, :, n] .* max(X3D[:, :, n],0) )) , origin = "lower",vmin=VMIN,vmax=VMAX)
        titlestring = @sprintf("%2.4f µm", wvl[n2] * 1e6)
        title(titlestring)
        xticks([])
        yticks([])

        if( n == (nw + 1 - SubRow + count_x) )
            xticks(collect(pos - 1), round(Int,indpix[pos] * 100000) / 100)
            xlabel("FOV (mas)")
            count_x += 1
        end
        if(n == (1 + count_y * SubRow))
            yticks(collect(pos - 1), round(Int,indpix[pos] * 100000) / 100)
            ylabel("FOV (mas)")
            count_y += 1
        end
    end

    figure("Projection")

    subplot(223); imshow(squeeze(mean(X3D,1),1).',aspect="auto");
    xticks(collect(pos - 1), round(Int,indpix[pos] * 100000) / 100)
    xlabel("FOV (mas)")
    ylabel("channels")
    title("spectrum 1")

    subplot(221); imshow(squeeze(mean(X3D,3),3), origin = "lower",aspect="auto");
    yticks(collect(pos - 1), round(Int,indpix[pos] * 100000) / 100)
    xticks([])
    ylabel("FOV (mas)")
    title("Gray")

    subplot(222); imshow( squeeze(mean(X3D,2) ,2) , origin = "lower",aspect="auto");
    yticks([])
    xlabel("channels")
    title("Spectrum 2")

end

    Folder="matiss"
    savepath = "matiss_bc2016.jld"

    dptype = "phase"
    dpprm =  0

    CountPlot = 5
    nbitermax = 100

    aff = true
    admm = true

    PlotFct = plotfunction

    FOV = 0.4
    nx = 64
    indwvl = []

    rho_y = 10
    rho_spat = .5
    rho_ps   = .05
    rho_spec = .05

    alpha = 1e6
    beta = 1e8

    lambda_spat =  1e-3
    lambda_spec = 1e-4
    lambda_L1 = 0.01

    epsilon = 1e-6

    eps1 = 1e-3
    eps2 = 1e-3

    xinit3D = []
    mask3D = PAINTER.mask(nx,round(Int, nx/2 - 2 ))

# initialize algorithm and run admm
    OIDATA, PDATA = PAINTER.painter(nbitermax = nbitermax, nx = nx, lambda_spat = lambda_spat,
                            lambda_spec = lambda_spec, rho_y = rho_y, rho_spat = rho_spat,
                            rho_spec = rho_spec, rho_ps = rho_ps, alpha = alpha, beta = beta,
                            eps1 = eps1, eps2 = eps2, FOV = FOV, indwvl = indwvl, admm = admm,
                            PlotFct = PlotFct, aff = aff, dptype = dptype,
                            dpprm = dpprm, Folder = Folder, lambda_L1=lambda_L1)

# save data struture in .jld files
    println("save results of gravity BC2016 data")
    PAINTER.paintersave(savepath,PDATA,OIDATA)

# load data struture in .jld files
    println("load results of gravity BC2016 data")
    PDATA, OIDATA = PAINTER.painterload(savepath)

    nothing


text:
using Paper
using Markdown

using Gadfly

Paper.rewire(Gadfly.Plot) do p
	Paper.addtochunk(convert(Paper.Tile, p))
end

Paper.launch()

###############################

Paper.chunk(:start)
Paper.title(3, "M-Splines")
Paper.chunk(:next)
Paper.title(3, "M-Splines")

"abcd" |> borderwidth(1px) |> bordercolor("#444") |> borderstyle(dashed)

Paper.chunk(:bottom)
"text"

function mspline(i, k, ts, x)
	if k == 1
		idx = searchsortedlast(ts, x)
		(idx == 0 | idx == length(ts)) && return 0.
		(idx != i) && return 0.
		return 1. / ( ts[ idx+1 ] - ts[ idx ] )
	else
		ti, tik  = ts[ max(i, 1) ], ts[ min(i+k, length(ts)) ]
		(ti == tik) && return 0
		num = (x-ti) * mspline(i, k-1, ts, x) + (tik-x) * mspline(i+1, k-1, ts, x)
		return k * num / ( (k-1) * (tik-ti) )
	end
end

tex("M_k(x) = \\frac{(x-t_i).M_{k-1}^{i}(x)+(t_ik-x).M_{k-1}^{i+1}(x)}{(k-1)(t_ik-t_i)}") |> pad(1em)

md"""
- exemple

```
	function mspline(i, k, ts, x)
		if k == 1
			idx = searchsortedlast(ts, x)
			(idx == 0 | idx == length(ts)) && return 0.
			(idx != i) && return 0.
			return 1. / ( ts[ idx+1 ] - ts[ idx ] )
		else
			ti, tik  = ts[ max(i, 1) ], ts[ min(i+k, length(ts)) ]
			(ti == tik) && return 0
			num = (x-ti) * mspline(i, k-1, ts, x) + (tik-x) * mspline(i+1, k-1, ts, x)
			return k * num / ( (k-1) * (tik-ti) )
		end
	end
```

- finalizer(x, function)
"""

vskip(1em)

Paper.chunk(:body)
title(1, "avec 4 courbes")

ts = [0., 0.3, 0.7, 1.]
pl = plot([ x -> mspline(i, 3, ts, x) for i in -1:3] , 0, 1) ;
vbox(pl, caption(md"This is a Graph of **M-Splines**")) |>
	packacross(center) |>
	borderstyle(solid)

Paper.chunk(:body2)
title(1, "avec plus de courbes")

ts = [0:0.1:1.]
pl = plot([ x -> mspline(i, 3, ts, x) for i in -1:length(ts)-1] , 0, 1) ;
vbox(pl,
	 caption(md"This is a Graph of **M-Splines** \\n 11 points")) |>
	packacross(center) |>
	borderstyle(solid)

5+6

"aze"
"xyz"

include("examples/mini.jl")


text:
"""
Abstract type for storing neuroimaging sensors.

Other types inherit from the Sensor type.
And common functions can be run on all sensors sub types.
All sensors have a label and coordinate.
Some sensors also store additional information.
For example, fNIRS sensors may hold wavelength information.

All Sensor types support the following functions:

* `label()`
* `labels()`
* `x()`
* `y()`
* `z()`

```julia
my_sensor = # Create a electrode, optode etc
label(my_sensor)  # Returns the sensor name
x(my_sensor)      # Returns the x coordinate of the sensor
```
"""
abstract type Sensor end

"""
Electrode sensor type used in EEG measurements.
"""
mutable struct Electrode <: Sensor
    label::AbstractString
    coordinate::Coordinate
    info::Dict
end

"""
Optode abstract sensor type used in fNIRS measrurements.
"""
abstract type Optode <: Sensor end

"""
Source optode sensor type used in fNIRS measurements.
"""
mutable struct Source <: Optode
    label::AbstractString
    coordinate::Coordinate
    info::Dict
end

"""
Detector optode sensor type used in fNIRS measurements.
"""
mutable struct Detector <: Optode
    label::AbstractString
    coordinate::Coordinate
    info::Dict
end

import Base.show
function show(s::S) where {S<:Sensor}
    println(
        "Sensor: $(s.label) $(typeof(s)) - ($(s.coordinate.x |> u"mm"), $(s.coordinate.y |> u"mm"), $(s.coordinate.z |> u"mm")) ($(typeof(s.coordinate)))",
    )
end

function show(s::Array{S}) where {S<:Sensor}
    println("$(length(s)) sensors: $(typeof(s[1])) ($(typeof(s[1].coordinate)))")
end

label(s::S) where {S<:Sensor} = s.label
label(s::Array{S,1}) where {S<:Sensor} = [si.label for si in s]
labels(s::S) where {S<:Sensor} = label(s)
labels(s::Array{S}) where {S<:Sensor} = label(s)

x(s::S) where {S<:Sensor} = s.coordinate.x
y(s::S) where {S<:Sensor} = s.coordinate.y
z(s::S) where {S<:Sensor} = s.coordinate.z
x(s::Array{S}) where {S<:Sensor} = [si.coordinate.x for si in s]
y(s::Array{S}) where {S<:Sensor} = [si.coordinate.y for si in s]
z(s::Array{S}) where {S<:Sensor} = [si.coordinate.z for si in s]


text:
module ReadStat

import StringEncodings

# Load libreadstat from our deps.jl
const depsjl_path = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
if !isfile(depsjl_path)
    error("ReadStat not installed properly, run Pkg.build(\"ReadStat\"), restart Julia and try again")
end
include(depsjl_path)

##############################################################################
##
## Import
##
##############################################################################

using DataValues: DataValueVector
using Dates

export ReadStatDataFrame, read_dta, read_sav, read_por, read_sas7bdat

##############################################################################
##
## Julia types that mirror C types
##
##############################################################################

const READSTAT_TYPE_STRING      = Cint(0)
const READSTAT_TYPE_CHAR        = Cint(1)
const READSTAT_TYPE_INT16       = Cint(2)
const READSTAT_TYPE_INT32       = Cint(3)
const READSTAT_TYPE_FLOAT       = Cint(4)
const READSTAT_TYPE_DOUBLE      = Cint(5)
const READSTAT_TYPE_LONG_STRING = Cint(6)

const READSTAT_ERROR_OPEN       = Cint(1)
const READSTAT_ERROR_READ       = Cint(2)
const READSTAT_ERROR_MALLOC     = Cint(3)
const READSTAT_ERROR_USER_ABORT = Cint(4)
const READSTAT_ERROR_PARSE      = Cint(5)

##############################################################################
##
## Pure Julia types
##
##############################################################################

struct ReadStatValue
    union::Int64
    readstat_types_t::Cint
    tag::Cchar
    @static if Sys.iswindows()
        bits::Cuint
    else
        bits::UInt8
    end
end

const Value = ReadStatValue

mutable struct ReadStatDataFrame
    data::Vector{Any}
    headers::Vector{Symbol}
    types::Vector{DataType}
    labels::Vector{String}
    formats::Vector{String}
    storagewidths::Vector{Csize_t}
    measures::Vector{Cint}
    alignments::Vector{Cint}
    val_label_keys::Vector{String}
    val_label_dict::Dict{String, Dict{Any,String}}
    rows::Int
    columns::Int
    filelabel::String
    timestamp::DateTime
    format::Clong

    ReadStatDataFrame() =
        new(Any[], Symbol[], DataType[], String[], String[], Csize_t[], Cint[], Cint[],
        String[], Dict{String, Dict{Any,String}}(), 0, 0, "", Dates.unix2datetime(0), 0)
end

include("C_interface.jl")

##############################################################################
##
## Julia functions
##
##############################################################################

function handle_info!(obs_count::Cint, var_count::Cint, ds_ptr::Ptr{ReadStatDataFrame})
    ds = unsafe_pointer_to_objref(ds_ptr)
    ds.rows = obs_count
    ds.columns = var_count
    return Cint(0)
end

function handle_metadata!(metadata::Ptr{Nothing}, ds_ptr::Ptr{ReadStatDataFrame})
    ds = unsafe_pointer_to_objref(ds_ptr)
    ds.filelabel = readstat_get_file_label(metadata)
    ds.timestamp = Dates.unix2datetime(readstat_get_modified_time(metadata))
    ds.format = readstat_get_file_format_version(metadata)
    ds.rows = readstat_get_row_count(metadata)
    ds.columns = readstat_get_var_count(metadata)
    return Cint(0)
end

get_name(variable::Ptr{Nothing}) = Symbol(readstat_variable_get_name(variable))

function get_label(var::Ptr{Nothing})
    ptr = ccall((:readstat_variable_get_label, libreadstat), Cstring, (Ptr{Nothing},), var)
    ptr == C_NULL ? "" : unsafe_string(ptr)
end

function get_format(var::Ptr{Nothing})
    ptr = ccall((:readstat_variable_get_format, libreadstat), Cstring, (Ptr{Nothing},), var)
    ptr == C_NULL ? "" : unsafe_string(ptr)
end

function get_type(variable::Ptr{Nothing})
    data_type = readstat_variable_get_type(variable)

    if data_type == READSTAT_TYPE_STRING
        return String
    elseif data_type == READSTAT_TYPE_CHAR
        return Int8
    elseif data_type == READSTAT_TYPE_INT16
        return Int16
    elseif data_type == READSTAT_TYPE_INT32
        return Int32
    elseif data_type == READSTAT_TYPE_FLOAT
        return Float32
    elseif data_type == READSTAT_TYPE_DOUBLE
        return Float64
    end
    return Nothing
end

get_storagewidth(variable::Ptr{Nothing}) = readstat_variable_get_storage_width(variable)

get_measure(variable::Ptr{Nothing}) = readstat_variable_get_measure(variable)

get_alignment(variable::Ptr{Nothing}) = readstat_variable_get_measure(variable)

function handle_variable!(var_index::Cint, variable::Ptr{Nothing},
                         val_label::Cstring,  ds_ptr::Ptr{ReadStatDataFrame})
    col = var_index + 1
    ds = unsafe_pointer_to_objref(ds_ptr)

    push!(ds.val_label_keys, (val_label == C_NULL ? "" : unsafe_string(val_label)))
    push!(ds.headers, get_name(variable))
    push!(ds.labels, get_label(variable))
    push!(ds.formats, get_format(variable))
    jtype = get_type(variable)
    push!(ds.types, jtype)
    push!(ds.data, DataValueVector{jtype}(ds.rows))
    push!(ds.storagewidths, get_storagewidth(variable))
    push!(ds.measures, get_measure(variable))
    push!(ds.alignments, get_alignment(variable))

    return Cint(0)
end

function get_type(val::Value)
    data_type = readstat_value_type(val)
    return [String, Int8, Int16, Int32, Float32, Float64, String][data_type + 1]
end

Base.convert(::Type{Int8}, val::Value) = ccall((:readstat_int8_value, libreadstat), Int8, (Value,), val)
Base.convert(::Type{Int16}, val::Value) = ccall((:readstat_int16_value, libreadstat), Int16, (Value,), val)
Base.convert(::Type{Int32}, val::Value) = ccall((:readstat_int32_value, libreadstat), Int32, (Value,), val)
Base.convert(::Type{Float32}, val::Value) = ccall((:readstat_float_value, libreadstat), Float32, (Value,), val)
Base.convert(::Type{Float64}, val::Value) = ccall((:readstat_double_value, libreadstat), Float64, (Value,), val)
function Base.convert(::Type{String}, val::Value)
    ptr = ccall((:readstat_string_value, libreadstat), Cstring, (Value,), val)
    ptr ≠ C_NULL ? unsafe_string(ptr) : ""
end
as_native(val::Value) = convert(get_type(val), val)

function handle_value!(obs_index::Cint, variable::Ptr{Nothing},
                       value::ReadStatValue, ds_ptr::Ptr{ReadStatDataFrame})
    ds = unsafe_pointer_to_objref(ds_ptr)
    var_index = readstat_variable_get_index(variable)
    if !readstat_value_is_missing(value, variable)
        readfield!(ds.data[var_index + 1], obs_index + 1, value)
    end

    return Cint(0)
end

function readfield!(dest::DataValueVector{String}, row, val::Value)
    ptr = ccall((:readstat_string_value, libreadstat), Cstring, (Value,), val)
    if ptr ≠ C_NULL
        @inbounds dest[row] = unsafe_string(ptr)
    end
end

function readfield!(dest::DataValueVector{Int8}, row, val::Value)
    @inbounds dest[row] = ccall((:readstat_int8_value, libreadstat), Int8, (Value,), val)
end

function readfield!(dest::DataValueVector{Int16}, row, val::Value)
    @inbounds dest[row] = ccall((:readstat_int16_value, libreadstat), Int16, (Value,), val)
end

function readfield!(dest::DataValueVector{Int32}, row, val::Value)
    @inbounds dest[row] = ccall((:readstat_int32_value, libreadstat), Int32, (Value,), val)
end

function readfield!(dest::DataValueVector{Float64}, row, val::Value)
    @inbounds dest[row] = ccall((:readstat_double_value, libreadstat), Float64, (Value,), val)
end

function readfield!(dest::DataValueVector{Float32}, row, val::Value)
    @inbounds dest[row] = ccall((:readstat_float_value, libreadstat), Float32, (Value,), val)
end

function handle_value_label!(val_labels::Cstring, value::Value, label::Cstring, ds_ptr::Ptr{ReadStatDataFrame})
    val_labels ≠ C_NULL || return Cint(0)
    ds = unsafe_pointer_to_objref(ds_ptr)
    dict = get!(ds.val_label_dict, unsafe_string(val_labels), Dict{Any,String}())
    dict[as_native(value)] = unsafe_string(label)

    return Cint(0)
end

function read_data_file(filename::AbstractString, filetype::Val)
    # initialize ds
    ds = ReadStatDataFrame()
    # initialize parser
    parser = Parser()
    # parse
    parse_data_file!(ds, parser, filename, filetype)
    # return dataframe instead of ReadStatDataFrame
    return ds
end

function Parser()
    parser = ccall((:readstat_parser_init, libreadstat), Ptr{Nothing}, ())
    info_fxn = @cfunction(handle_info!, Cint, (Cint, Cint, Ptr{ReadStatDataFrame}))
    meta_fxn = @cfunction(handle_metadata!, Cint, (Ptr{Nothing}, Ptr{ReadStatDataFrame}))
    var_fxn = @cfunction(handle_variable!, Cint, (Cint, Ptr{Nothing}, Cstring,  Ptr{ReadStatDataFrame}))
    val_fxn = @cfunction(handle_value!, Cint, (Cint, Ptr{Nothing}, ReadStatValue, Ptr{ReadStatDataFrame}))
    label_fxn = @cfunction(handle_value_label!, Cint, (Cstring, Value, Cstring, Ptr{ReadStatDataFrame}))
    ccall((:readstat_set_metadata_handler, libreadstat), Int, (Ptr{Nothing}, Ptr{Nothing}), parser, meta_fxn)
    ccall((:readstat_set_variable_handler, libreadstat), Int, (Ptr{Nothing}, Ptr{Nothing}), parser, var_fxn)
    ccall((:readstat_set_value_handler, libreadstat), Int, (Ptr{Nothing}, Ptr{Nothing}), parser, val_fxn)
    ccall((:readstat_set_value_label_handler, libreadstat), Int, (Ptr{Nothing}, Ptr{Nothing}), parser, label_fxn)
    return parser
end

function parse_data_file!(ds::ReadStatDataFrame, parser::Ptr{Nothing}, filename::AbstractString, filetype::Val)
    retval = readstat_parse(filename, filetype, parser, ds)
    readstat_parser_free(parser)
    retval == 0 ||  error("Error parsing $filename: $retval")
end

read_dta(filename::AbstractString) = read_data_file(filename, Val(:dta))
read_sav(filename::AbstractString) = read_data_file(filename, Val(:sav))
read_por(filename::AbstractString) = read_data_file(filename, Val(:por))
read_sas7bdat(filename::AbstractString) = read_data_file(filename, Val(:sas7bdat))

end #module ReadStat


text:
#! /usr/bin/julia

# Rosetta Code, Date format

ts = time()

println("The system time is (in ISO 8601 format):")
println(strftime("    %F %T %Z", ts))

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
@testset "minimal_nums_to_interval_test" begin
    @test Interval(-1.0,1.0) == interval(-1.0,1.0)
    @test Interval(-Inf,1.0) == interval(-Inf,1.0)
    @test Interval(-1.0,Inf) == interval(-1.0,Inf)
    @test Interval(-Inf,Inf) == interval(-Inf,Inf)
    @test_throws AssertionError Interval(NaN,NaN) == emptyinterval()
    @test_throws AssertionError Interval(1.0,-1.0) == emptyinterval()
    @test_throws AssertionError Interval(-Inf,-Inf) == emptyinterval()
    @test_throws AssertionError Interval(Inf,Inf) == emptyinterval()
end

@testset "minimal_text_to_interval_test" begin
end

@testset "minimal_interval_part_test" begin
end


text:
using LinearAlgebra

import LinearAlgebra: mul!, rmul!

lmul!(G::LinearAlgebra.Givens, ::Nothing) = nothing
rmul!(::Nothing, G::LinearAlgebra.Givens) = nothing

function svdvals2x2(d1, d2, e)
    d1sq = d1*d1
    d2sq = d2*d2
    esq  = e*e
    b = d1sq + d2sq + esq
    D = sqrt(abs2((d1 + d2)*(d1 - d2)) + esq*(2*(d1sq + d2sq) + esq))
    D2 = b + D
    λ1 = 2*d1sq*d2sq/D2
    λ2 = D2/2
    return minmax(sqrt(λ1), sqrt(λ2))
end

function svdIter!(B::Bidiagonal{T}, n1, n2, shift, U = nothing, Vᴴ = nothing) where T<:Real

    if B.uplo === 'U'

        d = B.dv
        e = B.ev

        G, r = givens(d[n1] - abs2(shift)/d[n1], e[n1], n1, n1 + 1)
        lmul!(G, Vᴴ)

        ditmp       = d[n1]
        ei1         = e[n1]
        di          = ditmp*G.c + ei1*G.s
        ei1         = -ditmp*G.s + ei1*G.c
        di1         = d[n1 + 1]
        bulge       = di1*G.s
        di1        *= G.c

        for i = n1:n2 - 2
            G, r      = givens(di, bulge, i, i + 1)
            rmul!(U, G')
            d[i]      = G.c*di + G.s*bulge
            ei        = G.c*ei1 + G.s*di1
            di1       = -G.s*ei1 + G.c*di1
            ei1       = e[i + 1]
            bulge     = G.s*ei1
            ei1      *= G.c

            G, r      = givens(ei, bulge, i + 1, i + 2)
            lmul!(G, Vᴴ)
            e[i]      = ei*G.c + bulge*G.s
            di        = di1*G.c + ei1*G.s
            ei1       = -di1*G.s + ei1*G.c
            di2       = d[i + 2]
            bulge     = di2*G.s
            di1       = di2*G.c
        end

        G, r      = givens(di, bulge, n2 - 1, n2)
        rmul!(U, G')
        d[n2 - 1] = G.c*di + G.s*bulge
        e[n2 - 1] = G.c*ei1 + G.s*di1
        d[n2]     = -G.s*ei1 + G.c*di1

    else
        throw(ArgumentError("lower bidiagonal version not implemented yet"))
    end

    return B
end

# See LAWN 3
function svdDemmelKahan!(B::Bidiagonal{T}, n1, n2, U = nothing, Vᴴ = nothing) where T<:Real

    if B.uplo === 'U'

        d = B.dv
        e = B.ev

        oldcs = one(T)
        G     = LinearAlgebra.Givens(1, 2, one(T), zero(T))
        Gold  = G

        for i = n1:n2 - 1
            G, r  = givens(d[i] * G.c, e[i], i, i + 1)
            lmul!(G, Vᴴ)
            if i != n1
                e[i - 1] = Gold.s * r
            end

            Gold, d[i] = givens(Gold.c * r, d[i + 1] * G.s, i, i + 1)
            rmul!(U, Gold')
        end

        h         = d[n2] * G.c
        e[n2 - 1] = h * Gold.s
        d[n2]     = h * Gold.c

    else
        throw(ArgumentError("lower bidiagonal version not implemented yet"))
    end

    return B
end

# Recurrence to estimate smallest singular value from LAWN3 Lemma 1
function estimate_σ⁻(dv::AbstractVector, ev::AbstractVector, n1::Integer, n2::Integer)
    λ  = abs(dv[n2])
    B∞ = λ
    for j in (n2 - 1):-1:n1
        λ = abs(dv[j])*(λ/(λ + abs(ev[j])))
        B∞ = min(B∞, λ)
    end

    μ  = abs(dv[n1])
    B1 = μ
    for j in n1:(n2 - 1)
        μ = abs(dv[j + 1])*(μ*(μ + abs(ev[j])))
        B1 = min(B1, μ)
    end

    return min(B∞, B1)
end

# The actual SVD solver routine
# Notice that this routine doesn't adjust the sign and sorts the values
function __svd!(B::Bidiagonal{T}, U = nothing, Vᴴ = nothing; tol = 100eps(T), debug = false) where T<:Real

    n = size(B, 1)
    n1 = 1
    n2 = n
    d = B.dv
    e = B.ev
    count = 0

    # See LAWN3 page 6 and 22
    σ⁻ = estimate_σ⁻(d, e, n1, n2)
    fudge = n
    thresh = tol*σ⁻

    if B.uplo === 'U'
        while true

            # Search for biggest index for non-zero off diagonal value in e
            for n2i = n2:-1:1
                if n2i == 1
                    # We are done!
                    return nothing
                else
                    if abs(e[n2i - 1]) > thresh
                        n2 = n2i
                        break
                    end
                end
            end

            # Search for largest sub-bidiagonal matrix ending at n2
            for _n1 = (n2 - 1):-1:1
                if _n1 == 1
                    n1 = _n1
                    break
                elseif abs(e[_n1 - 1]) < thresh
                    n1 = _n1
                    break
                end
            end

            debug && println("n1=", n1, ", n2=", n2, ", d[n1]=", d[n1], ", d[n2]=", d[n2], ", e[n1]=", e[n1], " e[n2-1]=", e[n2-1], " thresh=", thresh)

            # See LAWN 3 p 18 and
            # We approximate the smallest and largest singular values of the
            # current block to determine if it's safe to use shift or if
            # the zero shift algorithm is required to maintain high relative
            # accuracy
            σ⁻ = estimate_σ⁻(d, e, n1, n2)
            σ⁺ = max(maximum(view(d, n1:n2)), maximum(view(e, n1:(n2 - 1))))

            if fudge*tol*σ⁻ <= eps(σ⁺)
                svdDemmelKahan!(B, n1, n2, U, Vᴴ)
            else
                shift = svdvals2x2(d[n2 - 1], d[n2], e[n2 - 1])[1]
                if shift^2 < eps(B.dv[n1]^2)
                    # Shift is too small to affect the iteration so just use
                    # zero-shift algorithm anyway
                    svdDemmelKahan!(B, n1, n2, U, Vᴴ)
                else
                    svdIter!(B, n1, n2, shift, U, Vᴴ)
                end
            end

            count += 1
            debug && println("count=", count)
        end
    else
        # Just transpose the matrix
        error("SVD for lower triangular bidiagonal matrices isn't implemented yet.")
    end
end

function _svdvals!(B::Bidiagonal{T}; tol = eps(T), debug = false) where T<:Real
    __svd!(B, tol = tol, debug = debug)
    return sort(abs.(diag(B)), rev = true)
end

function _sort_and_adjust!(U, s, Vᴴ)
    n = length(s)

    # Adjust sign of singular values if necessary
    if any(!isposdef, s)
        rmul!(U, Diagonal(sign.(s)))
        map!(abs, s, s)
    end

    # Sort singular values if necessary
    # FixMe! Try to do this without (much) allocation
    p = sortperm(s, by = abs, rev = true)
    if p != 1:n
        U[:]  = U[:,p]
        s[:]  = s[p]
        Vᴴ[:] = Vᴴ[p,:]
    end

    return nothing
end

function _svd!(B::Bidiagonal{T}; tol = eps(T), debug = false) where T<:Real
    n  = size(B, 1)

    U  = Matrix{T}(I, n, n)
    Vᴴ = Matrix{T}(I, n, n)

    __svd!(B, U, Vᴴ, tol = tol, debug = debug)

    s = diag(B)

    # Sort and adjust the singular values. The algorithm used in __svd! might produce negative singular values and it will return the values in no particular ordering.
    _sort_and_adjust!(U, s, Vᴴ)

    return SVD(U, s, Vᴴ)
end

struct BidiagonalFactorization{T,S,U,V}
    bidiagonal::Bidiagonal{T,S}
    reflectors::U
    τl::V
    τr::V
end

function bidiagonalize!(A::AbstractMatrix)
    m, n = size(A)
    τl, τr = eltype(A)[], eltype(A)[]

    if m >= n
        # tall case: upper bidiagonal
        for i in 1:min(m, n)
            x = view(A, i:m, i)
            τi = LinearAlgebra.reflector!(x)
            push!(τl, τi)
            LinearAlgebra.reflectorApply!(x, τi, view(A, i:m, (i + 1):n))
            if i < n
                x = view(A, i, i + 1:n)
                conj!(x)
                τi = LinearAlgebra.reflector!(x)
                push!(τr, τi)
                LinearAlgebra.reflectorApply!(view(A, (i + 1):m, (i + 1):n), x, τi)
            end
        end

        bd = Bidiagonal(real(diag(A)), real(diag(A, 1)), :U)

        return BidiagonalFactorization{eltype(bd),typeof(bd.dv),typeof(A),typeof(τl)}(bd, A, τl, τr)
    else
        # wide case: lower bidiagonal
        for i in 1:min(m, n)
            x = view(A, i, i:n)
            conj!(x)
            τi = LinearAlgebra.reflector!(x)
            push!(τr, τi)
            LinearAlgebra.reflectorApply!(view(A, (i + 1):m, i:n), x, τi)
            if i < m
                x = view(A, i + 1:m, i)
                τi = LinearAlgebra.reflector!(x)
                push!(τl, τi)
                LinearAlgebra.reflectorApply!(x, τi, view(A, (i + 1):m, (i + 1):n))
            end
        end

        bd = Bidiagonal(real(diag(A)), real(diag(A, -1)), :L)

        return BidiagonalFactorization{eltype(bd),typeof(bd.dv),typeof(A),typeof(τl)}(bd, A, τl, τr)
    end
end

# This one is currently very type unstable but it is useful to be able to reuse
# existing elementary reflectors types such as QRPackedQ, LQPackedQ, and HessenbergQ
# such that multiplication methods are available for free.
function Base.getproperty(F::BidiagonalFactorization, s::Symbol)
    BD = getfield(F, :bidiagonal)
    R  = getfield(F, :reflectors)
    τl = getfield(F, :τl)
    τr = getfield(F, :τr)
    if BD.uplo === 'U'
        if s === :leftQ
            return LinearAlgebra.QRPackedQ(R, τl)
        elseif s === :rightQ
            if VERSION < v"1.3.0-DEV.243"
                return LinearAlgebra.HessenbergQ(copy(transpose(R[1:size(R,2),:])), τr)
            else
                factors = copy(transpose(R[1:size(R,2),:]))
                return LinearAlgebra.HessenbergQ{eltype(factors),typeof(factors),typeof(τr),false}('U', factors, τr)
            end
        else
            return getfield(F, s)
        end
    else
        if s === :leftQ
            if VERSION < v"1.3.0-DEV.243"
                return LinearAlgebra.HessenbergQ(R[:,1:size(R,1)], τl)
            else
                factors = R[:,1:size(R,1)]
                return LinearAlgebra.HessenbergQ{eltype(factors),typeof(factors),typeof(τr),false}('U', factors, τl)
            end
        elseif s === :rightQ
            # return transpose(LinearAlgebra.LQPackedQ(R, τr)) # FixMe! check that this shouldn't be adjoint
            LinearAlgebra.QRPackedQ(copy(transpose(R)), τr)
        else
            return getfield(F, s)
        end
    end
end

# For now, we define a generic lmul! and rmul! for HessenbergQ here
function lmul!(Q::LinearAlgebra.HessenbergQ, B::AbstractVecOrMat)

    m, n = size(B, 1), size(B, 2)

    if m != size(Q, 2)
        throw(DimensionMismatch(""))
    end

    for j in 1:n
        for l in length(Q.τ):-1:1
            τl = Q.τ[l]
            ṽ  = view(Q.factors, (l + 2):m, l)
            b  = view(B, (l + 2):m, j)
            vᴴb = B[l + 1, j]
            if length(ṽ) > 0
                vᴴb += ṽ'b
            end
            B[l + 1, j] -= τl*vᴴb
            b .-= ṽ .* τl .* vᴴb
        end
    end

    return B
end

function lmul!(adjQ::Adjoint{<:Any,<:LinearAlgebra.HessenbergQ}, B::AbstractVecOrMat)

    Q = parent(adjQ)
    m, n = size(B, 1), size(B, 2)

    if m != size(adjQ, 2)
        throw(DimensionMismatch(""))
    end

    for j in 1:n
        for l in 1:length(Q.τ)
            τl = Q.τ[l]
            ṽ  = view(Q.factors, (l + 2):m, l)
            b  = view(B, (l + 2):m, j)
            vᴴb = B[l + 1, j]
            if length(ṽ) > 0
                vᴴb += ṽ'b
            end
            B[l + 1, j] -= τl'*vᴴb
            b .-= ṽ .* τl' .* vᴴb
        end
    end

    return B
end

function rmul!(A::AbstractMatrix, Q::LinearAlgebra.HessenbergQ)

    m, n = size(A)

    if n != size(Q, 1)
        throw(DimensionMismatch(""))
    end

    for i in 1:m
        for l in 1:(n - 1)
            τl = Q.τ[l]
            ṽ  = view(Q.factors, (l + 2):n, l)
            aᵀ = transpose(view(A, i, (l + 2):n))
            aᵀv = A[i, l + 1]
            if length(ṽ) > 0
                aᵀv += aᵀ*ṽ
            end
            A[i, l + 1] -= aᵀv*τl
            aᵀ .-= aᵀv .* τl .* ṽ'
        end
    end

    return A
end

function rmul!(A::AbstractMatrix, adjQ::Adjoint{<:Any,<:LinearAlgebra.HessenbergQ})

    m, n = size(A)
    Q = parent(adjQ)

    if n != size(adjQ, 1)
        throw(DimensionMismatch(""))
    end

    for i in 1:m
        for l in (n - 1):-1:1
            τl = Q.τ[l]
            ṽ  = view(Q.factors, (l + 2):n, l)
            aᵀ = transpose(view(A, i, (l + 2):n))
            aᵀv = A[i, l + 1]
            if length(ṽ) > 0
                aᵀv += aᵀ*ṽ
            end
            A[i, l + 1] -= aᵀv*τl'
            aᵀ .-= aᵀv .* τl' .* ṽ'
        end
    end

    return A
end

function rmul!(A::AbstractMatrix, Q::LinearAlgebra.LQPackedQ)

    m, n = size(A)

    if n != size(Q, 1)
        throw(DimensionMismatch(""))
    end

    for i in 1:m
        for l in length(Q.τ):-1:1
            τl = Q.τ[l]
            ṽ  = view(Q.factors, l, (l + 1):n)
            aᵀ = transpose(view(A, i, (l + 1):n))
            aᵀv = A[i, l]
            if length(ṽ) > 0
                aᵀv += aᵀ*ṽ
            end
            A[i, l] -= aᵀv*τl
            aᵀ .-= aᵀv .* τl .* ṽ'
        end
    end

    return A
end

# Overload LinearAlgebra methods

LinearAlgebra.svdvals!(B::Bidiagonal{T}; tol = eps(T), debug = false) where T<:Real = _svdvals!(B, tol = tol, debug = debug)

"""
    svdvals!(A [, tol, debug])

Generic computation of singular values.
```jldoctest
julia> using LinearAlgebra, GenericLinearAlgebra, Quaternions

julia> n = 20;

julia> H = [big(1)/(i + j - 1) for i in 1:n, j in 1:n]; # The Hilbert matrix

julia> Float64(svdvals(H)[end]/svdvals(Float64.(H))[end] - 1) # The relative error of the LAPACK based solution in 64 bit floating point.
-0.9999999999447275

julia> A = qr([Quaternion(randn(4)...) for i in 1:3, j in 1:3]).Q *
           Diagonal([3, 2, 1]) *
           qr([Quaternion(randn(4)...) for i in 1:3, j in 1:3]).Q'; # A quaternion matrix with the singular value 1, 2, and 3.

julia> svdvals(A) ≈ [3, 2, 1]
true
```
"""
function LinearAlgebra.svdvals!(A::StridedMatrix; tol = eps(real(eltype(A))), debug = false)
    B = bidiagonalize!(A).bidiagonal
    # It doesn't matter that we transpose the bidiagonal matrix when we are only interested in the values
    return svdvals!(Bidiagonal(B.dv, B.ev, :U), tol = tol, debug = debug)
end

# FixMe! The full keyword is redundant for Bidiagonal and should be removed from Base
LinearAlgebra.svd!(B::Bidiagonal{T};
    tol = eps(T),
    full = false,
    # To avoid breaking on <Julia 1.3, the `alg` keyword doesn't do anything. Once we drop support for Julia 1.2
    # and below, we can make the keyword argument work correctly
    alg = nothing,
    debug = false) where T<:Real = _svd!(B, tol = tol, debug = debug)

"""
    svd!(A[, tol, full, debug])::SVD

A generic singular value decomposition (SVD). The implementation only uses Julia functions so the SVD can be computed for any element type provided that the necessary arithmetic operations are supported by the element type.

- `tol`: The relative tolerance for determining convergence. The default value is `eltype(T)` where `T` is the element type of the input matrix bidiagonal (i.e. after converting the matrix to bidiagonal form).

- `full`: Sepcifies if all the left and right singular vectors be returned or if only the vectors us to the number of rows and columns of the input matrix `A` should be returned (the default).

- `debug`: A Boolean flag to activate debug information during the executions of the algorithm. The default is false.

# Algorithm
...tomorrow

# Example

```jldoctest
julia> svd(big.([1 2; 3 4]))
SVD{BigFloat, BigFloat, Matrix{BigFloat}}
U factor:
2×2 Matrix{BigFloat}:
 -0.404554   0.914514
 -0.914514  -0.404554
singular values:
2-element Vector{BigFloat}:
 5.464985704219042650451188493284182533042584640492784181017488774646871847029449
 0.3659661906262578204229643842614005434788136943931877734325179702209382149672422
Vt factor:
2×2 Matrix{BigFloat}:
 -0.576048  -0.817416
 -0.817416   0.576048
```
"""
function LinearAlgebra.svd!(A::StridedMatrix{T};
    tol = eps(real(eltype(A))),
    full = false,
    # To avoid breaking on <Julia 1.3, the `alg` keyword doesn't do anything. Once we drop support for Julia 1.2
    # and below, we can make the keyword argument work correctly
    alg = nothing,
    debug = false) where T

    m, n = size(A)

    # Convert the input matrix A to bidiagonal form and return a BidiagonalFactorization object
    BF = bidiagonalize!(A)

    # The location of the super/sub-diagonal of the bidiagonal matrix depends on the aspect ratio of the input. For tall and square matrices, the bidiagonal matrix is upper whereas it is lower for wide input matrices as illustrated below. The 'O' entried indicate orthogonal (unitary) matrices.

    #      A      =   Q_l  *   B   * Q_r^H

    #     |x x x| = |O O O| |x x  | |O O O|
    #     |x x x|   |O O O|*|  x x|*|O O O|
    #     |x x x|   |O O O| |    x| |O O O|
    #     |x x x|   |O O O|
    #     |x x x|   |O O O|

    # |x x x x x| = |O O O| |x    | |O O O O O|
    # |x x x x x|   |O O O|*|x x  |*|O O O O O|
    # |x x x x x|   |O O O| |  x x| |O O O O O|

    _B = BF.bidiagonal
    B  = _B.uplo === 'U' ? _B : Bidiagonal(_B.dv, _B.ev, :U)

    # Compute the SVD of the bidiagonal matrix B
    F = _svd!(B, tol = tol, debug = debug)

    # Form the matrices U and Vᴴ by combining the singular vector matrices of the bidiagonal SVD with the Householder reflectors from the bidiagonal factorization.
    if _B.uplo === 'U'
        U  = Matrix{T}(I, m, full ? m : n)
        U[1:n,1:n] = F.U
        lmul!(BF.leftQ, U)
        Vᴴ = F.Vt*BF.rightQ'
    else
        U  = BF.leftQ*F.V
        Vᴴ = Matrix{T}(I, full ? n : m, n)
        Vᴴ[1:m,1:m] = F.U'
        rmul!(Vᴴ, BF.rightQ')
    end

    s = F.S

    return SVD(U, s, Vᴴ)
end