text:
export HIPError

struct HIPError <: Exception
    code::hipError_t
    msg::AbstractString
end
Base.show(io::IO, err::HIPError) = print(io, "HIPError(code $(err.code), $(err.msg))")

function HIPError(code::hipError_t)
    msg = status_message(code)
    return HIPError(code, msg)
end

function status_message(status)
    if status == hipSuccess
        return "the operation completed successfully"
    elseif status == hipErrorOutOfMemory
        return "out of memory"
    elseif status == hipErrorNotInitialized
        return "not initialized"
    elseif status == hipErrorDeinitialized
        return "already deinitialized"
    elseif status == hipErrorProfilerDisabled
        return "profiler disabled"
    elseif status == hipErrorProfilerNotInitialized
        return "profiler not initialized"
    elseif status == hipErrorProfilerAlreadyStarted
        return "profiler already started"
    elseif statps == hipErrorProfilerAlreadyStopped
        return "profiler already stopped"
    elseif status == hipErrorInsufficientDriver
        return "insufficient driver"
    elseif status == hipErrorInvalidImage
        return "invalid image"
    elseif status == hipErrorInvalidContext
        return "invalid context"
    elseif status == hipErrorContextAlreadyCurrent
        return "context is already current"
    elseif status == hipErrorMapFailed
        return "map failed"
    elseif status == hipErrorUnmapFailed
        return "unmap failed"
    elseif status == hipErrorArrayIsMapped
        return "array is already mapped"
    elseif status == hipErrorAlreadyMapped
        return "already mapped"
    elseif status == hipErrorNoBinaryForGpu
        return "no binary for GPU"
    elseif status == hipErrorAlreadyAcquired
        return "already acquired"
    elseif status == hipErrorNotMapped
        return "not mapped"
    elseif status == hipErrorNotMappedAsArray
        return "not mapped as array"
    elseif status == hipErrorNotMappedAsPointer
        return "not mapped as pointer"
    elseif status == hipErrorECCNotCorrectable
        return "ECC is not correctable"
    elseif status == hipErrorUnsupportedLimit
        return "unsupported limit"
    elseif status == hipErrorContextAlreadyInUse
        return "context is already in use"
    elseif status == hipErrorPeerAccessUnsupported
        return "peer access not supported"
    elseif status == hipErrorInvalidKernelFile # in CUDA CUDA_ERROR_INVALID_PTX
        return "invalid kernel file"
    elseif status == hipErrorInvalidGraphicsContext
        return "invalid graphics context"
    elseif status == hipErrorInvalidSource
        return "invalid source"
    elseif status == hipErrorFileNotFound
        return "file not found"
    elseif status == hipErrorSharedObjectSymbolNotFound
        return "shared object symbol not found"
    elseif status == hipErrorSharedObjectInitFailed
        return "shared object initialization failed"
    elseif status == hipErrorOperatingSystem
        return "operating system"
    elseif status == hipErrorSetOnActiveProcess
        return "set on active process"
    elseif status == hipErrorInvalidHandle
        return "invalid handle"
    elseif status == hipErrorNotFound
        return "not found"
    elseif status == hipErrorIllegalAddress
        return "illegal address"
    elseif status == hipErrorInvalidSymbol
        return "invalid symbol"
    # runtime errors:
    elseif status == hipErrorMissingConfiguration
        return "missing configuration"
    elseif status == hipErrorMemoryAllocation
        return "error allocating memory"
    elseif status == hipErrorInitializationError
        return "initialization error"
    elseif status == hipErrorLaunchFailure
        return "an exception occurred on the device while executing a kernel"
    elseif status == hipErrorPriorLaunchFailure
        return "failure before launch"
    elseif status == hipErrorLaunchTimeOut
        return "launch timed out"
    elseif status == hipErrorLaunchOutOfResources
        return "out of resources for launch"
    elseif status == hipErrorInvalidDeviceFunction
        return "invalid device function"
    elseif status == hipErrorInvalidConfiguration
        return "invalid configuration"
    elseif status == hipErrorInvalidDevice
        return "DeviceID must be in range 0...#compute-devices"
    elseif status == hipErrorInvalidValue
        return "one or more arguments outside of acceptable range"
    elseif status == hipErrorInvalidDevicePointer
        return "invalid device pointer"
    elseif status == hipErrorInvalidMemcpyDirection
        return "invalid memory copy direction"
    elseif hipErrorInvalidResourceHandle
        return "resource handle (hipEvent_t or hipStream_t) invalid"
    elseif status == hipErrorNotReady # not actually an error
        return "asynchronous operations not completed"
    elseif status == hipErrorNoDevice
        return "no device found"
    elseif status == hipErrorPeerAccessAlreadyEnabled
        return "peer access was already enabled from the current device"
    elseif status == hipErrorPeerAccessNotEnabled
        return "peer access was never enabled from the current device"
    elseif status == hipErrorRuntimeMemory
        return "HSA runtime memory call returned error"
    elseif status == hipErrorRuntimeOther
        return "HSA runtime error"
    elseif status == hipErrorHostMemoryAlreadyRegistered
        return "tried to lock a page-lockded memory"
    elseif status == hipErrorHostMemoryNotRegistered
        return "tried to unlock a non-page-locked memory"
    elseif status == hipErrorMapBufferObjectFailed
        return "IPC memory attach failed from ROCr"
    elseif status == hipAssert
        return "kernel calls produced assert"
    elseif status == hipErrorNotSupported
        return "hip API is not supported/implemented"
    else
        return "unknown error"
    end
end

macro check(hip_func)
    quote
        local err::hipError_t
        err = $(esc(hip_func::Expr))
        if err != hipSuccess && err != hipErrorNotReady
            throw(HIPError(err))
        end
        err
    end
end


text:
function snoop(package, tomlpath, snoopfile, outputfile, reuse = false, blacklist = [])

    command = """
    using Pkg, PackageCompiler
    """

    if tomlpath !== nothing
        command *= """
        empty!(Base.LOAD_PATH)
        # Take LOAD_PATH from parent process
        append!(Base.LOAD_PATH, $(repr(Base.LOAD_PATH)))
        Pkg.activate($(repr(tomlpath)))
        Pkg.instantiate()
        Pkg.resolve()
        Pkg.instantiate()
        """
    end

    command *= """
    # let's wrap the snoop file in a try catch...
    # This way we still do some snooping even if there is an error in the tests!
    try
        include($(repr(snoopfile)))
    catch e
        @warn("Snoop file errored. Precompile statements were recorded untill error!", exception = e)
    end
    """

    # let's use a file in the PackageCompiler dir,
    # so it doesn't get lost if later steps fail
    tmp_file = package_folder("precompile_tmp.jl")
    if !reuse
        run_julia(command, compile = "all", O = 0, g = 1, trace_compile = tmp_file)
    end
    used_packages = Set{String}() # e.g. from test/REQUIRE
    if package !== nothing
        push!(used_packages, string(package))
    end
    usings = ""
    if tomlpath !== nothing
        # add toml packages, in case extract_used_packages misses a package
        deps = get(TOML.parsefile(tomlpath), "deps", Dict{String, Any}())
        union!(used_packages, string.(keys(deps)))
    end
    if !isempty(used_packages)
        packages = join(setdiff(used_packages,string.(blacklist)), ", ")
        usings *= """
        using $packages
        for Mod in [$packages]
            isdefined(Mod, :__init__) && Mod.__init__()
        end
        """
    end

    line_idx = 0; missed = 0
    open(outputfile, "w") do io
        println(io, """
        # We need to use all used packages in the precompile file for maximum
        # usage of the precompile statements.
        # Since this can be any recursive dependency of the package we AOT compile,
        # we decided to just use them without installing them. An added
        # benefit is, that we can call __init__ this way more easily, since
        # incremental sysimage compilation won't call __init__ on `using`
        # https://github.com/JuliaLang/julia/issues/22910
        $usings
        # bring recursive dependencies of used packages and standard libraries into namespace
        for Mod in Base.loaded_modules_array()
            if !Core.isdefined(@__MODULE__, nameof(Mod))
                Core.eval(@__MODULE__, Expr(:const, Expr(:(=), nameof(Mod), Mod)))
            end
        end
        """)
        for line in eachline(tmp_file)
            line_idx += 1
            # replace function instances, which turn up as typeof(func)().
            # TODO why would they be represented in a way that doesn't actually work?
            line = replace(line, r"typeof\(([\u00A0-\uFFFF\w_!´\.]*@?[\u00A0-\uFFFF\w_!´]+)\)\(\)" => s"\1")
            # Is this ridicilous? Yes, it is! But we need a unique symbol to replace `_`,
            # which otherwise ends up as an uncatchable syntax error
            line = replace(line, r"\b_\b" => "🐃")
            try
                expr = Meta.parse(line, raise = true)
                if expr.head != :incomplete
                    # after all this, we still need to wrap into try catch,
                    # since some anonymous symbols won't be found...
                    println(io, "try;", line, "; catch e; @debug \"couldn't precompile statement $line_idx\" exception = e; end")
                else
                    missed += 1
                    @debug "Incomplete line in precompile file: $line"
                end
            catch e
                missed += 1
                @debug "Parse error in precompile file: $line" exception=e
            end
        end
    end
    @info "used $(line_idx - missed) out of $line_idx precompile statements"
    outputfile
end

function snoop_packages(packages::Symbol...; blacklist = [], file = package_folder("incremental_precompile.jl"))
    finaltoml = Dict{Any, Any}(
        "deps" => Dict(),
        "compat" => Dict(),
    )
    toml_path = package_folder("Project.toml")
    manifest_dict = Dict{String, Vector{Dict{String, Any}}}()
    open(file, "w") do compile_io
        println(compile_io, "# Precompile file for $(join(packages, " "))")
        # make sure we have all packages from toml installed
        println(compile_io, """
        using Pkg
        empty!(Base.LOAD_PATH)
        # Take LOAD_PATH from parent process
        append!(Base.LOAD_PATH, $(repr(Base.LOAD_PATH)))
        Pkg.activate($(repr(toml_path)))
        Pkg.instantiate()
        """)
        for package in packages
            precompiles = package_folder(string(package), "incremental_precompile.jl")
            toml, snoopfile = package_toml(package)
            snoop(package, toml, snoopfile, precompiles, false, blacklist)
            pkg_toml = TOML.parsefile(toml)
            manifest = joinpath(dirname(toml), "Manifest.toml")
            if isfile(manifest) # not all get a manifest (only if pkg operations are executed I suppose)
                pkg_manifest = TOML.parsefile(manifest)
                for (name, pkgs) in pkg_manifest
                    pkg_vec = get!(()-> Dict{String, Any}[], manifest_dict, name)
                    append!(pkg_vec, pkgs); unique!(pkg_vec)
                end
            end
            merge!(finaltoml["deps"], get(pkg_toml, "deps", Dict()))
            merge!(finaltoml["compat"], get(pkg_toml, "compat", Dict()))
            println(compile_io)
            write(compile_io, read(precompiles))
        end
    end
    finaltoml["name"] = "PackagesPrecompile"
    write_toml(toml_path, finaltoml)
    manifest_path = package_folder("Manifest.toml")
    # make sure we don't reuse old manifests
    isfile(manifest_path) && rm(manifest_path)
    if !isempty(manifest_dict)
        write_toml(manifest_path, manifest_dict)
    end
    return toml_path, file
end

ext: jl


text:
module DataDepsGenerators
using Gumbo, Cascadia, AbstractTrees
using JSON
using HTTP
using Dates

using InteractiveUtils: subtypes

export generate, citation_text, remove_cite_version
export UCI, GitHub, DataDryad, DataOneV2, CKAN, DataCite, Figshare, JSONLD_Web, JSONLD_DOI

abstract type DataRepo end

struct GeneratorError{T<:DataRepo} <: Exception
    repo::T
    message::String
end

GeneratorError(repo)=GeneratorError(repo, "")

Base.showerror(io::IO, e::GeneratorError) = print(io, e.repo, " generator was not suitable. $(e.message)")

include("utils.jl")

include("metadata.jl")
include("generate.jl")
include("format_codeblock.jl")

include("misc_extractors.jl")

include("APIs/UCI.jl")
include("APIs/GitHub.jl")
include("APIs/DataDryad.jl")
include("APIs/DataOneV2/DataOneV2.jl")
include("APIs/CKAN.jl")
include("APIs/DataCite.jl")
include("APIs/Figshare.jl")
include("APIs/JSONLD/JSONLD.jl")

end # module


text:
using Base.Broadcast: broadcasted
bstyle = Base.Broadcast.DefaultArrayStyle{1}()

@testset "broadcast" begin
    for (title, mr,br,sr) in (("OneTo", DynamicAxis(10),         OneTo(10),              static(OneTo(10))),
                              ("StepRangeLen", mutable(StepRangeLen(1., 1., 10)), StepRangeLen(1., 1., 10), static(StepRangeLen(1., 1., 10))))
        @testset "$title" begin
            @test broadcasted(bstyle, +, mr) == broadcast(+, br) == broadcasted(bstyle, +, sr)
            @test broadcasted(bstyle, -, mr) == broadcast(-, br) == broadcasted(bstyle, -, sr)
            @test broadcasted(bstyle, +, 1, mr) == broadcast(+, 1, br) == broadcasted(bstyle, +, 1, sr)
            @test broadcasted(bstyle, -, 1, mr) == broadcast(-, 1, br) == broadcasted(bstyle, -, 1, sr)
            @test broadcasted(bstyle, +, mr, 1) == broadcast(+, br, 1) == broadcasted(bstyle, +, sr, 1)
            @test broadcasted(bstyle, -, mr, 1) == broadcast(-, br, 1) == broadcasted(bstyle, -, sr, 1)
            @test broadcasted(bstyle, *, 1, mr) == broadcast(*, 1, br) == broadcasted(bstyle, *, 1, sr)
            @test broadcasted(bstyle, *, mr, 1) == broadcast(*, br, 1) == broadcasted(bstyle, *, sr, 1)
            @test broadcasted(bstyle, \, 1, mr) == broadcast(\, 1, br) == broadcasted(bstyle, \, 1, sr)
            @test broadcast(bstyle, /, mr, 2) == broadcast(bstyle, /, sr, 2) == broadcast(bstyle, /, br, 2)
            @test -(mr, (2 .* br)) == -(br, (2 .* br)) == -(sr, (2 .* br))
            @test -((2 .* br), mr) == -((2 .* br), br) == -((2 .* br), sr)
            @test +(mr, (1 .+ br)) == +(br, (1 .+ br)) == +(sr, (1 .+ br))
            @test +((1 .+ br), mr) == +((1 .+ br), br) == +((1 .+ br), sr)
            @test *(mr, 2) == *(sr, 2) == *(br, 2)
            @test *(2, mr) == *(2, sr) == *(2, br)
        end
    end
end


text:
using StatsBase: countmap

rand(DCHSBM_sampler(x->1, [1,1,1,2,2], fill(1, 5), 3), 10)
rand(DCHSBM_sampler(x->1/FBD.sorted_unique_count(x)^2, repeat(1:3, inner=10000000), fill(1, 30000000), 2))

Z = [1,1,1,2,2]
θ = 1:5
kmax = 3
m = 1_000_000
intensity = x -> hash(x)/typemax(UInt)+1 # determinitic random ∈ [1,2]
sampler = DCHSBM_sampler(intensity, Z, θ, kmax)

graph = rand(sampler, m)

nodes = []
for e in graph
    append!(nodes, e)
end
degrees = countmap(nodes)

@testset "shape" begin
    @test length(graph) == m
    @test all(length.(graph) .== kmax)
    @test Set(keys(degrees)) == Set(axes(Z, 1))
end

tolerance = 6
max_rtol = ceil(tolerance/sqrt(minimum(values(degrees)))*10_000)/10_000
@testset "degree distribution (rtol ≤ $max_rtol)" begin
    for group = [1:3, 4:5]
        for i in group
            for j in group
                rtol = tolerance/sqrt(min(degrees[i], degrees[j]))
                @test rtol ≤ max_rtol
                @test degrees[i]/degrees[j] ≈ θ[i]/θ[j] rtol=rtol
            end
        end
    end
end

max_stars_repo_name: LilithHafner/FBD

max_issues_repo_name: LilithHafner/FBD

max_forks_repo_name: LilithHafner/FBD


text:
export MagneticFieldMeasurement, saveMagneticFieldAsHDF5, loadMagneticFieldMeasurement, loadMagneticField, addMeasuredPosition

Base.@kwdef mutable struct MagneticFieldMeasurement
  "Description of the dataset."
  description::Union{String, Nothing} = nothing
  "Positions of the measured field values."
  positions::Union{Positions, Missing} = missing
  "Field values at the specific positions."
  fields::Union{Array{typeof(1.0u"T"), 2}, Missing} = missing
  "Error of the measured field values."
  fieldsError::Union{Array{typeof(1.0u"T"), 2}, Nothing} = nothing
  "Frequency of the measured field."
  fieldsFrequency::Union{Vector{typeof(1.0u"Hz")}, Nothing} = nothing
  "Applied current while measuring the matching position."
  currents::Union{Array{typeof(1.0u"A"), 2}, Nothing} = nothing
  "Timestamp of the matching position."
  timestamp::Union{Vector{DateTime}, Nothing} = nothing
  "Offset of the hall probe's active areas."
  sensorOffset::Union{Vector{typeof(1.0u"m")}, Nothing} = nothing
  "Temperature at the start of the measurement or, if defined as a vector, the temperature when measuring the matching position."
  temperature::Union{typeof(1.0u"°C"), Vector{typeof(1.0u"°C")}, Nothing} = nothing
end

function MagneticFieldMeasurement(filename::String)
  return loadMagneticFieldMeasurement(filename)
end

function MagneticFieldMeasurement(file::HDF5.File)
  return loadMagneticFieldMeasurement(file)
end

function saveMagneticFieldAsHDF5(measurement::MagneticFieldMeasurement, filename::String)
  h5open(filename, "w") do file
    saveMagneticFieldAsHDF5(measurement, file)
  end
end

function saveMagneticFieldAsHDF5(measurement::MagneticFieldMeasurement, file::HDF5.File)
  write(file, measurement.positions)

  if !isnothing(description(measurement))
    write(file, "/description", description(measurement))
  end

  if !ismissing(fields(measurement))
    write(file, "/fields", ustrip.(fields(measurement)))
  end

  if !isnothing(fieldsError(measurement))
    write(file, "/fieldsError", ustrip.(fieldsError(measurement)))
  end

  if !isnothing(fieldsFrequency(measurement))
    write(file, "/fieldsFrequency", ustrip.(fieldsFrequency(measurement)))
  end

  if !isnothing(currents(measurement))
    write(file, "/currents", ustrip.(currents(measurement)))
  end

  if !isnothing(timestamp(measurement))
    write(file, "/timestamp", string.(timestamp(measurement)))
  end

  if !isnothing(sensorOffset(measurement))
    write(file, "/sensorOffset", ustrip.(sensorOffset(measurement)))
  end

  if !isnothing(temperature(measurement))
    write(file, "/temperature", ustrip.(temperature(measurement)))
  end
end

function loadMagneticFieldMeasurement(filename::String)
  h5open(filename, "r") do file
    return loadMagneticFieldMeasurement(file)
  end
end

function loadMagneticFieldMeasurement(file::HDF5.File)
  splattingDict = Dict{Symbol, Any}()
  splattingDict[:positions] = Positions(file)

  if haskey(file, "fields")
    splattingDict[:fields] = read(file, "fields")u"T"
  else
    error("The HDF5 file for a magnetic field measurement must contain a field vector.")
  end

  if typeof(splattingDict[:positions]) == MeanderingGridPositions
    splattingDict[:fields] = splattingDict[:fields][:, getPermutation(splattingDict[:positions]), :]
    splattingDict[:positions] = splattingDict[:positions].grid
  end

  if haskey(file, "description")
    splattingDict[:description] = read(file, "description")
  end

  if haskey(file, "fieldsError")
    splattingDict[:fieldsError] = read(file, "fieldsError")u"T"
  end

  if haskey(file, "fieldsFrequency")
    splattingDict[:fieldsFrequency] = read(file, "fieldsFrequency")u"Hz"
  end

  if haskey(file, "currents")
    splattingDict[:currents] = read(file, "currents")u"A"
  end

  if haskey(file, "timestamp")
    splattingDict[:timestamp] = DateTime.(read(file, "timestamp"))
  end

  if haskey(file, "sensorOffset")
    splattingDict[:sensorOffset] = read(file, "sensorOffset")u"m"
  end

  if haskey(file, "temperature")
    splattingDict[:temperature] = read(file, "temperature")u"°C"
  end

  return MagneticFieldMeasurement(;splattingDict...)
end

# Alias for backwards compatibility
loadMagneticField(filename::String) = loadMagneticFieldMeasurement(filename)

# Create getter and setter for all fields of `MagneticFieldMeasurement`
for (fieldname, fieldtype) in zip(fieldnames(MagneticFieldMeasurement), fieldtypes(MagneticFieldMeasurement))
  fieldnameStr = string(fieldname)

  # At the moment, this should be a Union
  missingOrNothing = (fieldtype.b <: Union{Missing, Nothing}) ? fieldtype.b : fieldtype.a
  fieldtype = (fieldtype.b <: Union{Missing, Nothing}) ? fieldtype.a : fieldtype.b

  #@info "" fieldnameStr missingOrNothing fieldtype

  @eval begin
    export $fieldname

    # TODO: Add docstring from struct; I did not yet find a way to retrieve it
    function $(fieldname)(measurement::MagneticFieldMeasurement)::Union{$fieldtype, $missingOrNothing}
      return measurement.$fieldname
    end

    # TODO: Add docstring from struct; I did not yet find a way to retrieve it
    function $(fieldname)(measurement::MagneticFieldMeasurement, value::Union{$fieldtype, $missingOrNothing})
      measurement.$fieldname = value
    end
  end
end

function addMeasuredPosition(measurement::MagneticFieldMeasurement, pos::Vector; field=nothing, fieldError=nothing, fieldFrequency=nothing, current=nothing, timestamp=nothing, temperature=nothing)
  idx = posToLinIdx(measurement.positions, pos)

  if !isnothing(field)
    if ismissing(measurement.fields)
      measurement.fields = fill(0.0u"T", (length(measurement.positions), 3))
    end

    measurement.fields[idx, :] .= field
  end

  if !isnothing(fieldError)
    if isnothing(measurement.fieldsError)
      measurement.fieldsError = fill(0.0u"T", (length(measurement.positions), 3))
    end

    measurement.fieldsError[idx, :] .= fieldError
  end

  if !isnothing(fieldFrequency)
    if isnothing(measurement.fieldsFrequency)
      measurement.fieldsFrequency = fill(0.0u"Hz", length(measurement.positions))
    end

    measurement.fieldsFrequency[idx] = fieldFrequency
  end

  if !isnothing(current)
    if isnothing(measurement.currents)
      measurement.currents = fill(0.0u"A", (length(measurement.positions), length(current)))
    end

    measurement.currents[idx, :] .= current
  end

  if !isnothing(timestamp)
    if isnothing(measurement.timestamp)
      measurement.timestamp = fill(now(), length(measurement.positions))
    end

    measurement.timestamp[idx] = timestamp
  end

  if !isnothing(temperature)
    if isnothing(measurement.temperature)
      measurement.temperature = fill(-273.0u"°C", length(measurement.positions))
    end

    measurement.temperature[idx] = temperature
  end
end


text:
using SortingAlgorithms
using RCall, DataFrames, BenchmarkTools, IterTools, CSV
using Plots
#using PyCall
Plots.gr()

function string_sort_timing_cmp(N, K, idnumlen = 10)
    srand(1);
    svec1 = rand(["id"*dec(i,idnumlen) for i in 1:N÷K], N);
    julia_timing = @elapsed sort(svec1, alg = StringRadixSort);
    println(julia_timing)

    julia_timing_base = @elapsed sort(svec1);
    println(julia_timing_base)

    R"""
    # increase memory limit; only works on Windows
    if (Sys.info()["sysname"] == "Windows") memory.limit(32653)
    id3 = sample(sprintf("id%010d",1:($N/$K)), $N, TRUE) # small groups (char)
    pt = proc.time()
    sort(id3, method = "radix")
    r_timing = proc.time() - pt
    rm(id3)
    gc()
    """
    @rget r_timing

    println(r_timing[3])

    # py"""
    # import numpy as np
    # import timeit

    # # randChar is workaround for MemoryError in mtrand.RandomState.choice
    # # http://stackoverflow.com/questions/25627161/how-to-solve-memory-error-in-mtrand-randomstate-choice
    # def randChar(f, numGrp, N) :
    #     things = [f%x for x in range(numGrp)]
    #     return [things[x] for x in np.random.choice(numGrp, N)]

    # id3 = randChar("id%010d", $N//$K, $N)   # small groups (char)
    # tt = timeit.Timer("id3.sort()" ,"from __main__ import id3").timeit(1) # 6.8 seconds
    # """

    python_timing = py"tt"

    DataFrame(  lang = ["Julia", "Julia", "R", "Python"],
                alg  = ["radix", "base", "radix", "base"],
                timings = [julia_timing, julia_timing_base, r_timing[3], python_timing],
                N = repeat([N], inner = 4),
                K = repeat([K], inner = 4),
                idnumlen = repeat([idnumlen], inner = 4),
                grps=[N÷K for i=1:4])
end

# sorting a vector with no duplicates
res10m = string_sort_timing_cmp(10_000_000, 1, 10);
CSV.write("res10m.csv", res10m);
res10m = CSV.read("res10m.csv");

gb = bar(res10m[:lang].*"-".*res10m[:alg], res10m[:timings], label="Seconds")
title!(gb,"String sort performance: 10m unique id strings")
ylabel!(gb, "Seconds")
savefig(gb, "sort_perf_10m_u.png")

# sorting a vector with lots of duplicates
res10m100 = string_sort_timing_cmp(10_000_000, 100, 10);
CSV.write("res10m100.csv", res10m100);

res10m100 = CSV.read("res10m100.csv");
gb = bar(res10m100[:lang].*"-".*res10m100[:alg], res10m100[:timings], label="Seconds")
title!(gb,"String sort performance: 10m strings 100k unique values")
ylabel!(gb, "Seconds")
savefig(gb, "sort_perf_10m100_u.png")

string_sort_timing_cmp(1_000_000, 100, 3)

res = [string_sort_timing_cmp(n,k,idnumlen) for (n,k,idnumlen) = product(10.^(6:7),[1,100],[3,10])]
res = reduce(vcat,res)

sort!(res,cols = [:K,:N,:idnumlen,:timings])

res1m = string_sort_timing_cmp(1_000_000, 100,  3)
res10m = string_sort_timing_cmp(10_000_000, 100, 3)
res100m = string_sort_timing_cmp(100_000_000, 100, 3)
# res200m = string_sort_timing_cmp(200_000_000, 100, 3)
# res400m = string_sort_timing_cmp(400_000_000, 100, 3)
res = reduce(vcat,[res1m, res10m, res100m, res200m, res400m])

res1m = string_sort_timing_cmp(1_000_000, 100, 3)
res10m = string_sort_timing_cmp(10_000_000, 100,3)
res100m = string_sort_timing_cmp(100_000_000, 100,3)
# res200m = string_sort_timing_cmp(200_000_000, 100,3)
# res400m = string_sort_timing_cmp(400_000_000, 100,3)
res = reduce(vcat,[res1m, res10m, res100m, res200m, res400m])

res1m10 = string_sort_timing_cmp(1_000_000, 100)
res10m10 = string_sort_timing_cmp(10_000_000, 100)
res100m10 = string_sort_timing_cmp(100_000_000, 100)
# res200m10 = string_sort_timing_cmp(200_000_000, 100)
# res400m10 = string_sort_timing_cmp(400_000_000, 100)
res = reduce(vcat,[res1m, res10m, res100m, res200m, res400m])

N = 10_000_000
K = 100_000
x = rand(["id"*dec(i,3) for i in 1:N÷K], N);

srand(1);
N = 10_000_000
K = 1
svec1 = rand(["id"*dec(i,10) for i in 1:N÷K], N);
julia_timing = @elapsed sort(svec1, alg = StringRadixSort);
println(julia_timing)


text:
# Code from previous section
using PyPlot
using Statistics
A = imread("sample_photo.png")
B = mean(A, dims=3)[:,:,1]
function imshow_scale(A)
    # Like imshow(A) but scales the values to [0,1] and supports grayscale

    A .-= minimum(A)            # Scale and shift to [0,1]
    A ./= maximum(A)
    if ndims(A) < 3
        A = reshape(A, size(A,1), size(A,2), 1)
    end
    if size(A,3) == 1
        A = repeat(A, 1, 1, 3)  # Set R=G=B for grayscale
    end
    imshow(A)
    end;

function image_downsample(A)
    # For simplicity, ignore last row/column if size is odd
    sz = size(A)
    sz2 = sz .÷ 2
    B = (A[1:2:2sz2[1], 1:2:2sz2[2], :] .+
         A[2:2:2sz2[1], 1:2:2sz2[2], :] .+
         A[1:2:2sz2[1], 2:2:2sz2[2], :] .+
         A[2:2:2sz2[1], 2:2:2sz2[2], :]) ./ 4
end

C = copy(B);
for i = 1:4
    C = image_downsample(C);
    subplot(2,2,i), imshow_scale(C)
    title("$(i) levels downsampling")
end
tight_layout()

function image_upsample(A)
    sz = size(A)
    B = zeros(Float32, 2sz[1]-1, 2sz[2]-1, sz[3])
    B[1:2:end, 1:2:end, :] .= A
    B[2:2:end-1, 1:2:end, :] .= (A[1:end-1,:,:] .+ A[2:end,:,:]) ./ 2
    B[1:2:end, 2:2:end-1, :] .= (A[:,1:end-1,:] .+ A[:,2:end,:]) ./ 2
    B[2:2:end-1, 2:2:end-1, :] .= (A[1:end-1,1:end-1,:] .+ A[1:end-1,2:end,:]
                               .+  A[2:end,1:end-1,:] .+ A[2:end,2:end,:]) ./ 4
    return B
end

D = copy(C);
for i = 1:4
    D = image_upsample(D);
    subplot(2,2,i), imshow_scale(D)
    title("$(i) levels upsampling")
end
tight_layout()

subplot(1,2,1); imshow_scale(B); title("Original");
subplot(1,2,2); imshow_scale(D); title("Down+Up sampling");


text:
########################################################################
# Codes for the paper:
# An Anthropocene Map of Genetic Diversity
#
# Andreia Miraldo, Sen Li, Michael K. Borregaard, Alexander Floréz-Rodriguéz, Shyam Gopalakrishnan, Mirneza Risvanovic, Zhiheng Wang, Carsten Rahbek, Katharine A. Marske & David Nogués-Bravo
#
# Submitted to Science, 2016
# Code in this file by Michael K. Borregaard
#########################################################################

#Load the necessary library
using DataFrames, DataFramesMeta
include("GD_summary_functions.jl")

function boot_species(dat::DataFrame)
    by(dat, :cell) do df
        nrows = size(df, 1)
        rows = rand(1:nrows, nrows)
        mean(df[:x1][rows])
    end
end

function bootstrap_species(dat::DataFrame)
    boot_full = hcat([boot_species(dat)[:x1] for i in 1:100]...) #creates a matrix of cell values
    ret = by(df -> mean(df[:x1]), dat, :cell)
    names!(ret, [:site, :empirical])
    ret[:boot_sd] = vec(mapslices(std, boot_full, 2))
    ret[:boot_low] = vec(mapslices(x -> quantile(x,  [:0.025]), boot_full, 2))
    ret[:boot_high] = vec(mapslices(x -> quantile(x,  [:0.975]), boot_full, 2))
    ret
end

for typ in ["anthromes", "biomes", "latbands"]
    # Read in the matrices for mammals and amphibians and concatenate them
    am = ready_data("pairwise_$(typ)_amphibians_cytb.csv")
    ma = ready_data("pairwise_$(typ)_mammals_cytb.csv")
    macoi = ready_data("pairwise_$(typ)_mammals_coi.csv")
    dats = vcat(am, ma)
    # restrict to pairs with at least 0.5 overlap value

    amphibs, amphibsfull = calc_cellvalues(am)
    mamms, mammalfull  = calc_cellvalues(ma)
    macoi, macoifull = calc_cellvalues(macoi)
    total, full = calc_cellvalues(dats)

    bootfull = bootstrap_species(full)
    bootmams = bootstrap_species(mammalfull)
    bootamphibs = bootstrap_species(amphibsfull)
    bootmacoi = bootstrap_species(macoifull)

    boot_data = join(bootfull, bootmams, on = :site, kind = :left)
    boot_data = join(boot_data, bootamphibs, on = :site, kind = :left)
    boot_data = join(boot_data, bootmacoi, on = :site, kind = :left)
    names!(boot_data, [:site, :total_empirical, :total_SD, :total_0_025, :total_0_975, :mamm_empirical, :mamm_SD, :mamm_0_025, :mamm_0_975,:amph_empirical, :amph_SD, :amph_0_025, :amph_0_975,:macoi_empirical, :macoi_SD, :macoi_0_025, :macoi_0_975])

    writetable("$(typ)_bootstrap_numbers.csv", boot_data)
end

max_stars_repo_path: Miraldo et al. 2016/Calculate_Values_for_Bootstrapping_Anthromes.jl

max_stars_repo_name: mkborregaard/codes-for-published-papers

max_issues_repo_name: DalavanCloud/codes-for-published-papers

max_forks_repo_name: DalavanCloud/codes-for-published-papers


text:
function wrap_cell(ic, jc, l)

	iw = ic
	jw = jc

	if iw < 1
		iw = l
	elseif iw > l
		iw = 1
	end

	if jw < 1
		jw = l
	elseif jw > l
		jw = 1
	end

	return iw, jw
end


text:
export SpMatMul,addVectors

function SpMatMul(alpha::Union{Float64,Complex128,Complex64,Float32},AT::SparseCSCTypes,x::ArrayTypes,beta::Union{Float64,Complex128,Complex64},target::ArrayTypes,numCores::Int64)
# function calculates: target = beta*target + alpha*A*x
if hasParSpMatVec
	ParSpMatVec.Ac_mul_B!( alpha, AT, x, beta, target, numCores );
else
	Base.Ac_mul_B!(alpha,AT,x,beta,target);
end

return target;
end

function SpMatMul(AT::SparseCSCTypes,x::ArrayTypes,target::ArrayTypes,numCores::Int64)
alpha = one(eltype(x));
beta = zero(eltype(x));
if hasParSpMatVec
	ParSpMatVec.Ac_mul_B!( alpha, AT, x, beta, target, numCores );
else
	Base.Ac_mul_B!(alpha,AT,x,beta,target);
end
return target;
end

function addVectors(alpha::Union{Float64,Complex128,Complex64},x::ArrayTypes,target::ArrayTypes)
# target <- target + alpha*x;
BLAS.axpy!(length(x), alpha, x,1, target, 1);
return;
end


text:
module JuliaUtils

using DataFrames

export diff_functions_two_files

"""
```julia
function diff_functions_two_files(file_main::String, file_to_proof::String, col_name_main::String, col_name_to_proof::String)
```
    This function is supposed to get the function from file_main and check if the function name also appears in the second file.
    If this is not the case, these functions will be displayed in the result.
"""
function diff_functions_two_files(file_main::String, file_to_proof::String, col_name_main::String, col_name_to_proof::String)

    fileTextMain =open(f->read(f,String),file_main)
    fileTextToProof =open(f->read(f,String),file_to_proof)

    allMatches_Main = findall(r"(?i)(?s)function +(\w+\.|\w+)+?[a-zA-Z0-9_]*\((([^()])|(\([^()]*\)))*\)",fileTextMain)
    allMatches_to_proof = findall(r"(?i)(?s)function +(\w+\.|\w+)+?[a-zA-Z0-9_]*\((([^()])|(\([^()]*\)))*\)",fileTextToProof)

    eachFunctionPostgres = map(x->SubString(fileTextMain,x),allMatches_Main)
    eachFunctionOracle = map(x->SubString(fileTextToProof,x),allMatches_to_proof)

    main_FunctionCalls = [string(first(split(item,"("))) for item in eachFunctionOracle]
    to_proof_FunctionCalls= [string(first(split(item,"("))) for item in eachFunctionPostgres]

    overview = Matrix{Any}(missing,length(to_proof_FunctionCalls),2)

    for i in 1:length(to_proof_FunctionCalls)
        overview[i,1] = to_proof_FunctionCalls[i]
        if findfirst(x->to_proof_FunctionCalls[i] == x,main_FunctionCalls) !== nothing
            overview[i,2] = main_FunctionCalls[findfirst(x->to_proof_FunctionCalls[i] == x,main_FunctionCalls)]
        end
    end

    df = DataFrame(overview,[Symbol(col_name_main), Symbol(col_name_to_proof)])

    missedFunctions = filter(x->x[Symbol(col_name_main)] === missing,df)

end #end diff_functions_two_files

"""
```julia
function add_line_at_first_in_textfile(file_name, line)
```
This function add a text at first position at file
"""
function add_line_at_first_in_textfile(file_name, line)
    file_content = read(file_name,String)
    result_content = "$line \n" * file_content

    file_open = open(file_name,"w")
    write(file_open, result_content)
    close(file_open)
end

end


text:
tests = [
    "types",
    "fit",
    "discrete",
    "univariate",
    "truncate",
    "multinomial",
    "dirichlet",
    "multivariate",
    "mvnormal",
    "mvtdist",
    "kolmogorov",
    "edgeworth",
    "matrix",
    "vonmisesfisher",
    "compoundvariate",
    "conjugates",
    "conjugates_normal",
    "conjugates_mvnormal",
    "wishart",
    "mixture",
    "gradloglik"]

println("Running tests:")

for t in tests
    test_fn = "$t.jl"
    println(" * $test_fn")
    include(test_fn)
end


text:
function NeuralDensityEstimator(prob_layer_sizes, support)
    NeuralDensityEstimator(prob_layer_sizes, minimum(support), maximum(support))
end

function NeuralDensityEstimator(prob_layer_sizes, support_min::Real, support_max::Real)
    NeuralDensityEstimator(prob_layer_sizes, SVector(support_min), SVector(support_max))
end

"""
    NeuralDensityEstimator{N}(prob_layer_sizes, support_min::NTuple{N, Real}, support_max::NTuple{N,Real})

Defines a Density Estimator.

support_min and support_max define the boundries of the support.
They should be at least as large as the support (though for noncompact supports an approximation will be required)

`support_min` defined the "lower left corner"
and
`support_max` defines the "upper right corner"

Note the support boundries must be rectangular, if your support is not rectangular
you must make a domain transform on your observations before using this.

"""
function NeuralDensityEstimator{N}(prob_layer_sizes, support_min::StaticVector{N}, support_max::StaticVector{N})
    sess = Session(Graph())
    @tf begin
        t = placeholder(Float32, shape=[N, -1])
        @assert all(support_min .< support_max)
        smin = constant(reshape(collect(support_min), (N,1)))
        smax = constant(reshape(collect(support_max), (N,1)))

        layer_sizes= [N; prob_layer_sizes; 1]

        network_fun_stack = Function[Base.identity]

        for ii in 2:length(layer_sizes)
            below_size = layer_sizes[ii-1]
            above_size = layer_sizes[ii]

            Wii = get_variable("W_$ii", [above_size, below_size], Float32)
            Wii2  = Ops.mul(Wii, Wii; name = "W_$(ii)_squared")

            #bii = get_variable("b_$ii", [above_size, 1], Float32)
            #act_fun = z -> nn.sigmoid(Wii2*z .+ bii)

            act_fun = if ii!=length(layer_sizes)
                bii = get_variable("b_$ii", [above_size, 1], Float32)
                z -> nn.sigmoid(Wii2*z .+ bii)
            else
                z-> exp(Wii2*z)
            end
            push!(network_fun_stack, z->act_fun(network_fun_stack[ii-1](z)))
        end

        network = network_fun_stack[end]

        ysmin = TensorFlow.identity(network(smin))
        ysmax = TensorFlow.identity(network(smax))
        yt = TensorFlow.identity(network(t))

        denominator = reduce_prod(ysmax-ysmin) #area

        # Take a derivative with respect to each column, i.e. each variable
        #dytn = yt
        #for column_ii in 1:N
        #    t_col = gather(t, column_ii; name="t_$(column_ii)")
        #    grads = gradients(dytn, t_col)
        #    dytn = TensorFlow.identity; grads, name="dyt_$(column_ii)")
        #end
        dytn = gradients(yt, t)
        numerator = TensorFlow.identity(dytn)
        pdf = numerator/denominator
        @show pdf
#        @assert(get_shape(pdf,

        n_points = TensorFlow.shape(t)[2]
        loglikelihood = reduce_sum(log(numerator)) - n_points.*log(denominator)

        area_loss = (1f0.-denominator)^2
        working_loss = -1*loglikelihood + 0.1*area_loss

        optimizer = train.minimize(train.AdamOptimizer(), working_loss)

        # Conditioning
        # Make sure that ysmin~=1, and ysmax~=2
        condition_loss = (1f0 - ysmin)^2 + (2f0 - ysmax)^2
        condition_optimiser = train.minimize(train.AdamOptimizer(;name="adam_cond"), condition_loss)
    end

    run(sess, global_variables_initializer())

    NeuralDensityEstimator{N}(sess, optimizer, condition_optimiser, t, pdf)
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   International Geomagnetic Field Model.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
# ==============================================================================
#
#   [1] https://www.ngdc.noaa.gov/IAGA/vmod/igrf.html
#   [2] https://www.ngdc.noaa.gov/IAGA/vmod/igrf12.f
#   [3] https://www.mathworks.com/matlabcentral/fileexchange/34388-international-geomagnetic-reference-field--igrf--model
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export igrf, igrfd

################################################################################
#                                  Functions
################################################################################

"""
    igrfd(date::Number, [r,h]::Number, λ::Number, Ω::Number, T[, P, dP]; max_degree = 13, show_warns = true)

**IGRF Model**

*Current version: v13*

Compute the geomagnetic field vector [nT] at the date `date` [Year A.D.] and
position (`r`, `λ`, `Ω`).

The position representation is defined by `T`. If `T` is `Val(:geocentric)`,
then the input must be **geocentric** coordinates:

1. Distance from the Earth center `r` [m];
2. Geocentric latitude `λ` (-90°, +90°); and
3. Geocentric longitude `Ω` (-180°, +180°).

If `T` is `Val(:geodetic)`, then the input must be **geodetic** coordinates:

1 Altitude above the reference ellipsoid `h` (WGS-84) \\[m];
2. Geodetic latitude `λ` (-90°, +90°); and
3. Geodetic longitude `Ω` (-180°, +180°).

If `T` is omitted, then it defaults to `Val(:geocentric)`.

Notice that the output vector will be represented in the same reference system
selected by the parameter `T` (geocentric or geodetic). The Y-axis of the output
reference system always points East. In case of **geocentric coordinates**, the
Z-axis points toward the center of Earth and the X-axis completes a right-handed
coordinate system. In case of **geodetic coordinates**, the X-axis is tangent to
the ellipsoid at the selected location and points toward North, whereas the
Z-axis completes a right-hand coordinate system.

The optional arguments `P` and `dP` must be two matrices with at least
`max_degree + 1 × max_degree + 1` real numbers. If they are present, then they
will be used to store the Legendre coefficients and their derivatives. In this
case, no allocation will be performed when computing the magnetic field. If they
are not present, then 2 allocations will happen to create them.

# Keywords

- `max_degree::Int`: Maximum degree used in the spherical harmonics when
    computing the geomagnetic field. If it is higher than the available number
    of coefficients in the IGRF matrices, then it will be clamped. If it is
    equal of lower than 0, then it will be set to 1. (**Default** = 13)
- `show_warns::Bool`: Show warnings about the data (**Default** = `true`).

# Remarks

The `date` must be greater or equal to 1900 and less than or equal 2030. Notice
that a warning message is printed for dates greater than 2025.

# Disclaimer

This function is an independent implementation of the IGRF model. It contains a
more readable code than the original one in FORTRAN, because it uses features
available in Julia language.

"""
@inline function igrfd(
    date::Number,
    r::Number,
    λ::Number,
    Ω::Number;
    max_degree::Int = 13,
    show_warns::Bool = true
)
    return igrfd(
        date,
        r,
        λ,
        Ω,
        Val(:geocentric);
        max_degree = max_degree,
        show_warns = show_warns
    )
end

@inline function igrfd(
    date::Number,
    r::Number,
    λ::Number,
    Ω::Number,
    P::AbstractMatrix{T},
    dP::AbstractMatrix{T};
    max_degree::Int = 13,
    show_warns::Bool = true
) where T<:AbstractFloat
    return igrfd(
        date,
        r,
        λ,
        Ω,
        Val(:geocentric),
        P,
        dP;
        max_degree = max_degree,
        show_warns = show_warns
    )
end

@inline function igrfd(
    date::Number,
    rh::S1,
    λ::S2,
    Ω::S3,
    R::T;
    max_degree::Int = 13,
    show_warns::Bool = true
) where {
    S1<:Number,
    S2<:Number,
    S3<:Number,
    T<:Union{Val{:geocentric}, Val{:geodetic}}
}
    S = promote_type(float(S1), float(S2), float(S3))
    # Currently, the maximum degree allowed for IGRF is 13.
    max_degree = clamp(max_degree, 1, 13)

    P  = Matrix{S}(undef, max_degree + 1, max_degree + 1)
    dP = Matrix{S}(undef, max_degree + 1, max_degree + 1)

    return igrfd(
        date,
        rh,
        λ,
        Ω,
        R,
        P,
        dP;
        max_degree = max_degree,
        show_warns = show_warns
    )
end

@inline function igrfd(
    date::Number,
    rh::Number,
    λ::Number,
    Ω::Number,
    R::T,
    P::AbstractMatrix{S},
    dP::AbstractMatrix{S};
    max_degree::Int = 13,
    show_warns::Bool = true
) where {T<:Union{Val{:geocentric}, Val{:geodetic}}, S<:AbstractFloat}
    # Check if the latitude and longitude are valid.
    if ((λ < -90) || (λ > 90))
        error("The latitude must be between -90° and +90° rad.")
    end

    if ((Ω < -180) || (Ω > 180))
        error("The longitude must be between -180° and +180° rad.")
    end

    return igrf(
        date,
        rh,
        deg2rad(λ),
        deg2rad(Ω),
        R,
        P,
        dP;
        max_degree = max_degree,
        show_warns = show_warns
    )
end

"""
    igrf(date::Number, [r,h]::Number, λ::Number, Ω::Number, T[, P, dP]; max_degree = 13, show_warns = true)

**IGRF Model**

*Current version: v13*

Compute the geomagnetic field vector [nT] at the date `date` [Year A.D.] and
position (`r`, `λ`, `Ω`).

The position representation is defined by `T`. If `T` is `Val(:geocentric)`,
then the input must be **geocentric** coordinates:

1. Distance from the Earth center `r` [m];
2. Geocentric latitude `λ` (-π/2, +π/2) \\[rad]; and
3. Geocentric longitude `Ω` (-π, +π) \\[rad].

If `T` is `Val(:geodetic)`, then the input must be **geodetic** coordinates:

1. Altitude above the reference ellipsoid `h` (WGS-84) \\[m];
2. Geodetic latitude `λ` (-π/2, +π/2) \\[rad]; and
3. Geodetic longitude `Ω` (-π, +π) \\[rad].

If `T` is omitted, then it defaults to `Val(:geocentric)`.

Notice that the output vector will be represented in the same reference system
selected by the parameter `T` (geocentric or geodetic). The Y-axis of the output
reference system always points East. In case of **geocentric coordinates**, the
Z-axis points toward the center of Earth and the X-axis completes a right-handed
coordinate system. In case of **geodetic coordinates**, the X-axis is tangent to
the ellipsoid at the selected location and points toward North, whereas the
Z-axis completes a right-hand coordinate system.

The optional arguments `P` and `dP` must be two matrices with at least
`max_degree + 1 × max_degree + 1` real numbers. If they are present, then they
will be used to store the Legendre coefficients and their derivatives. In this
case, no allocation will be performed when computing the magnetic field. If they
are not present, then 2 allocations will happen to create them.

# Keywords

- `max_degree::Int`: Maximum degree used in the spherical harmonics when
    computing the geomagnetic field. If it is higher than the available number
    of coefficients in the IGRF matrices, then it will be clamped. If it is
    equal of lower than 0, then it will be set to 1. (**Default** = 13)
- `show_warns::Bool`: Show warnings about the data (**Default** = `true`).

# Remarks

The `date` must be greater or equal to 1900 and less than or equal 2030. Notice
that a warning message is printed for dates greater than 2025.

# Disclaimer

This function is an independent implementation of the IGRF model. It contains a
more readable code than the original one in FORTRAN, because it uses features
available in Julia language.

"""
function igrf(
    date::Number,
    r::Number,
    λ::Number,
    Ω::Number;
    max_degree::Int = 13,
    show_warns::Bool = true
)
    return igrf(date,
        r,
        λ,
        Ω,
        Val(:geocentric);
        max_degree = max_degree,
        show_warns = show_warns
    )
end

function igrf(
    date::Number,
    r::Number,
    λ::Number,
    Ω::Number,
    P::AbstractMatrix{T},
    dP::AbstractMatrix{T};
    max_degree::Int = 13,
    show_warns::Bool = true
) where T<:AbstractFloat
    return igrf(
        date,
        r,
        λ,
        Ω,
        Val(:geocentric),
        P,
        dP;
        max_degree = max_degree,
        show_warns = show_warns
    )
end

function igrf(
    date::Number,
    r::S1,
    λ::S2,
    Ω::S3,
    R::T;
    max_degree::Int = 13,
    show_warns::Bool = true
) where {
    S1<:Number,
    S2<:Number,
    S3<:Number,
    T<:Union{Val{:geocentric}, Val{:geodetic}}
}
    S = promote_type(float(S1), float(S2), float(S3))

    # Currently, the maximum degree allowed for IGRF is 13.
    max_degree = clamp(max_degree, 1, 13)

    P  = Matrix{S}(undef, max_degree + 1, max_degree + 1)
    dP = Matrix{S}(undef, max_degree + 1, max_degree + 1)

    return igrf(date,
        r,
        λ,
        Ω,
        R,
        P,
        dP;
        max_degree = max_degree,
        show_warns = show_warns
    )
end

function igrf(
    date::Number,
    r::Number,
    λ::Number,
    Ω::Number,
    ::Val{:geocentric},
    P::AbstractMatrix{T},
    dP::AbstractMatrix{T};
    max_degree::Int = 13,
    show_warns::Bool = true
) where T<:AbstractFloat

    # Model data
    # ==========

    max_year = 2030    # Maximum year in which we can compute.
    rel_year = 2025    # Maximum year with reliable data.
    dat_year = 2020    # Last year in which real measurements are available.

    # Input verification
    # ==================

    # Check the data, since this model is valid for years between 1900 and
    # `max_year`.
    if ((date < 1900) || (date > max_year))
        error("This IGRF version will not work for years outside the interval [1900, $max_year).")
    end

    # Check if the latitude and longitude are valid.
    if ((λ < -π/2) || (λ > π/2))
        error("The latitude must be between -π/2 and +π/2 rad.")
    end

    if ((Ω < -π) || (Ω > π))
        error("The longitude must be between -π and +π rad.")
    end

    # Warn the user that for dates after the year `rel_year` the accuracy maybe
    # reduced.
    if show_warns && (date > rel_year)
        @warn("The magnetic field computed with this IGRF version may be of reduced accuracy for years greater than $rel_year.")
    end

    # If the `max_degree` is equal or lower than 0, we must clamp it to 1.
    max_degree = max(max_degree, 1)

    # Input variables conversion
    # ==========================

    # Convert latitude / longitude to colatitude and east-longitude.
    θ = T(π / 2) - λ
    ϕ = (Ω >= 0) ? Ω : T(2π) + Ω

    # The input variable `r` is in [m], but all the algorithm requires it to be
    # in [km].
    r /= 1000

    # Preliminary setup
    # =================

    # Compute the epoch that will be used to obtain the coefficients. This is
    # necessary because the IGRF provides coefficients every 5 years. Between
    # two epochs, those coefficients must be interpolated.
    idx   = floor(Int, clamp((date - 1900) * 0.2 + 1, 0, (rel_year - 1900) / 5))
    epoch = 1900 + (idx - 1) * 5

    # Compute the fraction of time from the epoch of the coefficient selected by
    # `idx`.
    Δt = date - epoch

    # Compute the maximum spherical harmonic degree for the selected date.
    n_max = (epoch < 1995) ? 10 : 13

    # Check if the user wants a lower degree.
    n_max = min(max_degree, n_max)

    # Make sure we have the required amount of space in the matrices.
    rows, cols = size(P)

    if (rows < n_max + 1) || (cols < n_max + 1)
        error("Matrix `P` must have at least $(n_max + 1) rows and columns.")
    end

    rows, cols = size(dP)
    if ((rows < n_max + 1) || (cols < n_max + 1))
        error("Matrix `dP` must have at least $(n_max+1) rows and columns.")
    end

    # Compute the Schmidt quasi-normalized associated Legendre functions and
    # their first order derivative, neglecting the phase term.
    legendre!(Val(:schmidt), P, θ, false, n_max, n_max)
    dlegendre!(Val(:schmidt), dP, θ, P, false, n_max, n_max)

    # Parameters and auxiliary variables
    # ==================================

    # Auxiliary variables to select the IGRF coefficients.
    H = H_igrf
    G = G_igrf

    # Get the last index of the matrices G and H, which must be the same.
    ~, gend = size(G)
    ~, hend = size(H)

    (gend != hend) &&
    error("Internal error: the matrices G and H must have the same number of columns.")

    # Reference radius [km].
    a = T(6371.2)

    # Auxiliary variables to decrease the computational burden.
    sin_ϕ,  cos_ϕ  = sincos(1ϕ)
    sin_2ϕ, cos_2ϕ = sincos(2ϕ)

    ratio = a/r
    fact  = ratio

    # Initialization of variables
    # ===========================

    dVr = T(0)  # Derivative of the Geomagnetic potential w.r.t. r.
    dVθ = T(0)  # Derivative of the Geomagnetic potential w.r.t. θ.
    dVϕ = T(0)  # Derivative of the Geomagnetic potential w.r.t. ϕ.
    ΔG  = T(0)  # Auxiliary variable to interpolate the G coefficients.
    ΔH  = T(0)  # Auxiliary variable to interpolate the H coefficients.
    kg  = 1     # Index to obtain the values of the matrix `G`.
    kh  = 1     # Index to obtain the values of the matrix `H`.

    # Geomagnetic potential
    # =====================

    @inbounds for n in 1:n_max
        aux_dVr = T(0)
        aux_dVθ = T(0)
        aux_dVϕ = T(0)

        # Compute the contributions when `m = 0`
        # ======================================

        # Get the coefficients in the epoch and interpolate to the desired
        # time.
        Gnm_e0 = T(G[kg,idx+2])

        if date < dat_year
            Gnm_e1 = T(G[kg,idx+3])
            ΔG     = (Gnm_e1 - Gnm_e0) / 5
        else
            ΔG     = T(G[kg,gend])
        end

        Gnm  = Gnm_e0 + ΔG * Δt
        kg  += 1

        aux_dVr += -(n + 1) / r * Gnm * P[n+1,1]
        aux_dVθ += Gnm * dP[n+1,1]

        # Sine and cosine with m = 1
        # ==========================
        #
        # This values will be used to update recursively `sin(m*ϕ)` and
        # `cos(m*ϕ)`, reducing the computational burden.
        sin_mϕ   = +sin_ϕ    # sin( 1*λ_gc)
        sin_m_1ϕ = T(0)      # sin( 0*λ_gc)
        sin_m_2ϕ = -sin_ϕ    # sin(-1*λ_gc)
        cos_mϕ   = +cos_ϕ    # cos( 1*λ_gc)
        cos_m_1ϕ = T(1)      # cos( 0*λ_gc)
        cos_m_2ϕ = +cos_ϕ    # cos(-2*λ_gc)

        # Other auxiliary variables that depend only on `n`
        # =================================================

        fact_dVr = T(n + 1) / T(r)

        # Compute the contributions when `m ∈ [1,n]`
        # ==========================================

        for m = 1:n
            # Compute recursively `sin(m*ϕ)` and `cos(m*ϕ)`.
            sin_mϕ = 2cos_ϕ * sin_m_1ϕ - sin_m_2ϕ
            cos_mϕ = 2cos_ϕ * cos_m_1ϕ - cos_m_2ϕ

            # Compute the coefficients `G_nm` and `H_nm`
            # ==========================================

            # Get the coefficients in the epoch and interpolate to the
            # desired time.
            Gnm_e0 = T(G[kg,idx+2])
            Hnm_e0 = T(H[kh,idx+2])

            if date < dat_year
                Gnm_e1 = T(G[kg,idx+3])
                Hnm_e1 = T(H[kh,idx+3])
                ΔG     = (Gnm_e1 - Gnm_e0) / 5
                ΔH     = (Hnm_e1 - Hnm_e0) / 5
            else
                ΔG = T(G[kg,gend])
                ΔH = T(H[kh,hend])
            end

            Gnm  = Gnm_e0 + ΔG * Δt
            Hnm  = Hnm_e0 + ΔH * Δt
            kg  += 1
            kh  += 1

            GcHs_nm = Gnm * cos_mϕ + Hnm * sin_mϕ
            GsHc_nm = Gnm * sin_mϕ - Hnm * cos_mϕ

            # Compute the contributions for `m`
            # =================================

            aux_dVr += -fact_dVr * GcHs_nm * P[n+1,m+1]
            aux_dVθ += GcHs_nm * dP[n+1,m+1]
            aux_dVϕ += (θ == 0) ?
                -m * GsHc_nm * dP[n+1,m+1] :
                -m * GsHc_nm * P[n+1,m+1]

            # Update the values for the next step
            # ===================================

            sin_m_2ϕ = sin_m_1ϕ
            sin_m_1ϕ = sin_mϕ
            cos_m_2ϕ = cos_m_1ϕ
            cos_m_1ϕ = cos_mϕ
        end

        # Perform final computations related to the summation in `n`
        # ==========================================================

        # fact = (a/r)^(n+1)
        fact *= ratio

        # aux_<> *= (a/r)^(n+1)
        aux_dVr *= fact
        aux_dVϕ *= fact
        aux_dVθ *= fact

        dVr += aux_dVr
        dVϕ += aux_dVϕ
        dVθ += aux_dVθ
    end

    dVr *= a
    dVϕ *= a
    dVθ *= a

    # Compute the Geomagnetic field vector in the geocentric reference frame
    # ======================================================================

    x = +dVθ / r
    y = (θ == 0) ? -dVϕ / r : -dVϕ / (r * sin(θ))
    z = dVr

    B_gc = SVector{3, T}(x, y, z)

    return B_gc
end

function igrf(
    date::Number,
    h::Number,
    λ::Number,
    Ω::Number,
    ::Val{:geodetic},
    P::AbstractMatrix{T},
    dP::AbstractMatrix{T};
    max_degree::Int = 13,
    show_warns::Bool = true
) where T<:AbstractFloat

    # TODO: This method has a small error (≈ 0.01 nT) compared with the
    # `igrf12syn`. However, the result is exactly the same as the MATLAB
    # function in [3]. Hence, this does not seem to be an error in the
    # conversion from geodetic to geocentric coordinates. This is probably
    # caused by a numerical error. Further verification is necessary.

    # Convert the geodetic coordinates to geocentric coordinates.
    λ_gc, r = geodetic_to_geocentric(λ, h)

    # Compute the geomagnetic field in geocentric coordinates.
    B_gc = igrf(
        date,
        r,
        λ_gc,
        Ω,
        Val(:geocentric),
        P,
        dP;
        max_degree = max_degree,
        show_warns = show_warns
    )

    # Convert to geodetic coordinates.
    D_gd_gc = create_rotation_matrix(λ_gc - λ, :Y)
    B_gd    = D_gd_gc * B_gc

    return B_gd
end


text:
# # Registration Functions
#
#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/generated/example_placeholder.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/example_placeholder.ipynb)

# In this example page I will show how to work with the registration algorithms [`lap`](@ref), [`pflap`](@ref), [`sparse_lap`](@ref) and [`sparse_pflap`](@ref).

# I will generate some images to show the workings and function calls.

## using the package
using LAP_julia;

## default arguments:
img, imgw, flow = gen_init();

# See the [`Basic Interaction`](@ref basic_interaction) section or the [`Public Documentation`](@ref public_api) for more advanced input generation

# __These are the differences between the `target` (`img`) and `source` (`imgw`) images.__

imgoverlay(img, imgw, figtitle="Target vs Source")

# ## Full LAP Methods

# These functions are from the algorithms from the paper [Local All-Pass Geometric Deformations](https://www.researchgate.net/publication/320574179_Local_All-Pass_Geometric_Deformations).

# ### [`lap`](@ref) function

## set a filter half size, which is larger than the default max displacement: 10 used in gen_init
fhs = 15;
flow_est, source_reg = lap(img, imgw, fhs);

# This gives us a estimation of the flow and the registered source image

# __This is what the registered source image looks like:__
imgoverlay(img, source_reg, figtitle="lap: Target vs Registered Source")

# ### [`pflap`](@ref) function

flow_est, source_reg = pflap(img, imgw, display=false);

# This gives us a estimation of the flow and the registered source image

# __This is what the registered source image looks like:__
imgoverlay(img, source_reg, figtitle="pflap: Target vs Registered Source")

# ## Sparse LAP Methods

# These functions are insired by the [`lap`](@ref) and [`pflap`](@ref) functions.
# The difference being that the displacement vectors are calculated at points of high gradient and
# the result is fit into a global deformation.

# ### [`sparse_lap`](@ref) function

## set a filter half size, which is larger than the default max displacement: 10 used in gen_init
fhs = 15;
flow_est, source_reg = sparse_lap(img, imgw, fhs, display=false);

# This gives us a estimation of the flow and the registered source image

# __This is what the registered source image looks like:__
imgoverlay(img, source_reg, figtitle="sparse lap: Target vs Registered Source")

# ### [`sparse_pflap`](@ref) function

flow_est, source_reg = sparse_pflap(img, imgw, display=false);

# This gives us a estimation of the flow and the registered source image

# __This is what the registered source image looks like:__
imgoverlay(img, source_reg, figtitle="sparse pflap: Target vs Registered Source")


text:
import Base: convert, parse

abstract type Packet end

struct LiteralValue <: Packet
    version::Int
    value::Int
end

struct Operator <: Packet
    version::Int
    action::Function
    subpackets::Vector{Packet}
end

convert(::Type{Int}, bits::BitVector) = foldl((x, y) -> (x << 1) | y, bits)

popfirstn!(collection, n) = convert(typeof(collection), [popfirst!(collection) for _ in 1:n])

function parse(::Type{BitVector}, str::String)
    bitify(n) = collect(bitstring(n))[5:8] .== '1'
    reduce(vcat, bitify.(parse.(UInt8, collect(str), base = 16)))
end

function decode!(::Type{Packet}, bits::BitVector)
    type_id = convert(Int, bits[4:6])
    decode!(type_id == 4 ? LiteralValue : Operator, bits)
end

function decode!(::Type{LiteralValue}, bits::BitVector)
    version = convert(Int, popfirstn!(bits, 3))
    popfirstn!(bits, 3)  # skip type id as it is always 4
    value = empty(bits)
    more_groups = true
    while more_groups
        more_groups = popfirst!(bits)
        append!(value, popfirstn!(bits, 4))
    end
    LiteralValue(version, convert(Int, value .== 1))
end

function decode!(::Type{Operator}, bits::BitVector)
    version = convert(Int, popfirstn!(bits, 3))
    type_id = convert(Int, popfirstn!(bits, 3))
    action = type_id < 4 ? (sum, prod, minimum, maximum)[type_id + 1] : (>, <, ==)[type_id - 4]
    length_type_id = convert(Int, popfirst!(bits))
    subpackets = Packet[]
    if length_type_id == 1
        n_subpackets = convert(Int, popfirstn!(bits, 11))
        for _ in 1:n_subpackets
            push!(subpackets, decode!(Packet, bits))
        end
    else
        n_bits = convert(Int, popfirstn!(bits, 15))
        n_remaining_bits = length(bits) - n_bits
        while length(bits) > n_remaining_bits
            push!(subpackets, decode!(Packet, bits))
        end
    end
    Operator(version, action, subpackets)
end

versions(packet::Operator) = append!([packet.version], versions.(packet.subpackets)...)

versions(packet::LiteralValue) = packet.version

function execute(packet::Operator)
    packet.action in (sum, prod, minimum, maximum) ? (
        packet.action(execute.(packet.subpackets))
    ) : (
        packet.action(execute.(packet.subpackets)...)
    )
end

execute(packet::LiteralValue) = packet.value

bits = parse(BitVector, readline("day16/transmission.txt"))
packet = decode!(Packet, bits)

p1 = sum(versions(packet))

@show p1

p2 = execute(packet)

@show p2


text:
mutable struct SuperNodeTree
	snd::Array{Int64,1} #vertices of supernodes stored in one array
	snptr::Array{Int64,1} # vertices of supernode i are stored in snd[snptr[i]:snptr[i+1]-1]
	snd_par::Array{Int64,1}  # parent of supernode k is supernode j=snd_par[k]
	snd_post::Array{Int64,1} # post order of supernodal elimination tree
	post::Array{Int64} # post ordering of the vertices in elim tree σ(j) = v
	par::Array{Int64}
	cliques::Array{Int64,1} #vertices of cliques stored in one array
	chptr::Array{Int64,1} #points to the indizes where new clique starts in cliques array
	nBlk::Array{Int64,1} #sizes of submatrizes defined by each clique
	function SuperNodeTree(L)
		par = etree(L)
		child = child_from_par(par)
		post = post_order(par, child)
		# sorting of children from highest order one to lowest make sure that a node always gets
		# added to the supernode with the highest rep vertex
		sort_children!(child, post)
		# faster algorithm according to Vandenberghe p.308
		degrees = higher_degrees(L)
		snd, snptr, snd_par = find_supernodes(par, child, post, degrees)

		# TODO: amalgamate supernodes
		snd_child = child_from_par(snd_par)
		# # a post-ordering of the elimination tree is needed to make sure that in the
		# # psd completion step the matrix is completed from lower-right to top-left
		 snd_post = post_order(snd_par, snd_child)

		cliques, chptr, nBlk = find_cliques(L, snd, snptr, snd_par, post)

		new(snd, snptr, snd_par, snd_post, post, par, cliques, chptr, nBlk)
	end

end

function sort_children!(child, post)
	for v = 1:length(child)
		child[v] = sort(child[v], by = x -> post[x], rev = true)
	end
end
# -------------------------------------
# FUNCTION DEFINITIONS
# -------------------------------------
# given v=σ^-1(i) it returns i=σ(v)
function invert_order(sigma::Array{Int64,1})
	sigma_inv = zeros(Int64, length(sigma))
	for iii=1:length(sigma)
		sigma_inv[sigma[iii]] = iii
	end
	return sigma_inv
end

function num_cliques(sntree::SuperNodeTree)
	return length(sntree.chptr)
end

# elimination tree algorithm from H.Liu - A Compact Row Storage Scheme for Cholesky Factors Using Elimination Trees
# function etree_liu(g)
# 	N = length(g.adjacency_list)
# 	par = zeros(Int64,N)
# 	ancestor = zeros(Int64,N)

# 	elemSequence = g.reverse_ordering[collect(1:N)]
# 	for iii in elemSequence
# 		for vk in g.adjacency_list[iii]
# 			if g.ordering[vk] < g.ordering[iii]
# 				r = vk
# 				while (ancestor[r] != 0) && (ancestor[r] != iii)
# 					t = ancestor[r]
# 					ancestor[r] = iii
# 					r = t
# 				end
# 				if ancestor[r] == 0
# 					ancestor[r] = iii
# 					par[r] = iii
# 				end
# 			end
# 		end
# 	end

# 	return par
# end

# simplified version of my own elimination tree algorithm with simplified data structure (fastest)
function etree(L)
	N = size(L, 1)
	par = zeros(Int64, N)
	# loop over Vertices of graph
	for i=1:N
		value = i
		# number of i-neighbors with order higher than order of node i
		par_ = find_parent_direct(L, i)
		par[i] = par_
	end
	return par
end

# perform a depth-first-search to determine the post order of the tree defined by parent and children vectors
function post_order(par::Array{Int64,1}, child::Array{Array{Int64,1}})
	N = length(par)
	order = zeros(Int64, N)
	root = findall(x ->x == 0, par)[1]
	stack = Array{Int64}(undef, 0)
	iii = N
	push!(stack, root)
	while !isempty(stack)
		v = pop!(stack)
		order[v] = iii
		iii -= 1
		push!(stack, child[v]...)
	end
	post = collect(1:N)
	sort!(post, by = x-> order[x])
	return post
end

function child_from_par(par::Array{Int64,1})

	child = [Array{Int64}(undef, 0) for i = 1:length(par)]
	for i = 1:length(par)
		par_ = par[i]
		par_ != 0 && push!(child[par_], i)
	end
	return child
end

function get_snd(sntree::SuperNodeTree, ind::Int64)
	N = length(sntree.snptr)
	if ind == N
		return sntree.snd[sntree.snptr[ind]:end]
	else
		return sntree.snd[sntree.snptr[ind]:sntree.snptr[ind + 1] - 1]
	end
end

function get_clique(sntree::SuperNodeTree, ind::Int64)
	N = length(sntree.chptr)
	if ind == N
		return sntree.cliques[sntree.chptr[ind]:end]
	else
		return sntree.cliques[sntree.chptr[ind]:sntree.chptr[ind + 1] - 1]
	end
end

function print_cliques(sntree::SuperNodeTree)
	N = length(sntree.chptr)
	chptr = sntree.chptr
	println("Cliques of Graph:")
	for iii = 1:N
		if iii != N
			clique = sntree.cliques[chptr[iii]:chptr[iii + 1] - 1]
		else
			clique = sntree.cliques[chptr[iii]:end]
		end
		println("$(iii): $(clique)")
	end
end

function print_supernodes(sntree::SuperNodeTree)
	N = length(sntree.snptr)
	snptr = sntree.snptr
	println("Supernodes of Graph:")
	for iii = 1:N
		if iii != N
			sn = sntree.snd[snptr[iii]:snptr[iii + 1] - 1]
		else
			sn = sntree.snd[snptr[iii]:end]
		end
		println("$(iii): $(sn)")
	end
end

function check_degree_condition(v::Int64, w::Int64, degrees::Array{Int64,1})
	return degrees[v] > degrees[w] - 1
end

# Algorithm from A. Poten and C. Sun: Compact Clique Tree Data Structures in Sparse Matrix Factorizations (1989)
function pothen_sun(par::Array{Int64,1}, child::Array{Array{Int64,1}}, post::Array{Int64,1}, degrees::Array{Int64,1})
	N = length(par)
	sn_ind = -1 * ones(Int64, N) # if snInd[v] < 0 then v is a rep vertex, otherwise v ∈ supernode[snInd[v]]
	supernode_par = -1 * ones(Int64, N)

	# go through vertices in post_order
	for v in post
		child_ind = 0
		# check degree condition for all of v's childs from highest to lowest

		for (iii, w) in enumerate(child[v])
			# if not deg+(v) > deg+(w) - 1 for a certain w, set u to be w in snd(u), add v to snd(u)
			if !check_degree_condition(v, w, degrees)
				sn_ind[w] < 0 ? (u = w) : (u = sn_ind[w])
				sn_ind[v] = u
				child_ind = iii
				break
			end
		end

		# if v is still a rep vertex (i.e. above loop didnt find a child that fulfilled degree condition)
		if sn_ind[v] < 0
			u = v
		end

		for (iii, w) in enumerate(child[v])
			if iii != child_ind
				sn_ind[w] < 0 ? (x = w) : x = sn_ind[w]
				supernode_par[x] = u
			end
		end
	end

	# representative vertices
	reprv = findall(x-> x < 0, sn_ind)
	# vertices that are the parent of representative vertices
	repr_par = supernode_par[reprv]
	# take into account that all non-representative arrays are removed from the parent structure
	sn_par = zeros(Int64, length(reprv))

	for (iii, rp) in enumerate(repr_par)
		ind = findfirst(x -> x == rp, reprv)
		ind == nothing && (ind = 0)
		sn_par[iii] = ind
	end

	return sn_par, sn_ind
end

function find_supernodes(par::Array{Int64,1}, child::Array{Array{Int64,1}}, post::Array{Int64,1}, degrees::Array{Int64,1})
	supernode_par, snInd = pothen_sun(par, child, post, degrees)
	# number of vertices
	N = length(par)
	# number of representative vertices == number of supernodes
	Nrep = length(supernode_par)

	snode = zeros(Int64, N)
	snode_list = [Array{Int64}(undef, 0) for i = 1:N]
	snptr = zeros(Int64,Nrep)

	for iii in post
		f = snInd[iii]
		if f < 0
			push!(snode_list[iii], iii)

		else
			push!(snode_list[f], iii)
		end
	end

	kkk = 1
	jjj = 1
	for (iii, list) in enumerate(snode_list)
		len = length(list)
		if len > 0
			snptr[jjj] = kkk
			snode[kkk:kkk + len - 1] = list
			kkk += len
			jjj += 1
		end
	end
	return snode, snptr, supernode_par

end

function find_cliques(L, snodes::Array{Int64,1}, snptr::Array{Int64,1}, supernode_par::Array{Int64,1}, post::Array{Int64,1})
	postInv = invert_order(post)

	Nc = length(supernode_par)
	cliques = Array{Int64}(undef, 0)
	nBlk = zeros(Int64, Nc)
	chptr = zeros(Int64, Nc)
	jjj = 1

	for iii = 1:Nc
		if iii < Nc
			vRep = snodes[snptr[iii]:snptr[iii + 1] - 1][1]
		else
			vRep = snodes[snptr[iii]:end][1]
		end
		adjPlus = find_higher_order_neighbors(L, vRep)
		deg = length(adjPlus) + 1
		cliques = [cliques; vRep; adjPlus]
		chptr[iii] = jjj
		nBlk[iii] = Base.power_by_squaring(length(adjPlus) + 1, 2)
		jjj += deg
	end

	return cliques, chptr, nBlk

end

# Given a sparse lower triangular matrix L find the neighboring vertices of v
function find_neighbors(L::SparseMatrixCSC, v::Int64)
	L = Symmetric(L)
	col_ptr = L.colptr
	row_val = L.rowval
	return row_val[col_ptr[v]:col_ptr[v + 1] - 1]
end

function find_higher_order_neighbors(L::SparseMatrixCSC, v::Int64)
	v == size(L, 1) && return 0
	col_ptr = L.colptr
	row_val = L.rowval
	return row_val[col_ptr[v]:col_ptr[v + 1] - 1]
end

function find_parent_direct(L::SparseMatrixCSC, v::Int64)
	v == size(L, 1) && return 0
	col_ptr = L.colptr
	row_val = L.rowval
	return row_val[col_ptr[v]]
end

# findall the cardinality of adj+(v) for all v in V
function higher_degrees(L::SparseMatrixCSC)
	N = size(L, 1)
	degrees = zeros(Int64, N)
	col_ptr = L.colptr
	row_val = L.rowval

	for v = 1:N-1
		degrees[v] = col_ptr[v + 1] - col_ptr[v]
	end
	return degrees
end

# -------------------------------------
# FUNCTION DEFINITIONS
# -------------------------------------

# this assumes a sparse lower triangular matrix L
function connect_graph!(L::SparseMatrixCSC)
	# unconnected blocks don't have any entries below the diagonal in their right-most column
	m = size(L, 1)
	row_val = L.rowval
	col_ptr = L.colptr
	for j = 1:m-1
		connected = false
		for k in col_ptr[j]:col_ptr[j+1]-1
			if row_val[k] > j
				connected  = true
				break
			end
		end
		if !connected
			L[j+1, j] = 1
		end
	end
end

function find_graph!(ci, rows::Array{Int64, 1}, N::Int64, C::AbstractConvexSet)
	row_val, col_val = row_ind_to_matrix_indices(rows, N, C)
	F = QDLDL.qdldl(sparse(row_val, col_val, ones(length(row_val))), logical = true)
	# this takes care of the case that QDLDL returns an unconnected adjacency matrix L
	connect_graph!(F.L)
	ci.L = F.L
	return F.perm
end

# given an array [rows] that represent the nonzero entries of a vectorized NxN matrix,
# return the rows and columns of the nonzero entries of the original matrix
function row_ind_to_matrix_indices(rows::Array{Int64,1}, N::Int64, ::PsdCone{Float64})
	row_val = zeros(Int64, length(rows))
	col_val = zeros(Int64, length(rows))
	for (ind, r) in enumerate(rows)
		_rem = mod(r, N)
		fl = fld(r, N)
		if _rem == 0
			row_val[ind] = N
			col_val[ind] = fl
		else
			row_val[ind] = _rem
			col_val[ind] = fl + 1
		end
	end
	return row_val, col_val
end

# given an array [rows] that represent the nonzero entries of the vectorized upper triangular part of a NxN matrix,
# return the rows and columns of the nonzero entries of the original matrix
function row_ind_to_matrix_indices(rows::Array{Int64,1}, N::Int64, ::COSMO.PsdConeTriangle{Float64})
	#  allocate conservative here since we don't know how many diagonal entries are contained in row
	row_val = zeros(Int64, 2 * length(rows) + N)
	col_val = zeros(Int64, 2 * length(rows) + N)
	ind = 1
	_step = 1
	for (iii, r) in enumerate(rows)
		i, j = COSMO.svec_to_mat(r)

		row_val[ind] = i
		col_val[ind] = j
		ind += 1

		if i != j
			row_val[ind] = j
			col_val[ind] = i
			ind += 1
		end
	end
	return row_val[1:ind - 1], col_val[1:ind - 1]
end

# Given a linear index find the col of the corresponding upper triangular matrix
function svec_get_col(x::Int64)
	c = (sqrt(8 * x + 1) - 1) / 2
	if c % 1 != 0
		return Int((floor(c) + 1))
	else
		return Int(c)
	end
end

# Given a linear index find the row of the corresponding upper triangular matrix
function svec_get_row(x::Int64)
	c = get_col(x) - 1
	k = (c + 1) * c / 2
	return (x - k)::Int64
end

# Given a linear index find the row and column of the corresponding upper triangular matrix
function svec_to_mat(ind::Int64)
	c = svec_get_col(ind)
	k = (c - 1) * c / 2
	r = Int(ind - k)
	return r::Int64, c::Int64
end


text:
using Test
import ArchGDAL; const AG = ArchGDAL

@testset "Iterator interface Window Iterator" begin
    ds = AG.readraster("ospy/data4/aster.img")
    band = AG.getband(ds, 1)
    window = AG.windows(band)
    @test Base.IteratorSize(window) == Base.HasShape{2}()
    @test Base.IteratorEltype(window) == Base.HasEltype()
    @test eltype(window) == Tuple{UnitRange{Int},UnitRange{Int}}
    @test size(window) == (79, 89)
    @test length(window) == 7031
end


text:
include("datasstructures/rsparsebitset.jl")

"""
    TableConstraint

Efficient implentation of the constraint enforcing: ∃ j ∈ ⟦1,m⟧, such that ∀ i ∈ ⟦1,n⟧ xᵢ=table[i,j].

The datastructure and the functions using it are inspired by:
Demeulenaere J. et al. (2016) Compact-Table: Efficiently Filtering Table Constraints with Reversible Sparse Bit-Sets.
In: Rueher M. (eds) Principles and Practice of Constraint Programming. CP 2016. Lecture Notes in Computer Science, vol 9892. Springer, Cham.
https://doi.org/10.1007/978-3-319-44953-1_14

    TableConstraint(scope, active, table, currentTable, modifiedVariables, unfixedVariables, supports, residues)

TableConstraint default constructor.

# Arguments
- `scope::Vector{<:AbstractIntVar}`: the ordered variables present in the table.
- `table::Matrix{Int}`: the original table describing the constraint (impossible assignment are filtered).
- `currentTable::RSparseBitSet{UInt64}`: the reversible representation of the table.
- `modifiedVariables::Vector{Int}`: vector with the indexes of the variables modified since the last propagation.
- `unfixedVariables::Vector{Int}`: vector with the indexes of the variables which are not binding.
- `supports::Dict{Pair{Int,Int},BitVector}`: dictionnary which, for each pair (variable => value), gives the support of this pair.
- `residues::Dict{Pair{Int,Int},Int}`: dictionnary which, for each pair (variable => value), gives the residue of this pair.

    TableConstraint(variables, table, supports, trailer)

Create a CompactTable constraint from the `variables`, with values given in `table` and supports given in `supports`.

This constructor gives full control on both `table` and `supports`. The attributes are not duplicated and remains
linked to the variables given to the constructor. It should be used only to avoid duplication of `table` when using
the same table constraint many times with different variables. *WARNING*: all variables must have the same domain;
`table` and `supports` must be cleaned.

# Arguments
- `variables::Vector{<:AbstractIntVar}`: vector of variables of size (n, ).
- `table::Matrix{Int}`: matrix of the constraint of size (n, m).
- `supports::Dict{Pair{Int,Int},BitVector}`: for each assignment, a list of the tuples supporting this assignment.
- `trailer::Trailer`: the trailer of the model.
"""
struct TableConstraint <: Constraint
    scope::Vector{<:AbstractIntVar}
    active::StateObject{Bool}
    initialized::StateObject{Bool}
    table::Matrix{Int}
    currentTable::RSparseBitSet{UInt64}
    modifiedVariables::Set{Int}
    unfixedVariables::Set{Int}
    supports::Dict{Pair{Int,Int},BitVector}
    residues::Dict{Pair{Int,Int},Int}

    TableConstraint(scope, active, initialized, table, currentTable, modifiedVariables, unfixedVariables, supports, residues) = new(scope, active, initialized, table, currentTable, modifiedVariables, unfixedVariables, supports, residues)

    function TableConstraint(variables::Vector{<:AbstractIntVar}, table::Matrix{Int}, supports::Dict{Pair{Int,Int},BitVector}, trailer)
        @assert length(variables) == size(table, 1)

        nVariables, nTuples = size(table)

        constraint = new(
            variables,
            StateObject{Bool}(true, trailer),
            StateObject{Bool}(false, trailer),
            table,
            RSparseBitSet{UInt64}(nTuples, trailer),
            Set{Int}(),
            Set(collect(1:nVariables)),
            supports,
            buildResidues(supports)
        )

        for variable in variables
            addOnDomainChange!(variable, constraint)
        end

        return constraint
    end

end

"""
    TableConstraint(variables, table, trailer)

Create a CompactTable constraint from the `variables`, with values given in `table`.

This is the recommended constructor, as it safely builds the constraint and isolates it. If you choose any of the 2 other constructor,
beware how you create your attributes, and be sure not to modify them at any time.

# Arguments
- `variables::Vector{<:AbstractIntVar}`: vector of variables of size (n, ).
- `table::Matrix{Int}`: matrix of the constraint of size (n, m).
- `trailer::Trailer`: the trailer of the model.
"""
function TableConstraint(variables::Vector{<:AbstractIntVar}, table::Matrix{Int}, trailer)
    @assert length(variables) == size(table, 1)
    cleanedTable = cleanTable(variables, table)

    active = StateObject{Bool}(true, trailer)
    initialized = StateObject{Bool}(false, trailer)
    nVariables, nTuples = size(cleanedTable)
    currentTable = RSparseBitSet{UInt64}(nTuples, trailer)
    modifiedVariables = Set{Int}()
    unfixedVariables = Set(collect(1:nVariables))

    supports = buildSupport(variables, cleanedTable)
    cleanSupports!(supports, variables)

    residues = buildResidues(supports)
    constraint = TableConstraint(
        variables,
        active,
        initialized,
        cleanedTable,
        currentTable,
        modifiedVariables,
        unfixedVariables,
        supports,
        residues
    )

    for variable in variables
        addOnDomainChange!(variable, constraint)
    end

    return constraint
end

"""
    bitVectorToUInt64Vector(bitset)::Vector{UInt64}

Convert a Julia BitVector to a vector preformatted for the RSparseBitSet{UInt64}.

# Arguments
- `bitset::BitVector`: the BitVector to convert.
"""
function bitVectorToUInt64Vector(bitset::BitVector)::Vector{UInt64}
    return [bitreverse(chunk) for chunk in bitset.chunks]
end

"""
    cleanTable(variables, table)::Matrix{Int}

Return a table without the columns corresponding to impossible assignments.

# Arguments
- `variables::Vector{<:AbstractIntVar}`: the variables of the constraint.
- `table::Matrix{Int}`: the table of the constraint.
"""
function cleanTable(variables::Vector{<:AbstractIntVar}, table::Matrix{Int})::Matrix{Int}
    tuplesToDrop = Set{Int}()
    for (idx, col) in enumerate(eachcol(table))
        isValid = all([value in variable.domain for (variable, value) in zip(variables, col)])
        if !isValid
            push!(tuplesToDrop, idx)
        end
    end
    return table[:, setdiff(1:size(table, 2), tuplesToDrop)]
end

"""
    buildSupport(variables, table)::Dict{Pair{Int,Int},BitVector}

Return the support dictionnary, giving for each pair wich columns of the table uses it.

# Arguments
- `variables::Vector{<:AbstractIntVar}`: the variables of the constraint.
- `table::Matrix{Int}`: the table of the constraint.
"""
function buildSupport(variables::Vector{<:AbstractIntVar}, table::Matrix{Int})::Dict{Pair{Int,Int},BitVector}
    nVariables, nTuples = size(table)

    support = Dict{Pair{Int,Int},BitVector}()
    for (idx, var) in enumerate(variables), a in var.domain
        support[idx => a] = BitVector(fill(false, nTuples))
    end

    for (idx, col) in enumerate(eachcol(table)), (variable, value) in enumerate(col)
        support[variable => value][idx] = true
    end

    return support
end

"""
    cleanSupports!(supports, variables)

Remove the pair (variable => value) from `supports` and from the variables domain, when they aren't present in the table.

# Arguments
- `supports::Dict{Pair{Int,Int},BitVector}`: the previously generated supports of the constraint.
- `variables::Vector{<:AbstractIntVar}`: the variables of the constraint.
"""
function cleanSupports!(supports::Dict{Pair{Int,Int},BitVector}, variables::Vector{<:AbstractIntVar})::Nothing
    for ((variable, value), support) in supports
        if !any(support)
            delete!(supports, variable => value)
        end
    end
    return
end

"""
    buildResidues(variables, supports)::Dict{Pair{Int,Int},Int}

Return the residues from the variables and the cleaned supports.

# Arguments
- `variables::Vector{<:AbstractIntVar}`: the variables of the constraint.
- `supports::Dict{Pair{Int,Int},BitVector}`: the previously generated supports of the constraint.
"""
function buildResidues(supports::Dict{Pair{Int,Int},BitVector})::Dict{Pair{Int,Int},Int}
    residues = Dict{Pair{Int,Int},Int}()
    n = 64
    for (key, support) in supports
        residues[key] = Int(ceil(findfirst(support)/n))
    end
    return residues
end

"""
    initialPruning!(constraint, prunedValues)

Store all the assignments without support in `prunedValues` for them to be pruned.
"""
function initialPruning!(constraint::TableConstraint, prunedValues::Vector{Vector{Int}})
    for (idx, variable) in enumerate(constraint.scope), value in variable.domain
        if !((idx => value) in keys(constraint.supports))
            push!(prunedValues[idx], value)
        end
    end
    return
end

"""
    updateTable!(constraint, prunedValues)

Remove the outdated columns from the reversible table of the constraint.

This function is directly inspired by Demeulenaere J. et al. paper.
"""
function updateTable!(constraint::TableConstraint, prunedValues::Vector{Vector{Int}})::Bool
    for variable in constraint.modifiedVariables
        clearMask!(constraint.currentTable)
        if length(prunedValues[variable]) < length(constraint.scope[variable].domain)
            for value in prunedValues[variable]
                support = bitVectorToUInt64Vector(constraint.supports[variable => value])
                addToMask!(constraint.currentTable, support)
            end
            reverseMask!(constraint.currentTable)
        else
            for value in constraint.scope[variable].domain
                support = bitVectorToUInt64Vector(constraint.supports[variable => value])
                addToMask!(constraint.currentTable, support)
            end
        end
        intersectWithMask!(constraint.currentTable)
        if isempty(constraint.currentTable)
            return false
        end
    end
    return true
end

"""
    pruneDomains!(constraint, prunedValues)

Remove the impossible values from the variables domain update prunedValues.

This function is directly inspired by Demeulenaere J. et al. paper.
"""
function pruneDomains!(constraint::TableConstraint)::Vector{Vector{Int}}
    prunedValues = Vector{Vector{Int}}(undef, length(constraint.scope))
    for i = 1:length(constraint.scope)
        prunedValues[i] = Int[]
    end

    for variable in constraint.unfixedVariables, value in constraint.scope[variable].domain
        index = constraint.residues[variable => value]
        support = bitVectorToUInt64Vector(constraint.supports[variable => value])
        if constraint.currentTable.words[index].value & support[index] == UInt64(0)
            index = intersectIndex(constraint.currentTable, support)
            if index != -1
                constraint.residues[variable => value] = index
            else
                push!(prunedValues[variable], value)
            end
        end
    end
    return prunedValues
end

"""
    propagate!(constraint, toPropagate::Set{Constraint}, prunedDomains::CPModification)

`TableConstraint` propagation function. Implement the full procedure of the paper.
"""
function propagate!(constraint::TableConstraint, toPropagate::Set{Constraint}, prunedDomains::CPModification)
    if !constraint.active.value
        return true
    end

    initialPrunedValues = nothing
    if !constraint.initialized.value
        initialPrunedValues = Vector{Vector{Int}}(undef, length(constraint.scope))
        for i = 1:length(constraint.scope)
            initialPrunedValues[i] = Int[]
        end
        initialPruning!(constraint, initialPrunedValues)
        setValue!(constraint.initialized, true)
        for (variable, toRemove) in enumerate(initialPrunedValues), v in toRemove
            remove!(constraint.scope[variable].domain, v)
        end
    end

    # Store all the changes concerning the variables in `constraint.scope`
    prunedValues = Vector{Vector{Int}}(undef, length(constraint.scope))
    empty!(constraint.modifiedVariables)
    for (idx, variable) in enumerate(constraint.scope)
        if haskey(prunedDomains, variable.id)
            push!(constraint.modifiedVariables, idx)
            prunedValues[idx] = copy(prunedDomains[variable.id])
        end
    end

    # Sometime an assignment has already been pruned before the initialization eventhough
    # it doesn't have a support .In that case it shouldn't be used to update the table.
    for variable in constraint.modifiedVariables
        filter!(prunedValues[variable]) do value
            haskey(constraint.supports, variable => value)
        end
    end

    # Update unfixed variables
    empty!(constraint.unfixedVariables)
    union!(constraint.unfixedVariables, findall(var -> length(var.domain) > 1, constraint.scope))

    # Remove impossible supports from the table
    if !updateTable!(constraint, prunedValues)
        return false
    end

    # Restrict domains to the new supports
    newPrunedValues = pruneDomains!(constraint)
    if !isnothing(initialPrunedValues)
        for i = 1:length(constraint.scope)
            append!(newPrunedValues[i], initialPrunedValues[i])
        end
    end
    for (prunedVar, var) in zip(newPrunedValues, constraint.scope)
        if !isempty(prunedVar)
            for val in prunedVar
                remove!(var.domain, val)
            end
            triggerDomainChange!(toPropagate, var)
            addToPrunedDomains!(prunedDomains, var, prunedVar)
        end
    end

    if constraint in toPropagate
        pop!(toPropagate, constraint)
    end

    return true
end

function Base.show(io::IO, ::MIME"text/plain", con::TableConstraint)
    table = split(repr("text/plain", con.table), '\n')[2:end]
    maxlen = Base.maximum(x -> length(x.id), con.scope)

    println(io, string(typeof(con)), ": active = ", con.active)
    for i = 1:length(con.scope)
        println(io, "   ", rpad(con.scope[i].id, maxlen),"  in  ", table[i])
    end
    print(io, "\n   With domains:")
    for var in con.scope
        print(io, "\n      ", var)
    end
end

function Base.show(io::IO, con::TableConstraint)
    print(io, string(typeof(con)), ": (", join([x.id for x in con.scope], ", "), ") in ", con.table, ", active = ", con.active)
end

variablesArray(constraint::TableConstraint) = constraint.scope


text:
using Distributions, SpecialFunctions, Roots, LinearAlgebra
# import SpecialFunctions: lbeta
"""
GNas12: Modified algorithm of Nassar et al., Nat Neuro, 2012 for Gaussian prior
            for the Gaussian task (section 2.3.1)
"""
struct GNas12 <: GLearner
    χ_0::Float64
    ν_0::Float64
    mu_0::Float64
    sigma_0::Float64 # std of observations
    sigma::Float64 # std of observations
    changeprobability::Float64
    m::Float64
    ρ::Float64
    χ_n::Array{Float64, 1}
    ν_n::Array{Float64, 1}
    r_n::Array{Float64, 1}
    var_n::Array{Float64, 1}
    mu_n::Array{Float64, 1} #parameter we are interested in
    Sgm_n::Array{Float64, 1} # Bayes Factor Surprise
    γ_n::Array{Float64, 1}
    Ssh_n::Array{Float64, 1} # Shannon Surprise
    Py_0_n::Array{Float64, 1}
end
function GNas12(; mu_0 = 0., sigma_0 = 1., sigma = 0.1, changeprobability = .001)
    ρ = sigma^2 / sigma_0^2
    m = changeprobability/ (1. - changeprobability)
    ν_0 = calcν(sigma, sigma_0)
    χ_0 = calcχ(mu_0, ν_0, sigma)
    var_n_0 = sigma^2 / ν_0
    mu_n_0 = sigma * χ_0 / ν_0
    GNas12(χ_0, ν_0, mu_0, sigma_0, sigma, changeprobability, m, ρ, [χ_0], [ν_0],
                [1.], [var_n_0], [mu_n_0], [1.], [m/(1. + m)], [1.], [1.])
end
export GNas12
function update!(learnerG::Union{GNas12, GNas12original}, y)
    calcSgm!(learnerG, y)
    calcγ!(learnerG) # Ω = learnerG.m * learnerG.Sgm_n[1]/(1. + learnerG.m * learnerG.Sgm_n[1])
    calcShannon!(learnerG, y)
    Δrule_mu_n = learnerG.mu_n[1] + (y - learnerG.mu_n[1])/(learnerG.ρ + learnerG.r_n[1] + 1.)
    Δrule_mu_0 = learnerG.mu_0 + (y - learnerG.mu_0)/(learnerG.ρ + 1.)
    computemu_n!(learnerG, Δrule_mu_n, Δrule_mu_0)
    computevar_n!(learnerG, Δrule_mu_n, Δrule_mu_0)
    computer_n!(learnerG)
    transformto_χ_ν!(learnerG)
    # learnerG.Sgm_n[1] = copy(Sgm)
end
export update!
function computevar_n!(learnerG::Union{GNas12, GNas12original},
                    Δrule_mu_n, Δrule_mu_0)
    α = (learnerG.ρ + learnerG.γ_n[1]*learnerG.r_n[1] + 1.) / (learnerG.ρ + learnerG.r_n[1] + 1.)
    learnerG.var_n[1] = learnerG.sigma^2 * α / (learnerG.ρ + 1.)
    learnerG.var_n[1] += (1. - learnerG.γ_n[1]) * learnerG.γ_n[1] * (Δrule_mu_n - Δrule_mu_0)^2
end
function computer_n!(learnerG::Union{GNas12, GNas12original})
    learnerG.r_n[1] = learnerG.sigma^2/learnerG.var_n[1] - learnerG.sigma^2/learnerG.sigma_0^2
end


text:
using Atmosphere
using Test

@testset "Temperature and Pressure: Fit" begin

	#---Parameters, Inputs, etc.
	h = [0.0; 5004.0; 11019.0; 15035.0; 20063.0; 26107.0; 32162.0; 39241.0; 47350.0]#[0.0; 5.0; 11.0; 15.0; 20.0; 26.0; 32.0; 39.0; 47.0]*1e3

	#values from the 1976 US Standard Atmosphere Data https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19770009539.pdf, Table 1: Geopotential Altitdue, H
	Tcheck = [288.150; 255.650; 216.650; 216.650; 216.650; 222.650; 228.650; 248.250; 270.650] #K
	Pcheck = [101325.0; 54019.0; 22632.0; 12044.0; 5474.8; 2153.0; 868.01; 318.22; 110.90] #Pa

	#---Run Function/Tests
	for i=1:length(h)
		T,P = Atmosphere.temp_presfit(h[i])
		@test isapprox(T, Tcheck[i], atol=3.75)
		@test isapprox(P, Pcheck[i], atol=1.0825e2)
	end #for length h

end #Temp and Pres: Fit

@testset "Temperature and Pressure: Table" begin

	#---Parameters, Inputs, etc.
	h = [0.0; 5004.0; 11019.0; 15035.0; 20063.0; 26107.0; 32162.0; 39241.0; 47350.0; 49381.0; 51413.0; 61591.0; 71802.0; 78969.0; 85638.0]
	#values from the 1976 US Standard Atmosphere Data https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19770009539.pdf, Table 1: Geopotential Altitdue
	Tcheck = [288.150; 255.650; 216.650; 216.650; 216.650; 222.650; 228.650; 248.250; 270.650; 270.650; 270.650; 242.650; 214.650; 200.650; 187.650] #K
	Pcheck = [101325.0; 54019.0; 22632.0; 12044.0; 5474.8; 2153.0; 868.01; 318.22; 110.90; 86.162; 66.938; 17.660; 3.9564; 1.2501; 0.39814] #Pa

	#---Run Function/Tests
	T = zeros(length(h))
	P = zeros(length(h))
	for i=1:length(h)
		T, P = Atmosphere.temp_prestable(h[i])
		@test isapprox(T, Tcheck[i], atol=3.75e-3)
		@test isapprox(P, Pcheck[i], atol=1.155e1)
	end #for length h

end #Temp and Pres: Table

@testset "Ideal Gas: Density" begin

	#---Parameters, Inputs, etc.
	R = 8.31432e3 #N m kmol−1 K−1 gas constant used in standard atmosphere tables
	M = 28.9645 #g/mol mean molar mass for dry air
	Rs = R/M #Rspecific
	T = [288.150; 216.650; 216.650; 228.650; 270.650; 270.650; 214.650; 187.650] #K
	P = [101325.0; 22632.0; 5474.8; 868.01; 110.90; 66.938; 3.9564; .39814] #Pa

	#---Run Function/Tests
	rho = zeros(length(P))
	for i=1:length(P)
		rho = Atmosphere.density(T[i],P[i])
		@test isapprox(rho, P[i]/(Rs*T[i]), atol=1e-10)
	end #for length P

end #Ideal Gas Density

@testset "Ideal Gas: Speed of Sound" begin

	#---Parameters, Inputs, etc.
	gamma = 1.4  #ratio of specific heats (air)
	R = 8.31432e3 #N m kmol−1 K−1 gas constant used in standard atmosphere tables
	M = 28.9645 #g/mol mean molar mass for dry air
	Rs = R/M #Rspecific
	T = [288.150; 216.650; 216.650; 228.650; 270.650; 270.650; 214.650; 187.650] #K

	#---Initialize and Run Function
	sos = zeros(length(T))
	for i=1:length(T)
		sos = Atmosphere.speedofsound(T[i])
		@test isapprox(sos, sqrt(gamma*Rs*T[i]), atol=1e-10)
	end #for length T

end #Ideal Gas: SoS

@testset "Sutherland's Equation: Dynamic Viscosity" begin

	#---Parameters, Inputs, etc.
	Tsl = 288.15 #sea level temperature (K)
	musl = 0.0000181206 #sea level viscosity (N-s/m^2)
	Sc = 113.0 #constant in Sutherlands formula for air from Sutherland, W. (1893), "The viscosity of gases and molecular force"
	T = [288.150; 216.650; 216.650; 228.650; 270.650; 270.650; 214.650; 187.650] #K

	#---Initialize and Run Function
	for i=1:length(T)
		mu = Atmosphere.viscosity(T[i])
		@test isapprox(mu, musl*(T[i]/Tsl)^(3.0/2.0)*(Tsl+Sc)/(T[i]+Sc), atol=3.25e-8)
	end

end #Southerland's Eqn.: Viscosity


text:
using Test
using BinaryBuilderBase
using BinaryBuilderBase: platform_dlext, platform_exeext

@testset "Wrappers utilities" begin
    @test nbits(Platform("i686", "linux")) == 32
    @test nbits(Platform("x86_64", "linux"; march="avx")) == 64
    @test nbits(Platform("armv7l", "linux")) == 32
    @test nbits(Platform("aarch64", "linux"; cuda="10.1")) == 64
    @test nbits(Platform("powerpc64le", "linux")) == 64
    @test nbits(AnyPlatform()) == 64

    @test proc_family(Platform("i686", "linux")) == "intel"
    @test proc_family(Platform("x86_64", "linux"; march="avx")) == "intel"
    @test proc_family(Platform("armv7l", "linux")) == "arm"
    @test proc_family(Platform("aarch64", "linux"; cuda="10.1")) == "arm"
    @test proc_family(Platform("powerpc64le", "linux")) == "power"
    @test proc_family(AnyPlatform()) == "any"

    @test platform_exeext(Platform("i686", "linux")) == ""
    @test platform_exeext(Platform("x86_64", "freebsd")) == ""
    @test platform_exeext(Platform("x86_64", "macos")) == ""
    @test platform_exeext(Platform("i686", "windows")) == ".exe"
    @test platform_exeext(Platform("x86_64", "linux"; march="avx512")) == ""

    @test aatriplet(Platform("x86_64", "linux"; libc="glibc", libgfortran_version=v"4.0.0", march="avx", cuda="9.2")) == "x86_64-linux-gnu"
end

# Are we using docker? If so, test that the docker runner works...
@testset "Runner utilities" begin
    # Test that is_ecryptfs works for something we're certain isn't encrypted
    if isdir("/proc")
        isecfs = (false, "/proc/")
        @test_logs (:info, "Checking to see if /proc/ is encrypted...") @test BinaryBuilderBase.is_ecryptfs("/proc"; verbose=true) == isecfs
        @test_logs (:info, "Checking to see if /proc/ is encrypted...") @test BinaryBuilderBase.is_ecryptfs("/proc/"; verbose=true) == isecfs
        @test_logs (:info, "Checking to see if /proc/not_a_file is encrypted...") @test BinaryBuilderBase.is_ecryptfs("/proc/not_a_file"; verbose=true) == isecfs
    else
        @test_logs (:info, "Checking to see if /proc/ is encrypted...") @test BinaryBuilderBase.is_ecryptfs("/proc"; verbose=true) == (false, "/proc")
        @test_logs (:info, "Checking to see if /proc/ is encrypted...") @test BinaryBuilderBase.is_ecryptfs("/proc/"; verbose=true) == (false, "/proc/")
        @test_logs (:info, "Checking to see if /proc/not_a_file is encrypted...") @test BinaryBuilderBase.is_ecryptfs("/proc/not_a_file"; verbose=true) == (false, "/proc/not_a_file")
    end

    if isa(preferred_runner(), BinaryBuilderBase.DockerRunner)
        @testset "Docker image importing" begin
            # First, delete the docker image, in case it already existed
            BinaryBuilderBase.delete_docker_image()

            # Next, import it and ensure that doesn't throw
            rootfs = first(BinaryBuilderBase.choose_shards(platform))
            mktempdir() do dir
                @test BinaryBuilderBase.import_docker_image(rootfs, dir; verbose=true) === nothing
            end

            # Test that deleting the docker image suceeds, now that we know
            # it exists
            @test BinaryBuilderBase.delete_docker_image()
        end
    end

    @testset "hello world" begin
        mktempdir() do dir
            ur = preferred_runner()(dir; platform=Platform("x86_64", "linux"; libc="musl"))
            iobuff = IOBuffer()
            @test run(ur, `/bin/bash -c "echo test"`, iobuff)
            seek(iobuff, 0)
            # Test that we get the output we expect (e.g. the second line is `test`)
            @test split(String(read(iobuff)), "\n")[2] == "test"
        end
    end

    @testset "Compilation" begin
        mktempdir() do dir
            ur = preferred_runner()(dir; platform=Platform("x86_64", "linux"; libc="musl"))
            iobuff = IOBuffer()
            @test run(ur, `/bin/bash -c "echo 'int main() {return 0;}' | cc -x c -"`, iobuff; tee_stream=devnull)
            seekstart(iobuff)
            @test split(String(read(iobuff)), "\n")[2] == ""
        end
    end

    @testset "Locking microarchitecture" begin
        mktempdir() do dir
            platform = Platform("x86_64", "linux"; libc="musl")
            cmd = `/bin/bash -c "echo 'int main() {return 0;}' | cc -x c -march=native -"`
            ur = preferred_runner()(dir; platform=platform, lock_microarchitecture=true)
            iobuff = IOBuffer()
            @test !run(ur, cmd, iobuff; tee_stream=devnull)
            seekstart(iobuff)
            @test split(String(read(iobuff)), "\n")[2] == "Cannot force an architecture"

            ur = preferred_runner()(dir; platform=platform, lock_microarchitecture=false)
            iobuff = IOBuffer()
            @test run(ur, cmd, iobuff; tee_stream=devnull)
            seekstart(iobuff)
            @test split(String(read(iobuff)), "\n")[2] == ""
        end
    end

    @testset "Unsafe flags" begin
        mktempdir() do dir
            platform = Platform("x86_64", "linux"; libc="musl")
            cmd = `/bin/bash -c "echo 'int main() {return 0;}' | cc -x c -Ofast -"`
            ur = preferred_runner()(dir; platform=platform, allow_unsafe_flags=false)
            iobuff = IOBuffer()
            @test !run(ur, cmd, iobuff; tee_stream=devnull)
            seekstart(iobuff)
            lines = split(String(read(iobuff)), "\n")
            @test lines[2] == "You used one or more of the unsafe flags: -Ofast, -ffast-math, -funsafe-math-optimizations"
            @test lines[3] == "Please repent."

            ur = preferred_runner()(dir; platform=platform, allow_unsafe_flags=true)
            iobuff = IOBuffer()
            @test run(ur, cmd, iobuff; tee_stream=devnull)
            seekstart(iobuff)
            @test split(String(read(iobuff)), "\n")[2] == ""
        end
    end
end