text:
## Token and registration related API

"""
    register(code::String) -> String
    register() -> String

Create a new authentification token from a code obtained on https://my.remarkable.com/connect/desktop
If a code is not provided a tutorial is given.
"""
function register(code::String; path_to_token::String = "", kwargs...)
    data = Dict(
            "code" => code,
            "deviceDesc" => "desktop-windows",
            "deviceID" => string(uuid4())
        )
    @info "Registering device"
    response = HTTP.request_json(
                    "POST",
                    AUTH_API,
                    data;
                    kwargs...
    )
    token = String(response.body)
    @info "Token : $token"
    if isempty(path_to_token)
        path_to_token = pwd()
    end
    write(joinpath(path_to_token, ".token"), token)
    @info "Token saved at $(joinpath(path_to_token))"
    return token
end

function register()
    @warn """
        Hi, to use this API you first need to register this device.
        To do so go to https://my.remarkable.com/connect/desktop and ask for a new code!
        Once you have it run `register(code)` (where `code` is a `String`).
        It will return your authentification token and also save your token in a local file.
        """
end

## Items accessors

"""
    list_items(client::RemarkableClient) -> Collection

Return a `Collection` of all `Document`s and `Collection` present on the server
You can visualize them nicely via `print_tree`
"""
function list_items(client::RemarkableClient; kwargs...)
    @info "Listing all items"
    body = HTTP.request(client, "GET", STORAGE_API[] * ITEM_LIST; kwargs...)
    items = JSON.parse(String(body))
    docs = RemarkableObject[]
    for item in items
        item = Dict(Symbol(key)=>value for (key, value) in item)
        if item[:Type] == "DocumentType"
            doc = Document(;item...)
            push!(docs, doc)
        elseif item[:Type] == "CollectionType"
            collec = Collection(;item...)
            push!(docs, collec)
        end
    end
    return create_tree(docs)
end

"""
    get_item(client, id::String, download = false)
    get_item(client, id::RemarkableObject, download = false)

Return a `RemarkableObject` given an ID or an existing `RemarkableObject`,
using `download=true` will give the `BlobURLGet` to download the files
"""
function get_item(client::RemarkableClient, id::String, download::Bool = false; kwargs...)
    get_item(client, Document(ID = id), download; kwargs...)
end

function get_item(client::RemarkableClient, doc::RemarkableObject, download::Bool = false; kwargs...)
    query = Dict("doc" => doc.ID)
    if download
        query["withBlob"] = "true"
    end
    @info "Listing item"
    body = HTTP.request(
                    client,
                    "GET",
                    STORAGE_API[] * ITEM_LIST;
                    query = query,
                    kwargs...
                    )
    item = JSON.parse(String(body))
    return Document(first(item))
end

## Item modificators

"""
    delete_item!(client, id::String)
    delete_item!(client, obj::remarkable)

Delete the object from your collection (online)
"""
function delete_item!(client::RemarkableClient, id::String; kwargs...)
    delete_item!(client, Document(ID = id); kwargs...)
end

function delete_item!(client::RemarkableClient, obj::RemarkableObject; kwargs...)
    @info "Deleting item `$(title(obj))`"
    return storage_request(client, "PUT", "delete", obj_to_dict(obj); kwargs...);
end

"""
    update_metadata!(client, obj::RemarkableObject) ->

Update the metadata of an object, can be used to modify a file or create a
Collection
"""
function update_metadata!(client::RemarkableClient, obj::RemarkableObject; kwargs...)
    @info "Updating item metadata `$(title(obj))`"
    storage_request(client, "PUT", UPDATE_STATUS, obj_to_dict(obj); kwargs...)
end

"""
    create_folder!(client, name::String, parent::String="") ->

Create a folder (Collection) in `parent` (root by default)
"""
function create_folder!(client::RemarkableClient, name::String, parent::String = ""; kwargs...)
    item = Collection(
                    Parent = parent,
                    VissibleName = name,
                    )
    @info "Creating folder $name"
    res = update_metadata!(client, item; kwargs...)
    return Collection(res)
end

## Download files

"""
    download_document(client, id::String, [path_target::String]) -> ZipFile Body
    download_document(client, doc::Document, [path_target::String]) -> Zipfile Body

Download a document object with given id/doc.
The document is always given as a `ZipFile` and can be saved via `write(filepath, body)`
if `path_target` is given, the zip file is automatically written and named.
"""
function Base.download(client::RemarkableClient, id::String; kwargs...)
    download(client, Document(ID= id); kwargs...)
end

function Base.download(client::RemarkableClient, doc::Document; kwargs...)
    doc = get_item(client, doc.ID, true)
    @info "Downloading data"
    return HTTP.request(client, "GET", doc.BlobURLGet; kwargs...)
end

function Base.download(client::RemarkableClient, id::String, path_target::String; kwargs...)
    download_document(client, Document(ID=id), path_target; kwargs...)
end

function Base.download(client::RemarkableClient, doc::Document, path_target::String; kwargs...)
    file_name = isempty(title(doc)) ?
                    doc.ID :
                    (ispdf(doc) ?
                        title(doc)[1:end-4] :
                        title(doc))
    file_path = joinpath(path_target, file_name * ".zip")
    body = download(client, doc; kwargs...)
    write(file_path, body)
    @info "File downloaded at $(file_path)"
    return file_path
end

function download_pdf(client::RemarkableClient, doc::Document, path_target::String; kwargs...)
    file_path = download(client, doc, path_target; kwargs...)
    file_name = isempty(title(doc)) ?
                    doc.ID :
                    (ispdf(doc) ?
                        title(doc)[1:end-4] :
                        title(doc))
    z = ZipFile.Reader(file_path)
    for f in z.files
        if endswith(f, ".pdf")
            write(joinpath(path_target, file_name * ".pdf"), f)
            @info "Extracted $(file_name).pdf"
            return joinpath(path_target, file_name * ".pdf")
        end
    end
end

## Upload files

"""
    create_upload_request(client, doc::Document=Document()) -> Document

Create a request to upload a document with a given id.
"""
function create_upload_request(client::RemarkableClient, doc::Document=Document(); kwargs...)
    @info "Creating upload request"
    data = storage_request(client, "PUT", "upload/request", obj_to_dict(doc); kwargs...)
    if isempty(data["BlobURLPut"])
        error("Failed to get URL for upload")
    end
    return Document(doc; BlobURLPut = data["BlobURLPut"])
end

"""
    upload_document!(client, obj::RemarkableObject, zip) -> RemarkableObject

Upload `zip` file (actual zip file) with metadata from obj
"""
function upload_document!(client::RemarkableClient, obj::RemarkableObject, zipfile; kwargs...)
    obj = create_upload_request(client, obj; kwargs...)
    @info "Uploading data"
    body = HTTP.request(client, "PUT", obj.BlobURLPut, Dict(), zipfile; kwargs...)
    item = update_metadata!(client, obj)
    return item
end

"""
    upload_pdf!(client, pdf_path::String, pdf_name::String = basename(pdf_path), parent)

Create a document for the given pdf and upload it.
"""
function upload_pdf!(client::RemarkableClient, pdf_path::String, pdf_name::String = basename(pdf_path), parent::String = ""; kwargs...)
    upload_pdf!(client, read(pdf_path), pdf_name, parent; kwargs...)
end

function upload_pdf!(client::RemarkableClient, pdf, pdfname::String, parent::String=""; kwargs...)
    doc = Document(
        Parent = parent,
        VissibleName = pdfname
    )
    ## Create the ZIP file here with ZipFile
    tmpdir = mktempdir()
    tmpfile = joinpath(tmpdir, doc.ID * ".zip")
    zip = ZipFile.Writer(tmpfile)

    pdffile = ZipFile.addfile(zip, doc.ID * ".pdf")
    write(pdffile, pdf)

    pagedata = ZipFile.addfile(zip, doc.ID * ".pagedata")
    write(pagedata, "")

    content = ZipFile.addfile(zip, doc.ID * ".content")
    json_content = Dict(
        "extraMetaData" => Dict(),
        "fileType" => "pdf",
        "lastOpenedPage" => 0,
        "lineHeight" => -1,
        "margins" => 100,
        "pageCount" => 0,
        "textScale" => 1,
        "transform" => Dict()
    )
    write(content, JSON.json(json_content))

    close(zip)
    # zip = ZipFile.Reader(tmpfile)
    @info "Zip file temporarily saved at $(tmpfile)."
    upload_document!(client, doc, read(joinpath(tmpfile)); kwargs...)
end


text:
using GNSSTools
using Statistics
using ProgressMeter
using JLD
using LinearAlgebra
using Optim
using Distributions
using QuadGK
using PyPlot
pygui(true)
include("kalman_filter_gain_error_functions.jl")

function eval_tracking(CN0s, f_s;
                       iterations=iterations,
                       T=1e-3, M=5, p_fa=1e-7,
                       state_num=3,
                       sigtype=define_l1ca_code_type(),
                       channel="I",
                       h_parms=[1.9e-6, 2e-6, 9.7e-8, 9.1e-11, 9.9e-15] ./ 16.3676e6^2,
                       qₐ=1e-9,
                       σω=1,
                       acquisition_T=20e-3,
                       fine_acq_T=100e-3,
                       t_length=1,
                       dll_b=1,
                       fd_rate=0,
                       f_if=0,
                       P0=missing)
    prn = 1
    fd_range = 2/acquisition_T
    sig_freq = sigtype.sig_freq
    h₀ = h_parms[3]
    h₋₁ = h_parms[2]
    h₋₂ = h_parms[1]
    replica = definereplica(sigtype, f_s, T)
    signal = definesignal(sigtype, f_s, t_length;
                          prn=prn, fd_rate=fd_rate, f_if=f_if,
                          include_phase_noise=true,
                          include_thermal_noise=true,
                          include_adc=true, include_carrier=true,
                          receiver_h_parms=h_parms)
    σᵩ_theo = Array{Float64}(undef, length(CN0s))
    σᵩ_exp = Array{Float64}(undef, length(CN0s))
    σᵩ_exp_std = Array{Float64}(undef, length(CN0s))
    σ_unfilt = Array{Float64}(undef, length(CN0s))
    σ_unfilt_std = Array{Float64}(undef, length(CN0s))
    pᵩs = Array{Float64}(undef, length(CN0s))
    σᵩ²_unfilt_theo = Array{Float64}(undef, length(CN0s))
    pᵩs_steady_state = Array{Float64}(undef, length(CN0s))
    Bs = Array{Float64}(undef, length(CN0s))
    Ps = Array{Float64}(undef, state_num, length(CN0s))
    Ks = Array{Float64}(undef, state_num, length(CN0s))
    P_fulls = Array{Float64}(undef, length(CN0s), state_num, state_num)
    N_over_2 = floor(Int, t_length/T/2)
    p = Progress(iterations*length(CN0s), 1, "Processing...")
    for i in 1:length(CN0s)
        CN0 = CN0s[i]
        CN0_linear = 10^(CN0/10)
        σᵩ² = GNSSTools.phase_noise_variance(CN0, T)
        R = [σᵩ²]
        if ismissing(P0)
            if state_num == 3
                P0 = diagm([σᵩ², (2π/(2*T))^2, (2π*σω)^2])
            elseif state_num == 2
                P0 = diagm([σᵩ², (2π/(2*T))^2])
            else
                @error("Invalid value for state_num.")
            end
        end
        sigma_phis = Array{Float64}(undef, iterations)
        pᵩ = Array{Float64}(undef, iterations)
        P = Array{Float64}(undef, state_num, iterations)
        P_full = Array{Float64}(undef, iterations, state_num, state_num)
        K = Array{Float64}(undef, state_num, iterations)
        sigma_unfiltered = Array{Float64}(undef, iterations)
        for iteration in 1:iterations
            ϕ₀ = rand(0:0.0001:2π)
            f_d = rand(-fd_range:0.1:fd_range)
            if state_num == 3
                x = [ϕ₀, 2π*f_d, 2π*fd_rate] #.+ rand(MvNormal(P0))
                ϕ0_est = x[1]
                fd_est = x[2]/2π
                fd_rate_est = x[3]/2π
            else
                x = [ϕ₀, 2π*f_d] #.+ rand(MvNormal(P0))
                ϕ0_est = x[1]
                fd_est = x[2]/2π
            end
            # code_start_idx = rationalize(rand(1:0.25:replica.sample_num))
            code_start_idx = rand(1:1:replica.sample_num)
            # code_start_idx = 1
            definesignal!(signal; CN0=CN0, new_phase_noise=true,
                          new_thermal_noise=true, phi=ϕ₀, f_d=f_d,
                          code_start_idx=code_start_idx,
                          nADC=16, fd_rate=fd_rate)
            generatesignal!(signal)
            replica = definereplica(sigtype, f_s, T)
            trackresults = trackprn(signal, replica, prn, ϕ0_est,
                                    fd_est,
                                    code_start_idx,
                                    P0, R;
                                    DLL_B=dll_b,
                                    state_num=state_num,
                                    dynamickf=false,
                                    cov_mult=1,
                                    qₐ=qₐ,
                                    q_mult=1,
                                    h₀=h₀, h₋₂=h₋₂,
                                    R_mult=1,
                                    fd_rate=0)
            # acqresults, trackresults, corr_result, SNR_est = process(signal,
            #                             sigtype,
            #                             prn, channel;
            #                             fine_acq_method=:fft,
            #                             return_corrresult=true,
            #                             fd_center=f_d,
            #                             fd_range=fd_range,
            #                             h₀=h₀,
            #                             h₋₂=h₋₂,
            #                             q_mult=1,
            #                             cov_mult=1,
            #                             R_mult=1,
            #                             σᵩ²=σᵩ²,
            #                             σω=σω,
            #                             q_a=qₐ,
            #                             dynamickf=true,
            #                             state_num=state_num,
            #                             dll_b=dll_b,
            #                             acquisition_T=acquisition_T,
            #                             fine_acq_T=fine_acq_T,
            #                             tracking_T=T,
            #                             M=M,
            #                             show_plot=false)
            K[:,iteration] = mean(trackresults.K[:,N_over_2:end], dims=2)
            sigma_phis[iteration] = var(trackresults.dphi_filt[N_over_2:end])
            sigma_unfiltered[iteration] = var(trackresults.dphi_meas[N_over_2:end])
            P[:,iteration] = mean(trackresults.P[:,N_over_2:end], dims=2)
            P_full[iteration,:,:] = mean(trackresults.P_full[N_over_2:end,:,:], dims=1)
            C = calcC(T, state_num)
            pᵩ[iteration] = (C*P_full[iteration,:,:]*C')[1]
            next!(p)
        end
        pᵩs[i] = mean(pᵩ)
        σᵩ_exp[i] = sqrt(mean(sigma_phis))
        σᵩ_exp_std[i] = std(sqrt.(sigma_phis))
        σ_unfilt[i] = sqrt(mean(sigma_unfiltered))
        σ_unfilt_std[i] = std(sqrt.(sigma_unfiltered))
        Ps[:,i] = mean(P, dims=2)
        P_fulls[i,:,:] = mean(P_full, dims=1)
        K_mean = mean(K, dims=2)
        Ks[:,i] = K_mean
        B = get_B(init_H(K_mean...))
        A = calcA(T, state_num)
        C = calcC(T, state_num)
        Q = calcQ(T, h₀, h₋₂, qₐ, sig_freq, state_num)
        R = [σᵩ²]
        σᵩ²_unfilt_theo[i] = σᵩ²
        S = GNSSTools.dare(A', C', Q, GNSSTools.Diagonal(R)')
        # P_steady_state = A*S*A'
        P_steady_state = S
        pᵩ_steady_state = C*P_steady_state*C'
        pᵩs_steady_state[i] = pᵩ_steady_state[1]
        Bs[i] = B
        σᵩ_theo[i] = sqrt(phase_noise_variance(B, CN0, T, h₀, h₋₁, h₋₂))
    end
    return (σᵩ_theo, σᵩ_exp, Ps, Bs, σ_unfilt, Ks, pᵩs, P_fulls, pᵩs_steady_state,
            σᵩ²_unfilt_theo, σᵩ_exp_std, σ_unfilt_std)
end

h_parms = [1.9e-6, 2e-6, 9.7e-8, 9.1e-11, 9.9e-15] ./ 16.3676e6^2  # TCXO
# h_parms = [2.2e-9, 1.6e-10, 1.6e-10, 7.1e-13, 7.2e-16] ./ 10e6^2  # OCXO
h₀ = h_parms[3]
h₋₁ = h_parms[2]
h₋₂ = h_parms[1]

channel = "I"
sigtype = define_l1ca_code_type()
# Ts = [1e-3 2e-3 5e-3 10e-3]
Ts = [1e-3 10e-3]
f_s = 2.5e6
name = "l1ca"

# channel = "Q"
# sigtype = define_l5_code_type(;channel=channel)
# Ts = [1e-3 20e-3]
# f_s = 25e6
# name = "l5"

CN0s = Array(range(20, 50, step=1))
iterations = 25
p_fa = 1e-7
qₐ = 0.1
state_num = 2
σω = 1
acquisition_T=40e-3
fine_acq_T=100e-3
M=2
t_length=1
dll_b=1
fd_rate = 0

directory = "/home/sjbilardi/Dropbox/Apps/Overleaf/MastersThesis/"
file_name = "ch4_signal_tracking_performance_$(name)"
file = string(directory, file_name, ".jld")

σ_theoretical = []
σ_experimental = []
σ_unfiltered = []
Ks = []
Ps = []
pᵩs = []
Bs = []
P_fulls = []
pᵩs_steady_state = []
σᵩ²_unfilt_theo = []
σᵩ_exp_std = []
σ_unfilt_std = []

for i in 1:length(Ts)
    σ_theo, σ_exp, P, B, σ_unfilt, K, pᵩ, P_full,
    pᵩ_steady_state, σᵩ²_unfilt_theo_, σᵩ_exp_std_, σ_unfilt_std_  = eval_tracking(CN0s, f_s;
                                                     iterations=iterations,
                                                     T=Ts[i], p_fa=p_fa,
                                                     state_num=state_num,
                                                     sigtype=sigtype,
                                                     channel=channel,
                                                     qₐ=qₐ,
                                                     σω=σω,
                                                     acquisition_T=acquisition_T,
                                                     fine_acq_T=fine_acq_T,
                                                     M=M,
                                                     t_length=t_length,
                                                     dll_b=dll_b,
                                                     fd_rate=fd_rate)
    push!(σ_theoretical, σ_theo)
    push!(σ_experimental, σ_exp)
    push!(Ps, P)
    push!(Bs, B)
    push!(Ks, K)
    push!(σ_unfiltered, σ_unfilt)
    push!(pᵩs, pᵩ)
    push!(P_fulls, P_full)
    push!(pᵩs_steady_state, pᵩ_steady_state)
    push!(σᵩ²_unfilt_theo, σᵩ²_unfilt_theo_)
    push!(σᵩ_exp_std, σᵩ_exp_std_)
    push!(σ_unfilt_std, σ_unfilt_std_)
    save(file,
         "σ_theoretical", σ_theoretical,
         "σ_experimental", σ_experimental,
         "Ps", Ps,
         "Bs", Bs,
         "Ts", Ts,
         "CN0s", CN0s,
         "iterations", iterations,
         "p_fa", p_fa,
         "f_s", f_s,
         "M", M,
         "h_parms", h_parms,
         "sig_freq", sigtype.sig_freq,
         "q_a", qₐ,
         "state_num", state_num,
         "σω", σω,
         "acquisition_T", acquisition_T,
         "fine_acq_T", fine_acq_T,
         "t_length", t_length,
         "dll_b", dll_b,
         "Ks", Ks,
         "σ_unfiltered", σ_unfiltered,
         "pᵩs", pᵩs,
         "P_fulls", P_fulls,
         "pᵩs_steady_state", pᵩs_steady_state,
         "fd_rate", fd_rate,
         "σᵩ²_unfilt_theo", σᵩ²_unfilt_theo,
         "σᵩ_exp_std", σᵩ_exp_std,
         "σ_unfilt_std", σ_unfilt_std)
end

directory = "/home/sjbilardi/Dropbox/Apps/Overleaf/MastersThesis/"
file_names = [string(directory, "ch4_signal_tracking_performance_l1ca.jld"),
              string(directory, "ch4_signal_tracking_performance_l5.jld")]
cmap = get_cmap("viridis")

fig = figure(figsize=(7.5, 3))
for f in 1:1
# for f in 1:length(file_names)
    file_name = file_names[f]
    ax = fig.add_subplot(1, 2, f)
    CN0s,σ_theoretical,σ_experimental, Ts, σ_unfiltered, σᵩ_exp_std, σ_unfilt_std = load(file_name,
                                             "CN0s",
                                             "σ_theoretical",
                                             "σ_experimental",
                                             "Ts",
                                             "σ_unfiltered",
                                             "σᵩ_exp_std",
                                             "σ_unfilt_std")
    σs_to_plot = length(σ_theoretical)
    σ_to_plot = Array(1:σs_to_plot)
    for i in 1:σs_to_plot
        if i in σ_to_plot
            color = cmap(float(i)/σs_to_plot)
            # label = "T = $(floor(Int, Ts[i]*1e3))ms, f_s = $(round(f_s*1e-6, digits=1))MHz"
            label = "T = $(floor(Int, Ts[i]*1e3))ms"
            # ax.plot(CN0s, sqrt.(pᵩs_steady_state[i]).*(180/π), "-", color=color,
                    # label=label)
            # ax.plot(CN0s, σ_theoretical[i].*(180/π), ":", color=color)
            # ax.plot(CN0s, σ_experimental[i].*(180/π), ":", color=color)
            # ax.plot(CN0s, σ_experimental[i].*(180/π), "+", color=color)
            # ax.plot(CN0s, sqrt.(σᵩ²_unfilt_theo[i]).*(180/π), "-", color=color)
            # ax.plot(CN0s, sqrt.(Ps[i][1,:]).*(180/π), ":", color=color)
            # ax.plot(CN0s, sqrt.(Ps[i][1,:]).*(180/π), ".", color=color)
            # ax.plot(CN0s, sqrt.(pᵩs[i]).*(180/π), ".", color=color)
            ax.plot(CN0s, sqrt.(GNSSTools.phase_noise_variance.(CN0s, Ts[i])).*(180/π),
                    "-", color=color)
            ax.errorbar(CN0s, σ_unfiltered[i].*(180/π), yerr=σ_unfilt_std[i],
                        marker=".", color=color)
            ax.plot([], [], ".-", color=color,
                    label="T=$(floor(Int, Ts[i]*1000))ms")
        end
    end
    ax.set_yscale("log")
    # hlines(y=10log10(15), xmin=CN0s[1], xmax=CN0s[end], color="grey",
        #    linestyle=":", label=string(L"15^\circ", " Tracking Limit"))
    legend()
    xlabel("C/N₀ (dB⋅Hz)")
    ylim([0.0, 60])
    ylabel("Unfiltered σᵩ (degrees)")
    letter = ('a' + (f-1))
    title("($(letter))")
end
subplots_adjust(wspace=0.3, bottom=0.15, left=0.08, right=0.92, top=0.92)
# savefig(string(directory, "figures/ch4_signal_tracking_performance_unfiltered_sigma-phi.pdf"), dpi=300)

fig = figure(figsize=(7.5, 3))
for f in 1:1
# for f in 1:length(file_names)
    file_name = file_names[f]
    ax = fig.add_subplot(1, 2, f)
    CN0s,σ_theoretical,σ_experimental, Ts, σ_unfiltered,
    σᵩ_exp_std, σ_unfilt_std, pᵩs_steady_state = load(file_name,
                                             "CN0s",
                                             "σ_theoretical",
                                             "σ_experimental",
                                             "Ts",
                                             "σ_unfiltered",
                                             "σᵩ_exp_std",
                                             "σ_unfilt_std",
                                             "pᵩs_steady_state")
    σs_to_plot = length(σ_theoretical)
    σ_to_plot = Array(1:σs_to_plot)
    for i in 1:σs_to_plot
        if i in σ_to_plot
            color = cmap(float(i)/σs_to_plot)
            # label = "T = $(floor(Int, Ts[i]*1e3))ms, f_s = $(round(f_s*1e-6, digits=1))MHz"
            label = "T = $(floor(Int, Ts[i]*1e3))ms"
            # ax.plot(CN0s, sqrt.(pᵩs_steady_state[i]).*(180/π), "-", color=color,
                    # label=label)
            ax.plot(CN0s, sqrt.(pᵩs_steady_state[i]).*(180/π), "-", color=color)
            # ax.plot(CN0s, sqrt.(GNSSTools.phase_noise_variance.(CN0s, Ts[i])).*(180/π), ":", color=color)
            # ax.plot(CN0s, σ_experimental[i].*(180/π), ":", color=color)
            ax.errorbar(CN0s, σ_experimental[i].*(180/π), yerr=σᵩ_exp_std[i],
                        marker=".", color=color)
            # ax.plot(CN0s, σ_unfiltered[i].*(180/π), "x", color=color)
            # ax.plot(CN0s, sqrt.(σᵩ²_unfilt_theo[i]).*(180/π), "-", color=color)
            # ax.plot(CN0s, sqrt.(Ps[i][1,:]).*(180/π), ":", color=color)
            # ax.plot(CN0s, sqrt.(Ps[i][1,:]).*(180/π), ".", color=color)
            # ax.plot(CN0s, sqrt.(pᵩs[i]).*(180/π), ".", color=color)
            ax.plot([], [], ".-", color=color,
                    label="T=$(floor(Int, Ts[i]*1000))ms")
        end
    end
    ax.set_yscale("log")
    hlines(y=10log10(15), xmin=CN0s[1], xmax=CN0s[end], color="grey",
           linestyle=":", label=string(L"15^\circ", " Tracking Limit"))
    legend()
    xlabel("C/N₀ (dB⋅Hz)")
    ylim([0.0, 60])
    ylabel("Filtered σᵩ (degrees)")
    letter = ('a' + (f-1))
    title("($(letter))")
end
subplots_adjust(wspace=0.3, bottom=0.15, left=0.08, right=0.92, top=0.92)
# savefig(string(directory, "figures/ch4_signal_tracking_performance_filtered_sigma-phi.pdf"), dpi=300)

ext: jl


text:
# For vectors
function permute_dims_expr(perm, D, C, B, R)
    s = Expr(:tuple) # size
    x = Expr(:tuple) # stride
    o = Expr(:tuple) # offsets
    Rnew = Expr(:tuple) # rank
    Dnew = Expr(:tuple) # dense
    Cnew = -1
    Bnew = -1
    N = length(perm)
    for n ∈ 1:N
        p = perm[n]
        push!(s.args, Expr(:ref, :s, p))
        push!(x.args, Expr(:ref, :x, p))
        push!(o.args, Expr(:ref, :o, p))
        push!(Rnew.args, R[p])
        push!(Dnew.args, D[p])
        if C == p
            Cnew = n
        end
        if B == p
            Bnew = n
        end
    end
    Dnew, Cnew, Bnew, Rnew, s, x, o
end
@generated function Base.permutedims(A::PtrArray{S,D,T,N,C,B,R}, ::Val{P}) where {S,D,T,N,C,B,R,P}
    Dnew, Cnew, Bnew, Rnew, s, x, o = permute_dims_expr(P, D, C, B, R)
    quote
        $(Expr(:meta,:inline))
        s = A.size
        ptr = A.ptr
        x = ptr.strd
        o = ptr.offsets
        sp = StridedPointer{$T,$N,$Cnew,$Bnew,$Rnew}(ptr.p, $x, $o)
        PtrArray(sp, $s, Val{$Dnew}())
    end
end

@inline Base.adjoint(A::AbstractStrideMatrix) = permutedims(A, Val{(2,1)}())
@inline Base.transpose(A::AbstractStrideMatrix) = permutedims(A, Val{(2,1)}())

@generated function Base.adjoint(a::PtrArray{S,D,T,1,C,B,R,X,O}) where {S,D,T,C,B,R,X,O}
    s = Expr(:tuple, :(One()), Expr(:ref, :s, 1))
    x₁ = Expr(:ref, :x, 1)
    x = Expr(:tuple, x₁, x₁)
    o = Expr(:tuple, :(One()), Expr(:ref, :o, 1))
    R1 = R[1]
    Rnew = Expr(:tuple, R1+1, R1)
    Dnew = Expr(:tuple, true, D[1])
    Cnew = C == 1 ? 2 : C
    quote
        $(Expr(:meta,:inline))
        s = a.size
        ptr = a.ptr
        x = ptr.strd
        o = ptr.offsets
        sp = StridedPointer{$T,2,$Cnew,$B,$Rnew}(ptr.p, $x, $o)
        PtrArray(sp, $s, Val{$Dnew}())
    end
end

@inline Base.transpose(a::AbstractStrideVector) = adjoint(a)


text:
"""
    AbstractModulus

  Assuming α^n = p[1] + p[2]α + ... + p[n]α^(n-1), create expressions for α^k, k=n, ... ,2n-2 and powers 0,...,n-1of p's companion matrix.
  """
abstract AbstractModulus{I<:Integer}

function modulusfields{I}(p::Array{I,1})
  n=length(p)
  polynomials = Array(Array{I,1},n-1)
  polynomials[1]=p
  A = zeros(I,(n,n))
  A[:,n]=p
  if n>1
    for i in 2:(n-1)
      polynomials[i] = vcat([0],polynomials[i-1][1:(n-1)])+(polynomials[i-1][n])*p
      A[i,i-1]=1
    end
    A[n,n-1]=1
  end
  matrices = Array(Array{I,2},n)
  for i in 1:n
    matrices[i]=A^(i-1)
  end
  return polynomials, matrices
end

immutable Modulus{I<:Integer} <: AbstractModulus{I}
  polynomials::Array{Array{I,1},1}
  matrices::Array{Array{I,2},1}
  symbol::String

  function (::Type{Modulus}){I}(p::Array{I,1}; symbol::String="α")
    n = length(p)
    n > 1 || error("n==1, a trivial case, is not supported")
    polynomials,matrices = modulusfields(p)
    new{I}(polynomials,matrices,symbol)
  end

end

immutable KummerModulus{I<:Integer} <: AbstractModulus{I}
  polynomials::Array{Array{I,1},1}
  matrices::Array{Array{I,2},1}
  symbol::String

  function (::Type{KummerModulus}){I}(n::I; symbol::String="ζ")
    n > 1 || error("n==$n, a trivial case, is not supported")
    p = cyclotomic(n)
    polynomials,matrices = modulusfields(-p[1:(end-1)])
    new{I}(polynomials,matrices,symbol)
  end

end

immutable QIModulus{I<:Integer} <: AbstractModulus{I}
  polynomials::Array{Array{I,1},1}
  matrices::Array{Array{I,2},1}
  symbol::String
  d::I
  k::I

  function (::Type{QIModulus}){I}(d::I; symbol::String="α")
    issquarefree(d) || error("$d is not square free")
    if mod(d,4) == 1
      # adjoin (1+sqrt(d))/2
      k = I((d-1)/4)
      polynomials, matrices = modulusfields([k,d])
    else
      # mod(d,4) is either 2 or 3, adjoin sqrt(d)
      polynomials, matrices = modulusfields([d,I(0)])
      k = I(0)
    end
    new{I}(polynomials, matrices, symbol, d, k)
  end

end

immutable Element{I<:Integer}
  coefficients::Array{I,1}
  modulus::AbstractModulus{I}

  function (::Type{Element}){I}(coef::Array{I,1}, modulus::AbstractModulus{I})
     length(coef)==length(modulus.polynomials[1]) || error("unequal degrees ", length(coef), " and ", length(modulus.polynomials[1]),".")
    new{I}(coef,modulus)
  end

end

## Conversion

function Matrix{I}(z::Element{I})
  p = z.coefficients
  n = length(p)
  matrices = z.modulus.matrices
  ans = p[1]*matrices[1]
  for i in 2:n
    ans += p[i]*matrices[i]
  end
  ans
end

## Display

function helper{I<:Integer}(p::Array{I,1}, symbol::String)
  tmp=""
  for i in eachindex(p)
    if p[i] != 0
      sgn = p[i] < 0 ? "-" : "+"
      if i==1
        tmp=tmp*"$(p[1])"
      elseif i==2
        tmp=tmp*"$(sgn)$(abs(p[i]))$(symbol)"
      else
        tmp=tmp*"$(sgn)$(abs(p[i]))$(symbol)^$(i-1)"
      end
    end
  end
  return length(tmp) == 0 ? "0" : tmp
end

function string{I<:Integer}(z::Element{I})
  helper(z.coefficients, z.modulus.symbol)
end

function string{I<:Integer}(m::AbstractModulus{I})
  "$(m.symbol)^$(length(m.polynomials[1])) = $(helper(m.polynomials[1],m.symbol))"
end

function show{I<:Integer}(io::IO, z::Element{I})
  print(io, "$(typeof(z))\n")
  print(io, string(z))
end

function show{I<:Integer}(io::IO, m::AbstractModulus{I})
  print(io, "$(typeof(m))\n")
  print(io, string(m))
end

max_stars_repo_name: WilCrofter/ZAdjoin

max_issues_repo_name: WilCrofter/ZAdjoin

max_forks_repo_name: WilCrofter/ZAdjoin


text:
@testset "Flux reports" begin

@testset "function nodeFluxReport" begin
    @test isempty( nodeFluxReport( mpSim, [-6.0, -12.0], :in, "active",
        "A junior", "Career", "Blort" ) )
    @test isempty( nodeFluxReport( mpSim, [12.0, 8.0, 0.0], :bang, "active",
        "A junior", "Career", "Blort" ) )
    @test isempty( nodeFluxReport( mpSim, [12.0, 8.0, 0.0], :in, "Blort" ) )
    @test isempty( nodeFluxReport( mpSim, -6, :in, "active",
        "A junior", "Career", "Blort" ) )

    report = nodeFluxReport( mpSim, 12, :in, "active", "A junior", "Career",
        "Blort" )
    @test !haskey( report, "Blort" ) && all( haskey.( Ref( report ), ["active",
        "Career", "A junior"] ) )
    @test all( report["active"][:, "external => active"] .== 30 )
    @test all( report["A junior"][:, "other => A junior"] .== 10 )
    @test all( report["Career"][:, "other => Career"] .==
        vcat( 0, 0, 24, 24, fill( 16, 22 ) ) )

    report = nodeFluxReport( mpSim, 12, :out, "active", "A junior", "Career",
        "Blort" )
    @test all( report["active"][:, "active => external"] .==
        vcat( 0, 0, 6, 6, fill( 14, 6 ), 38, 38, fill( 30, 14 ) ) )
    @test all( report["A junior"][:, "A junior => other"] .==
        vcat( 0, 0, fill( 10, 24 ) ) )
    @test all( report["Career"][:, "Career => other"] .==
        vcat( fill( 0, 10 ), 24, 24, fill( 16, 14 ) ) )

    report = nodeFluxReport( mpSim, 12, :within, "active", "A junior", "Career",
        "Blort" )
    @test all( report["active"][:, "within active"] .== [0, 0,
        24, 24, 16, 26, 26, 21, 16, 16, 26, 26, 21, 16, 16, 26, 26, 21, 16, 16,
        26, 26, 21, 16, 16, 26] )
end  # @testset "function nodeFluxReport"

@testset "function transitionFluxReport" begin
    @test isempty( transitionFluxReport( mpSim, [-12.0, -8.0], "Promotion",
        ("A junior", "A senior"), ("Reserve", "Reserve junior", "B senior") ) )
    @test isempty( transitionFluxReport( mpSim, [12.0, 24.0], "Blook",
        ("Foo", "Bar"), ("A senior", "B junior"), ("Blook", "Foo", "Bar"),
        ("Reserve", "A junior", "A senior") ) )
    @test isempty( transitionFluxReport( mpSim, -12.0, "Promotion",
        ("A junior", "A senior"), ("Reserve", "Reserve junior", "B senior") ) )
    report = transitionFluxReport( mpSim, 12, "attrition", "Promotion", "EW",
        "B-" )
    @test all( report[:, "attrition"] .== 0 )
    @test all( report[:, "Promotion"] .== [0, 0, 16, 16, 16, 26, 26, 21, 16,
        16, 26, 26, 21, 16, 16, 26, 26, 21, 16, 16, 26, 26, 21, 16, 16, 26] )
    @test all( report[:, "EW"] .== 30 )
    @test all( report[:, "B-"] .==
        vcat( 0, 0, 6, 6, fill( 14, 22 ) ) )

    report = transitionFluxReport( mpSim, 12, ("external", "A junior"),
        ("B senior", ""), ("Reserve junior", "A senior") )
    @test all( report[:, "external => A junior"] .== 10 )
    @test all( report[:, "B senior => external"] .==
        vcat( fill( 0, 10 ), 7, 7, 8, 8, 8, 3, 3, 8, 8, 8, 3, 3, 8, 8, 8, 3 ) )
    @test all( report[:, "Reserve junior => A senior"] .==
        vcat( 0, 0, 4, 4, fill( 0, 22 ) ) )

    report = transitionFluxReport( mpSim, 12, ("attrition", "Master", "out"),
        ("PE", "A senior", ""), ("EW", "OUT", "A junior"),
        ("Reserve", "Reserve junior", "B senior") )
    @test all( report[:, "attrition: Master => external"] .== 0 )
    @test all( report[:, "PE: A senior => external"] .==
        vcat( fill( 0, 10 ), 7, 7, 3, 8, 8, 3, 3, 3, 8, 8, 3, 3, 3, 8, 8, 3 ) )
    @test all( report[:, "EW: external => A junior"] .== 10 )
    @test all( report[:, "Reserve: Reserve junior => B senior"] .==
        vcat( 0, 0, 4, 4, fill( 0, 22 ) ) )
end  # @testset "function transitionFluxReport"

end  # @testset "Flux reports"


text:
@testset "Monomial Vector" begin
    Mod.@polyvar x y
    @test x > y
    @test x^2 > y^2
    X = [x^2, x*y, y^2]
    @test isempty(@inferred monomials((x, y), 1:0))
    for (i, m) in enumerate(monomials((x, y), 2))
        @test m == X[i]
    end
    for (i, m) in enumerate(monomials((x, y), 2:2))
        @test m == X[i]
    end
    X = [x^2, y^2]
    for (i, m) in enumerate(monomials((x, y), 2, m -> m != x*y))
        @test m == X[i]
    end
    @test (@inferred monovectype([1, x])) <: AbstractArray{<:AbstractMonomial}
    @test (@inferred monovectype([x])) <: AbstractArray{<:AbstractMonomial}
    @test (@inferred monovec([1, x])) isa monovectype([1, x])
    @test (@inferred monovec([x])) isa monovectype([x])
    @test (@inferred monovec([1, 2], [1, x]))[2] isa AbstractArray{<:AbstractMonomial}
    @test (@inferred monovec([1], [x]))[2] isa AbstractArray{<:AbstractMonomial}
    @test length(monovec([y, x])) == 2
    X = monovec([x, 1, x*y])
    @test X == collect(X)
    @test nvariables(X) == 2
    @test variables(X)[1] == x
    @test variables(X)[2] == y
    @test X[2:3][1] == x
    @test X[2:3][2] == 1
    @test monovec(X[[3, 2]])[1] == x
    @test monovec(X[[3, 2]])[2] == 1
    # Documentation examples
    @test monovec([x*y, x, x*y, x^2*y, x]) == [x^2*y, x*y, x]
    @test monovectype([x*y, x, 1, x^2*y, x]) <: AbstractVector{typeof(x*y)}
    @test monovectype([x*y, x, x*y, x^2*y, x]) <: AbstractVector
    σ, smv = sortmonovec([x*y, x, x*y, x^2*y, x])
    @test smv == [x^2*y, x*y, x]
    @test σ[1] == 4
    @test σ[2] in (1, 3)
    @test σ[3] in (2, 5)
    @test mergemonovec([[x*y, x, x*y], [x^2*y, x]]) == [x^2*y, x*y, x]
    @test_throws ArgumentError monovec([1, 2], [x^2])
    σ, X = sortmonovec((y, x))
    @test σ == [2, 1]
    @test X == [x, y]
    @test monomialtype([x, y]) <: AbstractMonomial
    @test monomialtype([x^2, 1]) <: AbstractMonomial
    @test monomialtype([x*y, x+y]) <: AbstractMonomial

    @test monovec([x, x^2]) != monovec([x*y, x^2*y])
    @test monomials(x, 1:3) == monovec([x^3, x, x^2])
    @test monomials((x, y), 2) != monomials((x, y), 1)
end


text:
# MIT License
#
# Copyright (c) 2020 Joshua E Gentry

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
using .Pipeline

# =========================================================================================
"""
    unlv_copy(
        pipe::TessPipeline,
        filename::AbstractString
    )::Bool

Generate a UNLV test file from the pipeline and save it to the specified file.  Returns
`false` if there is a problem adding the UNLV generator to the output.

__Arguments:__

| T | Name     | Default | Description
|:--| :------- | :------ | :----------
| R | pipe     |         | The pipline to collect the UNLV text from.
| R | filename |         | The file to write the UNLV text to.

__Details:__

If the file exists it will be overwritten. The file will be in the Latin-1 encoding which
is what Tesseract generates by default.
"""
function unlv_copy(
            pipe::TessPipeline,
            filename::AbstractString,
        )
    local dir         = dirname(filename)
    local tmpfile, io = mktemp(dir;cleanup=false); close(io)
    local tmpname     = basename(tmpfile)

    # -------------------------------------------------------------------------------------
    # Create the text renderer.
    local ptr = pipeline_create_unlv(tmpfile)

    rm(tmpfile)

    if ptr == C_NULL
        return false
    end

    # -------------------------------------------------------------------------------------
    # Create the task to generate the output string.
    local task = pipeline_start_task(
        joinpath(dir, "$(tmpname).unlv"),
        filename
    )

    pipeline_add_renderer(pipe, ptr, task, RENDERER_UNLV)

    return true
end

# =========================================================================================
"""
    unlv_transcode(
        pipe::TessPipeline,
        filename::AbstractString
    )::Bool

Generate a UNLV test file from the pipeline and save it to the specified file.  Returns
`false` if there is a problem adding the UNLV generator to the output.

__Arguments:__

| T | Name     | Default | Description
|:--| :------- | :------ | :----------
| R | pipe     |         | The pipline to collect the UNLV text from.
| R | filename |         | The file to write the UNLV text to.

__Details:__

If the file exists it will be overwritten. This will transcode the file into UTF-8 to match
the rest of Julia.
"""
function unlv_transcode(
            pipe::TessPipeline,
            filename::AbstractString,
        )
    local path = mktempdir(;cleanup=false)

    # -------------------------------------------------------------------------------------
    # Create the text renderer.
    local ptr = pipeline_create_unlv(
        joinpath(path, "output")
    )

    if ptr == C_NULL
        rm(path;recursive = true)
        return false
    end

    # -------------------------------------------------------------------------------------
    # Create the task to generate the output string.
    local task = pipeline_start_task(path, "output.unlv", filename)

    pipeline_add_renderer(pipe, ptr, task, RENDERER_UNLV)

    return true
end

# =========================================================================================
"""
    tess_pipeline_unlv(
        pipe::TessPipeline,
        filename::AbstractString,
        utf8::Bool = true
    )::Bool

Generate a UNLV test file from the pipeline and save it to the specified file.  Returns
`false` if there is a problem adding the UNLV generator to the output.

__Arguments:__

| T | Name     | Default | Description
|:--| :------- | :------ | :----------
| R | pipe     |         | The pipline to collect the UNLV text from.
| R | filename |         | The file to write the UNLV text to.
| O | utf8     | `true`  | Should the output be transcoded to UTF-8?

__Details:__

If the file exists it will be overwritten.  The Tesseract library outputs a file in Latin-1
encoding (even through all other formats are UTF-8).  We can use Julia to convert the file
into UTF-8 to match the other encodings if desired.

__Examples:__

```jldoctest
using Tesseract

# Generate some pages to load.
write("page01.tiff", sample_tiff())
write("page02.tiff", sample_tiff())
write("page03.tiff", sample_tiff())

download_languages() # Make sure we have the data files.

instance = TessInst()
pipeline = TessPipeline(instance)

tess_pipeline_unlv(pipeline, "My Book.unlv")

tess_run_pipeline(pipeline, "My First Book") do add
    add(pix_read("page01.tiff"), 72)
    add(pix_read("page02.tiff"), 72)
    add(pix_read("page03.tiff"), 72)
end

for line in readlines("My Book.unlv")[1:10]
    println(line)
end

# output

No one would have believed in the last years of the
the nineteenth century that this world was being watched
watched keenly and closely by intelligences greater than
than man's and yet as mortal as his own; that as men busied
busied themselves about their various concerns they were
were scrutinised and studied, perhaps almost as narrowly as
as a man with a microscope might scrutinise the transient
transient creatures that swarm and multiply in a drop of
of water. With infinite complacency men went to and fro over
over this globe about their little affairs, serene in their
```

See also: [`tess_run_pipeline`](@ref), [`tess_pipeline_word_box`](@ref),
          [`tess_pipeline_lstm_box`](@ref), [`tess_pipeline_unlv_latin1`](@ref)
"""
function tess_pipeline_unlv(
            pipe::TessPipeline,
            filename::AbstractString,
            utf8::Bool = true
        )::Bool
    local result = false

    if utf8 == true
        result = unlv_transcode(pipe, filename)
    else
        result = unlv_copy(pipe, filename)
    end

    return result
end

# =========================================================================================
"""
    tess_pipeline_unlv(
        pipe::TessPipeline
    )::Union{TessOutput{String}, Nothing}

Generate an UNLV file from the pipeline and save it to a string.  Returns `nothing` if
there is a problem adding the UNLV generator to the output.

__Arguments:__

| T | Name | Default | Description
|:--| :--- | :------ | :----------
| R | pipe |         | The pipline to collect the UNLV text from.

__Details:__

The Tesseract generates Latin-1 text however this function will transcode it to UTF-8 to
interact with Julia correctly.

__Examples:__

```jldoctest
using Tesseract

# Generate some pages to load.
write("page01.tiff", sample_tiff())
write("page02.tiff", sample_tiff())
write("page03.tiff", sample_tiff())

download_languages() # Make sure we have the data files.

instance = TessInst()
pipeline = TessPipeline(instance)

book = tess_pipeline_unlv(pipeline)

tess_run_pipeline(pipeline, "My First Book") do add
    add(pix_read("page01.tiff"), 72)
    add(pix_read("page02.tiff"), 72)
    add(pix_read("page03.tiff"), 72)
end

for line in split(book[], "\\n")[1:10]
    println(line)
end

# output

No one would have believed in the last years of the
the nineteenth century that this world was being watched
watched keenly and closely by intelligences greater than
than man's and yet as mortal as his own; that as men busied
busied themselves about their various concerns they were
were scrutinised and studied, perhaps almost as narrowly as
as a man with a microscope might scrutinise the transient
transient creatures that swarm and multiply in a drop of
of water. With infinite complacency men went to and fro over
over this globe about their little affairs, serene in their
```

See also: [`tess_run_pipeline`](@ref), [`tess_pipeline_word_box`](@ref),
          [`tess_pipeline_lstm_box`](@ref), [`tess_pipeline_unlv_latin1`](@ref)
"""
function tess_pipeline_unlv(
            pipe::TessPipeline
        )::Union{TessOutput{String}, Nothing}
    local path = mktempdir(;cleanup=false)

    # -------------------------------------------------------------------------------------
    # Create the text renderer.
    local ptr = pipeline_create_unlv(
        joinpath(path, "output")
    )

    if ptr == C_NULL
        rm(path;recursive = true)
        return nothing
    end

    # -------------------------------------------------------------------------------------
    # Create the task to generate the output string.
    local out = TessOutput(String)
    local task = pipeline_start_task(path, "output.unlv", out, true)

    pipeline_add_renderer(pipe, ptr, task, RENDERER_UNLV)

    return out
end

# =========================================================================================
"""
    tess_pipeline_unlv_latin1(
        pipe::TessPipeline
    )::Union{TessOutput{Vector{UInt8}}, Nothing}

Generate an UNLV file from the pipeline and save it to a byte array.  Returns `nothing` if
there is a problem adding the UNLV generator to the output.

__Arguments:__

| T | Name | Default | Description
|:--| :--- | :------ | :----------
| R | pipe |         | The pipline to collect the UNLV text from.

__Details:__

The returned bytes will be in Latin-1 encoding, so a conversion would be required before
using the bytes as a string in Julia.  [`tess_pipeline_unlv`](@ref) will provide you with
the same data but already encoded in UTF-8 for Julia.

__Examples:__

```jldoctest
using Tesseract

# Generate some pages to load.
write("page01.tiff", sample_tiff())
write("page02.tiff", sample_tiff())
write("page03.tiff", sample_tiff())

download_languages() # Make sure we have the data files.

instance = TessInst()
pipeline = TessPipeline(instance)

bytes = tess_pipeline_unlv_latin1(pipeline)

tess_run_pipeline(pipeline, "My First Book") do add
    add(pix_read("page01.tiff"), 72)
    add(pix_read("page02.tiff"), 72)
    add(pix_read("page03.tiff"), 72)
end

println(string("Book generated: ", length(bytes[]), " bytes."))

# output

Book generated: 4410 bytes.
```

See also: [`tess_run_pipeline`](@ref), [`tess_pipeline_word_box`](@ref),
          [`tess_pipeline_lstm_box`](@ref), [`tess_pipeline_unlv`](@ref)
"""
function tess_pipeline_unlv_latin1(
            pipe::TessPipeline
        )::Union{TessOutput{Vector{UInt8}}, Nothing}
    local path = mktempdir(;cleanup=false)

    # -------------------------------------------------------------------------------------
    # Create the text renderer.
    local ptr = pipeline_create_unlv(
        joinpath(path, "output")
    )

    if ptr == C_NULL
        rm(path;recursive = true)
        return nothing
    end

    # -------------------------------------------------------------------------------------
    # Create the task to generate the output string.
    local out = TessOutput(Vector{UInt8})
    local task = pipeline_start_task(path, "output.unlv", out)

    pipeline_add_renderer(pipe, ptr, task, RENDERER_UNLV)

    return out
end

# =========================================================================================
"""
    tess_pipeline_unlv(
        func::Function,
        pipe::TessPipeline
    )::Bool

Generate an UNLV file from the pipeline and pass it back to the client via a callback.
Returns `false` if there is a problem adding the UNLV generator to the output.

__Arguments:__

| T | Name | Default | Description
|:--| :--- | :------ | :----------
| R | func |         | The function to call with the lines of text.
| R | pipe |         | The pipline to collect the text from.

__Details:__

The text will be passed to the caller one line at a time. The "\\n" line terminator will be
included with the text.  Tesseract generates UNLV text in Latin-1, this method will
transcode it to UTF-8 to work with Julia.

__Examples:__

```jldoctest
using Tesseract

# Generate some pages to load.
write("page01.tiff", sample_tiff())
write("page02.tiff", sample_tiff())
write("page03.tiff", sample_tiff())

download_languages() # Make sure we have the data files.

instance = TessInst()
pipeline = TessPipeline(instance)

count = 0
tess_pipeline_unlv(pipeline) do line
    global count
    if count < 10
        print(line)
    end
    count += 1
end

tess_run_pipeline(pipeline, "My First Book") do add
    add(pix_read("page01.tiff"), 72)
    add(pix_read("page02.tiff"), 72)
    add(pix_read("page03.tiff"), 72)
end

# output

No one would have believed in the last years of the
the nineteenth century that this world was being watched
watched keenly and closely by intelligences greater than
than man's and yet as mortal as his own; that as men busied
busied themselves about their various concerns they were
were scrutinised and studied, perhaps almost as narrowly as
as a man with a microscope might scrutinise the transient
transient creatures that swarm and multiply in a drop of
of water. With infinite complacency men went to and fro over
over this globe about their little affairs, serene in their
true
```

See also: [`tess_run_pipeline`](@ref), [`tess_pipeline_word_box`](@ref),
          [`tess_pipeline_lstm_box`](@ref), [`tess_pipeline_unlv_latin1`](@ref)
"""
function tess_pipeline_unlv(
            func::Function,
            pipe::TessPipeline
        )::Bool
    local path = mktempdir(;cleanup=false)

    # -------------------------------------------------------------------------------------
    # Create the text renderer.
    local ptr = pipeline_create_unlv(
        joinpath(path, "output")
    )

    if ptr == C_NULL
        rm(path;recursive = true)
        return false
    end

    # -------------------------------------------------------------------------------------
    # Create the task to generate the output string.
    local task = pipeline_start_task(path, "output.unlv", func, "", true)

    pipeline_add_renderer(pipe, ptr, task, RENDERER_UNLV)

    return true
end


text:
module SubHunt

# package code goes here

using Random
using Printf
using LinearAlgebra
using Parameters
using StaticArrays
using Distributions
using POMDPs
using POMDPModelTools
using RecipesBase
using ParticleFilters

import POMDPs: gen, support, discount, isterminal
import POMDPs: actions, actionindex, action, dimensions
import POMDPs: states, stateindex, transition
import POMDPs: observations, observation, obsindex
import POMDPs: initialstate, initialobs
import POMDPs: updater, update
import POMDPs: reward
import POMDPs: convert_s, convert_a, convert_o

export
    SubHuntPOMDP,
    SubState,
    SubVis,

    DSubHuntPOMDP,

    PingFirst

const Vec8 = SVector{8, Float64}
const Pos = SVector{2, Int}

struct SubState
    own::Pos
    target::Pos
    goal::Int
    aware::Bool
end

@with_kw struct SubHuntPOMDP <: POMDP{SubState, Int, Vec8}
    size::Int                       = 20
    ownspeed::Int                   = 3
    kill_radius::Float64            = 2.0
    p_aware_kill::Float64           = 0.6
    passive_detect_radius::Float64  = 3.0
    passive_std::Float64            = 5.0
    passive_detected_std::Float64   = 0.5
    active_std::Float64             = 0.5
    discount::Float64               = 0.99
end

POMDPs.discount(p::SubHuntPOMDP) = p.discount
POMDPs.isterminal(p::SubHuntPOMDP, s::SubState) = s == END_KILL || reached_goal(p, s)

function reached_goal(p::SubHuntPOMDP, s::SubState)
    dir = DIR[s.goal]
    progress = dot(dir, s.target)
    if sum(dir) == 1
        return progress == p.size
    else
        return progress == -1 # -1 is for case when dir is negative
    end
end

function POMDPs.reward(p::SubHuntPOMDP, s::SubState, a::Int, sp::SubState)
    if sp == END_KILL
        return 100.0
    else
        return 0.0
    end
end

include("enums.jl")
include("spaces_and_distributions.jl")
include("transition.jl")
include("observation.jl")

include("discrete.jl")

include("edge_cases.jl")
include("heuristics.jl")
include("visualization.jl")

end # module


text:
using ProgressMeter

"    timed_ljh_rewriter(src, dest, timeout_s, pulses_written, pulses_total, i)
Open the ljh file at path `src`. Copy it's header to a file at `dest`, overwiting any existing file with that name.
Then write records from `src` to `dest`, pausing between each write by the difference between successive timestamps.
The maximum pause time (in seconds) for a single record is given by `timeout_s`. After completion the files `src` and `dest`
will be identical. `pulses_written` is an array, this process will increment `pulses_written[i]` each time it writes a pulse.
`pulses_total` is an array, this process will write the number of pulses in src to `pulses_total[i]` "
function timed_ljh_rewriter(src, dest, timeout_s, fastforward, pulses_written, pulses_total, i)
  ljh0 = LJH.LJHFile(src)
  #write hearder to dest file
  f = open(dest,"w+")
  seekstart(ljh0.io)
  write(f, read(ljh0.io,ljh0.datastartpos))
  flush(f)
  seekstart(f)
  ljh1 = LJH.LJHFile(dest, f)
  pulses_total[i] = length(ljh0)

  if length(ljh0)==0
    return
  end
  record = ljh0[1]
  tlast = record.timestamp_usec
  t_start = time()
  tnext = t_start
  for record in ljh0
    tnext += clamp((record.timestamp_usec-tlast)*1e-6/fastforward,0,timeout_s)
    tlast = record.timestamp_usec
    to_sleep_s = max(0.0,tnext-time())
    if to_sleep_s>=0.001 #minimum sleep time is 0.001 s, don't bother sleeping
      sleep(to_sleep_s)
    end
    write(ljh1,record)
    pulses_written[i]+=1
  end
  close(ljh0)
  close(ljh1)
  return
end

function mattersimprogress(srcnames, destnames, tasks, channels, timeout_s, pulses_written, pulses_total)
  println("\nMatter Simulator Running on $(length(channels)) channels, with timeout = $timeout_s s")
  println("Channel numbers: $channels")
  println("First input file: $(first(srcnames))")
  println("First output file: $(first(destnames))")
  println("Process pid = $(getpid())")
  flush(STDOUT)
  # srcs_size = sum(stat(src).size for src in srcnames)
  # dests_size = sum(stat(dest).size for dest in destnames)
  ntotal = sum(pulses_total)
  nwritten = sum(pulses_written)
  p = Progress(ntotal,0.5,"Matter Simulator: ")
  i=0
  while nwritten < ntotal && !all(istaskdone.(tasks))
    ProgressMeter.update!(p,nwritten)
    sleep(0.5)
    i+=1
    if i==10
      i=0
      println("\n$nwritten pulses written in $(p.tlast-p.tfirst) seconds, $(nwritten/(p.tlast-p.tfirst)) cps avg")
    end
    nwritten = sum(pulses_written)
    # dests_size = sum(stat(dest).size for dest in destnames)
  end
  ProgressMeter.update!(p, nwritten) # make sure it reads 100% finished if it  is
  if nwritten == ntotal
    println("\nMatter Simulator finished.")
  else
    println("\nMatter Simulator finished without writing all data.")
  end
  println("Time elapsed: $(p.tlast-p.tfirst) seconds.")
  dests_size = sum(stat(dest).size for dest in destnames)
  println("Total bytes written $dests_size.")
end

"    mattersim(srcdir, destdir, timeout_s=0.01, maxchannels=240)"
function mattersim(srcdir, destdir, timeout_s=0.01, fastforward=1.0, maxchannels=240)
  maxchannels<=0 && error("maxchannels must be positive, was $maxchannels")
  ljhdict = LJH.allchannels(srcdir, maxchannels) # ordered dict mapping channel number to filename
  length(ljhdict)>0 || error("found no ljh files in $srcdir")
  isdir(destdir) || mkdir(destdir)
  destnames = LJH.fnames(destdir, keys(ljhdict))
  pulses_written = zeros(Int, length(ljhdict)) # used to record how many pulses have been written
  pulses_total = zeros(Int, length(ljhdict))
  # i is an index for writinginto pulses_written and pulses_total
  tasks = [@task timed_ljh_rewriter(src, dest, timeout_s, fastforward, pulses_written, pulses_total, i) for (i,(src,dest)) in enumerate(zip(values(ljhdict), destnames))]
  LJH.write_sentinel_file(destnames[1],true)
  println("Matter simulator wrote sentinel file: $(destnames[1]), true")
  progresstask = @async mattersimprogress(values(ljhdict), destnames, tasks, keys(ljhdict), timeout_s, pulses_written, pulses_total)
  schedule.(tasks)
  wait.(tasks)
  wait(progresstask)
  LJH.write_sentinel_file(destnames[1],false)
  println("Matter simulator wrote sentinel file: $(destnames[1]), false")
end


text:
#=

    This is a port of the SPPL model election.pynb

    This model:
    Summary Statistics
    parameters        mean        std   naive_se      mcse         ess      rhat   ess_per_sec
        Symbol     Float64    Float64    Float64   Float64     Float64   Float64       Float64

        param    299.7774    28.9717     0.2897    0.3653   6479.7518    0.9999      866.9724
            p      0.4814     0.0318     0.0003    0.0004   6508.3665    1.0002      870.8010
        votes   1925.2743   131.9134     1.3191    1.7207   6442.0098    1.0003      861.9226
          win      0.2902     0.4539     0.0045    0.0056   6597.9243    0.9999      882.7836

=#

using Turing, StatsPlots, Distributions, StatsBase
using CSV
include("jl_utils.jl")

#=
from sppl.sym_util import binspace
n = 4000
param ~= randint(low=250, high=350)
switch (param) cases (b in range(250, 350)):
    p ~= beta(a=277, b=b)
switch (p) cases (x in binspace(0, 1, 20)):
    votes ~= binom(n=n, p=(x.left + x.right)/2)
win ~= votes > 0.5*n

=#
@model function election(n=4000)
    param ~ DiscreteUniform(250,350)
    p ~ Beta(277,param)
    votes ~ Binomial(n,p)

    win ~ Dirac(votes > 0.5*n)

end

n = 4000
model = election(n)

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
chns = sample(model, PG(5), 10_000)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

# chns = sample(model, HMC(0.1,6), 1_000)
# chns = sample(model, NUTS(), 1_000)

display(chns)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
export safe_fragment

"""
    safe_fragment(object_dtags::Vector{NTuple{2,Int32}},
                  tool_dtags::Vector{NTuple{2,Int32}};
                  material_hierarchy::Vector{Material} = Material[])

The equivalent to performing `gmsh.model.occ.fragment(object_dtags, tool_dtags)`, but
preserving the physical groups of the highest dimensional entities.

Only the highest dimensional physical groups can be preserved due to the available
parent-child relationships from `gmsh.model.occ.fragment`.

In the event that two overlapping entities have material physical groups, the optional
material hierarchy is used to choose a single material for the resultant overlapping
entity/entities.
"""
function safe_fragment(object_dtags::Vector{NTuple{2,Int32}},
                       tool_dtags::Vector{NTuple{2,Int32}};
                       material_hierarchy::Vector{Material} = Material[])

    # Get all physical groups and their names
    dim = max(maximum(getindex.(object_dtags, 1)), maximum(getindex.(tool_dtags, 1)))
    groups = gmsh.model.get_physical_groups(dim)
    old_tags = [gmsh.model.get_entities_for_physical_group(dtag[1], dtag[2]) for dtag in groups]
    names = [gmsh.model.get_physical_name(dtag[1], dtag[2]) for dtag in groups]

    # Fragment
    nents = length(object_dtags) + length(tool_dtags)
    @info "... Fragmenting $nents entities"
    output_dtags, output_dtags_map = gmsh.model.occ.fragment(object_dtags, tool_dtags)

    # Update the physical groups
    # The parent at index i of input_dtags has children out_dtags_map[i]
    @info "... Updating physical groups"
    input_dtags = vcat(object_dtags, tool_dtags)
    new_tags = similar(old_tags)
    # For each physical group
    for (i, group) in enumerate(groups)
        gdim, gnum = group
        new_tags[i] = NTuple{2,Int32}[]
        # For each of the dim tags in the physical group
        for tag in old_tags[i]
            dtag = (gdim, tag)
            # If the dim tag was one of the entities in the fragment
            if dtag ∈ input_dtags
                # Get its children
                index = findfirst(x->x==dtag, input_dtags)
                children = output_dtags_map[index]
                # Add the children to the new physical group
                for child in children
                    if child[2] ∉ new_tags[i]
                        push!(new_tags[i], child[2])
                    end
                end
            else
                # If it wasn't in the fragment, no changes necessary.
                push!(new_tags[i], tag)
            end
        end
    end

    # Remove old groups and synchronize
    for name in names
        gmsh.model.remove_physical_name(name)
    end
    gmsh.model.occ.synchronize()

    # Process the material hierarchy, if it exists, so that each entity has one
    # or fewer materials
    if 0 < length(material_hierarchy)
        _process_material_hierarchy!(names, new_tags, material_hierarchy)
    end

    # Create the new physical groups
    for (i, group) in enumerate(groups)
        gdim, gnum = group
        ptag = gmsh.model.add_physical_group(gdim, new_tags[i])
        gmsh.model.set_physical_name(gdim, ptag, names[i])
    end

    # Apply material colors
    if 0 < length(material_hierarchy)
        color_material_physical_group_entities(material_hierarchy)
    end

    return output_dtags, output_dtags_map
end

function _process_material_hierarchy!(names::Vector{String},
                                      new_tags::Vector{Vector{Int32}},
                                      material_hierarchy::Vector{Material})
    material_indices = findall(x->startswith(x, "Material: "), names)
    material_names = names[material_indices]
    material_dict = Dict{String, BitSet}()
    for (i, name) in enumerate(material_names)
        material_dict[name] = BitSet(new_tags[material_indices[i]])
    end
    # Ensure each material group is present in the hierarchy and warn now if it's not.
    # Otherwise the error that occurs later is not as easy to figure out
    hierarchy_names = ["Material: "*mat.name for mat in material_hierarchy]
    for hierarchy_name in hierarchy_names
        if hierarchy_name ∉ material_names
            error("Physical groups does not contain: ", hierarchy_name)
        end
    end
    # Use the hierarchy to ensure that no entity has more than 1 material by removing
    # all shared elements from sets lower than the current set
    nmats = length(hierarchy_names)
    for i ∈ 1:nmats
        for j ∈ i+1:nmats
            setdiff!(material_dict[hierarchy_names[j]], material_dict[hierarchy_names[i]])
        end
    end
    for (mat_index, global_index) in enumerate(material_indices)
        new_tags[global_index] = collect(material_dict[material_names[mat_index]])
    end
    return material_dict
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
using LinearAlgebra
using Roots
K = [2 -1; -1 1]
M = [4 1;1 2]
a = (2*7e10)
b = (2700/12)
c(w) =  det((a*K)-((w^2)*(b*M)))
find_zeros(c,0,1000000)


text:
using Gadfly, Dates

dates = [Date("2017-08-07"), Date("2017-08-08"), Date("2017-08-09")]
date = [Date("2017-08-08")]
hstack(
    plot(y=dates, Geom.point, yintercept=date, Geom.hline),
    plot(x=dates, Geom.point, xintercept=date, Geom.vline) )


text:
# All functions related to reproduction.

"""
    reproduce!(patch)

Reproduction of individuals in a patch (default function).
"""
function reproduce!(patch::Patch) #TODO: refactor!
    #TODO This is one of the most compute-intensive functions in the model - optimise?
    for ind in patch.community
        ind.marked && continue # individual might not have established yet
        ind.size < ind.traits["repsize"] && continue
        metaboffs = setting("fertility") * ind.size^(-1/4) * exp(-act/(boltz*patch.temp))
        noffs = rand(Poisson(metaboffs))
        noffs < 1 && continue
        partners = findmate(patch.community, ind)
        if length(partners) < 1 && rand() < ind.traits["selfing"]
            partners = [ind]
        elseif length(partners) < 1
            continue
        end
        numpartners = Integer(round(ind.traits["numpollen"]))
        if numpartners == 0
            #XXX happens really often?
            # `numpollen` not handled specifically in `createtraits()`?
            @simlog("Individual cannot reproduce, `numpollen` too low.", 'd')
            continue
        end
        for ptn in 1:numpartners
            partner = rand(partners)
            # subtract offspring mass from parent #XXX is that sensible, does it make a difference?
            parentmass = ind.size - noffs * ind.traits["seedsize"]
            if parentmass <= 0
                break
            else
                ind.size = parentmass
            end
            append!(patch.seedbank, createoffspring(noffs, ind, partner))
        end
    end
    @simlog("Patch $(patch.id): $(length(patch.seedbank)) offspring", 'd')
end

"""
    greproduce!(patch)

Reproduction of individuals in a patch with global mating.
"""
function greproduce!(patch::Patch)
    for ind in patch.community
        ind.marked && continue # individual might not have established yet
        ind.size < ind.traits["repsize"] && continue
        metaboffs = setting("fertility") * ind.size^(-1/4) * exp(-act/(boltz*patch.temp))
        noffs = rand(Poisson(metaboffs))
        noffs < 1 && continue
        partners = findmate([(map(x -> x.community, world)...)...], ind)
        if length(partners) < 1 && rand() < ind.traits["selfing"]
            partners = [ind]
        elseif length(partners) < 1
            continue
        end
        numpartners = Integer(round(ind.traits["numpollen"]))
        for ptn in 1:numpartners
            partner = rand(partners, 1)[1]
            parentmass = ind.size - noffs * ind.traits["seedsize"] # subtract offspring mass from parent
            if parentmass <= 0
                continue
            else
                ind.size = parentmass
            end
            append!(patch.seedbank, createoffspring(noffs, ind, partner))
        end
    end
    @simlog("Patch $(patch.id): $(length(patch.seedbank)) offspring", 'd')
end

"""
    reproduce!(world)

Carry out reproduction on all patches.
"""
function reproduce!(world::Array{Patch,1})
    for patch in world
        if setting("globalmating")
            greproduce!(patch)
        else
            (patch.isisland || !setting("static")) && reproduce!(patch) # pmap(!,patch) ???
        end
    end
end

"""
    findmate(population, individual)

Find a reproduction partner for the given individual in the given population.
"""
function findmate(population::AbstractArray{Individual, 1}, ind::Individual)
    #XXX This function is pretty expensive
    indstate = ind.marked
    ind.marked = true
    mates = Individual[]
    startidx = rand(eachindex(population))
    mateidx = startidx
    while true
        mate = population[mateidx]
        if !mate.marked && iscompatible(mate, ind)
            push!(mates, mate)
            break
        end
        mateidx += 1
        mateidx > length(eachindex(population)) && (mateidx = 1)
        mateidx == startidx && break
    end
    ind.marked = indstate
    mates
end

"""
    createoffspring(noffs, individual, partner, dimorphism)

The main reproduction function. Take two organisms and create the given number
of offspring individuals. Returns an array of individuals.
"""
function createoffspring(noffs::Integer, ind::Individual, partner::Individual, dimorphism::Bool=false)
    #TODO This is a very compute-intensive function - optimise?
    offspring = Individual[]
    heterozygosity = setting("heterozygosity")
    for i in 1:noffs # pmap? this loop could be factorized!
        # offspring have different genomes due to recombination
        partnergenome = meiosis(partner.genome, false, heterozygosity)
        mothergenome = meiosis(ind.genome, true, heterozygosity)
        (isempty(partnergenome) || isempty(mothergenome)) && continue
        genome = vcat(partnergenome,mothergenome)
        if setting("degpleiotropy") == 0 && setting("linkage") == "none"
            traits = gettraitdictfast(genome, setting("traitnames"))
        else
            traits = gettraitdict(genome, setting("traitnames"))
        end
        marked = true
        fitness = 0.0
        newpartner = 0
        newsize = ind.traits["seedsize"]
        newid = rand(UInt64) #XXX Is this collision-safe?
        sex = hermaphrodite
        if dimorphism
            rand(Bool) ? sex = male : sex = female
        end
        lineage = ind.lineage
        if lineage != partner.lineage
            # hybrids are assigned the lineage of the parent they are phenotypically most similar to
            idiff = sum(t -> abs(traits[t]-ind.traits[t]), keys(traits))
            pdiff = sum(t -> abs(traits[t]-partner.traits[t]), keys(traits))
            (pdiff < idiff) && (lineage = partner.lineage)
        end
        push!(offspring, Individual(lineage, genome, traits, marked, fitness,
                                    fitness, newsize, sex, newpartner, newid))
    end
    offspring
end

max_stars_repo_name: CCTB-Ecomods/gemm

max_issues_repo_name: CCTB-Ecomods/gemm

max_forks_repo_name: CCTB-Ecomods/gemm


text:
function CoordinateDescentTM{T<:Signed, F<:AbstractFloat}(Y::Array{F}, N::T; MaxIterλ::T = 100)
    ## Check variable type
    #readtype = typeof(Y[1])
    #if readtype == Float32
    #    N = Int32(N)
    #    MaxIterλ = Int32(MaxIterλ)
    #else
    #    N = Int64(N)
    #    MaxIterλ = Int64(MaxIterλ)
    #end

    ## Build LASSO-path

    ## Number of available lambdas: Granularity of the regulation parameter
    # MaxIterλ = 100

    ## Potential Penalization Parameters Grid
    passo = Float64(1 / MaxIterλ)

    ## Linear grid
    grid_linear = Array{Float64}(collect(1+0.001:passo*10:10))

    ## Logarithmic Grid
    LASSO_Path = Array{Float64}(1 - log10(grid_linear))

    ## Compose Matrices
    X, GMat, arrayInnerY, array_σ2, λ_max = ComposeMatrices(N, Y)

    ## LASSO Path
    println("Initializing LASSO path.")

    ## Initializing parameters
    β_tilde = zeros(N,1)
    β_ols = zeros(N,1)
    convergenceFlag = 0
    ## Active set
    A_actual = Int[]

    ## Preallocating vector β_lasso for speed
    β_lasso = zeros(N,length(LASSO_Path))

    ## First for runs through the LASSO path
    for cont = 1:length(LASSO_Path)-2
        ## cont
        ## Computing next lambda
        λ = LASSO_Path[cont] * λ_max
        lastchange = 0
        it_while = 1
        while (convergenceFlag != 2)
            ##  Reference active set
            A_reference = A_actual

            ##  Cycle
            if it_while == 1
                j_range = collect(1:N)'
            else
                j_range = A_actual
            end

            for j = j_range

                aux, λEff, β_ols = CoordinateDescentCore(A_actual, GMat, j, β_tilde, β_ols, N, arrayInnerY, λ, array_σ2)

                A_actual, β_tilde, lastchange, breakThis = updateCandidates(λEff, β_ols, j, A_actual, lastchange, β_tilde)

                if breakThis == 1
                    break
                end

            end
            ##  Check if the active set has been altered in the last cycle
            convergenceFlag = convCheck(A_reference, A_actual)
            it_while = it_while + 1

        end
        it_while
        convergenceFlag = 0

        β_lasso[:,cont] = β_tilde
    end

    println("LASSO path finished.")

    ## BIC
    bestB = calcBestB(N, Y, X, β_lasso, array_σ2)

    return β_lasso, LASSO_Path, bestB, X, array_σ2
end


text:
"""
    _flatten(ex::Expr)

Remove the `:block` in the inner expression and put it at top level.
"""
function _flatten(ex::Expr)
    head = ex.head
    args = _flatten.(ex.args)
    if head == :block
        newargs = Any[]
        for e in args
            if e isa Expr && e.head == :block
                append!(newargs, e.args)
            else
                push!(newargs, e)
            end
        end
        return Expr(:block, newargs...)
    elseif head == :(=) && args[2] isa Expr && args[2].head == :block
        newargs = args[2].args
        newargs[end] = Expr(:(=), args[1], newargs[end])
        return Expr(:block, newargs...)
    elseif head == :call && args[1]==:scalar && args[2] isa Expr && args[2].head == :block
        newargs = args[2].args
        newargs[end] = Expr(:call, args[1], newargs[end])
        return Expr(:block, newargs...)
    else
        return Expr(head, args...)
    end
end
_flatten(e) = e

"""
    removelinenumbernode(ex::Expr)

Remove all the `LineNumberNode` in the `:block` structure.
"""
function removelinenumbernode(ex::Expr)
    if ex.head == :block
        return Expr(:block, (removelinenumbernode(e) for e in ex.args if !(e isa LineNumberNode))...)
    else
        return ex
    end
end
removelinenumbernode(ex) = ex

const tensoroperationsfunctions = (:similar_from_indices,
                                    :cached_similar_from_indices,
                                    :add!, :trace!, :contract!,
                                    :scalar, :IndexError)

"""
    addtensoroperations(ex::Expr)

Provide the path to the functions of the tensor operations that are used in instantiating
the tensor expressions.
"""
function addtensoroperations(ex::Expr)
    if ex.head == :call && ex.args[1] in tensoroperationsfunctions
        return Expr(ex.head, GlobalRef(TensorContractionsXD, ex.args[1]),
                        (addtensoroperations(ex.args[i]) for i in 2:length(ex.args))...)
    else
        return Expr(ex.head, (addtensoroperations(e) for e in ex.args)...)
    end
end
addtensoroperations(ex) = ex


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

include("rng_init.jl")


text:
module DataAPI

"""
    defaultarray(T, N)

For a given element type `T` and number of dimensions `N`, return the appropriate array
type.

The default definition returns `Array{T, N}`. This function is useful for custom types that
have a more efficient vectorized representation (usually using SOA optimizations).

This generic function is owned by DataAPI.jl itself, which is the sole provider of the
default definition.
"""
function defaultarray end
defaultarray(::Type{T}, N) where {T} = Array{T, N}

"""
    refarray(A::AbstractArray)

For a given array `A`, potentially return an optimized "ref array" representation of the
original array, which can allow for faster comparison and sorting.

The default definition just returns the input array. This function is useful for custom
array types which already store a "hashed"-like representation of elements where testing
equality or permuting elements in place can be much faster than the original scalar value,
like pooled arrays.

This generic function is owned by DataAPI.jl itself, which is the sole provider of the
default definition.
"""
function refarray end
refarray(A::AbstractArray) = A

"""
    refvalue(A, x)

For the *original* array `A`, and a "ref value" `x` taken from `refarray(A)`, return the
appropriate *original* value. `refvalue(A, refarray(A)[I...])` must be equal to `A[I...]`.

By default, `refvalue(A, x)` returns `x` (since `refarray(A)` returns `A` by default).
This allows recovering an original array element after operating on the "ref values".

This generic function is owned by DataAPI.jl itself, which is the sole provider of the
default definition.
"""
function refvalue end
refvalue(A::AbstractArray, x) = x

"""
    refpool(A)

Whenever available, return an indexable object `pool` such that, given the *original* array `A` and
a "ref value" `x` taken from `refarray(A)`, `pool[x]` is the appropriate *original* value. Return
`nothing` if such object is not available.

By default, `refpool(A)` returns `nothing`.

If `refpool(A)` is not `nothing`, then `refpool(A)[refarray(A)[I...]]`
must be equal to (according to `isequal`) and of the same type as `A[I...]`,
and the object returned by `refpool(A)` must implement the iteration and
indexing interfaces as well as the `length`, `eachindex`, `keys`, `values`, `pairs`,
`firstindex`, `lastindex`, and `eltype` functions
in accordance with the `AbstractArray` interface.

This generic function is owned by DataAPI.jl itself, which is the sole provider of the
default definition.
"""
function refpool end
refpool(A::AbstractArray) = nothing

"""
    invrefpool(A)

Whenever available, return an indexable object such that given an array `A`
for which `refpool(A)` is not `nothing`:

* for any valid index `x` into `refpool(A)`, `invrefpool(A)[refpool(A)[x]]` is equal to `x`
  (according to `isequal`) and of the same type as `x`;
* for any valid index `ix` into `invrefpool(A)` , `refpool(A)[invrefpool(A)[ix]]` is equal to `ix`
  (according to `isequal`) and of the same type as `ix`.

Additionally it is required that for `invrefpool(A)` the following methods are defined:

* `Base.haskey`: allowing to check if `ix` is a valid index into it.
* `Base.get`: allowing to get a value from it or a passed default value if it is not present.

By default, `invrefpool(A)` returns `nothing`.

If `invrefpool(A)` is not `nothing`, then `refpool(A)` also must not be `nothing`.

This generic function is owned by DataAPI.jl itself, which is the sole provider of the
default definition.
"""
function invrefpool end
invrefpool(A::AbstractArray) = nothing

"""
    describe(io::IO, x)

For an object `x`, print descriptive statistics to `io`.

This generic function is owned by StatsBase.jl, which is the sole provider of the default
definition.
"""
function describe end

"""
    levels(x)

Return a vector of unique values which occur or could occur in collection `x`,
omitting `missing` even if present. Values are returned in the preferred order
for the collection, with the result of [`sort`](@ref) as a default.
If the collection is not sortable then the order of levels is unspecified.

Contrary to [`unique`](@ref), this function may return values which do not
actually occur in the data, and does not preserve their order of appearance in `x`.
"""
function levels(x)
    T = Base.nonmissingtype(eltype(x))
    u = unique(x)
    # unique returns its input with copying for ranges
    # (and possibly for other types guaranteed to hold unique values)
    nmu = (u === x || Base.mightalias(u, x)) ? filter(!ismissing, u) : filter!(!ismissing, u)
    levs = convert(AbstractArray{T}, nmu)
    try
        sort!(levs)
    catch
    end
    levs
end

"""
    Between(first, last)

Select the columns between `first` and `last` from a table.
"""
struct Between{T1 <: Union{Int, Symbol}, T2 <: Union{Int, Symbol}}
    first::T1
    last::T2
end

Between(x::AbstractString, y::AbstractString) = Between(Symbol(x), Symbol(y))
Between(x::Union{Int, Symbol}, y::AbstractString) = Between(x, Symbol(y))
Between(x::AbstractString, y::Union{Int, Symbol}) = Between(Symbol(x), y)

"""
    All(cols...)

Select the union of the selections in `cols`. If `cols == ()`, select all columns.
"""
struct All{T<:Tuple}
    cols::T
    function All(args...)
        if !isempty(args)
            Base.depwarn("All(args...) is deprecated, use Cols(args...) instead", :All)
        end
        return new{typeof(args)}(args)
    end
end

"""
    Cols(cols...)

Select the union of the selections in `cols`. If `cols == ()`, select no columns.
"""
struct Cols{T<:Tuple}
    cols::T
    Cols(args...) = new{typeof(args)}(args)
end

"""
    unwrap(x)

For a given scalar argument `x`, potentially "unwrap" it to return the base wrapped value.
Useful as a generic API for wrapper types when the original value is needed.

The default definition just returns `x` itself, i.e. no unwrapping is performned.

This generic function is owned by DataAPI.jl itself, which is the sole provider of the
default definition.
"""
function unwrap end
unwrap(x) = x

end # module


text:
Threads.nthreads() == 20 ||
    error("Doc build on Noctua should be run with 20 Julia threads!")
println("--- :julia: Instantiating project")
using Pkg
Pkg.activate("..")
Pkg.instantiate()
Pkg.activate(".")
Pkg.instantiate()
push!(LOAD_PATH, joinpath(@__DIR__, ".."))
deleteat!(LOAD_PATH, 2)
println("+++ :julia: Building documentation")
include("make.jl")


text:
module Heading

    function main(main_head_text)
        println("\n\n==================================================================")
        println("+", " ", main_head_text)
        println("==================================================================\n")
    end

    function sub(txt)
        println(txt, "\n")
    end
end


text:
using SimpleHttpIO
using StreamReader

if !isdefined(:TEST_DIR)
    include("utils.jl")
end

facts("Test Request.jl") do
    context("parse headers") do

        context("parse headers from IO") do
            io = IOBuffer()
            sio = IOSocket(io)
            wclf(io)
            seekstart(io)

            frqe(;message="Empty first line.") do
                S.readfirstline(sio)
            end

            io = IOBuffer()
            sio = IOSocket(io)

            line = "first line"
            ws = wclf(io, line)
            seekstart(io)
            @fact S.readfirstline(sio) --> (ws, line)

            io = IOBuffer()
            sio = IOSocket(io)

            ws = wclf(io, "H1: H1 value")
            ws += wclf(io, "HM: HM value 1")
            ws += wclf(io, "HM: HM value 2")
            ws += wclf(io)

            seekstart(io)
            hs = sizeof(io.data)

            eh = ["H1"=>"H1 value","HM"=>String["HM value 1","HM value 2"]]

            s, lnum, h = S.readheaders(sio)
            @fact ws --> hs
            @fact s --> hs
            @fact lnum --> 4

            for (k, v) in h
                @fact h[k] --> eh[k] ("Mismatch " * k * " Header")
            end
        end

        context("parse info") do
            for line in String["GET /index.html", "GET /index.html HTTP/1.1",
                "GET /directory%20name/file.txt HTTP/1.1"]
                parts = split(line, " ")

                if length(parts) == 2
                    m, r = parts
                    pr = default_protocol = astuple(S.FAKE_PROTOCOL)
                else
                    m, r, pr = parts
                    pr = split(pr, "/")
                    default_protocol = N
                end

                method, resource, protocol = S.parseinfo(line, default_protocol)
                @fact method --> m
                @fact r --> resource
                @fact protocol --> pr
            end
        end

        context("read request info") do
            for line in String["GET /index.html", "GET /index.html HTTP/1.1",
                "GET /directory%20name/file.txt HTTP/1.1"]
                context("read request info at " * repr(line)) do
                    io = IOBuffer()
                    wclf(io, line)
                    seekstart(io)

                    parts = split(line, " ")

                    if length(parts) == 2
                        m, rs = parts
                        pr = default_protocol = astuple(S.FAKE_PROTOCOL)
                        protocol!(S.FAKE_PROTOCOL)
                    else
                        m, rs, pr = parts
                        pr = split(pr, "/")
                        default_protocol = N
                    end

                    r = S.Request(IOSocket(io))
                    S.readinfo(r, S.PROTOCOLS, default_protocol)

                    @fact r.method --> m
                    @fact r.resource --> rs

                    if pr != N
                        @fact r.protocol.name --> pr[1]
                        @fact r.protocol.version --> pr[2]
                    else
                        @fact r.protocol --> pr
                    end
                end
            end
        end

        context("post request") do
            types = ["application/x-www-form-urlencoded" => ("postvar=pvalue&postvar2=pvalue2",
                    ["postvar2"=>{"pvalue2"},"postvar"=>{"pvalue"}]),
                "application/json" => ("{\"postvar\":\"pvalue\",\"postvar2\":\"pvalue2\"}",
                    ["postvar2"=>"pvalue2","postvar"=>"pvalue"])]

            for (tk, (tdata, tdatae)) in types
                context("test $tk content-type") do
                    req_io = IOBuffer()

                    write(req_io, "POST /page?getvar=gvalue HTTP/1.1\r\n" *
                        "Host: 0.0.0.0:7000\r\n" *
                        "Content-Type: $tk\r\n" *
                        "Content-Length: $(sizeof(tdata))\r\n" *
                        "\r\n" *
                        tdata)

                    seekstart(req_io)

                    req = Request(IOSocket(req_io))
                    init(req, S.PROTOCOLS)
                    data_parsed = parsedata(req)

                    @fact data_parsed --> true
                    @fact req.protocol --> HTTP_1_1
                    @fact req.content_length --> sizeof(tdata)
                    @fact req.content_type --> tk
                    @fact req.get --> ["getvar"=>{"gvalue"}]
                    @fact req.post --> tdatae
                end
            end
        end
    end

    context("parse multipart/form-data") do
        sample_path = joinpath(TEST_DIR, "data", "multipart", "rawdata", "sample.txt")
        sio = IOSocket(open(sample_path, "r"))
        s, info = S.readinfo(sio, S.PROTOCOLS)
        size, linenum, headers = S.readheaders(sio)
        content_size = parseint(headers["Content-Length"])
        boundary = S.multipartboundary(headers["Content-Type"])
        s, d, files = S.parsempfd(sio, boundary)

        for f in files
            fpath = joinpath(TEST_DIR, "data", "multipart", "files", f.name)
            open(fpath, "r") do fio
                @fact f.io.data --> readbytes(fio, filesize(fpath)) "File Data mismatch '$(f.name)' -> '$(fpath)'"
            end
        end
    end

end


text:
using RadioAstroTools
using Test

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