text:
#####
##### utilities for testing sample correctness
#####

"""
$(SIGNATURES)

Run `K` chains of MCMC on `ℓ`, each for `N` samples, return a vector of position matrices
and EBFMI statistics as fields of a `NamedTuple`.

Keyword arguments are passed to `mcmc_with_warmup`.
"""
function run_chains(rng, ℓ, N, K; mcmc_args...)
    results = [mcmc_with_warmup(rng, ℓ, N; reporter = NoProgressReport(), mcmc_args...)
               for i in 1:K]
    (position_matrices = map(r -> position_matrix(r.chain), results),
     EBFMIs = map(r -> EBFMI(r.tree_statistics), results))
end

"""
$(SIGNATURES)

`R̂` (within/between variance) and `τ` (effective sample size coefficient) statistics for
position matrices, eg the output of `run_chains`.
"""
function mcmc_statistics(position_matrices)
    K = size(first(position_matrices), 1)
    R̂ = [potential_scale_reduction([mx[i, :] for mx in position_matrices]...) for i in 1:K]
    τ = mean([vec(mapslices(first ∘ ess_factor_estimate, mx; dims = 2))
              for mx in position_matrices])
    (R̂ = R̂, τ = τ)
end

"""
$(SIGNATURES)

Jitter in `(-ϵ, +ϵ)`. Useful for tiebreaking for the Kolmogorov-Smirnov tests.
"""
jitter(rng, len, ϵ = 64*eps()) = (2 * ϵ) .* (rand(rng, len) .-  0.5)

"""
$(SIGNATURES)

Run MCMC on `ℓ`, obtaining `N` samples from `K` independently adapted chains.

`R̂`, `τ`, Kolmogorov-Smirnov and Anderson-Darling `p`, and EBFMIs are obtained and compared
to thresholds that either *alert* or *fail*. The latter should be lax because of false
positives, the tests can be rather hair-trigger.

Output is sent to `io`. Specifically, `title` is printed for the first alert.

`mcmc_args` are passed down to `mcmc_with_warmup`.
"""
function NUTS_tests(rng, ℓ, title, N; K = 3, io = stdout, mcmc_args = NamedTuple(),
                    R̂_alert = 1.05, R̂_fail = 2 * (R̂_alert - 1) + 1,
                    τ_alert = 0.2, τ_fail = τ_alert * 0.5,
                    p_alert = 0.001, p_fail = p_alert * 0.1,
                    EBFMI_alert = 0.5, EBFMI_fail = EBFMI_alert / 2)
    @argcheck 1 < R̂_alert ≤ R̂_fail
    @argcheck 0 < τ_fail ≤ τ_alert
    @argcheck 0 < p_fail ≤ p_alert
    @argcheck 0 < EBFMI_fail < EBFMI_alert

    d = dimension(ℓ)
    _round(x) = round(x; sigdigits = 3) # for printing

    title_printed = false
    function _print_title_once()
        if !title_printed
            printstyled(io, "INFO while testing: $(title), dimension $(d)\n";
                        color = :blue, bold = true)
            title_printed = true
        end
    end
    @unpack position_matrices, EBFMIs = run_chains(RNG, ℓ, N, K; mcmc_args...)

    # mixing and autocorrelation diagnostics
    @unpack R̂, τ = mcmc_statistics(position_matrices)
    max_R̂ = maximum(R̂)
    if max_R̂ > R̂_alert
        _print_title_once()
        println(io, "ALERT some R̂ = $(_round.(R̂)) > $(_round(R̂_alert))" )
    end
    @test all(max_R̂ ≤ R̂_fail)

    min_τ = minimum(τ)
    if min_τ < τ_alert
        _print_title_once()
        println(io, "ALERT some τ = $(_round.(τ)) < $(_round(τ_alert))" )
    end
    @test all(min_τ ≥ τ_fail)

    min_EBFMI = minimum(EBFMIs)
    if min_EBFMI < EBFMI_alert
        _print_title_once()
        println(io, "ALERT some EBFMI = $(_round.(EBFMIs)) > $(_round(EBFMI_alert))" )
    end
    @test all(min_EBFMI ≥ EBFMI_fail)

    # distribution comparison tests
    Z = reduce(hcat, position_matrices)
    Z′ = samples(ℓ, 1000)
    pd_alert = p_alert / d      # a simple Bonferroni correction
    pd_fail = p_fail / d
    ϵ = jitter(rng, size(Z, 2))
    ϵ′ = jitter(rng, size(Z′, 2))
    for i in 1:d
        p_AD = pvalue(KSampleADTest(Z′[i, :], Z[i, :]))
        if p_AD ≤ pd_alert
            _print_title_once()
            println(io, "ALERT extreme Anderson-Darling ",
                    "p ≈ $(round(p_AD, sigdigits = 3)) < $(round(pd_alert))",
                    " for coordinate $(i) of $(d)")
        end
        @test p_AD ≥ pd_fail
    end
end


text:
import Zygote: Params, pullback

function inscribe_history(history_filepath, NN, loss)
    jldopen(history_filepath, "a") do file
        if !haskey(file, "loss")
            file["loss/1"] = loss
        else
            epochs = length(keys(file["loss"])) + 1
            file["loss/$epochs"] = loss
        end

        if !haskey(file, "neural_network")
            file["neural_network/1"] = NN
        else
            epochs = length(keys(file["neural_network"])) + 1
            file["neural_network/$epochs"] = NN
        end
    end
    return nothing
end

inscribe_history(::Nothing, args...) = nothing

function dense_spatial_causality_train!(loss, ps, data, opt; cb = () -> ())
    #= Breaks in Zygote due to mutation...  =#
    local training_loss
    ps = Params(ps)
    for d in data
        gs = gradient(ps) do
            training_loss = loss(d...)
            return training_loss
        end
        for i in 1:2:length(ps)
            nrows, ncols = size(ps[i])
            mask = [x >= y ? 1.0 : 0.0 for x in 1:nrows, y in 1:ncols]
            ps[i] .*= mask
        end
        Flux.update!(opt, ps, gs)
        cb()
    end

end

function train_neural_differential_equation!(NN, NDEType, algorithm, datasets, T_scaling, wT_scaling, iterations, opt, epochs; history_filepath=nothing, causal_penalty=nothing)

    ids = [id for id in keys(datasets)] |> sort

    nde_params = Dict(id => FreeConvectionNDEParameters(datasets[id], T_scaling, wT_scaling) for id in ids)
    T₀ = Dict(id => T_scaling.(interior(datasets[id]["T"])[1, 1, :, 1]) for id in ids)
    ndes = Dict(id => NDEType(NN, datasets[id]; iterations) for id in ids)

    true_sols = [T_scaling.(interior(datasets[id]["T"])[1, 1, :, iterations]) for id in ids]
    true_sols = cat(true_sols..., dims=2)

    function nde_loss()
        nde_sols = cat([solve_nde(ndes[id], NN, T₀[id], algorithm, nde_params[id]) |> Array for id in ids]..., dims=2)
        if !isnothing(causal_penalty)
            return Flux.mse(nde_sols, true_sols) + causal_penalty(NN)
        else
            return Flux.mse(nde_sols, true_sols)
        end
     end

    function nde_callback()
        mse_loss = nde_loss()
        @info @sprintf("Training free convection NDE... MSE loss: μ_loss::%s = %.12e", typeof(mse_loss), mse_loss)
        inscribe_history(history_filepath, NN, mse_loss)
        return nothing
    end

    Flux.train!(nde_loss, Flux.params(NN), Iterators.repeated((), epochs), opt, cb=nde_callback)

    return nothing
end

ext: jl


text:
algo = "GS"
endogenous_nbr = 6
exogenous_nbr = 3
exogenous_deterministic_nbr = 0
forward_indices = collect(1:2)
current_indices = collect(1:6)
backward_indices = collect(4:6)
both_indices = Vector{Int64}(undef, 0)
static_indices = [3]

lre_ws = LinearRationalExpectationsWs(algo,
                                  endogenous_nbr,
                                  exogenous_nbr,
                                  exogenous_deterministic_nbr,
                                  forward_indices,
                                  current_indices,
                                  backward_indices,
                                  both_indices,
                                  static_indices)

state_nbr = length(backward_indices)

ws = LinearRationalExpectations.VarianceWs(endogenous_nbr,
                                              state_nbr,
                                              exogenous_nbr,
                                              lre_ws)
nonstate_nbr = endogenous_nbr - state_nbr

function test_variance_matrices()
    A = zeros(endogenous_nbr, endogenous_nbr)
    nonstationary = true
    while nonstationary
        A = [zeros(endogenous_nbr, endogenous_nbr - state_nbr) randn(endogenous_nbr, state_nbr)]
        nonstationary = any(abs.(eigen(A[backward_indices, backward_indices]).values) .> 1.0)
    end
    A1 = A[backward_indices, backward_indices]
    A2 = A[lre_ws.non_backward_indices, lre_ws.backward_indices]
    B = rand(endogenous_nbr, exogenous_nbr)
    B1 = B[backward_indices, :]
    B2 = B[lre_ws.non_backward_indices, :]
    Σe = randn(exogenous_nbr, exogenous_nbr)
    Σe = transpose(Σe)*Σe
    Σy = zeros(endogenous_nbr, endogenous_nbr)
    LinearRationalExpectations.compute_variance!(Σy, A1, A2, B1, B2, Σe, ws)
    @test Σy[4:6, 4:6]  ≈ A[4:6, 4:6]*Σy[4:6, 4:6]*transpose(A[4:6, 4:6]) + B[4:6, :]*Σe*transpose(B[4:6, :])
    @test Σy[1:3, 1:3]  ≈ A[1:3, 4:6]*Σy[4:6, 4:6]*transpose(A[1:3, 4:6]) + B[1:3, :]*Σe*transpose(B[1:3, :])
    @test Σy[4:6, 1:3]  ≈ A[4:6, 4:6]*Σy[4:6, 4:6]*transpose(A[1:3, 4:6]) + B[4:6, :]*Σe*transpose(B[1:3, :])
    @test Σy ≈ A*Σy*transpose(A) + B*Σe*transpose(B)
end

@testset "Variance matrices" begin
    for i = 1:100
        test_variance_matrices()
    end
end

endogenous_nbr = 8
exogenous_nbr = 3

exogenous_deterministic_nbr = 0
forward_indices = collect(1:2)
current_indices = collect(1:8)
backward_indices = collect(4:8)
both_indices = Vector{Int64}(undef, 0)
static_indices = [3]

lre_ws = LinearRationalExpectationsWs(algo,
                                  endogenous_nbr,
                                  exogenous_nbr,
                                  exogenous_deterministic_nbr,
                                  forward_indices,
                                  current_indices,
                                  backward_indices,
                                  both_indices,
                                  static_indices)

state_nbr = length(backward_indices)

ws = LinearRationalExpectations.VarianceWs(endogenous_nbr,
                                              state_nbr,
                                              exogenous_nbr,
                                              lre_ws)
nonstate_nbr = endogenous_nbr - state_nbr

function test_variance_lre()
    A = zeros(endogenous_nbr, endogenous_nbr)
    nonstationary = true
    while nonstationary
        A1a = [1 0.5; 0 1]
        A1b = randn(2, 3)
        A2a = randn(3, 3)
        A2b = randn(3, 3)
        A = hcat(zeros(endogenous_nbr, endogenous_nbr - state_nbr),
                        vcat(hcat(zeros(3, 2),
                                  A2a),
                             hcat(A1a, A1b),
                             hcat(zeros(3, 2),
                                  A2b)))
        nonstationary = any(abs.(eigen(A[backward_indices, backward_indices]).values) .> 1.0)
    end

    A1 = A[backward_indices, backward_indices]
    A2 = A[lre_ws.non_backward_indices, backward_indices]
    B = rand(endogenous_nbr, exogenous_nbr)

    B1 = B[backward_indices, :]
    B2 = B[lre_ws.non_backward_indices, :]

    Σe = randn(exogenous_nbr, exogenous_nbr)
    Σe = transpose(Σe)*Σe
    Σyss = zeros(state_nbr, state_nbr)
    wsl = LinearRationalExpectations.LyapdWs(state_nbr)

    LinearRationalExpectations.extended_lyapd!(Σyss, A1, B1*Σe*transpose(B1), wsl)
    Σy = zeros(endogenous_nbr, endogenous_nbr)
    LinearRationalExpectations.compute_variance!(Σy, A1, A2, B1, B2, Σe, ws)
    sv = ws.stationary_variables
    sv = [1, 2, 3, 6, 7, 8]
    @test Σy[sv, sv] ≈ A[sv, sv]*Σy[sv, sv]*transpose(A[sv, sv]) + B[sv,:]*Σe*transpose(B[sv, :])
    lre_results = LinearRationalExpectationsResults(endogenous_nbr, exogenous_nbr, state_nbr)
    lre_results.gs1 .= A1
    lre_results.gns1 .= A2
    lre_results.hs1 .= B1
    lre_results.hns1 .= B2
    newΣy = similar(Σy)
    LinearRationalExpectations.compute_variance!(newΣy, lre_results, Σe, ws)
    @test newΣy[sv, sv] ≈ Σy[sv, sv]

    ws2 = LinearRationalExpectations.VarianceWs(endogenous_nbr,
                                               endogenous_nbr,
                                               exogenous_nbr,
                                               lre_ws)
end

@testset "Variance LRE" begin
    for i = 1:100
        test_variance_lre()
    end
end

@testset "Moments" begin

    n = 10
    Sigma_y = rand(n,n)
    Sigma_y = transpose(Sigma_y)*Sigma_y
    sd = sqrt.(diag(Sigma_y))
    C = similar(Sigma_y)
    LinearRationalExpectations.correlation!(C, Sigma_y, sd)
    @test C .* sd .* transpose(sd) ≈ Sigma_y

    C = LinearRationalExpectations.correlation(Sigma_y)
    @test C .* sd .* transpose(sd) ≈ Sigma_y

    order = 4
    backward_indices = [1, 3]
    ns = length(backward_indices)
    A0 = rand(n,ns)
    lre_results = LinearRationalExpectationsResults(n, 2, ns)
    lre_results.gs1 .= A0[[1, 3], :]
    lre_results.gns1 .= A0[[2, (4:n)...], :]
    lre_results.endogenous_variance .= Sigma_y
    stationary_variables = [true, true, true, true, false, false, false, false, false, false]
    nstat = count(stationary_variables)
    AA = [zeros(nstat,nstat) for i = 1:order]
    S1a = view(lre_results.endogenous_variance, backward_indices, :)
    LinearRationalExpectations.autocovariance!(AA, lre_results, zeros(ns, n), zeros(ns, n), zeros(n - ns, n), backward_indices, stationary_variables)
    A = A0[1:4, :]
    A1 = A[[1, 3], :]
    Sigma_ys = Sigma_y[[1, 3], 1:4]
    @test AA[1]  ≈ A*Sigma_ys
    @test AA[2]  ≈ A*A1*Sigma_ys
    @test AA[3]  ≈ A*A1*A1*Sigma_ys
    @test AA[4]  ≈ A*A1*A1*A1*Sigma_ys

    VV = [zeros(nstat) for i = 1:order]
    LinearRationalExpectations.autocovariance!(VV, lre_results, zeros(ns, n), zeros(ns, n), zeros(n - ns, n), backward_indices, stationary_variables)
    @test VV[1]  ≈ diag(AA[1])
    @test VV[2]  ≈ diag(AA[2])
    @test VV[3]  ≈ diag(AA[3])
    @test VV[4]  ≈ diag(AA[4])

end


text:
#=
surface_test:
- Julia version: 1.1.0
- Author: Jirik
- Date: 2019-03-19
=#
# using Revise
using Test
using Logging
using SparseArrays
using Distributed
if nprocs() == 1
    addprocs(3)
end

@everywhere using LarSurf
# using Plasm
# using LinearAlgebraicRepresentation
# Lar = LinearAlgebraicRepresentation

@testset "Get Fchar test" begin

    segmentation = LarSurf.data234()
    block_size = [2,2,2]
    ch = RemoteChannel(()->Channel{Int}(32));
    Flin2 = LarSurf.__grid_get_surface_Fchar_per_block_parallel_channel(segmentation, block_size, ch)
    # @test nnz(Flin2) == 6*4
    # println("get fchar test 5")

    # Flin1 = LarSurf.__grid_get_surface_Fchar(segmentation, block_size)
    # println("get fchar test 2")
    # Flin1 = LarSurf.__grid_get_surface_Fchar_per_block(segmentation, block_size)
    # @test nnz(Flin1) == 6*4
    #
    # println("get fchar test 3")
    # Flin2 = LarSurf.__grid_get_surface_Fchar_per_block_parallel_pmap(segmentation, block_size)
    # @test nnz(Flin2) == 6*4
    # println("get fchar test 4")
end


text:
#!/usr/bin/julia --color=yes

# THIS IS OLD USE COMBINED INSTEAD

using ArgParse

using JLD
using Multilane
using POMDPs
using MCTS
using DataFrames
using DataFramesMeta
using PGFPlots

function curve_plot(mp::AbstractDataFrame, test_key::AbstractString; label=test_key, style="")
    data = sort!(@where(mp, :test_key.==test_key), cols=[:lambda])
    Plots.Linear(data[:time_to_lane], data[:nb_brakes], legendentry=label, style=style)
end

function print_tex(p::TikzPictures.TikzPicture, name)
    println("""

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        % $name
        """)
    println(p.data)
    println("""

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        """)
end

function print_tex(p::Axis, name)
    print_tex(plot(p), name)
end

function show_svg(p::Axis)
    f = string(tempname(), ".svg")
    PGFPlots.save(f, plot(p))
    @spawn run(`display -alpha off $f`)
end

function ttl_at_brake(mp, test_key, nb_brakes)
    data = sort!(@where(mp, :test_key.==test_key), cols=[:lambda])
    i = findfirst(data[:nb_brakes] .<= nb_brakes)
    ttl = data[:time_to_lane]
    nb = data[:nb_brakes]
    return ttl[i-1] + (ttl[i]-ttl[i-1])*(nb_brakes-nb[i-1])/(nb[i]-nb[i-1])
end

function brake_at_ttl(mp, test_key, t)
    data = sort!(@where(mp, :test_key.==test_key), cols=[:lambda])
    ttl = data[:time_to_lane]
    i = findfirst(ttl .>= t)
    nb = data[:nb_brakes]
    return nb[i-1] + (nb[i]-nb[i-1])*(t-ttl[i-1])/(ttl[i]-ttl[i-1])
end

filenames = ["combined_results_Sep_2_11_31.jld", "combined_results_Sep_2_11_52.jld", "combined_results_Sep_8_21_57.jld"]
filenames = [joinpath(Pkg.dir("Multilane"), "data", f) for f in filenames]

if !isdefined(:results)
    results = load(first(filenames))

    for f in filenames[2:end]
        new_results = load(f)
        try
            results = merge_results!(results, new_results)
        catch ex
            warn("Careful merge of data failed, there may be inconsistent data")
            results = merge_results!(results, new_results, careful=false)
        end
    end
end

stats = results["stats"]

mp = by(stats, [:solver_key, :lambda, :problem_key, :solver_problem_key, :test_key]) do df
    DataFrame(steps_in_lane=mean(df[:steps_in_lane]),
              nb_brakes=mean(df[:nb_brakes]),
              time_to_lane=mean(df[:time_to_lane]),
              steps=mean(df[:steps]),
              )
end

@show unique(mp[:test_key])

# first plot = uniform
# a = Axis([
#     curve_plot(mp, "upper_bound_unif", label="Omniscience bound"),
#     curve_plot(mp, "assume_normal_unif", label="Assume all normal"),
# ], xmin=0,
# xmax=30,
# ymin=0,
# ymax=1.5,
# xlabel="Time to change lanes",
# ylabel="Braking actions per episode",
# style="grid=both", title="Behaviors Uniformly Distributed")
#
# # show_svg(a)
#
# print_tex(a, "UNIFORM")

# uniform with pomdp solutions
a = Axis([
    curve_plot(mp, "upper_bound_unif", label="Omniscient", style="dashed"),
    # Plots.Command("\\only<2->{"),
    curve_plot(mp, "assume_normal_unif", label="Single assumed behavior"),
    # Plots.Command("}"),
    curve_plot(mp, "mlmpc_unif", label="Most likely behavior"),
    # Plots.Command("\\only<3->{"),
    curve_plot(mp, "pomcp_unif", label="POMCP"),
    # Plots.Command("}")
], xmin=0,
xmax=30,
ymin=0,
ymax=1.5,
xlabel="Time to target lane",
ylabel="Braking actions per episode",
style="grid=both, legend cell align=left",
# title="Scenario 1: \$\\rho = 0\$"
)

@show ub = ttl_at_brake(mp, "upper_bound_unif", 0.5)
@show an = ttl_at_brake(mp, "assume_normal_unif", 0.5)
@show an - ub
@show (an - ub)/ub

@show bub = brake_at_ttl(mp, "upper_bound_unif", 10.0)
@show ban = brake_at_ttl(mp, "assume_normal_unif", 10.0)
@show ban - bub
@show (ban - bub)/bub

# print_tex(a, "UNIFORM SOLVERS")

PGFPlots.save("/home/zach/Devel/Behavior-aware_lane_changing/acc_17_paper/fig/uniform.tex", a, include_preamble=false)
# show_svg(a)

# correlate with pomdp solutions
a = Axis([
    curve_plot(mp, "upper_bound", label="Omniscient", style="dashed"),
    # Plots.Command("\\only<2->{"),
    curve_plot(mp, "assume_normal", label="Single assumed behavior"),
    # Plots.Command("}"),
    curve_plot(mp, "mlmpc", label="Most likely behavior"),
    # Plots.Command("\\only<3->{"),
    curve_plot(mp, "pomcp", label="POMCP"),
    # Plots.Command("}")
], xmin=0,
xmax=30,
ymin=0,
ymax=1.5,
xlabel="Time to target lane",
ylabel="Braking actions per episode",
style="grid=both, legend cell align=left",
# title="Scenario 2: \$\\rho = 1\$"
)

# calculate the ttl for 0.5 brakes
# between 4 and 5

# print_tex(a, "CORRELATED SOLVERS")

PGFPlots.save("/home/zach/Devel/Behavior-aware_lane_changing/acc_17_paper/fig/correlated.tex", a, include_preamble=false)
# show_svg(a)

# cors = [0.5, 0.75, 0.9]
cors = [0.75]
for cor in cors
    a = Axis([
        curve_plot(mp, @sprintf("upper_bound_%03d", 100*cor), label="Omniscient", style="dashed"),
        # Plots.Command("\\only<2->{"),
        curve_plot(mp, @sprintf("assume_normal_%03d", 100*cor), label="Single assumed behavior"),
        # Plots.Command("}"),
        curve_plot(mp, @sprintf("mlmpc_%03d", 100*cor), label="Most likely behavior"),
        # Plots.Command("\\only<3->{"),
        curve_plot(mp, @sprintf("pomcp_%03d", 100*cor), label="POMCP"),
        # Plots.Command("}")
    ], xmin=0,
    xmax=30,
    ymin=0,
    ymax=1.5,
    xlabel="Time to target lane",
    ylabel="Braking actions per episode",
    style="grid=both, legend cell align=left",
    # title="Scenario 3: \$\\rho = $cor\$"
    )

    # print_tex(a, "CORRELATED $cor")

    PGFPlots.save(@sprintf("/home/zach/Devel/Behavior-aware_lane_changing/acc_17_paper/fig/correlated_%03d.tex", 100*cor), a, include_preamble=false)

    # show_svg(a)
end


text:
@testset "construct by periods" begin
    @test NanoDate(yr,mn,dy,hr,mi,sc,ms,cs,ns) == nd
    @test NanoDate(Yr,Mn,Dy,Hr,Mi,Sc,Ms,Cs,Ns) == nd
    @test NanoDate(yr,mn,dy,hr,mi,sc,ms) == nd0
    @test NanoDate(Yr,Mn,Dy,Hr,Mi,Sc,Ms) == nd0
    @test NanoDate(yr,mn,dy) == nd00
    @test NanoDate(Yr,Mn,Dy) == nd00
end

@testset "construct by Dates types" begin
    @test NanoDate(dayt, tyme) == nd
    @test NanoDate(daytime, ananosecs) == nd
    @test NanoDate(daytime) == nd0
    @test NanoDate(dayt) == nd00
end


text:
## todo: parse and include
## https://www.unicode.org/Public/UNIDATA/Scripts.txt
## unicode_scripts

unicode_blocks = [
    ("Basic_Latin", 0x0000:0x007F),
    ("Latin-1_Supplement", 0x0080:0x00FF),
    ("Latin_Extended-A", 0x0100:0x017F),
    ("Latin_Extended-B", 0x0180:0x024F),
    ("IPA_Extensions", 0x0250:0x02AF),
    ("Spacing_Modifier_Letters", 0x02B0:0x02FF),
    ("Combining_Diacritical_Marks", 0x0300:0x036F),
    ("Greek_and_Coptic", 0x0370:0x03FF),
    ("Cyrillic", 0x0400:0x04FF),
    ("Cyrillic_Supplementary", 0x0500:0x052F),
    ("Armenian", 0x0530:0x058F),
    ("Hebrew", 0x0590:0x05FF),
    ("Arabic", 0x0600:0x06FF),
    ("Syriac", 0x0700:0x074F),
    ("Thaana", 0x0780:0x07BF),
    ("Devanagari", 0x0900:0x097F),
    ("Bengali", 0x0980:0x09FF),
    ("Gurmukhi", 0x0A00:0x0A7F),
    ("Gujarati", 0x0A80:0x0AFF),
    ("Oriya", 0x0B00:0x0B7F),
    ("Tamil", 0x0B80:0x0BFF),
    ("Telugu", 0x0C00:0x0C7F),
    ("Kannada", 0x0C80:0x0CFF),
    ("Malayalam", 0x0D00:0x0D7F),
    ("Sinhala", 0x0D80:0x0DFF),
    ("Thai", 0x0E00:0x0E7F),
    ("Lao", 0x0E80:0x0EFF),
    ("Tibetan", 0x0F00:0x0FFF),
    ("Myanmar", 0x1000:0x109F),
    ("Georgian", 0x10A0:0x10FF),
    ("Hangul_Jamo", 0x1100:0x11FF),
    ("Ethiopic", 0x1200:0x137F),
    ("Cherokee", 0x13A0:0x13FF),
    ("Unified_Canadian_Aboriginal_Syllabics", 0x1400:0x167F),
    ("Ogham", 0x1680:0x169F),
    ("Runic", 0x16A0:0x16FF),
    ("Tagalog", 0x1700:0x171F),
    ("Hanunoo", 0x1720:0x173F),
    ("Buhid", 0x1740:0x175F),
    ("Tagbanwa", 0x1760:0x177F),
    ("Khmer", 0x1780:0x17FF),
    ("Mongolian", 0x1800:0x18AF),
    ("Limbu", 0x1900:0x194F),
    ("Tai_Le", 0x1950:0x197F),
    ("Khmer_Symbols", 0x19E0:0x19FF),
    ("Phonetic_Extensions", 0x1D00:0x1D7F),
    ("Latin_Extended_Additional", 0x1E00:0x1EFF),
    ("Greek_Extended", 0x1F00:0x1FFF),
    ("General_Punctuation", 0x2000:0x206F),
    ("Superscripts_and_Subscripts", 0x2070:0x209F),
    ("Currency_Symbols", 0x20A0:0x20CF),
    ("Combining_Diacritical_Marks_for_Symbols", 0x20D0:0x20FF),
    ("Letterlike_Symbols", 0x2100:0x214F),
    ("Number_Forms", 0x2150:0x218F),
    ("Arrows", 0x2190:0x21FF),
    ("Mathematical_Operators", 0x2200:0x22FF),
    ("Miscellaneous_Technical", 0x2300:0x23FF),
    ("Control_Pictures", 0x2400:0x243F),
    ("Optical_Character_Recognition", 0x2440:0x245F),
    ("Enclosed_Alphanumerics", 0x2460:0x24FF),
    ("Box_Drawing", 0x2500:0x257F),
    ("Block_Elements", 0x2580:0x259F),
    ("Geometric_Shapes", 0x25A0:0x25FF),
    ("Miscellaneous_Symbols", 0x2600:0x26FF),
    ("Dingbats", 0x2700:0x27BF),
    ("Miscellaneous_Mathematical_Symbols-A", 0x27C0:0x27EF),
    ("Supplemental_Arrows-A", 0x27F0:0x27FF),
    ("Braille_Patterns", 0x2800:0x28FF),
    ("Supplemental_Arrows-B", 0x2900:0x297F),
    ("Miscellaneous_Mathematical_Symbols-B", 0x2980:0x29FF),
    ("Supplemental_Mathematical_Operators", 0x2A00:0x2AFF),
    ("Miscellaneous_Symbols_and_Arrows", 0x2B00:0x2BFF),
    ("CJK_Radicals_Supplement", 0x2E80:0x2EFF),
    ("Kangxi_Radicals", 0x2F00:0x2FDF),
    ("Ideographic_Description_Characters", 0x2FF0:0x2FFF),
    ("CJK_Symbols_and_Punctuation", 0x3000:0x303F),
    ("Hiragana", 0x3040:0x309F),
    ("Katakana", 0x30A0:0x30FF),
    ("Bopomofo", 0x3100:0x312F),
    ("Hangul_Compatibility_Jamo", 0x3130:0x318F),
    ("Kanbun", 0x3190:0x319F),
    ("Bopomofo_Extended", 0x31A0:0x31BF),
    ("Katakana_Phonetic_Extensions", 0x31F0:0x31FF),
    ("Enclosed_CJK_Letters_and_Months", 0x3200:0x32FF),
    ("CJK_Compatibility", 0x3300:0x33FF),
    ("CJK_Unified_Ideographs_Extension_A", 0x3400:0x4DBF),
    ("Yijing_Hexagram_Symbols", 0x4DC0:0x4DFF),
    ("CJK_Unified_Ideographs", 0x4E00:0x9FFF),
    ("Yi_Syllables", 0xA000:0xA48F),
    ("Yi_Radicals", 0xA490:0xA4CF),
    ("Hangul_Syllables", 0xAC00:0xD7AF),
    ("High_Surrogates", 0xD800:0xDB7F),
    ("High_Private_Use_Surrogates", 0xDB80:0xDBFF),
    ("Low_Surrogates", 0xDC00:0xDFFF),
    ("Private_Use_Area", 0xE000:0xF8FF),
    ("CJK_Compatibility_Ideographs", 0xF900:0xFAFF),
    ("Alphabetic_Presentation_Forms", 0xFB00:0xFB4F),
    ("Arabic_Presentation_Forms-A", 0xFB50:0xFDFF),
    ("Variation_Selectors", 0xFE00:0xFE0F),
    ("Combining_Half_Marks", 0xFE20:0xFE2F),
    ("CJK_Compatibility_Forms", 0xFE30:0xFE4F),
    ("Small_Form_Variants", 0xFE50:0xFE6F),
    ("Arabic_Presentation_Forms-B", 0xFE70:0xFEFF),
    ("Halfwidth_and_Fullwidth_Forms", 0xFF00:0xFFEF),
    ("Specials", 0xFFF0:0xFFFF)
]

unicode_class = [
    x[1] => x[2:4] for x in [
        (:L, "Letter", "any kind of letter from any language", Base.Unicode.UTF8PROC_CATEGORY_LU:Base.Unicode.UTF8PROC_CATEGORY_LO),
        (:Ll, "Lowercase Letter", "a lowercase letter that has an uppercase variant", Base.Unicode.UTF8PROC_CATEGORY_LL),
        (:Lu, "Uppercase Letter", "an uppercase letter that has a lowercase variant", Base.Unicode.UTF8PROC_CATEGORY_LU),
        (:Lt, "Titlecase Letter", "a letter that appears at the start of a word when only the first letter of the word is capitalized", Base.Unicode.UTF8PROC_CATEGORY_LT),
        (Symbol("L&"), "Cased Letter", "a letter that exists in lowercase and uppercase variants (combination of Ll, Lu and Lt)", Base.Unicode.UTF8PROC_CATEGORY_LL:Base.Unicode.UTF8PROC_CATEGORY_LT),
        (:Lm, "Modifier Letter", "a special character that is used like a letter", Base.Unicode.UTF8PROC_CATEGORY_LM),
        (:Lo, "Other Letter", "a letter or ideograph that does not have lowercase and uppercase variants", Base.Unicode.UTF8PROC_CATEGORY_LO),
        (:M, "Mark", "a character intended to be combined with another character (e.g. accents, umlauts, enclosing boxes, etc.)", Base.Unicode.UTF8PROC_CATEGORY_MN:Base.Unicode.UTF8PROC_CATEGORY_ME),
        (:Mn, "Non Spacing Mark", "a character intended to be combined with another character without taking up extra space (e.g. accents, umlauts, etc.)", Base.Unicode.UTF8PROC_CATEGORY_MN),
        (:Mc, "Spacing Combining Mark", "a character intended to be combined with another character that takes up extra space (vowel signs in many Eastern languages)", Base.Unicode.UTF8PROC_CATEGORY_MC),
        (:Me, "Enclosing Mark", "a character that encloses the character it is combined with (circle, square, keycap, etc.)", Base.Unicode.UTF8PROC_CATEGORY_ME),
        (:Z, "Separator", "any kind of whitespace or invisible separator", Base.Unicode.UTF8PROC_CATEGORY_ZS:Base.Unicode.UTF8PROC_CATEGORY_ZP),
        (:Zs, "Space Separator", "a whitespace character that is invisible, but does take up space", Base.Unicode.UTF8PROC_CATEGORY_ZS),
        (:Zl, "Line Separator", "line separator character U+2028", Base.Unicode.UTF8PROC_CATEGORY_ZL),
        (:Zp, "Paragraph Separator", "paragraph separator character U+2029", Base.Unicode.UTF8PROC_CATEGORY_ZP),
        (:S, "Symbol", "math symbols, currency signs, dingbats, box-drawing characters, etc.", Base.Unicode.UTF8PROC_CATEGORY_SM:Base.Unicode.UTF8PROC_CATEGORY_SO),
        (:Sm, "Math Symbol", "any mathematical symbol", Base.Unicode.UTF8PROC_CATEGORY_SM),
        (:Sc, "Currency Symbol", "any currency sign", Base.Unicode.UTF8PROC_CATEGORY_SC),
        (:Sk, "Modifier Symbol", "a combining character (mark) as a full character on its own", Base.Unicode.UTF8PROC_CATEGORY_SK),
        (:So, "Other Symbol", "various symbols that are not math symbols, currency signs, or combining characters", Base.Unicode.UTF8PROC_CATEGORY_SO),
        (:N, "Number", "any kind of numeric character in any script", Base.Unicode.UTF8PROC_CATEGORY_ND:Base.Unicode.UTF8PROC_CATEGORY_NO),
        (:Nd, "Decimal Digit Number", "a digit zero through nine in any script except ideographic scripts", Base.Unicode.UTF8PROC_CATEGORY_ND),
        (:Nl, "Letter Number", "a number that looks like a letter, such as a Roman numeral", Base.Unicode.UTF8PROC_CATEGORY_NL),
        (:No, "Other Number", "a superscript or subscript digit, or a number that is not a digit 0–9 (excluding numbers from ideographic scripts)", Base.Unicode.UTF8PROC_CATEGORY_NO),
        (:P, "Punctuation", "any kind of punctuation character", Base.Unicode.UTF8PROC_CATEGORY_PC:Base.Unicode.UTF8PROC_CATEGORY_PO),
        (:Pc, "Connector Punctuation", "a punctuation character such as an underscore that connects words", Base.Unicode.UTF8PROC_CATEGORY_PC),
        (:Pd, "Dash Punctuation", "any kind of hyphen or dash", Base.Unicode.UTF8PROC_CATEGORY_PD),
        (:Ps, "Open Punctuation", "any kind of opening bracket", Base.Unicode.UTF8PROC_CATEGORY_PS),
        (:Pe, "Close Punctuation", "any kind of closing bracket", Base.Unicode.UTF8PROC_CATEGORY_PE),
        (:Pi, "Initial Punctuation", "any kind of opening quote", Base.Unicode.UTF8PROC_CATEGORY_PI),
        (:Pf, "Final Punctuation", "any kind of closing quote", Base.Unicode.UTF8PROC_CATEGORY_PF),
        (:Po, "Other Punctuation", "any kind of punctuation character that is not a dash, bracket, quote or connector", Base.Unicode.UTF8PROC_CATEGORY_PO),
        (:C, "Other", "invisible control characters and unused code points", Base.Unicode.UTF8PROC_CATEGORY_CC:Base.Unicode.UTF8PROC_CATEGORY_CO),
        (:Cc, "Control", "an ASCII or Latin-1 control character: 0x00–0x1F and 0x7F–0x9F", Base.Unicode.UTF8PROC_CATEGORY_CC),
        (:Cf, "Format", "invisible formatting indicator", Base.Unicode.UTF8PROC_CATEGORY_CF),
        (:Cs, "Surrogate", "one half of a surrogate pair in UTF-16 encoding", Base.Unicode.UTF8PROC_CATEGORY_CS),
        (:Co, "Private Use", "any code point reserved for private use", Base.Unicode.UTF8PROC_CATEGORY_CO),
        (:Cn, "Unassigned", "any code point to which no character has been assigned", Base.Unicode.UTF8PROC_CATEGORY_CN)
        ]
        ]

unicode_abbrev = Dict([ v[3]=>k for (k,v) in CombinedParsers.unicode_class])

unicode_classes = Dict(unicode_class...)

function in_any(x,sets)
    for s in sets
        x in s && return true
    end
    return false
end

export UnicodeClass
struct UnicodeClass{I}
    class::I
end

"""
    UnicodeClass(unicode_category::Symbol...)

used in [`ValueIn`](@ref), [`ValueNotIn`](@ref)
and succeeds if char at cursor is in one of the unicode classes.

```jldoctest
julia> match(ValueIn(:L), "aB")
ParseMatch("a")

julia> match(ValueIn(:Lu), "aB")
ParseMatch("B")

julia> match(ValueIn(:N), "aA1")
ParseMatch("1")
```

# Supported Unicode classes
```jldoctest
julia> for (k,v) in CombinedParsers.unicode_class
         println(":",k, " is a ",v[1],", ", v[2],".")
       end
:L is a Letter, any kind of letter from any language.
:Ll is a Lowercase Letter, a lowercase letter that has an uppercase variant.
:Lu is a Uppercase Letter, an uppercase letter that has a lowercase variant.
:Lt is a Titlecase Letter, a letter that appears at the start of a word when only the first letter of the word is capitalized.
:L& is a Cased Letter, a letter that exists in lowercase and uppercase variants (combination of Ll, Lu and Lt).
:Lm is a Modifier Letter, a special character that is used like a letter.
:Lo is a Other Letter, a letter or ideograph that does not have lowercase and uppercase variants.
:M is a Mark, a character intended to be combined with another character (e.g. accents, umlauts, enclosing boxes, etc.).
:Mn is a Non Spacing Mark, a character intended to be combined with another character without taking up extra space (e.g. accents, umlauts, etc.).
:Mc is a Spacing Combining Mark, a character intended to be combined with another character that takes up extra space (vowel signs in many Eastern languages).
:Me is a Enclosing Mark, a character that encloses the character it is combined with (circle, square, keycap, etc.).
:Z is a Separator, any kind of whitespace or invisible separator.
:Zs is a Space Separator, a whitespace character that is invisible, but does take up space.
:Zl is a Line Separator, line separator character U+2028.
:Zp is a Paragraph Separator, paragraph separator character U+2029.
:S is a Symbol, math symbols, currency signs, dingbats, box-drawing characters, etc..
:Sm is a Math Symbol, any mathematical symbol.
:Sc is a Currency Symbol, any currency sign.
:Sk is a Modifier Symbol, a combining character (mark) as a full character on its own.
:So is a Other Symbol, various symbols that are not math symbols, currency signs, or combining characters.
:N is a Number, any kind of numeric character in any script.
:Nd is a Decimal Digit Number, a digit zero through nine in any script except ideographic scripts.
:Nl is a Letter Number, a number that looks like a letter, such as a Roman numeral.
:No is a Other Number, a superscript or subscript digit, or a number that is not a digit 0–9 (excluding numbers from ideographic scripts).
:P is a Punctuation, any kind of punctuation character.
:Pc is a Connector Punctuation, a punctuation character such as an underscore that connects words.
:Pd is a Dash Punctuation, any kind of hyphen or dash.
:Ps is a Open Punctuation, any kind of opening bracket.
:Pe is a Close Punctuation, any kind of closing bracket.
:Pi is a Initial Punctuation, any kind of opening quote.
:Pf is a Final Punctuation, any kind of closing quote.
:Po is a Other Punctuation, any kind of punctuation character that is not a dash, bracket, quote or connector.
:C is a Other, invisible control characters and unused code points.
:Cc is a Control, an ASCII or Latin-1 control character: 0x00–0x1F and 0x7F–0x9F.
:Cf is a Format, invisible formatting indicator.
:Cs is a Surrogate, one half of a surrogate pair in UTF-16 encoding.
:Co is a Private Use, any code point reserved for private use.
:Cn is a Unassigned, any code point to which no character has been assigned.
```
"""
UnicodeClass(abbrev::Symbol...) =
    UnicodeClass(tuple((unicode_classes[a][3] for a in abbrev)...))
UnicodeClass(abbrev::String...) =
    UnicodeClass((Symbol(a) for a in abbrev)...)
_ismatch(x::Char, set::UnicodeClass{<:Tuple})::Bool =
    in_any(Base.Unicode.category_code(x),set.class)
_ismatch(x::Char, set::UnicodeClass)::Bool =
    in(Base.Unicode.category_code(x),set.class)

_regex_string(x::UnicodeClass) =
    if haskey(unicode_abbrev, x.class)
        "\\p{$(unicode_abbrev[x.class])}"
    else
        join([_regex_string(UnicodeClass(s)) for s in x.class])
    end


text:
module TestData
    importall Base # so that we get warnings for conflicts
    using Base.Test
    using DataFrames
    using Compat

    #test_group("DataVector creation")
    dvint = @data([1, 2, NA, 4])
    dvint2 = data([5:8;])
    dvint3 = data(5:8)
    dvflt = @data([1.0, 2, NA, 4])
    dvstr = @data(["one", "two", NA, "four"])
    dvdict = DataArray(Dict, 4)    # for issue #199

    #test_group("constructors")
    df1 = DataFrame(Any[dvint, dvstr], [:Ints, :Strs])
    df2 = DataFrame(Any[dvint, dvstr])
    df3 = DataFrame(Any[dvint])
    df4 = convert(DataFrame, [1:4 1:4])
    df5 = DataFrame(Any[@data([1,2,3,4]), dvstr])
    df6 = DataFrame(Any[dvint, dvint, dvstr], [:A, :B, :C])
    df7 = DataFrame(x = dvint, y = dvstr)
    @test size(df7) == (4, 2)
    @test isequal(df7[:x], dvint)

    #test_group("description functions")
    @test size(df6, 1) == 4
    @test size(df6, 2) == 3
    @test names(df6) == [:A, :B, :C]
    @test names(df2) == [:x1, :x2]
    @test names(df7) == [:x, :y]

    #test_group("ref")
    @test df6[2, 3] == "two"
    @test isna(df6[3, 3])
    @test df6[2, :C] == "two"
    @test isequal(df6[:B], dvint)
    @test size(df6[[2,3]], 2) == 2
    @test size(df6[2,:], 1) == 1
    @test size(df6[[1, 3], [1, 3]]) == (2, 2)
    @test size(df6[1:2, 1:2]) == (2, 2)
    @test size(head(df6,2)) == (2, 3)
    # lots more to do

    #test_group("assign")
    df6[3] = @data(["un", "deux", "troix", "quatre"])
    @test df6[1, 3] == "un"
    df6[:B] = [4, 3, 2, 1]
    @test df6[1,2] == 4
    df6[:D] = [true, false, true, false]
    @test df6[1,4] == true
    delete!(df6, :D)
    @test names(df6) == [:A, :B, :C]
    @test size(df6, 2) == 3

    #test_group("NA handling")
    @test nrow(df5[complete_cases(df5), :]) == 3

    #test_context("SubDataFrames")

    #test_group("constructors")
    # single index is rows
    sdf6a = sub(df6, 1)
    sdf6b = sub(df6, 2:3)
    sdf6c = sub(df6, [true, false, true, false])
    @test size(sdf6a) == (1,3)
    sdf6d = sub(df6, [1,3], :B)
    @test size(sdf6d) == (2,1)

    #test_group("ref")
    @test sdf6a[1,2] == 4

    #test_context("Within")
    #test_group("Associative")

    #test_group("DataFrame")
    srand(1)
    N = 20
    #Cast to Int64 as rand() behavior differs between Int32/64
    d1 = pdata(rand(@compat(map(Int64, 1:2)), N))
    d2 = (@pdata ["A", "B", NA])[rand(@compat(map(Int64, 1:3)), N)]
    d3 = data(randn(N))
    d4 = data(randn(N))
    df7 = DataFrame(Any[d1, d2, d3], [:d1, :d2, :d3])

    #test_group("groupby")
    gd = groupby(df7, :d1)
    @test length(gd) == 2
    # @test isequal(gd[2]["d2"], PooledDataVector["A", "B", NA, "A", NA, NA, NA, NA])
    @test sum(gd[2][:d3]) == sum(df7[:d3][dropna(df7[:d1] .== 2)])

    g1 = groupby(df7, [:d1, :d2])
    g2 = groupby(df7, [:d2, :d1])
    @test sum(g1[1][:d3]) == sum(g2[1][:d3])

    res = 0.0
    for x in g1
        res += sum(x[:d1])
    end
    @test res == sum(df7[:d1])

    df8 = aggregate(df7[[1, 3]], sum)
    @test df8[1, :d1_sum] == sum(df7[:d1])

    df8 = aggregate(df7, :d2, [sum, length])
    @test size(df8, 1) == 3
    @test size(df8, 2) == 5
    @test df8[2, :d1_length] == 4
    @test isequal(df8, aggregate(groupby(df7, :d2), [sum, length]))

    df9 = df7 |> groupby([:d2]) |> [sum, length]
    @test isequal(df9, df8)
    df9 = aggregate(df7, :d2, [sum, length])
    @test isequal(df9, df8)

    df10 = DataFrame(
        Any[[1:4;], [2:5;], ["a", "a", "a", "b" ], ["c", "d", "c", "d"]],
        [:d1, :d2, :d3, :d4]
    )

    gd = groupby(df10, [:d3])
    ggd = groupby(gd[1], [:d3, :d4]) # make sure we can groupby subdataframes
    @test ggd[1][1, :d3] == "a"
    @test ggd[1][1, :d4] == "c"
    @test ggd[1][2, :d3] == "a"
    @test ggd[1][2, :d4] == "c"
    @test ggd[2][1, :d3] == "a"
    @test ggd[2][1, :d4] == "d"

    #test_group("reshape")
    d1 = DataFrame(a = repeat([1:3;], inner = [4]),
                   b = repeat([1:4;], inner = [3]),
                   c = randn(12),
                   d = randn(12),
                   e = map(string, 'a':'l'))

    stack(d1, :a)
    d1s = stack(d1, [:a, :b])
    d1s2 = stack(d1, [:c, :d])
    d1s3 = stack(d1)
    d1m = melt(d1, [:c, :d, :e])
    @test isequal(d1s[1:12, :c], d1[:c])
    @test isequal(d1s[13:24, :c], d1[:c])
    @test isequal(d1s2, d1s3)
    @test names(d1s) == [:variable, :value, :c, :d, :e]
    @test isequal(d1s, d1m)
    d1m = melt(d1[[1,3,4]], :a)
    @test names(d1m) == [:variable, :value, :a]

    stackdf(d1, :a)
    d1s = stackdf(d1, [:a, :b])
    d1s2 = stackdf(d1, [:c, :d])
    d1s3 = stackdf(d1)
    d1m = meltdf(d1, [:c, :d, :e])
    @test isequal(d1s[1:12, :c], d1[:c])
    @test isequal(d1s[13:24, :c], d1[:c])
    @test isequal(d1s2, d1s3)
    @test names(d1s) == [:variable, :value, :c, :d, :e]
    @test isequal(d1s, d1m)
    d1m = meltdf(d1[[1,3,4]], :a)
    @test names(d1m) == [:variable, :value, :a]

    d1s[:id] = [1:12; 1:12]
    d1s2[:id] = [1:12; 1:12]
    d1us = unstack(d1s, :id, :variable, :value)
    d1us2 = unstack(d1s2)
    d1us3 = unstack(d1s2, :variable, :value)
    @test isequal(d1us[:a], d1[:a])
    @test isequal(d1us2[:d], d1[:d])
    @test isequal(d1us2[:3], d1[:d])

    d2 = DataFrame(id1 = [:a, :a, :a, :b],
                   id2 = [:A, :B, :B, :B],
                   id3 = [:t, :f, :t, :f],
                   val = [.1, .2, .3, .4])

    #test_group("merge")

    srand(1)
    df1 = DataFrame(a = shuffle!([1:10;]),
                    b = [:A,:B][rand(1:2, 10)],
                    v1 = randn(10))

    df2 = DataFrame(a = shuffle!(reverse([1:5;])),
                    b2 = [:A,:B,:C][rand(1:3, 5)],
                    v2 = randn(5))

    m1 = join(df1, df2, on = :a)
    @test isequal(m1[:a], @data([1, 2, 3, 4, 5]))
    # TODO: Re-enable
    # m2 = join(df1, df2, on = :a, kind = :outer)
    # @test isequal(m2[:b2], DataVector["A", "B", "B", "B", "B", NA, NA, NA, NA, NA])
    # @test isequal(m2[:b2], DataVector["B", "B", "B", "C", "B", NA, NA, NA, NA, NA])

    df1 = DataFrame(a = [1, 2, 3],
                    b = ["America", "Europe", "Africa"])
    df2 = DataFrame(a = [1, 2, 4],
                    c = ["New World", "Old World", "New World"])

    m1 = join(df1, df2, on = :a, kind = :inner)
    @test isequal(m1[:a], @data([1, 2]))

    m2 = join(df1, df2, on = :a, kind = :left)
    @test isequal(m2[:a], @data([1, 2, 3]))

    m3 = join(df1, df2, on = :a, kind = :right)
    @test isequal(m3[:a], @data([1, 2, 4]))

    m4 = join(df1, df2, on = :a, kind = :outer)
    @test isequal(m4[:a], @data([1, 2, 3, 4]))

    # test with NAs (issue #185)
    df1 = DataFrame()
    df1[:A] = @data(["a", "b", "a", NA])
    df1[:B] = @data([1, 2, 1, 3])

    df2 = DataFrame()
    df2[:A] = @data(["a", NA, "c"])
    df2[:C] = @data([1, 2, 4])

    m1 = join(df1, df2, on = :A)
    @test size(m1) == (3,3)
    @test isequal(m1[:A], @data([NA,"a","a"]))

    m2 = join(df1, df2, on = :A, kind = :outer)
    @test size(m2) == (5,3)
    @test isequal(m2[:A], @data([NA,"a","a","b","c"]))

    srand(1)
    df1 = DataFrame(
        a = [:x,:y][rand(1:2, 10)],
        b = [:A,:B][rand(1:2, 10)],
        v1 = randn(10)
    )

    df2 = DataFrame(
        a = [:x,:y][[1,2,1,1,2]],
        b = [:A,:B,:C][[1,1,1,2,3]],
        v2 = randn(5)
    )
    df2[1,:a] = NA

    # # TODO: Restore this functionality
    # m1 = join(df1, df2, on = [:a,:b])
    # @test isequal(m1[:a], DataArray(["x", "x", "y", "y", fill("x", 5)]))
    # m2 = join(df1, df2, on = ["a","b"], kind = :outer)
    # @test isequal(m2[10,:v2], NA)
    # @test isequal(m2[:a], DataVector["x", "x", "y", "y", "x", "x", "x", "x", "x", "y", NA, "y"])

    srand(1)
    function spltdf(d)
        d[:x1] = map(x -> x[1], d[:a])
        d[:x2] = map(x -> x[2], d[:a])
        d[:x3] = map(x -> x[3], d[:a])
        d
    end
    df1 = DataFrame(
        a = ["abc", "abx", "axz", "def", "dfr"],
        v1 = randn(5)
    )
    df1 = spltdf(df1)
    df2 = DataFrame(
        a = ["def", "abc","abx", "axz", "xyz"],
        v2 = randn(5)
    )
    df2 = spltdf(df2)

    # m1 = join(df1, df2, on = :a)
    # m2 = join(df1, df2, on = [:x1, :x2, :x3])
    # @test isequal(sort(m1[:a]), sort(m2[:a]))

    #test_group("New DataVector constructors")
    dv = DataArray(Int, 5)
    @test all(isna(dv))
    dv = DataArray(Float64, 5)
    @test all(isna(dv))
    dv = @data(zeros(5))
    @test all(dv .== 0.0)
    dv = @data(ones(5))
    @test all(dv .== 1.0)

    # No more NA corruption
    dv = @data(ones(10_000))
    @test !any(isna(dv))

    PooledDataArray(falses(2), falses(2))
    PooledDataArray(falses(2), trues(2))

    # Test vectorized comparisons work for DataVector's and PooledDataVector's
    @data([1, 2, NA]) .== 1
    @pdata([1, 2, NA]) .== 1
    @data(["1", "2", NA]) .== "1"
    @pdata(["1", "2", NA]) .== "1"

    # Test unique()
    #test_group("unique()")
    # TODO: Restore this
    # dv = DataArray(1:4)
    # dv[4] = NA
    # @test (1 in unique(dv))
    # @test (2 in unique(dv))
    # @test (3 in unique(dv))
    # @test (NA in unique(dv))

    # test nonunique() with extra argument
    df1 = DataFrame(a = ["a", "b", "a", "b", "a", "b"], b = 1:6, c = [1:3;1:3])
    df = vcat(df1, df1)
    @test find(nonunique(df)) == collect(7:12)
    @test find(nonunique(df, :)) == collect(7:12)
    @test find(nonunique(df, Colon())) == collect(7:12)
    @test find(nonunique(df, :a)) == collect(3:12)
    @test find(nonunique(df, [:a, :c])) == collect(7:12)
    @test find(nonunique(df, [1, 3])) == collect(7:12)
    @test find(nonunique(df, 1)) == collect(3:12)

    # Test unique() with extra argument
    @test unique(df) == df1
    @test unique(df, :) == df1
    @test unique(df, Colon()) == df1
    @test unique(df, 2:3) == df1
    @test unique(df, 3) == df1[1:3,:]
    @test unique(df, [1, 3]) == df1
    @test unique(df, [:a, :c]) == df1
    @test unique(df, :a) == df1[1:2,:]

    #test unique!() with extra argument
    unique!(df, [1, 3])
    @test df == df1

    #test_group("find()")
    dv = DataArray([true, false, true])
    @test isequal(find(dv), [1, 3])

    pdv = PooledDataArray([true, false, true])
    @test isequal(find(pdv), [1, 3])

    dv[1] = NA
    @test isequal(find(dv), [3])

    pdv[1] = NA
    @test isequal(find(pdv), [3])
end


text:
## Loading packages.
using MicroHH
using Tullio
using Statistics

## Loading settings.
include("drycbl_settings.jl")

## Initialize the model.
n_domains = 2
m = Model("drycbl", n_domains, settings, float_type)

## Create the initials fields.
for i in 1:n_domains
    f = m.fields[i]; g = m.grid[i]
    s = @view f.s[g.is:g.ie, g.js:g.je, g.ks:g.ke]
    z = @view g.z[g.ks:g.ke]
    rand2d = rand(g.itot, g.jtot)
    rand2d .-= mean(rand2d)
    s[:, :, 1] .+= rand2d[:, :]
    f.s_gradbot[:, :] .= - 0.1 / settings[i]["fields"]["visc"]
    f.s_gradtop[:, :] .= 0.003
    @tullio s[i, j, k] += 0.003 * z[k]
end

## Save the restart files.
save_model(m)


text:
using Mosek

# Demonstrates how to solve a linear optimization problem using the
# MOSEK API and modify and re-optimize the problem.

function production(productNames     :: Vector{String},
                    processNames     :: Vector{String},
                    timeRequirements :: Array{Float64,2},
                    timeResources    :: Vector{Float64},
                    profit           :: Vector{Float64},
                    # Modify A entry
                    mod_i               :: Int,
                    mod_j               :: Int,
                    mod_cof             :: Float64,

                    # Add product Name
                    addProductName             :: String,
                    addProductTimeRequirements :: Vector{Float64},
                    addProfit                  :: Float64,
                    # Add product Name
                    addProcessName      :: String,
                    addTimeRequirements :: Vector{Float64},
                    addTimeResources    :: Float64 )
    n = length(productNames)
    m = length(processNames)
    maketask() do t
        appendvars(t,n)
        appendcons(t,m)

        for j in 1:n
            putvarname(t,Int32(j),productNames[j])
            putvarbound(t,Int32(j),MSK_BK_LO,0.0,Inf)
        end

        for i in 1:m
            putconname(t,Int32(i),processNames[i])
            putconbound(t,i, MSK_BK_UP, -Inf,timeResources[i])
            putarow(t,i,Int32[1:n...],timeRequirements[:,i])
        end

        putobjsense(t,MSK_OBJECTIVE_SENSE_MAXIMIZE)
        putclist(t,Int32[1:n...],profit)

        optimize(t)
        xx = getxx(t,MSK_SOL_BAS)
        println("Initial solution:")
        for i in 1:n
            @printf("    %-15s: %10.1f\n",productNames[i],xx[i])
        end

        ############### Make a change to the A matrix ###############
        putaij(t,Int32(mod_i),Int32(mod_j),mod_cof)

        optimize(t)
        xx = getxx(t,MSK_SOL_BAS)
        println("After modified coefficient:")
        for i in 1:n
            @printf("    %-15s: %10.1f\n",productNames[i],xx[i])
        end

        ############## Add a new variable / column ##################
        appendvars(t,Int32(1))
        putvarname(t,Int32(n+1),addProductName)
        putcj(t,Int32(n+1),addProfit)
        putacol(t,Int32(n+1),Int32[1:n+1...],addProductTimeRequirements)
        putvarbound(t,Int32(n+1),MSK_BK_LO,0.0,Inf)

        # reoptimize using simplex
        putintparam(t,MSK_IPAR_OPTIMIZER, MSK_OPTIMIZER_FREE_SIMPLEX)

        optimize(t)
        xx = getxx(t,MSK_SOL_BAS)
        println("After added product (column):")
        for i in 1:n
            @printf("    %-15s: %10.1f\n",productNames[i],xx[i])
        end
        @printf("    %-15s: %10.1f\n",addProductName,xx[n+1])

        ################ Add a new constraint / row #################
        appendcons(t,Int32(1))
        putconname(t,Int32(m+1),addProcessName)
        putarow(t,Int32(m+1),Int32[1:n...],addTimeRequirements)
        putconbound(t,Int32(n+1),MSK_BK_UP,-Inf,addTimeResources)

        optimize(t)
        xx = getxx(t,MSK_SOL_BAS)
        println("After added process (row):")
        for i in 1:n
            @printf("    %-15s: %10.1f\n",productNames[i],xx[i])
        end
        @printf("    %-15s: %10.1f\n",addProductName,xx[n+1])
    end
end

timeRequirements = Float64[ 2.0  3.0  2.0
                            4.0  2.0  3.0
                            3.0  3.0  2.0 ]
timeResources    = Float64[100000.0, 50000.0, 60000.0]
profit           = Float64[1.5, 2.5, 3.0]
processNames     = String["Assembly", "Polish","Packing" ]
productNames     = String["Chairs", "Tables", "Beds"]

production(productNames,processNames,timeRequirements,timeResources,profit,
           1,1,3.0, # modify a coefficient
           "Cutting_board", Float64[4.0,0.0,1.0], 1.0, # additional product
           "Varnish", Float64[1.0,2.0,1.0,1.0], 30000.0) # additional process


text:
"""
    Methods

List of methods available.

- `ExtensiveForm`: deterministic equivalent formulation solve
- `Benders`: Benders decomposition
- `DW`: Dantzig-Wolfe decomposition
- `Dual`: Dual decomposition
"""
@enum Methods begin
    ExtensiveForm
    Benders
    DW
    Dual
end

mutable struct DSPProblem
    p::Ptr{Cvoid}

    numRows::Dict{Int, Int}
    numCols::Dict{Int, Int}
    primVal::Float64
    dualVal::Float64
    colVal::Dict{Int, Vector{Float64}}
    rowVal::Vector{Float64}
    objective_sense::Float64
    status::Int
    solve_time::Float64

    # Number of blocks
    nblocks::Int

    # Array of block ids:
    # The size of array is not necessarily same as nblocks,
    # as block ids may be distributed to multiple processors.
    block_ids::Vector{Integer}

    is_stochastic::Bool

    dro

    # solve_type should be one of these:
    solve_type::Methods

    # MPI settings
    comm
    comm_size::Int
    comm_rank::Int

    # linear and quadratic constraints of subproblems
    is_qc::Dict{Int64,Bool}
    quadConstrs::Dict{Int64, Dict{Int64,AbstractConstraint}}
    linConstrs::Dict{Int64, Dict{Int64,AbstractConstraint}}

    function DSPProblem(fake::Bool = false)
        prob = new(
            C_NULL, # p
            Dict(), # numRows
            Dict(), # numCols
            NaN, # primVal
            NaN, # dualVal
            Dict(), # colVal
            [], # rowVal
            1., # objective_sense
            3998, # status
            0., # solve_time
            -1, # nblocks
            [], # block_ids
            false, # is_stochastic
            nothing, # dro
            DW, # solve_type
            nothing, # comm
            1, # comm_size
            0, # comm_rank
            Dict(), # is_qc
            Dict(), # quadratic constraints
            Dict() # linear constraints
        )
        if !fake
            prob.p = createEnv()
            finalizer(freeEnv, prob)
        end

        return prob
    end
end

###############################################################################
# Help functions
###############################################################################

macro dsp_ccall(func, args...)
    @static if !Sys.iswindows()
        return esc(quote
            ccall(($func, "libDsp"), $(args...))
        end)
    end
    @static if Sys.iswindows()
        return esc(quote
            ccall(($func, "libDsp"), stdcall, $(args...))
        end)
    end
end

createEnv() = @dsp_ccall("createEnv", Ptr{Cvoid}, ())

function freeEnv(dsp::DSPProblem)
    if dsp.p != C_NULL
        freeModel(dsp)
        @dsp_ccall("freeEnv", Cvoid, (Ref{Ptr{Cvoid}},), dsp.p)
        dsp.p = C_NULL
    end
    dsp.comm = nothing
    dsp.comm_size = 1
    dsp.comm_rank = 0
end

function freeModel(dsp::DSPProblem)
    @dsp_ccall("freeModel", Cvoid, (Ptr{Cvoid},), dsp.p)
    dsp.numRows = Dict()
    dsp.numCols = Dict()
    dsp.primVal = NaN
    dsp.dualVal = NaN
    dsp.colVal = Dict()
    dsp.rowVal = []
    dsp.objective_sense = 1.
    dsp.status = 3998
    dsp.solve_time = 0.
    dsp.nblocks = -1
    dsp.block_ids = []
    dsp.is_stochastic = false
    dsp.solve_type = Dual
    dsp.is_qc = Dict()
    dsp.quadConstrs = Dict()
    dsp.linConstrs = Dict()
end

###############################################################################
# Load problems
###############################################################################

readParamFile(dsp::DSPProblem, param_file::AbstractString) = @dsp_ccall(
    "readParamFile", Cvoid,
    (Ptr{Cvoid}, Ptr{UInt8}),
    dsp.p, param_file)

readSmps(dsp::DSPProblem, filename::AbstractString) = @dsp_ccall(
    "readSmps", Cvoid,
    (Ptr{Cvoid}, Ptr{UInt8}),
    dsp.p, filename)

loadFirstStage(dsp::DSPProblem, start, index, value, clbd, cubd, ctype, obj, rlbd, rubd) = @dsp_ccall(
    "loadFirstStage", Cvoid,
    (Ptr{Cvoid}, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{UInt8}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
    dsp.p, start, index, value, clbd, cubd, ctype, obj, rlbd, rubd)

loadQCQPFirstStage(dsp::DSPProblem, start, index, value, clbd, cubd, ctype, obj, qobjrowindex, qobjcolindex, qobjvalue, qobjnum, rlbd, rubd, nqrows::Int, linnzcnt::Vector{Int}, quadnzcnt::Vector{Int}, rhs::Vector{Float64}, sense::Vector{Int}, linstart::Vector{Int}, linind::Vector{Int}, linval::Vector{Float64}, quadstart::Vector{Int}, quadrow::Vector{Int}, quadcol::Vector{Int}, quadval::Vector{Float64}) = @dsp_ccall(
    "loadQCQPFirstStage", Cvoid,
    (Ptr{Cvoid}, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{UInt8}, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint},Ptr{Cint},Ptr{Cdouble}),
    dsp.p, start, index, value, clbd, cubd, ctype, obj, qobjrowindex, qobjcolindex, qobjvalue, qobjnum, rlbd, rubd, nqrows, convert(Vector{Cint}, linnzcnt), convert(Vector{Cint}, quadnzcnt), rhs, convert(Vector{Cint}, sense), convert(Vector{Cint}, linstart), convert(Vector{Cint}, linind), linval, convert(Vector{Cint}, quadstart), convert(Vector{Cint}, quadrow), convert(Vector{Cint}, quadcol), quadval)

loadSecondStage(dsp::DSPProblem, id, probability, start, index, value, clbd, cubd, ctype, obj, rlbd, rubd) = @dsp_ccall(
    "loadSecondStage", Cvoid,
    (Ptr{Cvoid}, Cint, Cdouble, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{UInt8}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
    dsp.p, id, probability, start, index, value, clbd, cubd, ctype, obj, rlbd, rubd)

loadQCQPSecondStage(dsp::DSPProblem, id, probability, start, index, value, clbd, cubd, ctype, obj, qobjrowindex, qobjcolindex, qobjvalue, qobjnum, rlbd, rubd, nqrows::Int, linnzcnt::Vector{Int}, quadnzcnt::Vector{Int}, rhs::Vector{Float64}, sense::Vector{Int}, linstart::Vector{Int}, linind::Vector{Int}, linval::Vector{Float64}, quadstart::Vector{Int}, quadrow::Vector{Int}, quadcol::Vector{Int}, quadval::Vector{Float64}) = @dsp_ccall(
    "loadQCQPSecondStage", Cvoid,
    (Ptr{Cvoid}, Cint, Cdouble, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{UInt8}, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint},Ptr{Cint},Ptr{Cdouble}),
    dsp.p, id, probability, start, index, value, clbd, cubd, ctype, obj, qobjrowindex, qobjcolindex, qobjvalue, qobjnum, rlbd, rubd, nqrows, convert(Vector{Cint}, linnzcnt), convert(Vector{Cint}, quadnzcnt), rhs, convert(Vector{Cint}, sense), convert(Vector{Cint}, linstart), convert(Vector{Cint}, linind), linval, convert(Vector{Cint}, quadstart), convert(Vector{Cint}, quadrow), convert(Vector{Cint}, quadcol), quadval)

loadBlockProblem(dsp::DSPProblem, id, ncols, nrows, numels, start, index, value, clbd, cubd, ctype, obj, rlbd, rubd) = @dsp_ccall(
    "loadBlockProblem", Cvoid, (
    Ptr{Cvoid}, Cint, Cint, Cint, Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{UInt8}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
    dsp.p, id, ncols, nrows, numels, start, index, value, clbd, cubd, ctype, obj, rlbd, rubd)

updateBlocks(dsp::DSPProblem) = @dsp_ccall("updateBlocks", Cvoid, (Ptr{Cvoid},), dsp.p)

###############################################################################
# solve functions
###############################################################################

for func in [:freeSolver,
             :solveDe,
             :solveBd,
             :solveDd,
             :solveDw]
    strfunc = string(func)
    @eval begin
        function $func(dsp::DSPProblem)
            return @dsp_ccall($strfunc, Cvoid, (Ptr{Cvoid},), dsp.p)
        end
    end
end

for func in [:solveBdMpi, :solveDdMpi, :solveDwMpi]
    strfunc = string(func)
    if @isdefined(MPI)
        @eval begin
            function $func(dsp::DSPProblem)
                return @dsp_ccall($strfunc, Cvoid, (Ptr{Cvoid}, MPI.MPI_Comm), dsp.p, dsp.comm)
	        end
    	end
    else
        @eval begin
            function $func(dsp::DSPProblem)
                error("MPI package is required to use this function.")
			end
		end
	end
end

###############################################################################
# Get functions
###############################################################################

for (func,rtn) in [(:getNumScenarios, Cint),
                   (:getNumSubproblems, Cint),
                   (:getTotalNumRows, Cint),
                   (:getTotalNumCols, Cint),
                   (:getStatus, Cint),
                   (:getNumIterations, Cint),
                   (:getNumNodes, Cint),
                   (:getWallTime, Cdouble),
                   (:getPrimalBound, Cdouble),
                   (:getDualBound, Cdouble),
                   (:getNumCouplingRows, Cint),
                   (:getVersionMajor, Cint),
                   (:getVersionMinor, Cint),
                   (:getVersionPatch, Cint)]
    strfunc = string(func)
    @eval begin
        function $func(dsp::DSPProblem)
            @assert(dsp.p != C_NULL)
            return @dsp_ccall($strfunc, $rtn, (Ptr{Cvoid},), dsp.p)
        end
    end
end

function getSolution(dsp::DSPProblem, num::Integer)
    sol = zeros(num)
    if dsp.comm_rank == 0
        @dsp_ccall("getPrimalSolution", Cvoid, (Ptr{Cvoid}, Cint, Ptr{Cdouble}), dsp.p, num, sol)
    end
    return sol
end
getSolution(dsp::DSPProblem) = getSolution(dsp, getTotalNumCols(dsp))

function getDualSolution(dsp::DSPProblem, num::Integer)
	sol = zeros(num)
    if dsp.comm_rank == 0
        @dsp_ccall("getDualSolution", Cvoid, (Ptr{Cvoid}, Cint, Ptr{Cdouble}), dsp.p, num, sol)
    end
    return sol
end
getDualSolution(dsp::DSPProblem) = getDualSolution(dsp, getNumCouplingRows(dsp))

function getVersion(dsp::DSPProblem)
    major = getVersionMajor(dsp)
    minor = getVersionMinor(dsp)
    patch = getVersionPatch(dsp)
    return "$(major).$(minor).$(patch)"
end

###############################################################################
# Set functions
###############################################################################

setNumberOfScenarios(dsp::DSPProblem, nscen::Int) = @dsp_ccall(
    "setNumberOfScenarios", Cvoid,
    (Ptr{Cvoid}, Cint),
    dsp.p, convert(Cint, nscen))

setDimensions(dsp::DSPProblem, ncols1::Int, nrows1::Int, ncols2::Int, nrows2::Int) = @dsp_ccall(
    "setDimensions", Cvoid,
    (Ptr{Cvoid}, Cint, Cint, Cint, Cint),
    dsp.p, convert(Cint, ncols1), convert(Cint, nrows1), convert(Cint, ncols2), convert(Cint, nrows2))

setIntPtrParam(dsp::DSPProblem, name::String, n::Int, v::Vector{Int}) = @dsp_ccall(
    "setIntPtrParam", Cvoid,
    (Ptr{Cvoid}, Ptr{UInt8}, Cint, Ptr{Cint}),
    dsp.p, name, convert(Cint, n), convert(Vector{Cint}, v))

###############################################################################
# Other functions
###############################################################################
writeMps!(dsp::DSPProblem, filename::AbstractString) = @dsp_ccall(
    "writeMps", Cvoid,
    (Ptr{Cvoid}, Ptr{UInt8}),
    dsp.p, "$filename")

# end # end of module


text:
"""
    definesignal(signal_type::SignalType, f_s, t_length; prn=1,
                 f_if=0., f_d=0., fd_rate=0., Tsys=535.,
                 CN0=45., phi=0., nADC=4, include_carrier=true,
                 include_adc=true, include_thermal_noise=true,
                 code_start_idx=1., include_databits_I=true,
                 include_databits_Q=true, include_phase_noise=true,
                 phase_noise_scaler=1/10, name="custom",
                 skip_noise_generation=false, allocate_noise_vectors=true,
                 receiver_h_parms=h_parms_tcxo[1])

Define properties of locally generated generic signal
based off its type, PRN, etc.

Required Arguments:

- `signal_type::SignalType`: user defined signal type
    * see `definecodetype` and `definesignaltype`
- `f_s`: signal sampling frequency in Hz
- `t_length`: signal length in seconds

Optional Arguments:

- `prn`: satellite PRN number `(default = 1)`
- `f_if`: IF frequency in Hz `(default = 0Hz)`
- `f_d`: signal Doppler frequency in Hz `(default = 0Hz)`
- `fd_rate`: signal Doppler rate in Hz/s `(default = 0Hz/s)`
- `Tsys`: receiver noise temperature in Kelvin `(default = 535K)`
- `CN0`: signal carrier-to-noise ratio (C/N₀) `(default = 45dB⋅Hz)`
- `phi`: initial carrier phase in rad `(default = 0rad)`
- `nADC`: receiver bit depth `(default = 8 bits)`
- `include_carrier`: flag for modulating codes onto carrier `(default = true)`
- `include_adc`: flag for performing ADC quantization on signal `(default = true)`
- `include_thermal_noise`: flag for including thermal noise in signal generation
                           `(default = true)`
- `code_start_idx`: starting index in `ReplicaSignal.data` where all codes
                    start `(default = 1)`
- `include_databits_I`: flag for including I channel databits in simulated signal
                        `(default = true)`
- `include_databits_Q`: flag for including Q channel databits in simulated signal
                        `(default = true)`
- `include_phase_noise`: flag for including phase noise `(default = true)`
- `phase_noise_scaler`: standard deviation of phase noise in rad
                        `(default = 1/10rad)`
- `name`: name of signal `(default = custom)`
- `skip_noise_generation`: flag for skipping thermal and phase noise generation
                           `(default = false)`
- `allocate_noise_vectors`: flag for allocating memory for thermal and phase
                            noise vectors `(default = true)`
    * if set to `true`, thermal and phase noise vectors have the same length as
      `ReplicaSignal.data` vector
    * if set to `false`, thermal and phase noise vectors have length set to 0
      and the `include_thermal_noise` and `include_phase_noise` flags are set
      to `false`
- `receiver_h_parms`: vector of oscillator h parameters for use in simulating
                      oscillator phase noise
    * format is `[h₋₂, h₋₁, h₀, h₁, h₂]`
- `start_t`: the start time of the signal in seconds `(default = 0 seconds)`

Returns:

- `ReplicaSignal` struct
"""
function definesignal(signal_type::SignalType, f_s, t_length; prn=1,
                      f_if=0., f_d=0., fd_rate=0., Tsys=535.,
                      CN0=45., phi=0., nADC=8, include_carrier=true,
                      include_adc=true, include_thermal_noise=true,
                      code_start_idx=1., include_databits_I=true,
                      include_databits_Q=true, include_phase_noise=true,
                      phase_noise_scaler=1/10, name="custom",
                      skip_noise_generation=false, allocate_noise_vectors=true,
                      receiver_h_parms=h_parms_tcxo[4], start_t=0)
    # Obtain sample number based off the sampling frequency and duration of
    # signal
    sample_num = floor(Int, f_s*t_length)
    # Calculate code chipping rates with Doppler applied for all codes on
    # each channel and their initial code phases
    sig_freq = signal_type.sig_freq
    I_codes = signal_type.I_codes
    Q_codes = signal_type.Q_codes
    # Get adjusted code chipping rates and chipping rate rates for I channel
    # codes
    if signal_type.include_I
        f_code_d_I = Array{Float64}(undef, I_codes.code_num)
        f_code_dd_I = Array{Float64}(undef, I_codes.code_num)
        init_code_phases_I = Array{Float64}(undef, I_codes.code_num)
        for i in 1:signal_type.I_codes.code_num
            # I channel
            f_code = I_codes.chipping_rates[i]
            code_length = I_codes.code_lengths[i]
            f_code_d, f_code_dd = calc_doppler_code_rate(f_code, sig_freq,
                                                         f_d, fd_rate)
            f_code_d_I[i] = f_code_d
            f_code_dd_I[i] = f_code_dd
            init_code_phases_I[i] = calcinitcodephase(code_length, f_code_d,
                                                      f_code_dd, f_s, code_start_idx)
        end
        if I_codes.databits
            if ~include_databits_I
                I_codes.include_codes[end] = false
            else
                I_codes.include_codes[end] = true
            end
        end
    else
        # Set these vectors to zero length if there are no I channel codes
        f_code_d_I = Array{Float64}(undef, 0)
        f_code_dd_I = Array{Float64}(undef, 0)
        init_code_phases_I = Array{Float64}(undef, 0)
    end
    # Get adjusted code chipping rates and chipping rate rates for Q channel
    # codes
    if signal_type.include_Q
        f_code_d_Q = Array{Float64}(undef, Q_codes.code_num)
        f_code_dd_Q = Array{Float64}(undef, Q_codes.code_num)
        init_code_phases_Q = Array{Float64}(undef, Q_codes.code_num)
        for i in 1:signal_type.Q_codes.code_num
            # Q channel
            f_code = Q_codes.chipping_rates[i]
            code_length = Q_codes.code_lengths[i]
            f_code_d, f_code_dd = calc_doppler_code_rate(f_code, sig_freq,
                                                         f_d, fd_rate)
            f_code_d_Q[i] = f_code_d
            f_code_dd_Q[i] = f_code_dd
            init_code_phases_Q[i] = calcinitcodephase(code_length, f_code_d,
                                                      f_code_dd, f_s,
                                                      code_start_idx)
        end
        if Q_codes.databits
            if ~include_databits_Q
                Q_codes.include_codes[end] = false
            else
                Q_codes.include_codes[end] = true
            end
        end
    else
        # Set these vectors to zero length if there are no Q channel codes
        f_code_d_Q = Array{Float64}(undef, 0)
        f_code_dd_Q = Array{Float64}(undef, 0)
        init_code_phases_Q = Array{Float64}(undef, 0)
    end
    # Allocate space for signal. No simulated data is stored in yet. After
    # define the signal, user can use `generatesignal!(signal)` to generate
    # the I/Q samples which will be stored in `ReplicaSignal.data`.
    data = Array{Complex{Float64}}(undef, sample_num)
    isreplica = false
    noexp = false
    # Generate thermal noise and phase noise
    if allocate_noise_vectors
        if skip_noise_generation
            # Space for thermal and phase noise vectors are allocated. User can
            # use `definesignal!(signal; new_thermal_noise=true, new_phase_noise=true)`
            # to generate the noise in the vectors.
            thermal_noise = Array{Complex{Float64}}(undef, sample_num)
            phase_noise = Array{Float64}(undef, sample_num)
        else
            # Generate thermal and phase noise. Thermal noise is unscaled
            # but is scaled in `generatesignal!` when the signal is generated.
            # Phase noise is scaled based off the value of `phase_noise_scaler`.
            thermal_noise = randn(Complex{Float64}, sample_num)
            phase_noise = generate_phase_noise(t_length, f_s,
                                               signal_type.sig_freq,
                                               receiver_h_parms)
        end
    else
        # Create 0 sized thermal and phase noise vectors and set the
        # `include_thermal_noise` and `include_phase_noise` flags to `false`.
        thermal_noise = Array{Complex{Float64}}(undef, 0)
        phase_noise = Array{Float64}(undef, 0)
        include_thermal_noise = false
        include_phase_noise = false
    end
    return ReplicaSignal(name, prn, f_s, t_length, f_if, f_d, fd_rate, Tsys,
                         CN0, phi, nADC, code_start_idx, init_code_phases_I,
                         init_code_phases_Q, data, include_carrier,
                         include_adc, include_thermal_noise, include_databits_I,
                         include_databits_Q, include_phase_noise, f_code_d_I,
                         f_code_dd_I, f_code_d_Q, f_code_dd_Q, sample_num,
                         isreplica, noexp, thermal_noise, phase_noise,
                         signal_type, receiver_h_parms, start_t)
end

"""
    definesignal!(signal::ReplicaSignal;
                  prn=signal.prn, f_if=signal.f_if, f_d=signal.f_d,
                  fd_rate=signal.fd_rate, Tsys=signal.Tsys,
                  CN0=signal.CN0, phi=signal.phi, nADC=signal.nADC,
                  include_carrier=signal.include_carrier,
                  include_adc=signal.include_adc,
                  include_thermal_noise=signal.include_thermal_noise,
                  code_start_idx=signal.code_start_idx,
                  include_databits_I=signal.include_databits_I,
                  include_databits_Q=signal.include_databits_Q,
                  include_phase_noise=signal.include_phase_noise,
                  phase_noise_scaler=1/10, name=signal.name,
                  new_thermal_noise=false, new_phase_noise=false,
                  isreplica=signal.isreplica, noexp=signal.noexp,
                  new_databits=false)

Redefine properties of locally generated generic signal
based off its type, PRN, etc.

Required Arguments:

- `signal::ReplicaSignal`: the signal already defined by the user using
                            `definesignal`

Optional Arguments with Defaults Equal to `signal` Field Values:

- `prn`: satellite PRN number
- `f_if`: IF frequency in Hz
- `f_d`: signal Doppler frequency in Hz
- `fd_rate`: signal Doppler rate in Hz/s
- `Tsys`: receiver noise temperature in Kelvin
- `CN0`: signal carrier-to-noise ratio (C/N₀)
- `phi`: initial carrier phase in rad
- `nADC`: receiver bit depth
- `include_carrier`: flag for modulating codes onto carrier
- `include_adc`: flag for performing ADC quantization on signal
- `include_thermal_noise`: flag for including thermal noise in signal generation
- `code_start_idx`: starting index in `ReplicaSignal.data` where all codes
                    start
- `include_databits_I`: flag for including I channel databits in simulated signal
- `include_databits_Q`: flag for including Q channel databits in simulated signal
- `include_phase_noise`: flag for including phase noise
- `name`: name of signal
- `isreplica`: `[DEPRICATED]` flag to identify whether signal is being for processing
    * replica signals being used for processing simulated or real signals
      should not include any noise source
    * this is usually set to `true` when being used by processing functions
      such as `courseacquisition!`, `fineacquisition`, and `trackprn`, only if
      this signal structure is being used as the replica signal, and is NOT
      the signal being processed
    * if `true`, a `ReplicaSignal` struct will not have noise added onto it
      and will not undergo ADC quantization
    * signal generation will be done using the second method of `generatesignal!`,
      `generatesignal!(signal::ReplicaSignal, isreplica::Bool)`
- `noexp`: used only if second method of `generatesignal!`, discussed
           above, is used
    * does not modulate codes onto carrier

Other Optional Arguments:

- `phase_noise_init_phase`: initial phase for phase noise `(default = 0 rad)`
- `new_thermal_noise`: flag to generate new phase noise for signal
                       `(default = false)`
- `new_phase_noise`: flag to generate new phase noise for signal
                     `(default = false)`
- `new_databits`: flag to generate new databits for signal `(default = false)`
- `start_t`: the start time of the signal in seconds `(default = 0 seconds)`

Modifies and Returns:

- `signal::ReplicaSignal`
"""
function definesignal!(signal::ReplicaSignal;
                       prn=signal.prn, f_if=signal.f_if, f_d=signal.f_d,
                       fd_rate=signal.fd_rate, Tsys=signal.Tsys,
                       CN0=signal.CN0, phi=signal.phi, nADC=signal.nADC,
                       include_carrier=signal.include_carrier,
                       include_adc=signal.include_adc,
                       include_thermal_noise=signal.include_thermal_noise,
                       code_start_idx=signal.code_start_idx,
                       include_databits_I=signal.include_databits_I,
                       include_databits_Q=signal.include_databits_Q,
                       include_phase_noise=signal.include_phase_noise,
                       phase_noise_init_phase=0, name=signal.name,
                       new_thermal_noise=false, new_phase_noise=false,
                       isreplica=signal.isreplica, noexp=signal.noexp,
                       receiver_h_parms=signal.receiver_h_parms,
                       new_databits=false, start_t=signal.start_t)
    # Calculate code chipping rates with Doppler applied for all codes on
    # each channel and their initial code phases
    f_s = signal.f_s
    signal_type = signal.signal_type
    sample_num = signal.sample_num
    sig_freq = signal_type.sig_freq
    I_codes = signal_type.I_codes
    Q_codes = signal_type.Q_codes
    # Get adjusted code chipping rates and chipping rate rates for I channel
    # codes
    if signal_type.include_I
        for i in 1:signal_type.I_codes.code_num
            # I channel
            f_code = I_codes.chipping_rates[i]
            code_length = I_codes.code_lengths[i]
            f_code_d, f_code_dd = calc_doppler_code_rate(f_code, sig_freq,
                                                         f_d, fd_rate)
            signal.f_code_d_I[i] = f_code_d
            signal.f_code_dd_I[i] = f_code_dd
            signal.init_code_phases_I[i] = calcinitcodephase(code_length,
                                                             f_code_d,
                                                             f_code_dd, f_s,
                                                             code_start_idx)
        end
        if I_codes.databits
            if ~include_databits_I
                signal.signal_type.I_codes.include_codes[end] = false
            else
                signal.signal_type.I_codes.include_codes[end] = true
            end
            # Generate new databits for I channel if `new_databits` flag is
            # `true`
            if new_databits
                prn_keys = collect(keys(I_codes.codes[1]))

                if I_codes.similar_databits
                    # If one array of databits was used for all PRNs, then the
                    # databit arrays for all PRNs are references to the same
                    # one. Therefore, changing one, will change them all.
                    rand!(signal.signal_type.I_codes.codes[end][prn_keys[1]], 0:1)
                else
                    # If a dictionary of databits was given instead, then each
                    # databit array is different. They are then separately
                    # regenerated.
                    for key in prn_keys
                        rand!(signal.signal_type.I_codes.codes[end][key], 0:1)
                    end
                end
            end
        end
    end
    # Get adjusted code chipping rates and chipping rate rates for Q channel
    # codes
    if signal_type.include_Q
        for i in 1:signal_type.Q_codes.code_num
            # Q channel
            f_code = Q_codes.chipping_rates[i]
            code_length = Q_codes.code_lengths[i]
            f_code_d, f_code_dd = calc_doppler_code_rate(f_code, sig_freq,
                                                         f_d, fd_rate)
            signal.f_code_d_Q[i] = f_code_d
            signal.f_code_dd_Q[i] = f_code_dd
            signal.init_code_phases_Q[i] = calcinitcodephase(code_length,
                                                             f_code_d,
                                                             f_code_dd, f_s,
                                                             code_start_idx)
        end
        # Generate new databits for Q channel if `new_databits` flag is
        # `true`
        if Q_codes.databits
            if ~include_databits_Q
                signal.signal_type.Q_codes.include_codes[end] = false
            else
                signal.signal_type.Q_codes.include_codes[end] = true
            end
            if new_databits
                prn_keys = collect(keys(Q_codes.codes[1]))
                if Q_codes.similar_databits
                    # If one array of databits was used for all PRNs, then the
                    # databit arrays for all PRNs are references to the same
                    # one. Therefore, changing one, will change them all.
                    rand!(signal.signal_type.Q_codes.codes[end][prn_keys[1]], 0:1)
                else
                    for key in prn_keys
                        # If a dictionary of databits was given instead, then each
                        # databit array is different. They are then separately
                        # regenerated.
                        rand!(signal.signal_type.Q_codes.codes[end][key], 0:1)
                    end
                end
            end
        end
    end
    # Generate new thermal noise and phase noise vectors if `new_thermal_noise`
    # and/or `new_phase_noise` flags are set to `true`
    if new_thermal_noise
        randn!(signal.thermal_noise)
    end
    if new_phase_noise
        generate_phase_noise!(signal.phase_noise, signal.t_length,
                              signal.signal_type.sig_freq,
                              receiver_h_parms;
                              phi_init=phase_noise_init_phase)
    end
    signal.name = name
    signal.prn = prn
    signal.f_if = f_if
    signal.f_d = f_d
    signal.fd_rate = fd_rate
    signal.Tsys = Tsys
    signal.CN0 = CN0
    signal.phi = phi
    signal.nADC = nADC
    signal.code_start_idx = code_start_idx
    signal.include_carrier = include_carrier
    signal.include_adc = include_adc
    signal.include_thermal_noise = include_thermal_noise
    signal.include_databits_I = include_databits_I
    signal.include_databits_Q = include_databits_Q
    signal.include_phase_noise = include_phase_noise
    signal.isreplica = isreplica
    signal.noexp = noexp
    signal.start_t = start_t
    return signal
end

"""
    definereplica(signal_type::SignalType, f_s, t_length; prn=1,
                  f_if=0., f_d=0., fd_rate=0., phi=0.,
                  include_carrier=true, code_start_idx=1.,
                  include_databits_I=true, include_databits_Q=true,
                  name="custom",)

Define properties of locally generated generic replica signal
based off its type, PRN, etc.

This is used for processing `GNSSSignal` structs. It does not
contain noise souces in it. Do not use for GNSS signal simulation.
Use definesignal and definesignal! methods for simulating GNSS
signals with noise.

Required Arguments:

- `signal_type::SignalType`: user defined signal type
    * see `definecodetype` and `definesignaltype`
- `f_s`: signal sampling frequency in Hz
- `t_length`: signal length in seconds

Optional Arguments:

- `prn`: satellite PRN number `(default = 1)`
- `f_if`: IF frequency in Hz `(default = 0Hz)`
- `f_d`: signal Doppler frequency in Hz `(default = 0Hz)`
- `fd_rate`: signal Doppler rate in Hz/s `(default = 0Hz/s)`
- `phi`: initial carrier phase in rad `(default = 0rad)`
- `include_carrier`: flag for modulating codes onto carrier `(default = true)`
- `code_start_idx`: starting index in `ReplicaSignal.data` where all codes
                    start `(default = 1)`
- `include_databits_I`: flag for including I channel databits in simulated signal
                        `(default = true)`
- `include_databits_Q`: flag for including Q channel databits in simulated signal
                        `(default = true)`
- `name`: name of signal `(default = custom)`
- `start_t`: the start time of the signal in seconds `(default = 0 seconds)`

Returns:

- `ReplicaSignal` struct
"""
function definereplica(signal_type::SignalType, f_s, t_length; prn=1,
                       f_if=0., f_d=0., fd_rate=0., phi=0.,
                       include_carrier=true, code_start_idx=1.,
                       include_databits_I=true, include_databits_Q=true,
                       name="custom", start_t=0)
    return definesignal(signal_type, f_s, t_length; prn=prn,
                        f_if=f_if, f_d=f_d, fd_rate=fd_rate, phi=phi,
                        include_carrier=include_carrier, include_adc=false,
                        include_thermal_noise=false,
                        code_start_idx=code_start_idx,
                        include_databits_I=include_databits_I,
                        include_databits_Q=include_databits_Q,
                        include_phase_noise=false, name=name,
                        skip_noise_generation=true,
                        allocate_noise_vectors=false,
                        start_t=start_t)
end

"""
    definereplica!(signal::ReplicaSignal;
                   prn=signal.prn, f_if=signal.f_if, f_d=signal.f_d,
                   fd_rate=signal.fd_rate, phi=signal.phi,
                   include_carrier=signal.include_carrier,
                   code_start_idx=signal.code_start_idx,
                   include_databits_I=signal.include_databits_I,
                   include_databits_Q=signal.include_databits_Q,
                   name=signal.name)

Redefine properties of locally generated generic replica signal
based off its type, PRN, etc.

This is used for processing `GNSSSignal` structs. It does not
contain noise souces in it. Do not use for GNSS signal simulation.
Use definesignal and definesignal! methods for simulating GNSS
signals with noise.

Required Arguments:

- `signal::ReplicaSignal`: the signal already defined by the user using
                            `definesignal`

Optional Arguments with Defaults Equal to `signal` Field Values:

- `prn`: satellite PRN number
- `f_if`: IF frequency in Hz
- `f_d`: signal Doppler frequency in Hz
- `fd_rate`: signal Doppler rate in Hz/s
- `phi`: initial carrier phase in rad
- `include_carrier`: flag for modulating codes onto carrier
- `code_start_idx`: starting index in `ReplicaSignal.data` where all codes
                    start
- `include_databits_I`: flag for including I channel databits in simulated signal
- `include_databits_Q`: flag for including Q channel databits in simulated signal
- `name`: name of signal
- `start_t`: the start time of the signal in seconds `(default = 0 seconds)`

Modifies and Returns:

- `signal::ReplicaSignal`
"""
function definereplica!(signal::ReplicaSignal;
                        prn=signal.prn, f_if=signal.f_if, f_d=signal.f_d,
                        fd_rate=signal.fd_rate, phi=signal.phi,
                        include_carrier=signal.include_carrier,
                        code_start_idx=signal.code_start_idx,
                        include_databits_I=signal.include_databits_I,
                        include_databits_Q=signal.include_databits_Q,
                        name=signal.name, start_t=signal.start_t)
return definesignal!(signal;
                     prn=prn, f_if=f_if, f_d=f_d,
                     fd_rate=fd_rate, phi=phi,
                     include_carrier=include_carrier,
                     include_adc=false, include_thermal_noise=false,
                     code_start_idx=code_start_idx,
                     include_databits_I=include_databits_I,
                     include_databits_Q=include_databits_Q,
                     include_phase_noise=false, name=name,
                     new_thermal_noise=false, new_phase_noise=false,
                     new_databits=false,
                     start_t=signal.start_t)
end

"""
    check_length(a, N)

Check that an array, `a`, has length `N` or 1.
If it has length 1, return `a`, if its length
is `N`, return `fill(a, N)`.

Required Arguments:

- `a`: array whose length to check
- `N`: the length to compare with `length(a)`

Returns:

- array that is either the same as `a` or `fill(a, N)`
"""
function check_length(a, N)
    @assert (length(a) == N) || (length(a) == 1)
    if length(a) > 1
        @assert length(a) == N
        return a
    else
        return fill(a, N)
    end
end

"""
    definesignal(prn::Vector, signal_type::SignalType, f_s, t_length;
                 f_if=0., f_d=0., fd_rate=0., Tsys=535.,
                 CN0=45., phi=0., nADC=4, include_carrier=true,
                 include_adc=true, include_thermal_noise=true,
                 code_start_idx=1., include_databits_I=true,
                 include_databits_Q=true, include_phase_noise=true,
                 phase_noise_scaler=1/10, name="custom",
                 receiver_h_parms=h_parms_tcxo[1], start_t=0)
"""
function definesignal(prn::Vector{Int}, signal_type, f_s, t_length;
                      f_if=0., f_d=0., fd_rate=0., Tsys=535.,
                      CN0=45., phi=0., nADC=4, include_carrier=true,
                      include_adc=true, include_thermal_noise=true,
                      code_start_idx=1., include_databits_I=true,
                      include_databits_Q=true, include_phase_noise=true,
                      phase_noise_scaler=1/10, name="custom",
                      receiver_h_parms=h_parms_tcxo[1], start_t=0)
    # Check that parameters are either length of 1 or N. If they are length 1,
    # then they are globally asigned to all signals defined.
    N = length(prn)
    f_d = check_length(f_d, N)
    fd_rate = check_length(fd_rate, N)
    CN0 = check_length(CN0, N)
    phi = check_length(phi, N)
    code_start_idx = check_length(code_start_idx, N)
    include_databits_I = check_length(include_databits_I, N)
    include_databits_Q = check_length(include_databits_Q, N)
    # Check that `signal_type` is either an array of signal types or singular
    if isa(signal_type, Array{eltype(signal_type)})
       @assert length(signal_type) == N
    elseif isa(signal_type, SignalType)
        signal_type = fill(signal_type, N)
    else
        error("Invalid `signal_type` specified.")
    end
    # Define vector of `length(prn)` `ReplicaSignal` structs with 0 `t_length`
    replica_signals = Array{ReplicaSignal}(undef, N)
    B_Is = []
    B_Qs = []

    for i in 1:N
        push!(B_Is, signal_type[i].B_I)
        push!(B_Qs, signal_type[i].B_Q)
        replica_signals[i] = definesignal(signal_type[i], f_s, 0.; prn=prn[i],
                                          f_if=f_if, f_d=f_d[i],
                                          fd_rate=fd_rate[i], Tsys=Tsys[i],
                                          CN0=CN0[i], phi=phi[i], nADC=nADC,
                                          code_start_idx=code_start_idx[i],
                                          include_databits_I=include_databits_I[i],
                                          include_databits_Q=include_databits_Q[i],
                                          include_thermal_noise=false,
                                          include_phase_noise=false,
                                          include_carrier=include_carrier,
                                          include_adc=false,
                                          skip_noise_generation=true,
                                          allocate_noise_vectors=false)
    end
    B_I = maximuma(xor.(1, ismissing.(B_Is)))
    B_Q = maximuma(xor.(1, ismissing.(B_Qs)))
    B = max(B_I, B_Q)
    sample_num = floor(Int, f_s*t_length)
    data = Array{Complex{Float64}}(undef, sample_num)
    if include_thermal_noise
        thermal_noise = randn(Complex{Float64}, sample_num)
    else
        thermal_noise = Array{Complex{Float64}}(undef, sample_num)
    end
    if include_phase_noise
        phase_noise = generate_phase_noise(t_length, f_s,
                                           signal_type[1].sig_freq,
                                           receiver_h_parms)
    else
        phase_noise = Array{Float64}(undef, sample_num)
    end
    return ReplicaSignals(name, replica_signals, data, t_length, f_s, f_if,
                          B, Tsys, sample_num, nADC, include_carrier,
                          include_adc,  include_thermal_noise,
                          include_phase_noise, thermal_noise, phaser_noise,
                          receiver_h_parms, start_t)
end

"""
    definesignal!(signal::ReplicaSignals, prn::Vector{Int}, signal_type;
                  f_if=signal.f_if, f_d=0., fd_rate=0., Tsys=535.,
                  CN0=45., phi=0., nADC=4,
                  include_carrier=signal.include_carrier,
                  include_adc=signal.include_adc,
                  include_thermal_noise=signal.include_thermal_noise,
                  code_start_idx=1., include_databits_I=true,
                  include_databits_Q=true,
                  include_phase_noise=signal.include_phase_noise,
                  name="custom", new_thermal_noise=false,
                  new_phase_noise=false,
                  receiver_h_parms=signal.receiver_h_parms,
                  start_t=signal.start_t)
"""
function definesignal!(signal::ReplicaSignals, prn::Vector{Int}, signal_type;
                       f_if=signal.f_if, f_d=0., fd_rate=0., Tsys=535.,
                       CN0=45., phi=0., nADC=4,
                       include_carrier=signal.include_carrier,
                       include_adc=signal.include_adc,
                       include_thermal_noise=signal.include_thermal_noise,
                       code_start_idx=1., include_databits_I=true,
                       include_databits_Q=true,
                       include_phase_noise=signal.include_phase_noise,
                       name="custom", new_thermal_noise=false,
                       new_phase_noise=false,
                       receiver_h_parms=signal.receiver_h_parms,
                       start_t=signal.start_t)
    # Check that parameters are either length of 1 or N. If they are length 1,
    # then they are globally asigned to all signals defined.
    f_s = signal.f_s
    t_length = signal.t_length
    include_carrier = signal.include_carrier
    N = length(prn)
    f_d = check_length(f_d, N)
    fd_rate = check_length(fd_rate, N)
    CN0 = check_length(CN0, N)
    phi = check_length(phi, N)
    code_start_idx = check_length(code_start_idx, N)
    include_databits_I = check_length(include_databits_I, N)
    include_databits_Q = check_length(include_databits_Q, N)
    signal.include_thermal_noise = include_thermal_noise
    signal.include_phase_noise = include_phase_noise
    signal.Tsys = Tsys
    signal.start_t = start_t
    # Check that `signal_type` is either an array of signal types or singular
    if isa(signal_type, Array{eltype(signal_type)})
       @assert length(signal_type) == N
    elseif isa(signal_type, SignalType)
        signal_type = fill(signal_type, N)
    else
        error("Invalid `signal_type` specified.")
    end
    # Define vector of `length(prn)` `ReplicaSignal` structs with 0 `t_length`
    signal.replica_signals = Array{ReplicaSignal}(undef, N)
    B_Is = []
    B_Qs = []
    for i in 1:N
        push!(B_Is, signal_type[i].B_I)
        push!(B_Qs, signal_type[i].B_Q)
        signal.replica_signals[i] = definesignal(signal_type[i], f_s, 0.;
                                                 prn=prn[i],
                                                 f_if=f_if, f_d=f_d[i],
                                                 fd_rate=fd_rate[i], Tsys=Tsys[i],
                                                 CN0=CN0[i], phi=phi[i], nADC=nADC,
                                                 code_start_idx=code_start_idx[i],
                                                 include_databits_I=include_databits_I[i],
                                                 include_databits_Q=include_databits_Q[i],
                                                 include_thermal_noise=false,
                                                 include_phase_noise=false,
                                                 include_carrier=include_carrier,
                                                 include_adc=false,
                                                 skip_noise_generation=true,
                                                 allocate_noise_vectors=false)
    end
    B_I = maximuma(xor.(1, ismissing.(B_Is)))
    B_Q = maximuma(xor.(1, ismissing.(B_Qs)))
    B = max(B_I, B_Q)
    signal.B = B
    if new_thermal_noise
        randn!(signal.thermal_noise)
    end
    if new_phase_noise
        generate_phase_noise!(signal.phase_noise, t_length,
                              signal.replica_signals[1].signal_type.sig_freq,
                              receiver_h_parms)
    end
    return signal
end


text:
# This file was generated, do not modify it. # hide
Xs = source()
ys = source()

DecisionTreeRegressor = @load DecisionTreeRegressor pkg=DecisionTree
atom = DecisionTreeRegressor()

machines = (machine(atom, Xs, ys) for i in 1:100)

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials