text:
##################################################
##### COVID Patient Re-Distribution -- Texas #####
##################################################
printstyled("COVID Patient Re-Distribution -- Texas\n", bold=true)

#############
### Setup ###
#############
println("Setup")

## Includes ##
println("Includes")
using Dates
using Distributions
using JuMP
using Serialization
ENV["COLUMNS"] = 1000;

projectbasepath = "../../";
include(normpath(projectbasepath, "src/models/PatientAllocation.jl"));
include(normpath(projectbasepath, "src/util/PatientAllocationFigures.jl"));

include("tx_data.jl");

## Configuation ##
println("Config")
results_basepath = abspath(joinpath(projectbasepath, "results/"));
paperfigures_basepath = abspath(joinpath(projectbasepath, "figures/"));

if !isdir(results_basepath) mkpath(results_basepath) end;
if !isdir(paperfigures_basepath) mkpath(paperfigures_basepath) end;

shared_config = (
    rundate = today(),
    region = "Texas",
    region_abbrev = "tx",
    alloc_level = "tsa",
    results_basepath = results_basepath,
    paperfigures_basepath = paperfigures_basepath,
);

saves_path = "$(shared_config.results_basepath)/$(shared_config.region_abbrev)/$(shared_config.rundate)/saves/"
if !isdir(saves_path) mkpath(saves_path) end;

start_date = Date(2020, 6, 15);
end_date   = Date(2020, 8, 15);
date_range = collect(start_date : Day(1) : end_date);
T = length(date_range);

locations_limit = nothing;
focus_locations = [
    "Austin",
    "Dallas/Ft. Worth",
    "El Paso",
    "Galveston",
    "Houston",
    "San Antonio",
];

pct_beds_available = (
    icu = 0.5,
    ward = 0.35,
    allpat = 0.4,
);
los_dist = (
    icu = Weibull(1.58, 13.32),
    ward = Weibull(1.38, 12.88),
    allpat = Weibull(1.38, 12.88),
);

## Helper Functions ##
println("Helper functions")

function save_results(experiment_name, _sent)
    serialize(joinpath(saves_path, "sent_$(experiment_name).jldata"), _sent)
    return
end;

function print_solve_metrics(_model)
    println("termination status:       ", termination_status(_model))
    println("solve time:               ", round(solve_time(_model), digits=3), "s")
    println("objective function value: ", round(objective_value(_model), digits=3))
    return
end;

function make_data_figures(_data)
    _config = merge(shared_config, (experiment = "data",))
    _sent = zeros(Int, _data.N, _data.N, _data.T)
    _results = PatientAllocationFigures.compute_results(_config, _data, _sent, use_rounding=true)

    PatientAllocationFigures.plot_active_total(_config, _data, _results, display=false)

    return
end;

function make_all_figures_base(experiment_name, _data, _sent)
    _config = merge(shared_config, (experiment = experiment_name,))

    _results = PatientAllocationFigures.compute_results(_config, _data, _sent, use_rounding=true)

    PatientAllocationFigures.plot_metrics(_config, _data, _results, display=false)

    PatientAllocationFigures.plot_estimates_total(_config, _data, display=false)

    PatientAllocationFigures.plot_overflow_distribution(_config, _data, _results, display=false)
    PatientAllocationFigures.plot_load_distribution(_config, _data, _results, display=false)
    PatientAllocationFigures.plot_maxload_distribution(_config, _data, _results, display=false)

    PatientAllocationFigures.plot_sent_total(_config, _data, _results, display=false)
    PatientAllocationFigures.plot_transfers(_config, _data, _results, display=false)
    PatientAllocationFigures.plot_transfers(_config, _data, _results, display=false, subset=_data.focus_locations)

    PatientAllocationFigures.plot_load(_config, _data, _results, display=false)
    PatientAllocationFigures.plot_load(_config, _data, _results, display=false, subset=_data.focus_locations)

    PatientAllocationFigures.plot_active(_config, _data, _results, display=false)
    PatientAllocationFigures.plot_active(_config, _data, _results, display=false, subset=_data.focus_locations)

    PatientAllocationFigures.plot_figures_list(_config)

    return _results
end;

function make_all_figures_block(experiment_name, _data, _sent)
    _config = merge(shared_config, (experiment = experiment_name,))

    _results = PatientAllocationFigures.compute_results_block(_config, _data, _sent, use_rounding=true);

    PatientAllocationFigures.plot_active_total(_config, _data.icu, _results.icu, bedtype=:icu, display=false);
    PatientAllocationFigures.plot_active_total(_config, _data.ward, _results.ward, bedtype=:ward, display=false);

    PatientAllocationFigures.plot_metrics(_config, _data.icu, _results.icu, bedtype=:icu, display=false);
    PatientAllocationFigures.plot_metrics(_config, _data.ward, _results.ward, bedtype=:ward, display=false);

    PatientAllocationFigures.plot_estimates_total(_config, _data.icu, bedtype=:icu, display=false);
    PatientAllocationFigures.plot_estimates_total(_config, _data.ward, bedtype=:ward, display=false);

    PatientAllocationFigures.plot_overflow_distribution(_config, _data.icu, _results.icu, bedtype=:icu, display=false);
    PatientAllocationFigures.plot_load_distribution(_config, _data.icu, _results.icu, bedtype=:icu, display=false);
    PatientAllocationFigures.plot_maxload_distribution(_config, _data.icu, _results.icu, bedtype=:icu, display=false);

    PatientAllocationFigures.plot_overflow_distribution(_config, _data.ward, _results.ward, bedtype=:ward, display=false);
    PatientAllocationFigures.plot_load_distribution(_config, _data.ward, _results.ward, bedtype=:ward, display=false);
    PatientAllocationFigures.plot_maxload_distribution(_config, _data.ward, _results.ward, bedtype=:ward, display=false);

    PatientAllocationFigures.plot_sent_total(_config, _data.icu, _results.icu, bedtype=:icu, display=false);
    PatientAllocationFigures.plot_transfers(_config, _data.icu, _results.icu, bedtype=:icu, display=false);
    PatientAllocationFigures.plot_transfers(_config, _data.icu, _results.icu, bedtype=:icu, display=false, subset=_data.focus_locations);

    PatientAllocationFigures.plot_sent_total(_config, _data.ward, _results.ward, bedtype=:ward, display=false);
    PatientAllocationFigures.plot_transfers(_config, _data.ward, _results.ward, bedtype=:ward, display=false);
    PatientAllocationFigures.plot_transfers(_config, _data.ward, _results.ward, bedtype=:ward, display=false, subset=_data.focus_locations);

    PatientAllocationFigures.plot_load(_config, _data.icu, _results.icu, bedtype=:icu, display=false);
    PatientAllocationFigures.plot_load(_config, _data.icu, _results.icu, bedtype=:icu, display=false, subset=_data.focus_locations);

    PatientAllocationFigures.plot_load(_config, _data.ward, _results.ward, bedtype=:ward, display=false);
    PatientAllocationFigures.plot_load(_config, _data.ward, _results.ward, bedtype=:ward, display=false, subset=_data.focus_locations);

    PatientAllocationFigures.plot_active(_config, _data.icu, _results.icu, add_title=true, bedtype=:icu, display=false);
    PatientAllocationFigures.plot_active(_config, _data.icu, _results.icu, add_title=true, bedtype=:icu, display=false, subset=_data.focus_locations);

    PatientAllocationFigures.plot_active(_config, _data.ward, _results.ward, add_title=true, bedtype=:ward, display=false);
    PatientAllocationFigures.plot_active(_config, _data.ward, _results.ward, add_title=true, bedtype=:ward, display=false, subset=_data.focus_locations);

    return _results
end;

function make_all_figures_robust(experiment_name, _data, _sent, _sent_base)
    _config = merge(shared_config, (experiment = experiment_name,))

    _results = PatientAllocationFigures.compute_results(_config, _data, _sent, use_rounding=true)
    _results = merge(_results, (sent_robust = _sent, sent_nonrobust = _sent_base));

    PatientAllocationFigures.plot_active_robust(_config, _data, _results, display=false)
    PatientAllocationFigures.plot_active_robust(_config, _data, _results, display=false, subset=_data.focus_locations)

    PatientAllocationFigures.plot_active_samples_notransfers(_config, _data, _results, display=false)
    PatientAllocationFigures.plot_active_samples_notransfers(_config, _data, _results, display=false, subset=_data.focus_locations)

    PatientAllocationFigures.plot_robust_overflow_distribution(_config, _data, _results, debug=false, display=false)

    PatientAllocationFigures.plot_figures_list(_config)

    return _results
end;

############
### Data ###
############
println("Data")

println("Load data")
data = TexasData.load_data_tx(
    date_range,
    los_dist, pct_beds_available,
    focus_locations, locations_limit,
    use_rounding=true,
);

println("Save data")
serialize(joinpath(saves_path, "data.jldata"), data);
make_data_figures(data.icu);

##################
### Run models ###
##################
println("Run models")

## No Transfer Model ##
println("No Transfer Model")
sent_null = zeros(Float64, data.N, data.N, data.T);
results_null = make_all_figures_base("no_transfers", data.icu, sent_null);
println("===========================================================\n")

## Base Model ##
println("Base Model")
model_base = PatientAllocation.patient_allocation(
    data.icu.beds,
    data.icu.initial,
    data.icu.discharged,
    data.icu.admitted,
    data.icu.adj,
    los=data.icu.los_dist,
    verbose=true,
)
sent_base = value.(model_base[:sent])
save_results("base", sent_base)
print_solve_metrics(model_base)
results_base = make_all_figures_base("base", data.icu, sent_base);
println("===========================================================\n")

## Operational Model ##
println("Operational Model")
model_operational = PatientAllocation.patient_allocation(
    data.icu.beds,
    data.icu.initial,
    data.icu.discharged,
    data.icu.admitted,
    data.icu.adj,
    los=data.icu.los_dist,

    smoothness_penalty = 0.01,
    sent_penalty = 0.01,
    no_artificial_overflow = true,
    no_worse_overflow = true,
    capacity_cushion = 0.05,

    verbose=true,
)
sent_operational = value.(model_operational[:sent])
save_results("operational", sent_operational)
print_solve_metrics(model_operational)
results_operational = make_all_figures_base("operational", data.icu, sent_operational);
println("===========================================================\n")

## Robust Base Model ##
println("Robust Base Model")
model_robust = PatientAllocation.patient_allocation_robust(
    data.icu.beds,
    data.icu.initial,
    data.icu.discharged,
    data.icu.admitted,
    data.icu.admitted_uncertainty,
    data.icu.adj,
    los=data.icu.los_dist,

    Γ = 7,

    verbose=true,
)
sent_robust = value.(model_robust[:sent])
save_results("robust_base", sent_robust)
print_solve_metrics(model_robust)
make_all_figures_base("robust_base", data.icu, sent_robust);
results_robust = make_all_figures_robust("robust_base", data.icu, sent_robust, sent_base);
println("===========================================================\n")

## Robust Operational Model ##
println("Robust Operational Model")
model_robust_operational = PatientAllocation.patient_allocation_robust(
    data.icu.beds,
    data.icu.initial,
    data.icu.discharged,
    data.icu.admitted,
    data.icu.admitted_uncertainty,
    data.icu.adj,
    los=data.icu.los_dist,

    Γ = 7,

    smoothness_penalty = 0.01,
    sent_penalty = 0.01,
    no_artificial_overflow = true,
    no_worse_overflow = true,
    capacity_cushion = 0.05,

    verbose=true,
)
sent_robust_operational = value.(model_robust_operational[:sent])
save_results("robust_operational", sent_robust_operational)
print_solve_metrics(model_robust_operational)
make_all_figures_base("robust_operational", data.icu, sent_robust_operational);
results_robust_operational = make_all_figures_robust("robust_operational", data.icu, sent_robust_operational, sent_operational);
println("===========================================================\n")

## Block Base Model ##
println("Block Base Model")
model_block = PatientAllocation.patient_block_allocation(
    data.carepaths.beds,
    data.carepaths.initial,
    data.carepaths.discharged,
    data.carepaths.admitted,
    data.carepaths.los_bygroup,
    data.carepaths.adj,
    data.carepaths.group_transfer_graph,
    data.carepaths.bedtype_bygroup,
    verbose=true,
)
sent_block = value.(model_block[:sent])
save_results("block_base", sent_block)
print_solve_metrics(model_block)
results_block = make_all_figures_block("block_base", data.carepaths, sent_block);
println("===========================================================\n")

## Block Operational Model ##
println("Block Operational Model")
model_block_operational = PatientAllocation.patient_block_allocation(
    data.carepaths.beds,
    data.carepaths.initial,
    data.carepaths.discharged,
    data.carepaths.admitted,
    data.carepaths.los_bygroup,
    data.carepaths.adj,
    data.carepaths.group_transfer_graph,
    data.carepaths.bedtype_bygroup,

    smoothness_penalty = 0.01,
    sent_penalty = 0.01,
    no_artificial_overflow = true,
    no_worse_overflow = true,
    capacity_cushion = 0.05,

    verbose=true,
)
sent_block_operational = value.(model_block_operational[:sent])
save_results("block_operational", sent_block_operational)
print_solve_metrics(model_block_operational)
results_block_operational = make_all_figures_block("block_operational", data.carepaths, sent_block_operational);
println("===========================================================\n")

## Compare Results ##
println("Compare Results")
PatientAllocationFigures.plot_metrics_compare(shared_config, data.icu, [
    "no_transfer" => results_null,
    "base" => results_base,
    "operational" => results_operational,
    "robust_base" => results_robust,
    "robust_operational" => results_robust_operational,
], display_table=true)

##################
#### Complete ####
##################
println("\nComplete")

max_stars_repo_name: flixpar/covid-resource-allocation

max_issues_repo_name: flixpar/covid-resource-allocation

max_forks_repo_name: flixpar/covid-resource-allocation


text:
### Single precision interfaces (only direct calls implemented)

## 1D

"""
    nufft1d1!(xj      :: Array{Float32},
               cj      :: Array{ComplexF32},
               iflag   :: Integer,
               eps     :: Float32,
               fk      :: Array{ComplexF32}
               [, opts :: nufft_opts]
             )

Single precision type-1 1D complex nonuniform FFT. Output stored in fk.
"""
function nufft1d1!(xj      :: Array{Float32},
                    cj      :: Array{ComplexF32},
                    iflag   :: Integer,
                    eps     :: Float32,
                    fk      :: Array{ComplexF32},
                    opts    :: nufft_opts = finufft_default_opts())
    nj = length(xj)
    @assert length(cj)==nj
    ms = length(fk)
    ret = ccall( (:finufftf1d1, libfinufft),
                 Cint,
                 (BIGINT,
                  Ref{Float32},
                  Ref{ComplexF32},
                  Cint,
                  Float32,
                  BIGINT,
                  Ref{ComplexF32},
                  Ref{nufft_opts}),
                 nj, xj, cj, iflag, eps, ms, fk, opts
                 )
    check_ret(ret)
end

"""
    nufft1d2!(xj      :: Array{Float32},
               cj      :: Array{ComplexF32},
               iflag   :: Integer,
               eps     :: Float32,
               fk      :: Array{ComplexF32}
               [, opts :: nufft_opts]
             )

Single precision type-2 1D complex nonuniform FFT. Output stored in cj.
"""
function nufft1d2!(xj      :: Array{Float32},
                    cj      :: Array{ComplexF32},
                    iflag   :: Integer,
                    eps     :: Float32,
                    fk      :: Array{ComplexF32},
                    opts    :: nufft_opts = finufft_default_opts())
    nj = length(xj)
    @assert length(cj)==nj
    ms = length(fk)
    ret = ccall( (:finufftf1d2, libfinufft),
                 Cint,
                 (BIGINT,
                  Ref{Float32},
                  Ref{ComplexF32},
                  Cint,
                  Float32,
                  BIGINT,
                  Ref{ComplexF32},
                  Ref{nufft_opts}),
                 nj, xj, cj, iflag, eps, ms, fk, opts
                 )
    check_ret(ret)
end

"""
    nufft1d3!(xj      :: Array{Float32},
               cj      :: Array{ComplexF32},
               iflag   :: Integer,
               eps     :: Float32,
               sk      :: Array{Float32},
               fk      :: Array{ComplexF32},
               [, opts :: nufft_opts]
              )

Single precision type-3 1D complex nonuniform FFT. Output stored in fk.
"""
function nufft1d3!(xj      :: Array{Float32},
                    cj      :: Array{ComplexF32},
                    iflag   :: Integer,
                    eps     :: Float32,
                    sk      :: Array{Float32},
                    fk      :: Array{ComplexF32},
                    opts    :: nufft_opts = finufft_default_opts())
    nj = length(xj)
    @assert length(cj)==nj
    nk = length(sk)
    @assert length(fk)==nk
    ret = ccall( (:finufftf1d3, libfinufft),
                 Cint,
                 (BIGINT,
                  Ref{Float32},
                  Ref{ComplexF32},
                  Cint,
                  Float32,
                  BIGINT,
                  Ref{Float32},
                  Ref{ComplexF32},
                  Ref{nufft_opts}),
                 nj, xj, cj, iflag, eps, nk, sk, fk, opts
                 )
    check_ret(ret)
end

## 2D

"""
    nufft2d1!(xj      :: Array{Float32},
               yj      :: Array{Float32},
               cj      :: Array{ComplexF32},
               iflag   :: Integer,
               eps     :: Float32,
               fk      :: Array{ComplexF32}
               [, opts :: nufft_opts]
             )

Single precision type-1 2D complex nonuniform FFT. Output stored in fk.
"""
function nufft2d1!(xj      :: Array{Float32},
                    yj      :: Array{Float32},
                    cj      :: Array{ComplexF32},
                    iflag   :: Integer,
                    eps     :: Float32,
                    fk      :: Array{ComplexF32},
                    opts    :: nufft_opts = finufft_default_opts())
    nj = length(xj)
    @assert length(yj)==nj
    @assert length(cj)==nj
    ms, mt = size(fk)
    ret = ccall( (:finufftf2d1, libfinufft),
                 Cint,
                 (BIGINT,
                  Ref{Float32},
                  Ref{Float32},
                  Ref{ComplexF32},
                  Cint,
                  Float32,
                  BIGINT,
                  BIGINT,
                  Ref{ComplexF32},
                  Ref{nufft_opts}),
                 nj, xj, yj, cj, iflag, eps, ms, mt, fk, opts
                 )
    check_ret(ret)
end

"""
    nufft2d2!(xj      :: Array{Float32},
               yj      :: Array{Float32},
               cj      :: Array{ComplexF32},
               iflag   :: Integer,
               eps     :: Float32,
               fk      :: Array{ComplexF32}
               [, opts :: nufft_opts]
             )

Single precision type-2 2D complex nonuniform FFT. Output stored in cj.
"""
function nufft2d2!(xj      :: Array{Float32},
                    yj      :: Array{Float32},
                    cj      :: Array{ComplexF32},
                    iflag   :: Integer,
                    eps     :: Float32,
                    fk      :: Array{ComplexF32},
                    opts    :: nufft_opts = finufft_default_opts())
    nj = length(xj)
    @assert length(yj)==nj
    @assert length(cj)==nj
    ms, mt = size(fk)
    ret = ccall( (:finufftf2d2, libfinufft),
                 Cint,
                 (BIGINT,
                  Ref{Float32},
                  Ref{Float32},
                  Ref{ComplexF32},
                  Cint,
                  Float32,
                  BIGINT,
                  BIGINT,
                  Ref{ComplexF32},
                  Ref{nufft_opts}),
                 nj, xj, yj, cj, iflag, eps, ms, mt, fk, opts
                 )
    check_ret(ret)
end

"""
    nufft2d3!(xj      :: Array{Float32},
               yj      :: Array{Float32},
               cj      :: Array{ComplexF32},
               iflag   :: Integer,
               eps     :: Float32,
               sk      :: Array{Float32},
               tk      :: Array{Float32},
               fk      :: Array{ComplexF32}
               [, opts :: nufft_opts]
              )

Single precision type-3 2D complex nonuniform FFT. Output stored in fk.
    """
function nufft2d3!(xj      :: Array{Float32},
                    yj      :: Array{Float32},
                    cj      :: Array{ComplexF32},
                    iflag   :: Integer,
                    eps     :: Float32,
                    sk      :: Array{Float32},
                    tk      :: Array{Float32},
                    fk      :: Array{ComplexF32},
                    opts    :: nufft_opts = finufft_default_opts())
    nj = length(xj)
    @assert length(yj)==nj
    @assert length(cj)==nj
    nk = length(sk)
    @assert length(tk)==nk
    @assert length(fk)==nk
    ret = ccall( (:finufftf2d3, libfinufft),
                 Cint,
                 (BIGINT,
                  Ref{Float32},
                  Ref{Float32},
                  Ref{ComplexF32},
                  Cint,
                  Float32,
                  BIGINT,
                  Ref{Float32},
                  Ref{Float32},
                  Ref{ComplexF32},
                  Ref{nufft_opts}),
                 nj, xj, yj, cj, iflag, eps, nk, sk, tk, fk, opts
                 )
    check_ret(ret)
end

## 3D

"""
    nufft3d1!(xj      :: Array{Float32},
               yj      :: Array{Float32},
               zj      :: Array{Float32},
               cj      :: Array{ComplexF32},
               iflag   :: Integer,
               eps     :: Float32,
               fk      :: Array{ComplexF32}
               [, opts :: nufft_opts]
             )

Single precision type-1 3D complex nonuniform FFT. Output stored in fk.
"""
function nufft3d1!(xj      :: Array{Float32},
                    yj      :: Array{Float32},
                    zj      :: Array{Float32},
                    cj      :: Array{ComplexF32},
                    iflag   :: Integer,
                    eps     :: Float32,
                    fk      :: Array{ComplexF32},
                    opts    :: nufft_opts = finufft_default_opts())
    nj = length(xj)
    @assert length(yj)==nj
    @assert length(zj)==nj
    @assert length(cj)==nj
    ms, mt, mu = size(fk)
    ret = ccall( (:finufftf3d1, libfinufft),
                 Cint,
                 (BIGINT,
                  Ref{Float32},
                  Ref{Float32},
                  Ref{Float32},
                  Ref{ComplexF32},
                  Cint,
                  Float32,
                  BIGINT,
                  BIGINT,
                  BIGINT,
                  Ref{ComplexF32},
                  Ref{nufft_opts}),
                 nj, xj, yj, zj, cj, iflag, eps, ms, mt, mu, fk, opts
                 )
    check_ret(ret)
end

"""
    nufft3d2!(xj      :: Array{Float32},
               yj      :: Array{Float32},
               zj      :: Array{Float32},
               cj      :: Array{ComplexF32},
               iflag   :: Integer,
               eps     :: Float32,
               fk      :: Array{ComplexF32}
               [, opts :: nufft_opts]
             )

Single precision type-2 3D complex nonuniform FFT. Output stored in cj.
"""
function nufft3d2!(xj      :: Array{Float32},
                    yj      :: Array{Float32},
                    zj      :: Array{Float32},
                    cj      :: Array{ComplexF32},
                    iflag   :: Integer,
                    eps     :: Float32,
                    fk      :: Array{ComplexF32},
                    opts    :: nufft_opts = finufft_default_opts())
    nj = length(xj)
    @assert length(yj)==nj
    @assert length(zj)==nj
    @assert length(cj)==nj
    ms, mt, mu = size(fk)
    ret = ccall( (:finufftf3d2, libfinufft),
                 Cint,
                 (BIGINT,
                  Ref{Float32},
                  Ref{Float32},
                  Ref{Float32},
                  Ref{ComplexF32},
                  Cint,
                  Float32,
                  BIGINT,
                  BIGINT,
                  BIGINT,
                  Ref{ComplexF32},
                  Ref{nufft_opts}),
                 nj, xj, yj, zj, cj, iflag, eps, ms, mt, mu, fk, opts
                 )
    check_ret(ret)
end

"""
    nufft3d3!(xj      :: Array{Float32},
               yj      :: Array{Float32},
               zj      :: Array{Float32},
               cj      :: Array{ComplexF32},
               iflag   :: Integer,
               eps     :: Float32,
               sk      :: Array{Float32},
               tk      :: Array{Float32},
               uk      :: Array{Float32},
               fk      :: Array{ComplexF32}
               [, opts :: nufft_opts]
              )

Single precision type-3 3D complex nonuniform FFT. Output stored in fk.
"""
function nufft3d3!(xj      :: Array{Float32},
                    yj      :: Array{Float32},
                    zj      :: Array{Float32},
                    cj      :: Array{ComplexF32},
                    iflag   :: Integer,
                    eps     :: Float32,
                    sk      :: Array{Float32},
                    tk      :: Array{Float32},
                    uk      :: Array{Float32},
                    fk      :: Array{ComplexF32},
                    opts    :: nufft_opts = finufft_default_opts())
    nj = length(xj)
    @assert length(yj)==nj
    @assert length(zj)==nj
    @assert length(cj)==nj
    nk = length(sk)
    @assert length(tk)==nk
    @assert length(uk)==nk
    @assert length(fk)==nk
    ret = ccall( (:finufftf3d3, libfinufft),
                 Cint,
                 (BIGINT,
                  Ref{Float32},
                  Ref{Float32},
                  Ref{Float32},
                  Ref{ComplexF32},
                  Cint,
                  Float32,
                  BIGINT,
                  Ref{Float32},
                  Ref{Float32},
                  Ref{Float32},
                  Ref{ComplexF32},
                  Ref{nufft_opts}),
                 nj, xj, yj, zj, cj, iflag, eps, nk, sk, tk, uk, fk, opts
                 )
    check_ret(ret)
end

# keep nufftf calls for compability
const nufftf1d1 = nufft1d1
const nufftf1d2 = nufft1d2
const nufftf1d3 = nufft1d3
const nufftf2d1 = nufft2d1
const nufftf2d2 = nufft2d2
const nufftf2d3 = nufft2d3
const nufftf3d1 = nufft3d1
const nufftf3d2 = nufft3d2
const nufftf3d3 = nufft3d3

const nufftf1d1! = nufft1d1!
const nufftf1d2! = nufft1d2!
const nufftf1d3! = nufft1d3!
const nufftf2d1! = nufft2d1!
const nufftf2d2! = nufft2d2!
const nufftf2d3! = nufft2d3!
const nufftf3d1! = nufft3d1!
const nufftf3d2! = nufft3d2!
const nufftf3d3! = nufft3d3!

ext: jl


text:
using Gnuplot, Random
Random.seed!(123)
n = 20
x, y = rand(n), rand(n)
@gp x y "w l lc 'black' t 'path'"
save(term ="pngcairo size 600,400", output ="path.png")

max_stars_repo_name: lazarusA/gnuplot-examples

max_issues_repo_name: lazarusA/gnuplot-examples

max_forks_repo_name: lazarusA/gnuplot-examples


text:
using Winston

println("-----------------------------------")
println("  Welcome to Lambda Plot.jl")
println("-----------------------------------")

mt_int = ARGS[1]
mt_float = ARGS[2]
xi = ARGS[3]

Data = readcsv("Data/Gauge_$(mt_int)_$(mt_float)_$(xi).csv")

t = Data[:,1];
λ = Data[:,2];
# gauge = vcat(Data[:,3], Data[:,4], Data[:,5], Data[:,6]); # yt, g1, g2, g3
# G = Data[:,7];

# Gauge Plot
p = FramedPlot(
    title="Lambda Plots",
    xlabel="t",
    ylabel="\\lambda");
C = Curve(t, λ)
setattr(C, "label", "\\lambda")
lgnd = Legend(.9, .9, [C])
add(p, C, lgnd)
savefig(p, "Fig/Lambda_$(mt_int)_$(mt_float)_$(xi).svg", (1000, 600))
run(`inkscape -z Fig/Lambda_$(mt_int)_$(mt_float)_$(xi).svg -e Fig/Lambda_$(mt_int)_$(mt_float)_$(xi).png -d 300 --export-background=WHITE`)
run(`rm Fig/Lambda_$(mt_int)_$(mt_float)_$(xi).svg`)

max_stars_repo_name: Axect/RGE

max_issues_repo_name: Axect/RGE

max_forks_repo_name: Axect/RGE


text:
import Base: put!, wait, isready, take!, fetch
mutable struct BufferChannel <: AbstractChannel{Any}
    buffer::Dict{Int, GameHistory}
    cond_take::Condition    # waiting for data to become available
    BufferChannel() = new(Dict(), Condition())
end

"""
Save key and value to remote buffer
"""
function put!(buffer_channel::BufferChannel, k::Int, v::GameHistory)
    buffer_channel.buffer[k] = v
    notify(buffer_channel.cond_take)
    return buffer_channel
end

"""
Pop the key and its associated value from remote buffer
"""
function take!(buffer_channel::BufferChannel, k::Int)
    v=fetch(buffer_channel,k)
    delete!(buffer_channel.buffer, k)
    return v
end

isready(buffer_channel::BufferChannel) = length(buffer_channel.buffer) ≥ 1
isready(buffer_channel::BufferChannel, k::Int) = haskey(buffer_channel.buffer,k)

"""
Obtain Value from the remote buffer for a given Key
"""
function fetch(buffer_channel::BufferChannel, k::Int)
    wait(buffer_channel,k)
    return buffer_channel.buffer[k]
end

function wait(buffer_channel::BufferChannel, k::Int)
    while !isready(buffer_channel, k)
        wait(buffer_channel.cond_take)
    end
end

function wait(buffer_channel::BufferChannel)
    while !isready(buffer_channel)
        wait(buffer_channel.cond_take)
    end
end

"""
Returns the remote buffer as a Dict
"""
function fetch(buffer_channel::BufferChannel)
    wait(buffer_channel)
    return buffer_channel.buffer
end


text:
using D91PVInversion
using IterativeSolvers
using Random
using Printf
using JLD2
import PyPlot; const plt = PyPlot

# Define functional form of background PV anomaly
function bgq′fun(x, y, z, A, L)
    return A*exp(-(x^2 + y^2)/L^2)
end

# Define functional form of PV anomaly
function q′fun(x, y, z, A, L, H, p)
    πc = p.π0 - 0.5
    return A*exp(-(x^2 + y^2)/L^2)*exp(-(z-πc)^2/H^2)
end

# Set up inversion domain
params = Params(Float64)
domain = Domain(params, size = (96, 96, 16), x = (-3, 3), y = (-3, 3))

# Compute idealized background fields
B = 10.0
L = 0.5
inv = ColumnarVortex(domain = domain, params = params)
initialize!(inv, bgq′fun, B, L)
Random.seed!(4321)
solve!(inv; verbose = true)
ψ0 = inv.ψ
ϕ0 = inv.ϕ
q0 = inv.q

# Compute diagnostics
u0 = diagnose_umag(ψ0, domain)
v0 = diagnose_v(ψ0, domain)
ζ0 = diagnose_ζ(ψ0, domain)
θ0 = diagnose_θ(ϕ0, domain)
nx, ny, nz = size(domain)
x = domain.xc[1:nx]
z = domain.zc[1:nz]
v0mid = 0.5*(v0[1:nx,ny÷2,1:nz] + v0[1:nx,ny÷2+1,1:nz])
θ0mid = 0.5*(θ0[1:nx,ny÷2,1:nz] + θ0[1:nx,ny÷2+1,1:nz])
ζ0mid = 0.5*(ζ0[1:nx,ny÷2,1:nz] + ζ0[1:nx,ny÷2+1,1:nz])
q0mid = 0.5*(q0[1:nx,ny÷2,1:nz] + q0[1:nx,ny÷2+1,1:nz])

# Calculate vertical penetration
u0bot = u0[1:nx,1:ny,1]
u0max = 0.5*(u0[1:nx,1:ny,nz÷2] + u0[1:nx,1:ny,nz÷2+1])
vpen = maximum(u0bot)/maximum(u0max)
print("Vertical penetration: ", vpen, "\n")

# Plot solution
fig, axes = plt.subplots(
    figsize = (9.5, 3), nrows = 1, ncols = 4,
    sharey = true, sharex = true, constrained_layout = true
)
axes[1].set_xlim([-2, 2])
axes[1].invert_yaxis()
c = axes[1].contour(x, z, q0mid', colors = "black")
Δc = length(c.levels) > 1 ? c.levels[2] - c.levels[1] : NaN
axes[1].set_title(@sprintf("q₀\nmax %.1e\ncontour interval %.1e", maximum(q0mid), Δc))
fig.canvas.draw()
plt.clabel(c, inline = true, inline_spacing = 0, manual = true)
c = axes[2].contour(x, z, v0mid', colors = "black")
Δc = length(c.levels) > 1 ? c.levels[2] - c.levels[1] : NaN
axes[2].set_title(@sprintf("v₀\nmax %.1e\ncontour interval %.1e", maximum(v0mid), Δc))
fig.canvas.draw()
plt.clabel(c, inline = true, inline_spacing = 0, manual = true)
c = axes[3].contour(x, z, ζ0mid', colors = "black")
Δc = length(c.levels) > 1 ? c.levels[2] - c.levels[1] : NaN
axes[3].set_title(@sprintf("ζ₀\nmax %.1e\ncontour interval %.1e", maximum(ζ0mid), Δc))
fig.canvas.draw()
plt.clabel(c, inline = true, inline_spacing = 0, manual = true)
c = axes[4].contour(x, z, θ0mid', colors = "black")
Δc = length(c.levels) > 1 ? c.levels[2] - c.levels[1] : NaN
axes[4].set_title(@sprintf("θ₀\nmax %.1e\ncontour interval %.1e", maximum(θ0mid), Δc))
fig.canvas.draw()
plt.clabel(c, inline = true, inline_spacing = 0, manual = true)

# Create inversion problem
inv = LinearizedInversion(
    ψ0 = ψ0, ϕ0 = ϕ0,
    domain = domain, params = params
)

# Initialize problem
A = 1.0
L = 0.1
H = 0.1
initialize!(inv, q′fun, A, L, H, params)

# Solve
solve!(inv; Pl = IterativeSolvers.Identity(), verbose = true)

# Compute diagnostics
ψ′ = inv.ψ′
ϕ′ = inv.ϕ′
q′ = inv.q′
u′ = diagnose_umag(ψ′, domain)
v′ = diagnose_v(ψ′, domain)
ζ′ = diagnose_ζ(ψ′, domain)
θ′ = diagnose_θ(ϕ′, domain)
nx, ny, nz = size(domain)
x = domain.xc[1:nx]
z = domain.zc[1:nz]
v′mid = 0.5*(v′[1:nx,ny÷2,1:nz] + v′[1:nx,ny÷2+1,1:nz])
θ′mid = 0.5*(θ′[1:nx,ny÷2,1:nz] + θ′[1:nx,ny÷2+1,1:nz])
ζ′mid = 0.5*(ζ′[1:nx,ny÷2,1:nz] + ζ′[1:nx,ny÷2+1,1:nz])
q′mid = 0.5*(q′[1:nx,ny÷2,1:nz] + q′[1:nx,ny÷2+1,1:nz])

# Calculate vertical penetration
u′bot = u′[1:nx,1:ny,1]
u′max = 0.5*(u′[1:nx,1:ny,nz÷2] + u′[1:nx,1:ny,nz÷2+1])
vpen = maximum(u′bot)/maximum(u′max)
print("Vertical penetration: ", vpen, "\n")

# Plot solution
fig, axes = plt.subplots(
    figsize = (9.5, 3), nrows = 1, ncols = 4,
    sharey = true, sharex = true, constrained_layout = true
)
c = axes[1].contour(x, z, q′mid', colors = "black")
Δc = length(c.levels) > 1 ? c.levels[2] - c.levels[1] : NaN
axes[1].set_title(@sprintf("q′\nmax %.1e\ncontour interval %.1e", maximum(q′mid), Δc))
c = axes[2].contour(x, z, v′mid', colors = "black")
Δc = length(c.levels) > 1 ? c.levels[2] - c.levels[1] : NaN
axes[2].set_title(@sprintf("v′\nmax %.1e\ncontour interval %.1e", maximum(v′mid), Δc))
c = axes[3].contour(x, z, ζ′mid', colors = "black")
Δc = length(c.levels) > 1 ? c.levels[2] - c.levels[1] : NaN
axes[3].set_title(@sprintf("ζ′\nmax %.1e\ncontour interval %.1e", maximum(ζ′mid), Δc))
c = axes[4].contour(x, z, θ′mid', colors = "black")
Δc = length(c.levels) > 1 ? c.levels[2] - c.levels[1] : NaN
axes[4].set_title(@sprintf("θ′\nmax %.1e\ncontour interval %.1e", maximum(θ′mid), Δc))
axes[1].set_xlim([-1, 1])
axes[1].invert_yaxis()
plt.show()

# Save output
save_inversion_results(@sprintf("output/invert_linear_2_A%.1f_B%.1f.jld2", A, B), inv)


text:
names = [
    "test_sortedindices.jl",
    "test_polynomials.jl",

    "test_bases.jl",
    "test_states.jl",

    "test_operators.jl",
    "test_operators_dense.jl",
    "test_sparsematrix.jl",
    "test_operators_sparse.jl",
    "test_operators_lazytensor.jl",
    "test_operators_lazysum.jl",
    "test_operators_lazyproduct.jl",

    "test_fock.jl",
    "test_spin.jl",
    "test_particle.jl",
    "test_manybody.jl",
    "test_nlevel.jl",
    "test_subspace.jl",
    "test_state_definitions.jl",

    "test_transformations.jl",

    "test_metrics.jl",
    "test_embed.jl",

    "test_superoperators.jl",

    "test_pauli.jl",
    "test_spinors.jl",

    "test_abstractdata.jl",

    "test_printing.jl"
]

detected_tests = filter(
    name->startswith(name, "test_") && endswith(name, ".jl"),
    readdir("."))

unused_tests = setdiff(detected_tests, names)
if length(unused_tests) != 0
    error("The following tests are not used:\n", join(unused_tests, "\n"))
end

unavailable_tests = setdiff(names, detected_tests)
if length(unavailable_tests) != 0
    error("The following tests could not be found:\n", join(unavailable_tests, "\n"))
end

for name=names
    if startswith(name, "test_") && endswith(name, ".jl")
        include(name)
    end
end


text:
#!/usr/bin/env julia --project
using Test
using ClimateMachine
using ClimateMachine.GenericCallbacks
using ClimateMachine.ODESolvers
using ClimateMachine.Mesh.Filters
using ClimateMachine.VariableTemplates
using ClimateMachine.Mesh.Grids: polynomialorder
using ClimateMachine.Ocean
using ClimateMachine.Ocean.HydrostaticBoussinesq
using ClimateMachine.Ocean.ShallowWater
using ClimateMachine.Ocean.SplitExplicit: VerticalIntegralModel
using ClimateMachine.Ocean.SplitExplicit01
using ClimateMachine.Ocean.OceanProblems

using ClimateMachine.Mesh.Topologies
using ClimateMachine.Mesh.Grids
using ClimateMachine.BalanceLaws: vars_state, Prognostic, Auxiliary
using ClimateMachine.DGMethods
using ClimateMachine.DGMethods.NumericalFluxes
using ClimateMachine.MPIStateArrays
using ClimateMachine.VTK
using ClimateMachine.Checkpoint
using ClimateMachine.SystemSolvers

using MPI
using LinearAlgebra
using StaticArrays
using Logging, Printf, Dates

using CLIMAParameters
using CLIMAParameters.Planet: grav
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

struct SplitConfig{N, D3, D2, S, M, AT}
    name::N
    dg_3D::D3
    dg_2D::D2
    solver::S
    mpicomm::M
    ArrayType::AT
end

function run_split_explicit(
    config::SplitConfig,
    timespan;
    dt_fast = 300,
    dt_slow = 300,
    refDat = (),
    restart = 0,
    analytic_solution = false,
)
    Q_3D, Q_2D, t0 = init_states(config, Val(restart))

    tout, timeend = timespan

    # @show dt_fast = floor(Int, 1 / (2 * sqrt(gravity * H) / minΔx)) # / 4
    # @show dt_slow = floor(Int, minΔx / 15) # / 4

    nout = ceil(Int64, tout / dt_slow)
    dt_slow = tout / nout

    timeendlocal = timeend + t0

    lsrk_3D = LSRK54CarpenterKennedy(config.dg_3D, Q_3D, dt = dt_slow, t0 = t0)
    lsrk_2D = LSRK54CarpenterKennedy(config.dg_2D, Q_2D, dt = dt_fast, t0 = t0)

    odesolver = config.solver(lsrk_3D, lsrk_2D;)

    vtkstep = [restart, restart, restart + 1, restart + 1]
    cbvector = make_callbacks(
        abspath(joinpath(ClimateMachine.Settings.output_dir, config.name)),
        vtkstep,
        nout,
        config.mpicomm,
        odesolver,
        config.dg_3D,
        config.dg_3D.balance_law,
        Q_3D,
        config.dg_2D,
        config.dg_2D.balance_law,
        Q_2D,
        timeendlocal,
    )

    eng0 = norm(Q_3D)
    @info @sprintf """Starting
    norm(Q₀) = %.16e
    ArrayType = %s""" eng0 config.ArrayType

    # slow fast state tuple
    Qvec = (slow = Q_3D, fast = Q_2D)
    solve!(Q_3D, odesolver; timeend = timeendlocal, callbacks = cbvector)

    if analytic_solution
        Qe_3D = init_ode_state(config.dg_3D, timeendlocal, init_on_cpu = true)
        Qe_2D = init_ode_state(config.dg_2D, timeendlocal, init_on_cpu = true)

        error_3D = euclidean_distance(Q_3D, Qe_3D) / norm(Qe_3D)
        error_2D = euclidean_distance(Q_2D, Qe_2D) / norm(Qe_2D)

        println("3D error = ", error_3D)
        println("2D error = ", error_2D)
        @test isapprox(error_3D, FT(0.0); atol = 0.005)
        @test isapprox(error_2D, FT(0.0); atol = 0.005)
    end

    ## Check results against reference
    ClimateMachine.StateCheck.scprintref(cbvector[end])
    if length(refDat) > 0
        @test ClimateMachine.StateCheck.scdocheck(cbvector[end], refDat)
    end

    return nothing
end

function init_states(config, ::Val{0})
    Q_3D = init_ode_state(config.dg_3D, FT(0); init_on_cpu = true)
    Q_2D = config.dg_3D.modeldata.Q_2D

    return Q_3D, Q_2D, 0
end

function init_states(config, ::Val{restart}) where {restart}
    Q_3D, A_3D, t0 = read_checkpoint(
        abspath(joinpath(ClimateMachine.Settings.output_dir, config.name)),
        "baroclinic",
        config.ArrayType,
        config.mpicomm,
        restart,
    )
    Q_2D_restart, A_2D, _ = read_checkpoint(
        abspath(joinpath(ClimateMachine.Settings.output_dir, config.name)),
        "barotropic",
        config.ArrayType,
        config.mpicomm,
        restart,
    )

    direction = EveryDirection()
    A_3D = restart_auxiliary_state(
        config.dg_3D.balance_law,
        config.dg_3D.grid,
        A_3D,
        direction,
    )
    A_2D = restart_auxiliary_state(
        config.dg_2D.balance_law,
        config.dg_2D.grid,
        A_2D,
        direction,
    )

    config.dg_3D.state_auxiliary .= A_3D
    config.dg_2D.state_auxiliary .= A_2D

    Q_3D = restart_ode_state(config.dg_3D, Q_3D; init_on_cpu = true)
    Q_2D = config.dg_3D.modeldata.Q_2D
    Q_2D .= Q_2D_restart

    return Q_3D, Q_2D, t0
end

function make_callbacks(
    vtkpath,
    vtkstep,
    nout,
    mpicomm,
    odesolver,
    dg_slow,
    model_slow,
    Q_slow,
    dg_fast,
    model_fast,
    Q_fast,
    timeend,
)
    mkpath(vtkpath)
    mkpath(vtkpath * "/slow")
    mkpath(vtkpath * "/fast")

    A_slow = dg_slow.state_auxiliary
    A_fast = dg_fast.state_auxiliary

    function do_output(span, vtkstep, model, dg, Q, A)
        outprefix = @sprintf(
            "%s/%s/mpirank%04d_step%04d",
            vtkpath,
            span,
            MPI.Comm_rank(mpicomm),
            vtkstep
        )
        @info "doing VTK output" outprefix
        statenames = flattenednames(vars_state(model, Prognostic(), eltype(Q)))
        auxnames = flattenednames(vars_state(model, Auxiliary(), eltype(Q)))
        writevtk(outprefix, Q, dg, statenames, A, auxnames)
    end

    do_output("slow", vtkstep[1], model_slow, dg_slow, Q_slow, A_slow)
    vtkstep[1] += 1
    cbvtk_slow = GenericCallbacks.EveryXSimulationSteps(nout) do (init = false)
        do_output("slow", vtkstep[1], model_slow, dg_slow, Q_slow, A_slow)
        vtkstep[1] += 1
        nothing
    end

    do_output("fast", vtkstep[2], model_fast, dg_fast, Q_fast, A_fast)
    vtkstep[2] += 1
    cbvtk_fast = GenericCallbacks.EveryXSimulationSteps(nout) do (init = false)
        do_output("fast", vtkstep[2], model_fast, dg_fast, Q_fast, A_fast)
        vtkstep[2] += 1
        nothing
    end

    starttime = Ref(now())
    cbinfo = GenericCallbacks.EveryXWallTimeSeconds(60, mpicomm) do (s = false)
        if s
            starttime[] = now()
        else
            energy = norm(Q_slow)
            @info @sprintf(
                """Update
                simtime = %8.2f / %8.2f
                runtime = %s
                norm(Q) = %.16e""",
                ODESolvers.gettime(odesolver),
                timeend,
                Dates.format(
                    convert(Dates.DateTime, Dates.now() - starttime[]),
                    Dates.dateformat"HH:MM:SS",
                ),
                energy
            )
        end
    end

    cbcs_dg = ClimateMachine.StateCheck.sccreate(
        [
            (Q_slow, "3D state"),
            (A_slow, "3D aux"),
            (Q_fast, "2D state"),
            (A_fast, "2D aux"),
        ],
        nout;
        prec = 12,
    )

    cb_checkpoint = GenericCallbacks.EveryXSimulationSteps(nout) do
        write_checkpoint(
            Q_slow,
            A_slow,
            odesolver,
            vtkpath,
            "baroclinic",
            mpicomm,
            vtkstep[3],
        )

        write_checkpoint(
            Q_fast,
            A_fast,
            odesolver,
            vtkpath,
            "barotropic",
            mpicomm,
            vtkstep[4],
        )

        rm_checkpoint(vtkpath, "baroclinic", mpicomm, vtkstep[3] - 1)

        rm_checkpoint(vtkpath, "barotropic", mpicomm, vtkstep[4] - 1)

        vtkstep[3] += 1
        vtkstep[4] += 1
        nothing
    end

    return (cbvtk_slow, cbvtk_fast, cbinfo, cb_checkpoint, cbcs_dg)
end

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
import Images
using Color, FixedPointNumbers, Base.Test

if !isdefined(:workdir)
    const workdir = joinpath(tempdir(), "Images")
end

if !isdir(workdir)
    mkdir(workdir)
end

a = rand(2,2)
aa = convert(Array{Ufixed8}, a)
fn = joinpath(workdir, "2by2.png")
Images.imwrite(a, fn)
b = Images.imread(fn)
@test convert(Array, b) == aa
Images.imwrite(aa, fn)
b = Images.imread(fn)
@test convert(Array, b) == aa
aaimg = Images.grayim(aa)
open(fn, "w") do file
    writemime(file, MIME("image/png"), aaimg, minpixels=0)
end
b = Images.imread(fn)
@test b == aaimg
aa = convert(Array{Ufixed16}, a)
Images.imwrite(aa, fn)
b = Images.imread(fn)
@test convert(Array, b) == aa
aa = Ufixed12[0.6 0.2;
              1.4 0.8]
open(fn, "w") do file
    writemime(file, MIME("image/png"), Images.grayim(aa), minpixels=0)
end
b = Images.imread(fn)
@test Images.data(b) == Ufixed8[0.6 0.2;
                                1.0 0.8]

img = Images.colorim(rand(3,2,2))
img24 = convert(Images.Image{RGB24}, img)
Images.imwrite(img24, fn)
b = Images.imread(fn)
imgrgb8 = convert(Images.Image{RGB{Ufixed8}}, img)
@test Images.data(imgrgb8) == Images.data(b)

# test writemime's use of restrict
abig = Images.grayim(rand(Uint8, 1024, 1023))
fn = joinpath(workdir, "big.png")
open(fn, "w") do file
    writemime(file, MIME("image/png"), abig, maxpixels=10^6)
end
b = Images.imread(fn)
@test Images.data(b) == convert(Array{Ufixed8,2}, Images.data(Images.restrict(abig, (1,2))))

# More writemime tests
a = Images.colorim(rand(Uint8, 3, 2, 2))
fn = joinpath(workdir, "2by2.png")
open(fn, "w") do file
    writemime(file, MIME("image/png"), a, minpixels=0)
end
b = Images.imread(fn)
@test Images.data(b) == Images.data(a)

abig = Images.colorim(rand(Uint8, 3, 1021, 1026))
fn = joinpath(workdir, "big.png")
open(fn, "w") do file
    writemime(file, MIME("image/png"), abig, maxpixels=10^6)
end
b = Images.imread(fn)
@test Images.data(b) == convert(Array{RGB{Ufixed8},2}, Images.data(Images.restrict(abig, (1,2))))

using Color
datafloat = reshape(linspace(0.5, 1.5, 6), 2, 3)
dataint = round(Uint8, 254*(datafloat .- 0.5) .+ 1)  # ranges from 1 to 255
# build our colormap
b = RGB(0,0,1)
w = RGB(1,1,1)
r = RGB(1,0,0)
cmaprgb = Array(RGB{Float64}, 255)
f = linspace(0,1,128)
cmaprgb[1:128] = [(1-x)*b + x*w for x in f]
cmaprgb[129:end] = [(1-x)*w + x*r for x in f[2:end]]
img = Images.ImageCmap(dataint, cmaprgb)
Images.imwrite(img,joinpath(workdir,"cmap.jpg"))

c = reinterpret(Images.BGRA{Ufixed8}, [0xf0884422]'')
fn = joinpath(workdir, "alpha.png")
Images.imwrite(c, fn)
C = Images.imread(fn)
# @test C[1] == c[1]  # disabled because Travis has a weird, old copy of ImageMagick for which this fails (see #261)
Images.imwrite(reinterpret(ARGB32, [0xf0884422]''), fn)
D = Images.imread(fn)
# @test D[1] == c[1]


text:
using HyperdimensionalComputing
using Documenter

DocMeta.setdocmeta!(HyperdimensionalComputing, :DocTestSetup, :(using HyperdimensionalComputing); recursive=true)

makedocs(;
    modules=[HyperdimensionalComputing],
    authors="Steff Taelman, Dimi Boeckaerts",
    repo="https://github.com/dimiboeckaerts/HyperdimensionalComputing.jl/blob/{commit}{path}#{line}",
    sitename="HyperdimensionalComputing.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://dimiboeckaerts.github.io/HyperdimensionalComputing.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/dimiboeckaerts/HyperdimensionalComputing.jl",
)


text:
################################################################################
#
#  nmod_poly.jl : Flint nmod_poly (polynomials over Z/nZ, small modulus)
#
################################################################################

export NmodPolyRing, nmod_poly, parent, base_ring, elem_type, length, zero,
       one, gen, isgen, iszero, var, deepcopy, show, truncate, mullow, reverse,
       shift_left, shift_right, divexact, divrem, rem, gcd, xgcd, resultant,
       evaluate, derivative, compose, interpolate, inflate, deflate, lift,
       isirreducible, issquarefree, factor, factor_squarefree,
       factor_distinct_deg, factor_shape, setcoeff!, canonical_unit,
       add!, sub!, mul!, call, PolynomialRing, check_parent, gcdx, mod,
       invmod, gcdinv, mulmod, powmod, zero!, one!

################################################################################
#
#  Type and parent object methods
#
################################################################################

parent(a::nmod_poly) = a.parent

base_ring(R::NmodPolyRing) = R.base_ring

base_ring(a::nmod_poly) = base_ring(parent(a))

parent_type(::Type{nmod_poly}) = NmodPolyRing

elem_type(::Type{nmod_poly}) = nmod_poly

elem_type(::Type{NmodPolyRing}) = nmod_poly

function check_parent(x::nmod_poly, y::nmod_poly)
  parent(x) != parent(y) && error("Parents must coincide")
  nothing
end

################################################################################
#
#   Basic helper
#
################################################################################

function lead_isunit(a::nmod_poly)
  d = degree(a)
  u = ccall((:nmod_poly_get_coeff_ui, :libflint), UInt, (Ref{nmod_poly}, Int), a, d)
  n = ccall((:n_gcd, :libflint), UInt, (UInt, UInt), u, modulus(a))
  return n==1
end

function Base.hash(a::nmod_poly, h::UInt)
   b = 0x53dd43cd511044d1%UInt
   for i in 0:length(a) - 1
      u = ccall((:nmod_poly_get_coeff_ui, :libflint), UInt, (Ref{nmod_poly}, Int), a, i)
      b = xor(b, xor(hash(u, h), h))
      b = (b << 1) | (b >> (sizeof(Int)*8 - 1))
   end
   return b
end

################################################################################
#
#  Basic manipulation
#
################################################################################

length(x::nmod_poly) = ccall((:nmod_poly_length, :libflint), Int,
                               (Ref{nmod_poly}, ), x)

degree(x::nmod_poly) = ccall((:nmod_poly_degree, :libflint), Int,
                               (Ref{nmod_poly}, ), x)

function coeff(x::nmod_poly, n::Int)
  n < 0 && throw(DomainError())
  return base_ring(x)(ccall((:nmod_poly_get_coeff_ui, :libflint), UInt,
          (Ref{nmod_poly}, Int), x, n))
end

function coeff_raw(x::nmod_poly, n::Int)
  return ccall((:nmod_poly_get_coeff_ui, :libflint), UInt,
                (Ref{nmod_poly}, Int), x, n)
end

zero(R::NmodPolyRing) = R(UInt(0))

one(R::NmodPolyRing) = R(UInt(1))

gen(R::NmodPolyRing) = R([zero(base_ring(R)), one(base_ring(R))])

isgen(a::nmod_poly) = (degree(a) == 1 &&
                              iszero(coeff(a,0)) && isone(coeff(a,1)))

iszero(a::nmod_poly) = Bool(ccall((:nmod_poly_is_zero, :libflint), Int32,
                              (Ref{nmod_poly}, ), a))

modulus(a::nmod_poly) = a.parent.n

modulus(R::NmodPolyRing) = R.n

var(R::NmodPolyRing) = R.S

function deepcopy_internal(a::nmod_poly, dict::ObjectIdDict)
  z = nmod_poly(modulus(a), a)
  z.parent = a.parent
  return z
end

################################################################################
#
#  AbstractString I/O
#
################################################################################

function show(io::IO, x::nmod_poly)
  if length(x) == 0
    print(io, "0")
  else
    cstr = ccall((:nmod_poly_get_str_pretty, :libflint), Ref{UInt8},
            (Ref{nmod_poly}, Ref{UInt8}), x, string(var(parent(x))))
    print(io, unsafe_string(cstr))
    ccall((:flint_free, :libflint), Void, (Ref{UInt8}, ), cstr)
  end
end

function show(io::IO, R::NmodPolyRing)
  print(io, "Univariate Polynomial Ring in ")
  print(io, string(var(R)))
  print(io, " over ")
  print(io, base_ring(R))
end

show_minus_one(::Type{nmod_poly}) = true

################################################################################
#
#  Canonicalization
#
################################################################################

canonical_unit(a::nmod_poly) = canonical_unit(lead(a))

################################################################################
#
#  Unary operations
#
################################################################################

function -(x::nmod_poly)
  z = parent(x)()
  ccall((:nmod_poly_neg, :libflint), Void,
          (Ref{nmod_poly}, Ref{nmod_poly}), z, x)
  return z
end

################################################################################
#
#   Binary operations
#
################################################################################

function +(x::nmod_poly, y::nmod_poly)
  check_parent(x,y)
  z = parent(x)()
  ccall((:nmod_poly_add, :libflint), Void,
          (Ref{nmod_poly}, Ref{nmod_poly},  Ref{nmod_poly}), z, x, y)
  return z
end

function -(x::nmod_poly, y::nmod_poly)
  check_parent(x,y)
  z = parent(x)()
  ccall((:nmod_poly_sub, :libflint), Void,
          (Ref{nmod_poly}, Ref{nmod_poly},  Ref{nmod_poly}), z, x, y)
  return z
end

function *(x::nmod_poly, y::nmod_poly)
  check_parent(x,y)
  z = parent(x)()
  ccall((:nmod_poly_mul, :libflint), Void,
          (Ref{nmod_poly}, Ref{nmod_poly},  Ref{nmod_poly}), z, x, y)
  return z
end

###############################################################################
#
#  Ad hoc binary operations
#
###############################################################################

function *(x::nmod_poly, y::UInt)
  z = parent(x)()
  ccall((:nmod_poly_scalar_mul_nmod, :libflint), Void,
          (Ref{nmod_poly}, Ref{nmod_poly}, UInt), z, x, y)
  return z
end

*(x::UInt, y::nmod_poly) = y*x

function *(x::nmod_poly, y::fmpz)
  z = parent(x)()
  t = fmpz()
  tt = UInt(0)
  ccall((:fmpz_mod_ui, :libflint), UInt,
                (Ref{fmpz}, Ref{fmpz}, UInt), t, y, parent(x).n)
  tt = ccall((:fmpz_get_ui, :libflint), UInt, (Ref{fmpz}, ), t)
  return x*tt
end

*(x::fmpz, y::nmod_poly) = y*x

*(x::nmod_poly, y::Integer) = x*fmpz(y)

*(x::Integer, y::nmod_poly) = y*x

function *(x::nmod_poly, y::nmod)
  (base_ring(x) != parent(y)) && error("Must have same parent")
  return x*y.data
end

*(x::nmod, y::nmod_poly) = y*x

function +(x::nmod_poly, y::UInt)
  z = parent(x)()
  ccall((:nmod_poly_add_ui, :libflint), Void,
    (Ref{nmod_poly}, Ref{nmod_poly}, UInt), z, x, y)
  return z
end

+(x::UInt, y::nmod_poly) = y + x

function +(x::nmod_poly, y::fmpz)
  z = parent(x)()
  t = fmpz()
  tt = UInt(0)
  ccall((:fmpz_mod_ui, :libflint), UInt,
                (Ref{fmpz}, Ref{fmpz}, UInt), t, y, parent(x).n)
  tt = ccall((:fmpz_get_ui, :libflint), UInt, (Ref{fmpz}, ), t)
  return +(x,tt)
end

+(x::fmpz, y::nmod_poly) = y + x

+(x::nmod_poly, y::Integer) = x + fmpz(y)

+(x::Integer, y::nmod_poly) = y + x

function +(x::nmod_poly, y::nmod)
  (base_ring(x) != parent(y)) && error("Elements must have same parent")
  return +(x, y.data)
end

+(x::nmod, y::nmod_poly) = y + x

function -(x::nmod_poly, y::UInt)
  z = parent(x)()
  ccall((:nmod_poly_sub_ui, :libflint), Void,
    (Ref{nmod_poly}, Ref{nmod_poly}, UInt), z, x, y)
  return z
end

-(x::UInt, y::nmod_poly) = -(y - x)

function -(x::nmod_poly, y::fmpz)
  z = parent(x)()
  t = fmpz()
  tt = UInt(0)
  ccall((:fmpz_mod_ui, :libflint), UInt,
                (Ref{fmpz}, Ref{fmpz}, UInt), t, y, parent(x).n)
  tt = ccall((:fmpz_get_ui, :libflint), UInt, (Ref{fmpz}, ), t)
  return -(x,tt)
end

-(x::fmpz, y::nmod_poly) = -(y - x)

-(x::nmod_poly, y::Integer) = x - fmpz(y)

-(x::Integer, y::nmod_poly) = -(y - x)

function -(x::nmod_poly, y::nmod)
  (base_ring(x) != parent(y)) && error("Elements must have same parent")
  return -(x,y.data)
end

-(x::nmod, y::nmod_poly) = -(y - x)

################################################################################
#
#  Powering
#
################################################################################

function ^(x::nmod_poly, y::Int)
  y < 0 && throw(DomainError())
  z = parent(x)()
  ccall((:nmod_poly_pow, :libflint), Void,
          (Ref{nmod_poly}, Ref{nmod_poly}, Int), z, x, y)
  return z
end

################################################################################
#
#  Comparison
#
################################################################################

function ==(x::nmod_poly, y::nmod_poly)
  check_parent(x, y)
  return Bool(ccall((:nmod_poly_equal, :libflint), Int32,
          (Ref{nmod_poly}, Ref{nmod_poly}), x, y))
end

isequal(x::nmod_poly, y::nmod_poly) = x == y

################################################################################
#
#  Ad hoc comparisons
#
################################################################################

function ==(x::nmod_poly, y::nmod)
  base_ring(x) != parent(y) && error("Incompatible base rings in comparison")
  if length(x) > 1
    return false
  elseif length(x) == 1
    u = ccall((:nmod_poly_get_coeff_ui, :libflint), UInt,
            (Ref{nmod_poly}, Int), x, 0)
    return u == y
  else
    return iszero(y)
  end
end

==(x::nmod, y::nmod_poly) = y == x

################################################################################
#
#  Truncation
#
################################################################################

function truncate(a::nmod_poly, n::Int)
  n < 0 && throw(DomainError())
  z = deepcopy(a)
  if length(z) <= n
    return z
  end
  ccall((:nmod_poly_truncate, :libflint), Void,
          (Ref{nmod_poly}, Int), z, n)
  return z
end

function mullow(x::nmod_poly, y::nmod_poly, n::Int)
  check_parent(x, y)
  n < 0 && throw(DomainError())
  z = parent(x)()
  ccall((:nmod_poly_mullow, :libflint), Void,
          (Ref{nmod_poly}, Ref{nmod_poly}, Ref{nmod_poly}, Int), z, x, y, n)
  return z
end

###############################################################################
#
#   Reversal
#
###############################################################################

function reverse(x::nmod_poly, len::Int)
  len < 0 && throw(DomainError())
  z = parent(x)()
  ccall((:nmod_poly_reverse, :libflint), Void,
          (Ref{nmod_poly}, Ref{nmod_poly}, Int), z, x, len)
  return z
end

###############################################################################
#
#   Shifting
#
###############################################################################

function shift_left(x::nmod_poly, len::Int)
  len < 0 && throw(DomainError())
  z = parent(x)()
  ccall((:nmod_poly_shift_left, :libflint), Void,
          (Ref{nmod_poly}, Ref{nmod_poly}, Int), z, x, len)
  return z
end

function shift_right(x::nmod_poly, len::Int)
  len < 0 && throw(DomainError())
  z = parent(x)()
  ccall((:nmod_poly_shift_right, :libflint), Void,
            (Ref{nmod_poly}, Ref{nmod_poly}, Int), z, x, len)
  return z
end

###############################################################################
#
#   Exact division
#
###############################################################################

function divexact(x::nmod_poly, y::nmod_poly)
  check_parent(x, y)
  iszero(y) && throw(DivideError())
  !lead_isunit(y) && error("Impossible inverse in divexact")
  z = parent(x)()
  ccall((:nmod_poly_div, :libflint), Void,
          (Ref{nmod_poly}, Ref{nmod_poly}, Ref{nmod_poly}), z, x, y)
  return z
end

################################################################################
#
#  Ad hoc exact division
#
################################################################################

function divexact(x::nmod_poly, y::nmod)
  base_ring(x) != parent(y) && error("Elements must have same parent")
  iszero(y) && throw(DivideError())
  return divexact(x, parent(x)(y))
end

function divexact(x::nmod_poly, y::fmpz)
  iszero(y) && throw(DivideError())
  return divexact(x, parent(x)(y))
end

div(x::nmod_poly, y::nmod_poly) = divexact(x,y)

function divexact(x::nmod_poly, y::Int)
  y == 0 && throw(DivideError())
  return divexact(x, parent(x)(y))
end

################################################################################
#
#  Division with remainder
#
################################################################################

function divrem(x::nmod_poly, y::nmod_poly)
  check_parent(x,y)
  iszero(y) && throw(DivideError())
  !lead_isunit(y) && error("Impossible inverse in divrem")
  q = parent(x)()
  r = parent(x)()
  ccall((:nmod_poly_divrem, :libflint), Void,
          (Ref{nmod_poly}, Ref{nmod_poly}, Ref{nmod_poly}, Ref{nmod_poly}),
          q, r, x, y)
  return q, r
end

################################################################################
#
#  Remainder
#
################################################################################

function rem(x::nmod_poly, y::nmod_poly)
  check_parent(x,y)
  iszero(y) && throw(DivideError())
  !lead_isunit(y) && error("Impossible inverse in rem")
  z = parent(x)()
  ccall((:nmod_poly_rem, :libflint), Void,
          (Ref{nmod_poly}, Ref{nmod_poly}, Ref{nmod_poly}), z, x, y)
  return z
end

mod(x::nmod_poly, y::nmod_poly) = rem(x,y)

################################################################################
#
#  GCD
#
################################################################################

function gcd(x::nmod_poly, y::nmod_poly)
  check_parent(x,y)
  !is_prime(modulus(x)) && error("Modulus not prime in gcd")
  z = parent(x)()
  ccall((:nmod_poly_gcd, :libflint), Void,
          (Ref{nmod_poly}, Ref{nmod_poly}, Ref{nmod_poly}), z, x, y)
  return z
end

function gcdx(x::nmod_poly, y::nmod_poly)
  check_parent(x,y)
  !is_prime(modulus(x)) && error("Modulus not prime in gcdx")
  g = parent(x)()
  s = parent(x)()
  t = parent(x)()
  ccall((:nmod_poly_xgcd, :libflint), Void,
          (Ref{nmod_poly}, Ref{nmod_poly}, Ref{nmod_poly}, Ref{nmod_poly},
           Ref{nmod_poly}), g, s, t, x, y)
  return g,s,t
end

function gcdinv(x::nmod_poly, y::nmod_poly)
  check_parent(x,y)
  !is_prime(modulus(x)) && error("Modulus not prime in gcdinv")
  length(y) <= 1 && error("Length of second argument must be >= 2")
  g = parent(x)()
  s = parent(x)()
  ccall((:nmod_poly_gcdinv, :libflint), Void,
          (Ref{nmod_poly}, Ref{nmod_poly}, Ref{nmod_poly}, Ref{nmod_poly}),
          g, s, x, y)
  return g,s
end

################################################################################
#
#  Modular arithmetic
#
################################################################################

function invmod(x::nmod_poly, y::nmod_poly)
  length(y) == 0 && error("Second argument must not be 0")
  check_parent(x,y)
  if length(y) == 1
    return parent(x)(inv(eval(x,coeff(y,0))))
  end
  z = parent(x)()
  r = ccall((:nmod_poly_invmod, :libflint), Int32,
          (Ref{nmod_poly}, Ref{nmod_poly}, Ref{nmod_poly}), z, x, y)
  r == 0 ? error("Impossible inverse in invmod") : return z
end

function mulmod(x::nmod_poly, y::nmod_poly, z::nmod_poly)
  check_parent(x,y)
  check_parent(y,z)
  w = parent(x)()
  ccall((:nmod_poly_mulmod, :libflint), Void,
        (Ref{nmod_poly}, Ref{nmod_poly}, Ref{nmod_poly}, Ref{nmod_poly}),
        w, x, y, z)
  return w
end

function powmod(x::nmod_poly, e::Int, y::nmod_poly)
  check_parent(x,y)
  z = parent(x)()

  if e < 0
    g, x = gcdinv(x, y)
    if g != 1
      error("Element not invertible")
    end
    e = -e
  end

  ccall((:nmod_poly_powmod_ui_binexp, :libflint), Void,
        (Ref{nmod_poly}, Ref{nmod_poly}, Int, Ref{nmod_poly}), z, x, e, y)

  return z
end

################################################################################
#
#  Resultant
#
################################################################################

function resultant(x::nmod_poly, y::nmod_poly,  check::Bool = true)
  if check
    check_parent(x,y)
    !is_prime(modulus(x)) && error("Modulus not prime in resultant")
  end
  r = ccall((:nmod_poly_resultant, :libflint), UInt,
          (Ref{nmod_poly}, Ref{nmod_poly}), x, y)
  return base_ring(x)(r)
end

################################################################################
#
#  Evaluation
#
################################################################################

function evaluate(x::nmod_poly, y::nmod)
  base_ring(x) != parent(y) && error("Elements must have same parent")
  z = ccall((:nmod_poly_evaluate_nmod, :libflint), UInt,
              (Ref{nmod_poly}, UInt), x, y.data)
  return parent(y)(z)
end

################################################################################
#
#  Derivative
#
################################################################################

function derivative(x::nmod_poly)
  z = parent(x)()
  ccall((:nmod_poly_derivative, :libflint), Void,
        (Ref{nmod_poly}, Ref{nmod_poly}), z, x)
  return z
end

################################################################################
#
#  Integral
#
################################################################################

function integral(x::nmod_poly)
  z = parent(x)()
  ccall((:nmod_poly_integral, :libflint), Void,
        (Ref{nmod_poly}, Ref{nmod_poly}), z, x)
  return z
end

################################################################################
#
#  Composition
#
################################################################################

function compose(x::nmod_poly, y::nmod_poly)
  check_parent(x,y)
  z = parent(x)()
  ccall((:nmod_poly_compose, :libflint), Void,
          (Ref{nmod_poly}, Ref{nmod_poly}, Ref{nmod_poly}), z, x, y)
  return z
end

################################################################################
#
#  Interpolation
#
################################################################################

function interpolate(R::NmodPolyRing, x::Array{nmod, 1},
                                      y::Array{nmod, 1})
  z = R()

  ax = Array{UInt}(length(x))
  ay = Array{UInt}(length(y))

  for i in 1:length(x)
    ax[i] = x[i].data

    ay[i] = y[i].data
  end
  ccall((:nmod_poly_interpolate_nmod_vec, :libflint), Void,
          (Ref{nmod_poly}, Ref{UInt}, Ref{UInt}, Int),
          z, ax, ay, length(x))
  return z
end

################################################################################
#
#  Inflation and Deflation
#
################################################################################

function inflate(x::nmod_poly, n::Int)
  n < 0 && throw(DomainError())
  z = parent(x)()
  ccall((:nmod_poly_inflate, :libflint), Void,
          (Ref{nmod_poly}, Ref{nmod_poly}, UInt), z, x, UInt(n))
  return z
end

function deflate(x::nmod_poly, n::Int)
  n < 0 && throw(DomainError())
  z = parent(x)()
  ccall((:nmod_poly_deflate, :libflint), Void,
          (Ref{nmod_poly}, Ref{nmod_poly}, UInt), z, x, UInt(n))
  return z
end

################################################################################
#
#  Lifting
#
################################################################################

doc"""
    function lift(R::FmpzPolyRing, y::nmod_poly)
> Lift from a polynomial over $\mathbb{Z}/n\mathbb{Z}$ to a polynomial over
> $\mathbb{Z}$ with minimal reduced nonnegative coefficients. The ring `R`
> specifies the ring to lift into.
"""
function lift(R::FmpzPolyRing, y::nmod_poly)
  z = fmpz_poly()
  ccall((:fmpz_poly_set_nmod_poly, :libflint), Void,
          (Ref{fmpz_poly}, Ref{nmod_poly}), z, y)
  z.parent = R
  return z
end

################################################################################
#
#  Irreducibility
#
################################################################################

doc"""
    isirreducible(x::nmod_poly)
> Return `true` if $x$ is irreducible, otherwise return `false`.
"""
function isirreducible(x::nmod_poly)
  !is_prime(modulus(x)) && error("Modulus not prime in isirreducible")
  return Bool(ccall((:nmod_poly_is_irreducible, :libflint), Int32,
          (Ref{nmod_poly}, ), x))
end

################################################################################
#
#  Squarefree testing
#
################################################################################

doc"""
    issquarefree(x::nmod_poly)
> Return `true` if $x$ is squarefree, otherwise return `false`.
"""
function issquarefree(x::nmod_poly)
   !is_prime(modulus(x)) && error("Modulus not prime in issquarefree")
   return Bool(ccall((:nmod_poly_is_squarefree, :libflint), Int32,
       (Ref{nmod_poly}, ), x))
end

################################################################################
#
#  Factorization
#
################################################################################

doc"""
    factor(x::nmod_poly)
> Return the factorisation of $x$.
"""
function factor(x::nmod_poly)
  fac, z = _factor(x)
  return Fac(parent(x)(z), fac)
end

function _factor(x::nmod_poly)
  !is_prime(modulus(x)) && error("Modulus not prime in factor")
  fac = nmod_poly_factor(x.mod_n)
  z = ccall((:nmod_poly_factor, :libflint), UInt,
          (Ref{nmod_poly_factor}, Ref{nmod_poly}), fac, x)
  res = Dict{nmod_poly,Int}()
  for i in 1:fac.num
    f = parent(x)()
    ccall((:nmod_poly_factor_get_nmod_poly, :libflint), Void,
            (Ref{nmod_poly}, Ref{nmod_poly_factor}, Int), f, fac, i-1)
    e = unsafe_load(fac.exp,i)
    res[f] = e
  end
  return res, base_ring(x)(z)
end

doc"""
    factor_squarefree(x::nmod_poly)
> Return the squarefree factorisation of $x$.
"""
function factor_squarefree(x::nmod_poly)
  !is_prime(modulus(x)) && error("Modulus not prime in factor_squarefree")
  return Fac(parent(x)(lead(x)), _factor_squarefree(x))
end

function _factor_squarefree(x::nmod_poly)
  fac = nmod_poly_factor(x.mod_n)
  ccall((:nmod_poly_factor_squarefree, :libflint), UInt,
          (Ref{nmod_poly_factor}, Ref{nmod_poly}), fac, x)
  res = Dict{nmod_poly,Int}()
  for i in 1:fac.num
    f = parent(x)()
    ccall((:nmod_poly_factor_get_nmod_poly, :libflint), Void,
            (Ref{nmod_poly}, Ref{nmod_poly_factor}, Int), f, fac, i-1)
    e = unsafe_load(fac.exp,i)
    res[f] = e
  end
  return res
end

doc"""
    factor_distinct_deg(x::nmod_poly)
> Return the distinct degree factorisation of a squarefree polynomial $x$.
"""
function factor_distinct_deg(x::nmod_poly)
  !issquarefree(x) && error("Polynomial must be squarefree")
  !is_prime(modulus(x)) && error("Modulus not prime in factor_distinct_deg")
  degs = Array{Int}(degree(x))
  degss = [ pointer(degs) ]
  fac = nmod_poly_factor(x.mod_n)
  ccall((:nmod_poly_factor_distinct_deg, :libflint), UInt,
          (Ref{nmod_poly_factor}, Ref{nmod_poly}, Ptr{Ref{Int}}),
          fac, x, degss)
  res = Dict{Int,nmod_poly}()
  for i in 1:fac.num
    f = parent(x)()
    ccall((:nmod_poly_factor_get_nmod_poly, :libflint), Void,
            (Ref{nmod_poly}, Ref{nmod_poly_factor}, Int), f, fac, i-1)
    res[degs[i]] = f
  end
  return res
end

function factor_shape(x::PolyElem{T}) where {T <: RingElem}
  res = Dict{Int, Int}()
  square_fac = factor_squarefree(x)
  for (f, i) in square_fac
    discdeg = factor_distinct_deg(f)
    for (j,g) in discdeg
      num = div(degree(g), j)*i
      if haskey(res, j)
        res[j] += num
      else
        res[j] = num
      end
    end
  end
  return res
end

################################################################################
#
#   Remove and valuation
#
################################################################################

doc"""
    remove(z::nmod_poly, p::nmod_poly)
> Computes the valuation of $z$ at $p$, that is, the largest $k$ such that
> $p^k$ divides $z$. Additionally, $z/p^k$ is returned as well.
>
> See also `valuation`, which only returns the valuation.
"""
function remove(z::nmod_poly, p::nmod_poly)
   check_parent(z,p)
   iszero(z) && error("Not yet implemented")
   z = deepcopy(z)
   v = ccall((:nmod_poly_remove, :libflint), Int,
               (Ref{nmod_poly}, Ref{nmod_poly}), z,  p)
   return v, z
end

function divides(z::nmod_poly, x::nmod_poly)
   q, r = divrem(z, x)
   return iszero(r), q
end

################################################################################
#
#  Speedups for rings over nmod_poly
#
################################################################################

function det(M::Generic.Mat{nmod_poly})
   rows(M) != cols(M) && error("Not a square matrix in det")

   if is_prime(modulus(base_ring(M)))
     return det_popov(M)
   end

   try
      return det_fflu(M)
   catch
      return det_df(M)
   end
end

################################################################################
#
#  Unsafe functions
#
################################################################################

function zero!(x::nmod_poly)
  ccall((:nmod_poly_zero, :libflint), Void,
                   (Ref{nmod_poly},), x)
  return x
end

function one!(a::nmod_poly)
  ccall((:nmod_poly_one, :libflint), Void, (Ref{nmod_poly}, ), a)
  return a
end

function fit!(x::nmod_poly, n::Int)
  ccall((:nmod_poly_fit_length, :libflint), Void,
                   (Ref{nmod_poly}, Int), x, n)
  return nothing
end

function setcoeff!(x::nmod_poly, n::Int, y::UInt)
  ccall((:nmod_poly_set_coeff_ui, :libflint), Void,
                   (Ref{nmod_poly}, Int, UInt), x, n, y)
  return x
end

function setcoeff!(x::nmod_poly, n::Int, y::Int)
  ccall((:nmod_poly_set_coeff_ui, :libflint), Void,
                   (Ref{nmod_poly}, Int, UInt), x, n, mod(y, x.mod_n))
  return x
end

function setcoeff!(x::nmod_poly, n::Int, y::fmpz)
  r = ccall((:fmpz_fdiv_ui, :libflint), UInt, (Ref{fmpz}, UInt), y, x.mod_n)
  ccall((:nmod_poly_set_coeff_ui, :libflint), Void,
                   (Ref{nmod_poly}, Int, UInt), x, n, r)
  return x
end

setcoeff!(x::nmod_poly, n::Int, y::Integer) = setcoeff!(x, n, fmpz(y))

setcoeff!(x::nmod_poly, n::Int, y::nmod) = setcoeff!(x, n, y.data)

function add!(z::nmod_poly, x::nmod_poly, y::nmod_poly)
  ccall((:nmod_poly_add, :libflint), Void,
          (Ref{nmod_poly}, Ref{nmod_poly},  Ref{nmod_poly}), z, x, y)
  return z
end

function addeq!(z::nmod_poly, y::nmod_poly)
  ccall((:nmod_poly_add, :libflint), Void,
          (Ref{nmod_poly}, Ref{nmod_poly},  Ref{nmod_poly}), z, z, y)
  return z
end

function sub!(z::nmod_poly, x::nmod_poly, y::nmod_poly)
  ccall((:nmod_poly_sub, :libflint), Void,
          (Ref{nmod_poly}, Ref{nmod_poly},  Ref{nmod_poly}), z, x, y)
  return z
end

function mul!(z::nmod_poly, x::nmod_poly, y::nmod_poly)
  ccall((:nmod_poly_mul, :libflint), Void,
          (Ref{nmod_poly}, Ref{nmod_poly},  Ref{nmod_poly}), z, x, y)
  return z
end

function mul!(z::nmod_poly, x::nmod_poly, y::UInt)
  ccall((:nmod_poly_scalar_mul_nmod, :libflint), Void,
            (Ref{nmod_poly}, Ref{nmod_poly}, UInt), z, x, y)
  return z
end

################################################################################
#
#  Promotion rules
#
################################################################################

promote_rule(::Type{nmod_poly}, ::Type{V}) where {V <: Integer} = nmod_poly

promote_rule(::Type{nmod_poly}, ::Type{fmpz}) = nmod_poly

promote_rule(::Type{nmod_poly}, ::Type{nmod}) = nmod_poly

###############################################################################
#
#   Polynomial substitution
#
###############################################################################

function (f::nmod_poly)(a::nmod)
   if parent(a) != base_ring(f)
      return subst(f, a)
   end
   return evaluate(f, a)
end

################################################################################
#
#  Parent object call overloads
#
################################################################################

function (R::NmodPolyRing)()
  z = nmod_poly(R.n)
  z.parent = R
  return z
end

function (R::NmodPolyRing)(x::fmpz)
  r = ccall((:fmpz_fdiv_ui, :libflint), UInt, (Ref{fmpz}, UInt), x, R.n)
  z = nmod_poly(R.n, r)
  z.parent = R
  return z
end

function (R::NmodPolyRing)(x::UInt)
  z = nmod_poly(R.n, x)
  z.parent = R
  return z
end

function (R::NmodPolyRing)(x::Integer)
  z = nmod_poly(R.n, x)
  z.parent = R
  return z
end

function (R::NmodPolyRing)(x::nmod_poly)
   R != parent(x) && error("Wrong parents")
   return x
end

function (R::NmodPolyRing)(x::nmod)
  base_ring(R) != parent(x) && error("Wrong parents")
  z = nmod_poly(R.n, x.data)
  z.parent = R
  return z
end

function (R::NmodPolyRing)(arr::Array{fmpz, 1})
  z = nmod_poly(R.n, arr)
  z.parent = R
  return z
end

function (R::NmodPolyRing)(arr::Array{UInt, 1})
  z = nmod_poly(R.n, arr)
  z.parent = R
  return z
end

(R::NmodPolyRing)(arr::Array{T, 1}) where {T <: Integer} = R(map(base_ring(R), arr))

function (R::NmodPolyRing)(arr::Array{nmod, 1})
  if length(arr) > 0
     (base_ring(R) != parent(arr[1])) && error("Wrong parents")
  end
  z = nmod_poly(R.n, arr)
  z.parent = R
  return z
end

function (R::NmodPolyRing)(x::fmpz_poly)
  z = nmod_poly(R.n, x)
  z.parent = R
  return z
end

################################################################################
#
#  Polynomial ring constructor
#
################################################################################

function PolynomialRing(R::NmodRing, s::AbstractString; cached=true)
   parent_obj = NmodPolyRing(R, Symbol(s), cached)

   return parent_obj, parent_obj([R(0), R(1)])
end


text:
module PreliminarySelection
    using GLMNet
    using ..ModelSelection

    export lasso!, lasso, lassoselection, PRELIMINARYSELECTION_EXTRAKEY, _lasso!, _lasso

    include("const.jl")
    include("strings.jl")
    include("utils.jl")
    include("core.jl")
end


text:
@testset "Vals" begin

    @test Vals() === Vals{Void}
    @test Vals(Int64) === Vals{MetaPair{Int64,Void}}
    @test Vals(Int64,Float64,Bool) === Vals{MetaPair{Int64,MetaPair{Float64,MetaPair{Bool,Void}}}}

    @test (v = Vals(); v = push!(v, Int64); v === Vals(Int64))
    @test (v = Vals(Int64); v = push!(v, Float64); v === Vals(Int64,Float64))

    @test (v = Vals(Int64); (p,v) = pop!(v); v === Vals() && p === Int64)
    @test (v = Vals(Int64,Float64); (p,v) = pop!(v); v === Vals(Int64) && p === Float64)
    @test (v = Vals(Int64,Float64,Bool); (p,v) = pop!(v); v === Vals(Int64,Float64) && p === Bool)

    @test (v = Vals(); v = unshift!(v, Int64); v === Vals(Int64))
    @test (v = Vals(Int64); v = unshift!(v, Float64); v === Vals(Float64,Int64))

    @test (v = Vals(Int64); (p,v) = shift!(v); v === Vals() && p === Int64)
    @test (v = Vals(Int64,Float64); (p,v) = shift!(v); v === Vals(Float64) && p === Int64)
    @test (v = Vals(Int64,Float64,Bool); (p,v) = shift!(v); v === Vals(Float64,Bool) && p === Int64)

end


text:
# Debug Macros

for (mac,fn,lvl) in ((:debug,    :(Logging.debug),    Logging.DEBUG),
                     (:info,     :(Logging.info),     Logging.INFO),
                     (:warn,     :(Logging.warn),     Logging.WARNING),
                     (:err,      :(Logging.err),      Logging.ERROR),
                     (:error,    :(Logging.err),      Logging.ERROR),
                     (:critical, :(Logging.critical), Logging.CRITICAL))

    @eval macro $mac(msg...)
        if length(msg) > 0 && isa(msg[1], Logging.Logger)
            level = msg[1].level
        else
            level = Logging._root.level
        end

        if $lvl > level
            :nothing
        else
            Expr(:call, $fn, msg...)
        end
    end

end


text:
using ITensors
using ITensorTTN

N = 3
z = 3
χ = 5

@show N, z, χ

s = siteinds("S=1/2", N; conserve_qns=true)
ψ = InfiniteTTN(s, n -> isodd(n) ? "↑" : "↓"; z=z, linkdims=χ)

@assert degree(ψ) == z

for g in 0:(N + 2)
  println()
  for b in 1:(z - 1)
    @show g, b
    @show commonind(ψ[g, 1], ψ[g + 1, b])
  end
end

# Transfer matrix at generation 2 transforming
# from branch 1 to 2 (towards the center)
t12 = TransferMatrix(ψ, 2, 1 => 2)

@show input_inds(t12)
@show output_inds(t12)

v1 = randomITensor(QN(), input_inds(t12))
v2 = randomITensor(QN(), output_inds(t12))
@show (dag(v2) * (t12 * v1))[]

# Transfer matrix at generation 1 transforming
# from branch 2 to 1 (towards the center)
t21 = TransferMatrix(ψ, 1, 2 => 1)

@show input_inds(t21)
@show output_inds(t21)

# Composition of transfer matrices
t_tot = t21 * t12

@show input_inds(t_tot)
@show output_inds(t_tot)

v1 = randomITensor(QN(), input_inds(t_tot))
@show inds(t_tot * v1)


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

julia_version = v"1.5.3"

# The version of this JLL is decoupled from the upstream version.
# Whenever we package a new upstream release, we initially map its
# version X.Y.Z to X00.Y00.Z00 (i.e., multiply each component by 100).
# So for example version 2.6.3 would become 200.600.300.
#
# Together, this allows to increment the patch level of the JLL for minor tweaks.
# If a rebuild of the JLL is needed which keeps the upstream version identical
# but breaks ABI compatibility for any reason, one can increment the minor or major
# version (depending on whether package using this JLL use `~` or `^` compat entries)
# e.g. go from 200.600.300 to 200.601.300 or 201.600.300
# Similar tricks can also be used to package prerelease versions; e.g. one might
# map a prerelease of 2.7.0 to 200.690.000.

name = "SDPA"
version = v"700.300.800"
upstream_version = v"7.3.8"

# Collection of sources required to build SDPABuilder
sources = [
    ArchiveSource("https://sourceforge.net/projects/sdpa/files/sdpa/sdpa_$(upstream_version).tar.gz",
                  "c7541333da2f0bb2d18e90dbf758ac7cc099f3f7da3f256b284b0725f96d4117")
    DirectorySource("./bundled")
]

MUMPS_seq_version = v"4.10.0"
MUMPS_seq_packagespec = PackageSpec(; name = "MUMPS_seq_jll",
                                    uuid = "d7ed1dd3-d0ae-5e8e-bfb4-87a502085b8d",
                                    version = MUMPS_seq_version)

METIS_version = v"4.0.3"
METIS_packagespec = PackageSpec(; name = "METIS_jll",
                                uuid = "d00139f3-1899-568f-a2f0-47f597d42d70",
                                version = METIS_version)

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/sdpa-*

# Remove misleading libtool files
rm -f ${prefix}/lib/*.la
update_configure_scripts

# Apply patches
atomic_patch -p1 $WORKSPACE/srcdir/patches/shared.diff
mv configure.in configure.ac
atomic_patch -p1 $WORKSPACE/srcdir/patches/lt_init.diff
autoreconf -vi

export CPPFLAGS="${CPPFLAGS} -I${prefix}/include -I$prefix/include/coin"
export CXXFLAGS="${CXXFLAGS} -std=c++11"
if [[ ${target} == *mingw* ]]; then
    export LDFLAGS="-L$prefix/bin"
elif [[ ${target} == *linux* ]]; then
    export LDFLAGS="-ldl -lrt"
fi

./configure --prefix=$prefix --with-pic --disable-pkg-config --build=${MACHTYPE} --host=${target} \
--enable-shared lt_cv_deplibs_check_method=pass_all \
--with-blas="-lopenblas" --with-lapack="-lopenblas" \
--with-mumps-libs="-L${prefix}/lib -ldmumps -lzmumps -lcmumps -lsmumps -lmumps_common -lmpiseq -lpord -lmetis -lopenblas -lgfortran -lpthread" \
--with-mumps-include="-I${prefix}/include/mumps_seq"

make -j${nproc}
make install

## Then build the libcxxwrap-julia wrapper
cd $WORKSPACE/srcdir/sdpawrap

mkdir build
cd build/

if [[ $target == i686-* ]] || [[ $target == arm-* ]]; then
    export processor=pentium4
else
    export processor=x86-64
fi

# Override compiler ID to silence the horrible "No features found" cmake error
if [[ $target == *"apple-darwin"* ]]; then
  macos_extra_flags="-DCMAKE_CXX_COMPILER_ID=AppleClang -DCMAKE_CXX_COMPILER_VERSION=10.0.0 -DCMAKE_CXX_STANDARD_COMPUTED_DEFAULT=11"
fi

cmake $macos_extra_flags \
      -DCMAKE_FIND_ROOT_PATH=${prefix} \
      -DCMAKE_INSTALL_PREFIX=${prefix} \
      -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
      -DJulia_PREFIX=${prefix} \
      -DSDPA_DIR=$prefix \
      -DMUMPS_INCLUDE_DIR="../../destdir/include/mumps_seq" \
      -DSDPA_LIBRARY="-lsdpa" \
      -D_GLIBCXX_USE_CXX11_ABI=1 \
      ..
cmake --build . --config Release --target install

if [[ $target == *w64-mingw32* ]] ; then
    cp $WORKSPACE/destdir/lib/libsdpawrap.dll ${libdir}
fi

install_license $WORKSPACE/srcdir/sdpa-*/COPYING
"""

# The products that we will ensure are always built
products = [
    ExecutableProduct("sdpa", :sdpa),
    LibraryProduct("libsdpa", :libsdpa),
    LibraryProduct("libsdpawrap", :libsdpawrap)
]

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
include("../../L/libjulia/common.jl")
platforms = libjulia_platforms(julia_version)
platforms = expand_cxxstring_abis(platforms)
platforms = expand_gfortran_versions(platforms)
filter!(p -> libgfortran_version(p) >= v"4", platforms)

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("libcxxwrap_julia_jll"),
    Dependency("OpenBLAS32_jll"),
    Dependency("CompilerSupportLibraries_jll"),
    BuildDependency(PackageSpec(name="libjulia_jll", version=julia_version)),
    BuildDependency(MUMPS_seq_packagespec),
    BuildDependency(METIS_packagespec),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
    preferred_gcc_version=v"8",
    julia_compat = "$(julia_version.major).$(julia_version.minor)")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
typealias NullableSubArray{T,N,P<:NullableArray,IV,LD} SubArray{T,N,P,IV,LD}

@inline function Base.isnull(V::NullableSubArray, I::Int...)
    @boundscheck checkbounds(V, I...)
    @inbounds return V.parent.isnull[Base.reindex(V, V.indexes, I)...]
end

@inline function Base.values(V::NullableSubArray, I::Int...)
    @boundscheck checkbounds(V, I...)
    @inbounds return V.parent.values[Base.reindex(V, V.indexes, I)...]
end

typealias FastNullableSubArray{T,N,P<:NullableArray,IV} SubArray{T,N,P,IV,true}

@inline function Base.isnull(V::FastNullableSubArray, i::Int)
    @boundscheck checkbounds(V, i)
    @inbounds return V.parent.isnull[V.first_index + V.stride1*i-1]
end

@inline function Base.values(V::FastNullableSubArray, i::Int)
    @boundscheck checkbounds(V, i)
    @inbounds return V.parent.values[V.first_index + V.stride1*i-1]
end

# We can avoid a multiplication if the first parent index is a Colon or UnitRange
typealias FastNullableContiguousSubArray{T,N,P<:NullableArray,I<:Tuple{Union{Colon, UnitRange}, Vararg{Any}}} SubArray{T,N,P,I,true}

@inline function Base.isnull(V::FastNullableContiguousSubArray, i::Int)
    @boundscheck checkbounds(V, i)
    @inbounds return V.parent.isnull[V.first_index + i - 1]
end

@inline function Base.values(V::FastNullableContiguousSubArray, i::Int)
    @boundscheck checkbounds(V, i)
    @inbounds return V.parent.values[V.first_index + i - 1]
end


text:
"""
```
B, flags = pmf(A; ωₘ, slope, dhₘ, dh₀, cellsize)
```
Applies the progressive morphological filter by *Zhang et al. (2003)* [^zhang2003] to `A`.

# Output
- `B::Array{T,2}` Maximum allowable values
- `flags::Array{Float64,2}` A sized array with window sizes if filtered, zero if not filtered.

Afterwards, one can retrieve the resulting mask for `A` by `A .<= B` or `flags .== 0.`.

# Arguments
- `A::Array{T,2}` Input Array
- `ωₘ::Float64=20.` Maximum window size [m]
- `slope::Float64=0.01` Terrain slope [m/m]
- `dhₘ::Float64=2.5` Maximum elevation threshold [m]
- `dh₀::Float64=0.2` Initial elevation threshold [m]
- `cellsize::Float64=1.` Cellsize in [m]

[^zhang2003]: Zhang, Keqi, Shu-Ching Chen, Dean Whitman, Mei-Ling Shyu, Jianhua Yan, and Chengcui Zhang. “A Progressive Morphological Filter for Removing Nonground Measurements from Airborne LIDAR Data.” IEEE Transactions on Geoscience and Remote Sensing 41, no. 4 (2003): 872–82. <https://doi.org/10.1109/TGRS.2003.810682>.
"""
function pmf(A::AbstractMatrix{T};
    ωₘ::Float64 = 20.0,
    slope::Float64 = 0.01,
    dhₘ::Float64 = 2.5,
    dh₀::Float64 = 0.2,
    cellsize::Float64 = 1.0,
    circular = false) where {T<:Real}

    # Compute windowsizes and thresholds
    ωₘ = round(Int, ωₘ / cellsize)
    κ_max = floor(Int, log2(ωₘ - 1))  # determine # iterations based on exp growth
    windowsizes = Int.(exp2.(1:κ_max)) .+ 1

    # Compute tresholds
    dwindows = vcat(windowsizes[1], windowsizes)  # prepend first element so we get 0 as diff
    window_diffs = [dwindows[i] - dwindows[i-1] for i = 2:length(dwindows)]
    height_tresholds = [min(dhₘ, slope * window_diff * cellsize + dh₀) for window_diff in window_diffs]
    @info "Using the following thresholds: $height_tresholds for the following windows: $windowsizes"

    # Set up arrays
    Af = copy(A)  # array to be morphed
    nan_mask = isnan.(Af)
    Af[nan_mask] .= Inf  # Replace NaN with Inf, as to always filter these

    B = copy(A)  # max_elevation raster
    out = copy(A)  # max_elevation raster

    flags = zeros(size(A))  # 0 = ground, other values indicate window size
    flags[nan_mask] .= NaN

    mask = falses(size(A))

    # Iterate over window sizes and height tresholds
    for (ωₖ, dhₜ) in zip(windowsizes, height_tresholds)
        if circular
            opening_circ!(Af, ωₖ, out)
        else
            opening!(Af, ωₖ, out)
        end
        mask .= (A .- Af) .> dhₜ
        for I in eachindex(flags)
            if mask[I] && flags[I] == 0
                flags[I] = ωₖ
            end
        end
        B .= min.(B, Af .+ dhₜ)
    end

    B, flags
end


text:
# NOTE: The interface has been modeled as far as possible on the EzXML.jl XML
# parser. This is so EzXML can be a drop in replacement. The benefits of this
# parser over EzXML is that it has no dependencies. It is all pure Julia code.
# The downside is that it just supports the most common XML features
import Base: setindex!, getindex, haskey, findfirst

export  Node, Document, ElementNode, TextNode, AttributeNode,
        nodename, iselement, istext, isattribute, hasroot, hasnode,
        nodecontent,
        countnodes, countattributes,
        nodes, elements, textnodes, attributes, eachattribute,
        root, setroot!,
        addchild!, addchildren!, addelement!,
        # XPath Query API
        locatefirst

"All nodes in XML DOM is some type of Node."
abstract type Node end

"Top level of an XML DOM"
mutable struct Document
    rootnode::Union{Node, Nothing}
end

"Top level of an XML DOM"
function Document()
    Document(nothing)
end

"""
Represents an attribute inside a tag.

# Example
Here `class` and `name` are examples of attributs belonging to parent node
`widget`.

    <widget class="QCheckBox" name="checkBox">
"""
mutable struct AttributeNode <: Node
    # parent::Node
    name::String
    value::String
end

"XML Node which can contain attributes and child nodes"
mutable struct ElementNode <: Node
    # parent::Node
    name::String
    attributes::Vector{AttributeNode}
    children::Vector{Node}
end

function ElementNode(name::AbstractString)
    ElementNode(name, AttributeNode[], Node[])
end

"Represents the text found between two tags. E.g. in `<foo>bar</foo>` bar is the `TextNode`"
mutable struct TextNode <: Node
    # parent::Node
    content::String
end

"Creates an element node named `name` with a text node containing text `value`"
function ElementNode(name::AbstractString, value::AbstractString)
    ElementNode(name, AttributeNode[], Node[TextNode(value)])
end

"Element node with children `nodes`"
function ElementNode(name::AbstractString, nodes::Array{T}) where T <: Node
    ElementNode(name, AttributeNode[], nodes)
end

"""Element node with attributes given like `ElementNode("widget", ["class"=>class, "name"=>name])`"""
function ElementNode(name::AbstractString, attributes::Vector{Pair{String, String}})
    ElementNode(name, [AttributeNode(name, value) for (name, value) in attributes], Node[])
end

###### Public API #############

function getindex(n::ElementNode, key::AbstractString)
    for m in n.attributes
        if m.name == key
            return m.value
        end
    end
    error("No attribute with key $key exist")
end

function setindex!(n::ElementNode, value::AbstractString, key::AbstractString)
    ii = findall(m->m.name == key, n.attributes)
    if isempty(ii)
        push!(n.attributes, AttributeNode(key, value))
    else
        n.attributes[ii[1]] = AttributeNode(key, value)
    end
end

getindex(n::ElementNode, i::Integer) = n.children[i]

"Get all child nodes under node `n`"
nodes(n::Node) = Node[]
nodes(n::ElementNode) = n.children

"Get all elements under node `n`"
elements(n::Node) = ElementNode[]
elements(n::ElementNode) = filter(iselement, nodes(n))

textnodes(n::Node) = TextNode[]
textnodes(n::ElementNode) = filter(istext, nodes(n))

"Get an array of attributes under node `n`"
attributes(n::Node) = AttributeNode[]
attributes(n::ElementNode) = n.attributes

"Gets a dictionary of attributes meant to use in a for loop for iteration"
eachattribute(n::Node) = AttributeNode[]
eachattribute(n::ElementNode) = n.attributes

"For an XML tag looking like `<foo>bar</foo>` the `nodename` would be foo"
nodename(n::Node) = ""
nodename(n::TextNode) = "text"
nodename(n::ElementNode) = n.name

"Check if node is an element node. Element nodes can contain child nodes"
iselement(n::Node) = false
iselement(n::ElementNode) = true

"Check if node is a text node. Text nodes represents the text you find between XML tags."
istext(n::Node) = false
istext(n::TextNode) = true

"Check if node is an attribute node. Attributes are on the form `name = \"value\"`"
isattribute(n::Node) = false
isattribute(n::AttributeNode) = true

"Number of child nodes"
countnodes(n::Node) = 0
countnodes(n::ElementNode) = length(n.children)

"Number of attributes. Typically onlye element nodes have attributes"
countattributes(n::Node) = 0
countattributes(n::ElementNode) = length(n.attributes)

"Check if a root node has been set of XML document"
hasroot(doc::Document) = doc.rootnode != nothing

"Get root node. Make sure you check if it exists first with `hasroot(n)`"
root(n::Document) = n.rootnode
setroot!(n::Document) = n.rootnode = n

"Get content of all text nodes under `n`"
nodecontent(n::TextNode) = n.content
nodecontent(n::Node) = join(map(nodecontent, nodes(n)))

"""
    hasnode(node)
Return if `node` has a child node.
"""
hasnode(n::Node) = false
hasnode(n::ElementNode) = !isempty(n.children)

function addelement!(parent::Node, name::AbstractString)
    child = ElementNode(name)
    addchild!(parent, child)
    child
end

"Add `child` node to `parent` node"
function addchild!(parent::Node, child::Node)
    error("Can't add children to nodes of type $(typeof(p))")
end

addchild!(parent::ElementNode, child::Node) = push!(parent.children, child)

function addchildren!(p::Node, children)
    error("Can't add children to nodes of type $(typeof(p))")
end

"""
    addchildren!(parent, children::Vector{Pair{String, String}})

A convenience function for easily adding child elements to a parent node `p`.

# Examples

    addchildren!(node, ["x" => "10", "y" => "20"])
"""
function addchildren!(p::ElementNode, children::Vector{Pair{String, String}})
    for child in children
        addchild!(p, ElementNode(first(child), last(child)))
    end
end

"""
    addchildren!(parent, children::Vector{Node})

For easily adding multiple children to a parent node.

# Examples

    addchildren!(node, [ElementNode("foo"), ElementNode("bar")])
"""
function addchildren!(p::ElementNode, children::Vector{T}) where T <: Node
    append!(p.children, children)
end

###### XPath API ######

"""
    locatefirst(xpath, node)
Locate first node with given path. This is simplified xpath syntax.

    vertibrates/warm-blooded/mammals/humans

"""
function locatefirst(xpath::AbstractString, node::Node)
    locatefirst(split(xpath, '/'), node)
end

"""
    locatefirst(xpath_array, node)
Locate first node with given path. This is simplified xpath syntax.

    locatefirst(["vertibrates", "warm-blooded", "mammals", "humans"], node)

"""
function locatefirst(xpath::Array{T}, node::Node) where T <: AbstractString
    if isempty(xpath)
        return node
    end
    tag = xpath[1]
    for n in nodes(node)
        if nodename(n) == tag
            return locatefirst(xpath[2:end], n)
        end
    end
    nothing
end

function  locatefirst(xpath::Array{T}, doc::Document) where T <: AbstractString
    if isempty(xpath) ||  !hasroot(doc)
        return nothing
    end

    r = root(doc)
    tag = xpath[1]

    if tag != nodename(r)
        return nothing
    end
    locatefirst(xpath[2:end], r)
end

function locatefirst(xpath::AbstractString, doc::Document)
    locatefirst(split(xpath, '/'), doc)
end

"""
    locatefirst(name, attribute, value, node) -> Node
Finds first node with `name` which has an `attribute` with `value`. E.g.
to locate a node `<egg foobar="spam"/>` you could write:
    locatefirst("egg", "foobar", "spam", parent_node)
"""
function locatefirst(xpath::AbstractString, attribute::AbstractString, value::AbstractString, node::Node)
    locatefirst(split(xpath, '/'), attribute, value, node)
end

function locatefirst(xpath::AbstractString, attribute::AbstractString, value::AbstractString, doc::Document)
    locatefirst(split(xpath, '/'), attribute, value, doc)
end

function locatefirst(xpath::Array{T}, attribute::AbstractString, value::AbstractString, doc::Document) where T <: AbstractString
    if isempty(xpath) ||  !hasroot(doc)
        return nothing
    end

    r = root(doc)
    tag = xpath[1]

    if tag != nodename(r)
        return nothing
    end
    locatefirst(xpath[2:end], attribute, value, r)
end

function locatefirst(xpath::Array{T}, attribute::AbstractString, value::AbstractString, node::Node) where T <: AbstractString
    if isempty(xpath) return nothing end
    name = xpath[end]
    m = locatefirst(xpath[1:end-1], node)
    if m == nothing
        return nothing
    end
    for n in nodes(m)
       if nodename(n) == name && haskey(n, attribute) && n[attribute] == value
           return n
       end
    end
    return nothing
end

"""
    haskey(node, attribute) -> Bool
Check if XML node has a particular attribute. E.g.`haskey(n, "foobar")`
would return `true` for `<egg foobar="spam"/>` but `false` for `<foobar egg="spam"/>`
"""
function haskey(n::Node, key::AbstractString)
    any(m->m.name == key, n.attributes)
end