text:
using JudiLing
using CSV
using Test
using DataFrames

@testset "make prelinguistic semantic matrix for utterance" begin
    utterance = DataFrame(CSV.File(joinpath(
        @__DIR__,
        "data",
        "utterance_mini.csv",
    )))
    s_obj_train = JudiLing.make_pS_matrix(utterance)

    utterance_val = utterance[101:end, :]
    s_obj_val = JudiLing.make_pS_matrix(utterance_val, s_obj_train)
end

@testset "make semantic matrix" begin
    data = DataFrame(
        X = ["A", "B", "C", "D", "E"],
        Y = ["M", "F", "F", "M", missing],
        Z = ["P", "P", "S", missing, "P"],
    )

    n_features = 3
    seed = 314

    L = JudiLing.make_L_matrix(
        data[1:5, :],
        ["Y"],
        ["Z"],
        ncol = n_features,
        seed = seed,
        isdeep = true,
    )

    S_train, S_val = JudiLing.make_S_matrix(
        data[1:2, :],
        data[3:5, :],
        ["Y"],
        ["Z"],
        L,
        add_noise = false,
    )

    @test S_train[1, :] == L.L[L.f2i["M"], :] + L.L[L.f2i["P"], :]
    @test S_train[2, :] == L.L[L.f2i["F"], :] + L.L[L.f2i["P"], :]
    @test S_val[1, :] == L.L[L.f2i["F"], :] + L.L[L.f2i["S"], :]
    @test S_val[2, :] == L.L[L.f2i["M"], :]
    @test S_val[3, :] == L.L[L.f2i["P"], :]

    S_train, S_val = JudiLing.make_S_matrix(
        data[1:2, :],
        data[3:5, :],
        ["Y"],
        L,
        add_noise = false,
    )

    @test S_train[1, :] == L.L[L.f2i["M"], :]
    @test S_train[2, :] == L.L[L.f2i["F"], :]
    @test S_val[1, :] == L.L[L.f2i["F"], :]
    @test S_val[2, :] == L.L[L.f2i["M"], :]
    @test S_val[3, :] == zeros(Float64, 3)

    S_train = JudiLing.make_S_matrix(
        data[1:5, :],
        ["Y"],
        ["Z"],
        L,
        add_noise = false,
    )

    @test S_train[1, :] == L.L[L.f2i["M"], :] + L.L[L.f2i["P"], :]
    @test S_train[2, :] == L.L[L.f2i["F"], :] + L.L[L.f2i["P"], :]
    @test S_train[3, :] == L.L[L.f2i["F"], :] + L.L[L.f2i["S"], :]
    @test S_train[4, :] == L.L[L.f2i["M"], :]
    @test S_train[5, :] == L.L[L.f2i["P"], :]

    S_train =
        JudiLing.make_S_matrix(data[1:5, :], ["Y"], L, add_noise = false)

    @test S_train[1, :] == L.L[L.f2i["M"], :]
    @test S_train[2, :] == L.L[L.f2i["F"], :]
    @test S_train[3, :] == L.L[L.f2i["F"], :]
    @test S_train[4, :] == L.L[L.f2i["M"], :]
    @test S_train[5, :] == zeros(Float64, 3)

    S_train, S_val = JudiLing.make_S_matrix(
        data[1:5, :],
        data[3:5, :],
        ["Y"],
        ["Z"],
        ncol = n_features,
        seed = seed,
        isdeep = true,
        add_noise = false,
    )

    @test S_train[1, :] == L.L[L.f2i["M"], :] + L.L[L.f2i["P"], :]
    @test S_train[2, :] == L.L[L.f2i["F"], :] + L.L[L.f2i["P"], :]
    @test S_val[1, :] == L.L[L.f2i["F"], :] + L.L[L.f2i["S"], :]
    @test S_val[2, :] == L.L[L.f2i["M"], :]
    @test S_val[3, :] == L.L[L.f2i["P"], :]

    S_train = JudiLing.make_S_matrix(
        data[1:5, :],
        ["Y"],
        ["Z"],
        ncol = n_features,
        seed = seed,
        isdeep = true,
        add_noise = false,
    )

    S_train, S_val = JudiLing.make_S_matrix(
        data[1:5, :],
        data[3:5, :],
        ["Y"],
        ncol = n_features,
        seed = seed,
        isdeep = true,
        add_noise = false,
    )

    S_train = JudiLing.make_S_matrix(
        data[1:5, :],
        ["Y"],
        ncol = n_features,
        seed = seed,
        isdeep = true,
        add_noise = false,
    )
end

@testset "make combined semantic matrix" begin
    data = DataFrame(
        X = ["A", "B", "C", "D", "E"],
        Y = ["M", "F", "F", "M", missing],
        Z = ["P", "P", "S", missing, "P"],
    )

    n_features = 3
    seed = 314

    L = JudiLing.make_L_matrix(
        data[1:5, :],
        ["Y"],
        ["Z"],
        ncol = n_features,
        seed = seed,
        isdeep = true,
    )

    S_train, S_val = JudiLing.make_combined_S_matrix(
        data[1:2, :],
        data[3:5, :],
        ["Y"],
        ["Z"],
        ncol = n_features,
        seed = seed,
        isdeep = true,
        add_noise = false,
    )

    @test S_train[1, :] == L.L[L.f2i["M"], :] + L.L[L.f2i["P"], :]
    @test S_train[2, :] == L.L[L.f2i["F"], :] + L.L[L.f2i["P"], :]
    @test S_val[1, :] == L.L[L.f2i["F"], :] + L.L[L.f2i["S"], :]
    @test S_val[2, :] == L.L[L.f2i["M"], :]
    @test S_val[3, :] == L.L[L.f2i["P"], :]

    S_train, S_val = JudiLing.make_combined_S_matrix(
        data[1:2, :],
        data[3:5, :],
        ["Y"],
        ncol = n_features,
        seed = seed,
        isdeep = true,
        add_noise = false,
    )

    @test S_train[1, :] == L.L[L.f2i["M"], :]
    @test S_train[2, :] == L.L[L.f2i["F"], :]
    @test S_val[1, :] == L.L[L.f2i["F"], :]
    @test S_val[2, :] == L.L[L.f2i["M"], :]
    @test S_val[3, :] == zeros(Float64, n_features)
end

@testset "make L matrix" begin
    data = DataFrame(
        X = ["A", "B", "C", "D", "E"],
        Y = ["M", "F", "F", "M", missing],
        Z = ["P", "P", "S", missing, "P"],
    )

    n_features = 3
    seed = 314

    L1 = JudiLing.make_L_matrix(
        data[1:5, :],
        ["Y"],
        ["Z"],
        ncol = n_features,
        seed = seed,
        isdeep = true,
    )

    @test size(L1.L) == (4, 3)
    @test L1.f2i == Dict("S" => 4, "M" => 1, "P" => 3, "F" => 2)
    @test L1.i2f == ["M", "F", "P", "S"]
    @test L1.ncol == 3

    L2 = JudiLing.make_L_matrix(
        data[1:5, :],
        ["Y"],
        ncol = n_features,
        seed = seed,
        isdeep = true,
    )

    @test size(L2.L) == (2, 3)
    @test L2.f2i == Dict("M" => 1, "F" => 2)
    @test L2.i2f == ["M", "F"]
    @test L2.ncol == 3

    L3 = JudiLing.make_L_matrix(
        data[1:5, :],
        ["Y"],
        ["Z"],
        ncol = n_features,
        seed = seed,
        isdeep = false,
    )

    @test size(L3.L) == (4, 3)
    @test L3.f2i == Dict("S" => 4, "M" => 1, "P" => 3, "F" => 2)
    @test L3.i2f == ["M", "F", "P", "S"]
    @test L3.ncol == 3

    L4 = JudiLing.make_L_matrix(
        data[1:5, :],
        ["Y"],
        ncol = n_features,
        seed = seed,
        isdeep = false,
    )

    @test size(L4.L) == (2, 3)
    @test L4.f2i == Dict("M" => 1, "F" => 2)
    @test L4.i2f == ["M", "F"]
    @test L4.ncol == 3

    L5 = JudiLing.make_combined_L_matrix(
        data[1:2, :],
        data[3:5, :],
        ["Y"],
        ["Z"],
        ncol = n_features,
        seed = seed,
        isdeep = false,
    )

    @test size(L5.L) == (4, 3)
    @test L5.f2i == Dict("S" => 4, "M" => 1, "P" => 3, "F" => 2)
    @test L5.i2f == ["M", "F", "P", "S"]
    @test L5.ncol == 3

    L6 = JudiLing.make_combined_L_matrix(
        data[1:2, :],
        data[3:5, :],
        ["Y"],
        ncol = n_features,
        seed = seed,
        isdeep = true,
    )

    @test size(L6.L) == (2, 3)
    @test L6.f2i == Dict("M" => 1, "F" => 2)
    @test L6.i2f == ["M", "F"]
    @test L6.ncol == 3
end


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

name = "qr_mumps"
version = v"3.0.1"

# Collection of sources required to complete build
sources = [
    ArchiveSource("http://buttari.perso.enseeiht.fr/qr_mumps/releases/qr_mumps-$version.tgz",
                  "451d3a25beff3c8e8ac3c03249333341b871aa4f55ff37d91837e2397497246d")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd qr_mumps*
mkdir build
cd build
cmake .. -DARITH="d;s;c;z" -DBUILD_SHARED_LIBS=ON -DCMAKE_INSTALL_PREFIX=$prefix -DQRM_ORDERING_AMD=OFF -DQRM_ORDERING_METIS=OFF -DQRM_ORDERING_SCOTCH=OFF -DQRM_WITH_STARPU=OFF -DQRM_WITH_CUDA=OFF -DCMAKE_BUILD_TYPE=Release -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN}
make -j${nproc}
make install
exit
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_gfortran_versions(supported_platforms())
filter!(p -> libgfortran_version(p) != v"3", platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libqrm_common", :libqrm_common),
    LibraryProduct("libcqrm", :libcqrm),
    LibraryProduct("libsqrm", :libsqrm),
    LibraryProduct("libdqrm", :libdqrm),
    LibraryProduct("libzqrm", :libzqrm)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="OpenBLAS32_jll", uuid="656ef2d0-ae68-5445-9ca0-591084a874a2"))
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae"))
]

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

ext: jl

max_stars_repo_name: bcoconni/Yggdrasil

max_issues_repo_name: bcoconni/Yggdrasil

max_forks_repo_name: bcoconni/Yggdrasil


text:
module IPPCore

    const libippcore = "libippcore"

    export  IPPInt,     ippint,
            IPP8u,      ipp8u,
            IPP8s,      ipp8s,
            IPP16u,     ipp16u,
            IPP16s,     ipp16s,
            IPP32u,     ipp32u,
            IPP32s,     ipp32s,
            IPP16f,     ipp16f,
            IPP32f,     ipp32f,
            IPP64f,     ipp64f,
            IPP16uc,    ipp16uc,
            IPP16sc,    ipp16sc,
            IPP32uc,    ipp32uc,
            IPP32sc,    ipp32sc,
            IPP16fc,    ipp16fc,
            IPP32fc,    ipp32fc,
            IPP64fc,    ipp64fc,
            IPPSuffix,  IPPTypeSignature,
            IppVersion, ippversion,
            IppStatus,  ippstatus_string,
            ippcputype, ippcpudescr, ippcpufreq

    # common types

    typealias IPPInt    Cint
    typealias IPP8u     Uint8
    typealias IPP8s     Int8
    typealias IPP16u    Uint16
    typealias IPP16s    Int16
    typealias IPP32u    Uint32
    typealias IPP32s    Int32
    typealias IPP16f    Float16
    typealias IPP32f    Float32
    typealias IPP64f    Float64
    typealias IPP16uc   Complex{Uint16}
    typealias IPP16sc   Complex{Int16}
    typealias IPP32uc   Complex{Uint32}
    typealias IPP32sc   Complex{Int32}
    typealias IPP16fc   Complex{Float16}
    typealias IPP32fc   Complex{Float32}
    typealias IPP64fc   Complex{Float64}

    ippint( x )  = convert( IPPInt, x )
    ipp8u( x )   = convert( IPP8u, x )
    ipp8s( x )   = convert( IPP8s, x )
    ipp16u( x )  = convert( IPP16u, x )
    ipp16s( x )  = convert( IPP16s, x )
    ipp32u( x )  = convert( IPP32u, x )
    ipp32s( x )  = convert( IPP32s, x )
    ipp16f( x )  = convert( IPP16f, x )
    ipp32f( x )  = convert( IPP32f, x )
    ipp64f( x )  = convert( IPP64f, x )
    ipp16uc( x ) = convert( IPP16uc, x )
    ipp16sc( x ) = convert( IPP16sc, x )
    ipp32uc( x ) = convert( IPP32uc, x )
    ipp32sc( x ) = convert( IPP32sc, x )
    ipp16fc( x ) = convert( IPP16fc, x )
    ipp32fc( x ) = convert( IPP32fc, x )
    ipp64fc( x ) = convert( IPP64fc, x )

    ## Get version

    immutable IppVersion
        major::Int
        minor::Int
        majorBuild::Int
        build::Int
    end

    function Base.show(io::IO, v::IppVersion)
        print(io, "IPP Version $(v.major).$(v.minor).$(v.majorBuild) (build=$(v.build))")
    end

    function ippversion()
        p = ccall((:ippGetLibVersion, libippcore), Ptr{Uint8}, ())
        nbytes = sizeof(IPPInt) * 4 + 4 + sizeof(Ptr{Uint}) * 4
        a = pointer_to_array(p, nbytes)
        buf = IOBuffer(a)

        major = int(read(buf, IPPInt))
        minor = int(read(buf, IPPInt))
        majorBuild = int(read(buf, IPPInt))
        build = int(read(buf, IPPInt))

        return IppVersion(major, minor, majorBuild, build)
    end

    ## Status string

    typealias IppStatus IPPInt
    function ippstatus_string(s::Integer)
        p = ccall((:ippGetStatusString, libippcore), Ptr{Uint8}, (IppStatus,), s)
        return bytestring(p)
    end

    ## CPU info

    ippcputype() = uint8(ccall((:ippGetCpuType, libippcore), IPPInt, ()))

    const cpudescrmap = (Uint8=>ASCIIString)[
        0x00 => "Unknown CPU",
        0x01 => "Intel(R) Pentium(R) processor",
        0x02 => "Intel(R) Pentium(R) processor with MMX(TM) technology",
        0x03 => "Intel(R) Pentium(R) Pro processor",
        0x04 => "Intel(R) Pentium(R) II processor",
        0x05 => "Intel(R) Pentium(R) III processor and Pentium(R) III Xeon(R) processor",
        0x06 => "Intel(R) Pentium(R) 4 processor and Intel(R) Xeon(R) processor",
        0x07 => "Intel(R) Pentium(R) 4 Processor with HT Technology",
        0x08 => "Intel(R) Pentium(R) 4 processor with Intel(R) SSE 3",
        0x09 => "Intel(R) Centrino(R) mobile technology",
        0x0a => "Intel(R) Core(TM) Solo processor",
        0x0b => "Intel(R) Core(TM) Duo processor",
        0x10 => "Intel(R) Itanium(R) processor",
        0x11 => "Intel(R) Itanium(R) 2 processor",
        0x20 => "Intel(R) 64 Instruction Set Architecture (ISA)",
        0x21 => "Intel(R) Core(TM) 2 Duo processor",
        0x22 => "Intel(R) Core(TM) 2 Quad processor",
        0x23 => "Intel(R) Core(TM) 2 processor with Intel(R) SSE 4.1",
        0x24 => "Intel(R) Atom(TM) processor",
        0x25 => "Intel(R) Core(TM) i7 processor",
        0x26 => "Next Generation CPU",
        0x40 => "Processor supports Intel(R) SSE instruction set",
        0x41 => "Processor supports Intel(R) SSE 2 instruction set",
        0x42 => "Processor supports Intel(R) SSE 3 instruction set",
        0x43 => "Processor supports Intel(R) SSSE 3 instruction set",
        0x44 => "Processor supports Intel(R) SSE 4.1 instruction set",
        0x45 => "Processor supports Intel(R) SSE 4.2 instruction set",
        0x46 => "Processor supports Intel(R) AVX instruction set",
        0x47 => "Processor supports Intel(R) AES New Instructions",
        0x49 => "Processor supports RDRRAND & Float16 instructions",
        0x4a => "Processor supports Intel(R) AVX 2 instruction set",
        0x4b => "Processor supports ADCX and ADOX instructions",
        0x60 => "Processor supports 64 bit extension"]

    function ippcpudescr(i::Uint8)
        if haskey(cpudescrmap, i)
            return cpudescrmap[i]
        else
            return @sprintf("Unknown Processor (code = %02x)", i)
        end
    end

    ippcpudescr(i::Integer) = ippcpudescr(uint8(i))
    ippcpudescr() = ippcpudescr(ippcputype())

    function ippcpufreq()  # in terms of Mhz
        r = IPPInt[0]
        ccall((:ippGetCpuFreqMhz, libippcore), IppStatus, (Ptr{Int64},), pointer(r))
        return int(r[1])
    end

end


text:
using DSP, Test

@testset "esprit" begin
    # create a sum of sinusoids in noise, and estimate their frequencies
    Fs = 10000.0           # sampling frequency in Hz
    duration = 1           # length of signal, in seconds
    n = Int(Fs * duration) # number of samples
    t = collect((1:n)/Fs)  # time vector
    frequencies = [1000.0 1250.0]
    amplitudes = [2.0 1.5]
    phases = [0.7 -1.0]
    x = exp.( 1im*2*π*t*frequencies .+ phases) * amplitudes'
    noise = randn(n, 2)*[1;1im]
    sigma = 0.1
    noise *= sigma
    x += noise
    M = 300
    p = 2                  # number of sinusoids to estimate
    frequencies_estimated = sort(esprit(x, M, p, Fs))
    @test isapprox(frequencies', frequencies_estimated; atol = 1e-2)
end


text:
using TextAnalysis, Flux
using LinearAlgebra
# include("../src/ner.jl")

function try_outs(ner, x_in, y_in)
    unique_labels = unique(ner.model.labels)
    num_labels = length(unique_labels)
    confusion_matrix = zeros(Int, (num_labels, num_labels))

    for (x_seq, y_seq) in zip(x_in, y_in)
        preds = ner(x_seq)

        for (pred, logit) in zip(preds, y_seq)
            confusion_matrix[findfirst(x -> x==pred, unique_labels), findfirst(x -> x==logit, unique_labels)] += 1
        end
    end

    # print(confusion_matrix)
    s1 = sum(confusion_matrix, dims=2)
    s2 = sum(confusion_matrix, dims=1)
    dg = diag(confusion_matrix)

    a = sum(dg ./ s1) / num_labels
    b = sum(dg ./ s2') / num_labels

    println("Precision and recall are:", a, " ", b)
    println("F1 is:", (2 * a * b) / (a + b))
end


text:
Base.:(+)(x::XFloat16, y::XFloat16) = reinterpret(XFloat16, (+)(reinterpret(Float32, x), reinterpret(Float32, y)))
Base.:(-)(x::XFloat16, y::XFloat16) = reinterpret(XFloat16, (-)(reinterpret(Float32, x), reinterpret(Float32, y)))
Base.:(*)(x::XFloat16, y::XFloat16) = reinterpret(XFloat16, (*)(reinterpret(Float32, x), reinterpret(Float32, y)))
Base.:(/)(x::XFloat16, y::XFloat16) = reinterpret(XFloat16, (/)(reinterpret(Float32, x), reinterpret(Float32, y)))
Base.:(\)(x::XFloat16, y::XFloat16) = reinterpret(XFloat16, (\)(reinterpret(Float32, x), reinterpret(Float32, y)))

Base.:(+)(x::XFloat32, y::XFloat32) = reinterpret(XFloat32, (+)(reinterpret(Float64, x), reinterpret(Float64, y)))
Base.:(-)(x::XFloat32, y::XFloat32) = reinterpret(XFloat32, (-)(reinterpret(Float64, x), reinterpret(Float64, y)))
Base.:(*)(x::XFloat32, y::XFloat32) = reinterpret(XFloat32, (*)(reinterpret(Float64, x), reinterpret(Float64, y)))
Base.:(/)(x::XFloat32, y::XFloat32) = reinterpret(XFloat32, (/)(reinterpret(Float64, x), reinterpret(Float64, y)))
Base.:(\)(x::XFloat32, y::XFloat32) = reinterpret(XFloat32, (\)(reinterpret(Float64, x), reinterpret(Float64, y)))


text:
include("imports.jl")

@testset ExtendedTestSet "showitems" begin
    image = Image(rand(RGB, 100, 100))
    image2 = Image(rand(RGB, 200, 500))
    kps = Keypoints([SVector(1., 2.), SVector(150., 200.)], getbounds(image))
    poly = Polygon([SVector(1., 2.), SVector(150., 200.)], getbounds(image))
    bb = BoundingBox([SVector(1., 2.), SVector(150., 200.)], getbounds(image))

    items = [image, image2, kps, poly, bb]

    @test_nowarn showitems(items)
    @test_nowarn showitems(items, showbounds = true)

end


text:
"Default maximum depth for trees."
const DEFAULT_MAX_TREE_DEPTH = 10

"""
The Tree data structure manages memory while sampling. In particular, preallocation is tricky while sampling from trees. Exploring a tree of depth N
will require evaluating 2^N positions, and calculating gradients and momentums for each, as well as turning statistics.
The tree allocates enough space for the maximum possible number of live variables (depth + 1), which is much smaller than 2^N for N > 2, and keeps
track of which are still live or have been freed and can thus be reused.

Supports a treedepth of up to 30. Code can be modified to support up to 62, but if you're anywhere near 30 you've probably got severe sampling issues.
The default maximum depth is 10. Runtime is exponential as a function of the depth you're hitting; a depth of 30 will take about a million times longer than 10.

The limit comes from using 32 bits to indicate whether each vector is allocated or not, allowing space for depth+2 vectors.
A 1 bit indicates the corresponding space is available; 0 means occupied.
"""
struct Tree{D,T,L}
    root::Ptr{T}
    depth::Int
    sptr::StackPointer
end
aligned_offset(::Tree{D,T,L}) where {D,T,L} = L*sizeof(T)
# struct FlaggedPhasePoint{D,T,L}
    # z::PhasePoint{D,T,L}
    # flag::UInt32
# end
struct FlaggedVector{D,T,L} <: PaddedMatrices.AbstractMutableFixedSizeVector{D,T,L}
    v::PtrVector{D,T,1,0,0,false}
    flag::UInt32
end
@inline function FlaggedVector{D,T,L}(ptr::Ptr{T}, flag::UInt32) where {D,T,L}
    FlaggedVector{D,T,L}(PtrVector{D,T,1,0,0,false}(ptr), flag)
end
@inline Base.pointer(v::FlaggedVector) = v.v.ptr

function Tree{D,T,L}(sptr::StackPointer, depth::Int = DEFAULT_MAX_TREE_DEPTH) where {D,T,L}
    root = pointer(sptr, T)
    depth₊ = depth + 3
    # set roots to zero so that all loads can properly be interpreted as bools without further processing on future accesses.
    tree = Tree{D,T,L}( root, depth₊, sptr + VectorizationBase.REGISTER_SIZE + 6L*sizeof(T)*depth₊ )
    clear!(tree)
    tree
end
clear!(tree::Tree) = SIMDPirates.vstore!(reinterpret(Ptr{UInt64}, tree.root), (VE(0xffffffffffffffff),VE(0xffffffffffffffff)))
clear_all_but_z!(tree::Tree, flag::UInt32) = SIMDPirates.vstore!(reinterpret(Ptr{UInt32}, tree.root), (VE(0xffffffff ⊻ flag),VE(0xffffffff),VE(0xffffffff),VE(0xffffffff)))

@generated Tree{D,T}(sptr::StackPointer, depth::Int = DEFAULT_MAX_TREE_DEPTH) where {D,T} = :(Tree{$D,$T,$(VectorizationBase.align(D,T))}(sptr, depth))
Tree(sptr::StackPointer, depth::Int, ::PtrVector{D,T,L}) where {D,T,L} = Tree{D,T,L}(sptr, dept)

"""
To allocate, we find the first unoccupied slot (corresponding to a bit of 1), and then set that bit to 0.
"""
function allocate!(root::Ptr)
    root32 = reinterpret(Ptr{UInt32}, root)
    allocations = VectorizationBase.load(root32)
    first_unallocated = leading_zeros(allocations)
    flag = 0x80000000 >> first_unallocated
    VectorizationBase.store!(root32, allocations ⊻ flag)
    first_unallocated, flag
end
function allocate!(tree::Tree, flag::UInt32, i=0)
    root32 = reinterpret(Ptr{UInt32}, tree.root) + i
    allocations = VectorizationBase.load(root32)
    VectorizationBase.store!(root32, allocations ⊻ flag)
end
function isallocated(tree::Tree, flag::UInt32, i=0)
    allocations = VectorizationBase.load(reinterpret(Ptr{UInt32}, tree.root))
    (allocations | flag) != allocations
end
function undefined_z(tree::Tree{D,T,L}) where {D,T,L}
    @unpack root = tree
    LT = aligned_offset(tree)
    stump = root + VectorizationBase.REGISTER_SIZE
    first_unallocated, flag = allocate!(root)
    # @show first_unallocated, flag
    # display(flag)
    # print("Defining z at $first_unallocated with flag "); display(flag)
    # println("Defining z at $first_unallocated with flag $(bitstring(flag))")
    # first_unallocated < tree.depth || println("Warning, results invalid: Defining z at $first_unallocated")
    # first_unallocated < tree.depth || display(stacktrace())
    @assert first_unallocated < tree.depth "Don't have space to allocate z!"
    stump + 3first_unallocated*LT, flag
end
function undefined_ρ♯(tree::Tree{D,T,L}) where {D,T,L}
    @unpack root, depth = tree
    LT = aligned_offset(tree)
    stump = root + VectorizationBase.REGISTER_SIZE + 4LT*depth
    first_unallocated, flag = allocate!(root + 8)
    # @assert first_unallocated < depth
    # first_unallocated < tree.depth || println("Warning, results invalid: Defining ρ♯ at $first_unallocated")
    # println("Defining ρ♯ at $first_unallocated with flag $(bitstring(flag))")
    # first_unallocated < tree.depth || display(stacktrace())
    @assert first_unallocated < (depth<<1) "Don't have space to allocate ρ♯!"
    FlaggedVector{D,T,L}( stump + first_unallocated*LT, flag )
end
function undefined_Σρ(tree::Tree{D,T,L}) where {D,T,L}
    @unpack root, depth = tree
    LT = aligned_offset(tree)
    stump = root + VectorizationBase.REGISTER_SIZE + 3LT*depth
    first_unallocated, flag = allocate!(root + 4)
    # @assert first_unallocated < depth
    # first_unallocated < tree.depth || println("Warning, results invalid: Defining Σρ at $first_unallocated")
    # first_unallocated < tree.depth || display(stacktrace())
    @assert first_unallocated < depth "Don't have space to allocate Σρ!"
    FlaggedVector{D,T,L}( stump + first_unallocated*LT, flag )
end
"""
To free, we simply set the bit back to 1 so that it may be allocated again.
"""
function free!(tree::Tree, flag::UInt32, offset::Int)
    @unpack root = tree
    # println("Freeing flag $(bitstring(flag)) at offset $offset")
    # print("Freeing flag "); display(flag); println(" at offset $offset")
    root32 = reinterpret(Ptr{UInt32}, root) + offset
    allocated = VectorizationBase.load(root32)
    flag != 0x00000000 && @assert (allocated | flag) != allocated "Double free!!!!"
    VectorizationBase.store!(root32, VectorizationBase.load(root32) | flag)
    nothing
end
free_z!(tree::Tree, flag::UInt32) = free!(tree, flag, 0)
free_ρ♯!(tree::Tree, flag::UInt32) = free!(tree, flag, 8)
free_Σρ!(tree::Tree, flag::UInt32) = free!(tree, flag, 4)

#####
##### Abstract tree/trajectory interface
#####

####
#### Directions
####

"Maximum number of iterations [`next_direction`](@ref) supports."
const MAX_DIRECTIONS_DEPTH = 32

"""
Internal type implementing random directions. Draw a new value with `rand`, see
[`next_direction`](@ref).
Serves two purposes: a fixed value of `Directions` is useful for unit testing, and drawing a
single bit flag collection economizes on the RNG cost.
"""
struct Directions
    flags::UInt32
end

Base.rand(rng::AbstractRNG, ::Type{Directions}) = Directions(rand(rng, UInt32))
Base.rand(rng::VectorizedRNG.AbstractPCG, ::Type{Directions}) = Directions(rand(rng, UInt32))

"""
$(SIGNATURES)
Return the next direction flag and the new state of directions. Results are undefined for
more than [`MAX_DIRECTIONS_DEPTH`](@ref) updates.
"""
function next_direction(directions::Directions)
    @unpack flags = directions
    Bool(flags & 0x01), Directions(flags >>> 1)
end

####
#### Trajectory interface
####

"""
    $(FUNCTIONNAME)(trajectory, z, is_forward)
Move along the trajectory in the specified direction. Return the new position.
"""
function move end

"""
    $(FUNCTIONNAME)(trajectory, τ)
Test if the turn statistics `τ` indicate that the corresponding tree is turning.
Will only be called on nontrivial trees (at least two nodes).
"""
function is_turning end

"""
    $(FUNCTIONNAME)(trajectory, τ₁, τ₂)
Combine turn statistics on trajectory. Implementation can assume that the trees that
correspond to the turn statistics have the same ordering.
"""
function combine_turn_statistics end

"""
    $(FUNCTIONNAME)(trajectory, v₁, v₂)
Combine visited node statistics for adjacent trees trajectory. Implementation should be
invariant to the ordering of `v₁` and `v₂` (ie the operation is commutative).
"""
function combine_visited_statistics end

"""
    $(FUNCTIONNAME)(trajectory, is_doubling::Bool, ω₁, ω₂, ω)
Calculate the log probability if selecting the subtree corresponding to `ω₂`. Being the log
of a probability, it is always `≤ 0`, but implementations are allowed to return and accept
values `> 0` and treat them as `0`.
When `is_doubling`, the tree corresponding to `ω₂` was obtained from a doubling step (this
can be relevant eg for biased progressive sampling).
The value `ω = logaddexp(ω₁, ω₂)` is provided for avoiding redundant calculations.
See [`biased_progressive_logprob2`](@ref) for an implementation.
"""
function calculate_logprob2 end

"""
    $(FUNCTIONNAME)(rng, trajectory, ζ₁, ζ₂, logprob2::Real, is_forward::Bool)
Combine two proposals `ζ₁, ζ₂` on `trajectory`, with log probability `logprob2` for
selecting `ζ₂`.
 `ζ₁` is before `ζ₂` iff `is_forward`.
"""
function combine_proposals end

"""
    ζωτ_or_nothing, v = $(FUNCTIONNAME)(trajectory, z, is_initial)
Information for a tree made of a single node. When `is_initial == true`, this is the first
node.
The first value is either
1. `nothing` for a divergent node,
2. a tuple containing the proposal `ζ`, the log weight (probability) of the node `ω`, the
turn statistics `τ` (never tested as with `is_turning` for leafs).
The second value is the visited node information.
"""
function leaf end

####
#### utilities
####

"""
$(SIGNATURES)
Combine turn statistics with the given direction. When `is_forward`, `τ₁` is before `τ₂`,
otherwise after.
Internal helper function.
"""
function combine_turn_statistics_in_direction(tree::Tree, trajectory, τ₁, τ₂, is_forward::Bool)
    if is_forward
        combine_turn_statistics(tree, trajectory, τ₁, τ₂)
    else
        combine_turn_statistics(tree, trajectory, τ₂, τ₁)
    end
end

function combine_proposals_and_logweights(
    rng, tree::Tree, trajectory, ζ₁, ζ₂, ω₁::Real, ω₂::Real, is_forward::Bool, is_doubling::Bool
)
    ω = logaddexp(ω₁, ω₂)
    logprob2 = calculate_logprob2(trajectory, is_doubling, ω₁, ω₂, ω)
    ζ = combine_proposals(rng, tree, trajectory, ζ₁, ζ₂, logprob2, is_forward)#, is_doubling)
    ζ, ω
end
# function combine_proposals_and_logweights(
    # rng, tree::Tree, trajectory, ζ₁, ζ₂, ω₁::Real, ω₂::Real, is_forward::Bool
# )
    # ω = logaddexp(ω₁, ω₂)
    # logprob2 = calculate_logprob2(trajectory, false, ω₁, ω₂, ω)
    # ζ = combine_proposals(rng, tree, trajectory, ζ₁, ζ₂, logprob2, is_forward)
    # ζ, ω
# end

"""
$(SIGNATURES)
Given (relative) log probabilities `ω₁` and `ω₂`, return the log probabiliy of
drawing a sample from the second (`logprob2`).
When `bias`, biases towards the second argument, introducing anti-correlations.
"""
function biased_progressive_logprob2(bias::Bool, ω₁::Real, ω₂::Real, ω = logaddexp(ω₁, ω₂))
    ω₂ - (bias ? ω₁ : ω)
end

####
#### abstract trajectory interface
####

"""
$(SIGNATURES)
Information about an invalid (sub)tree, using positions relative to the starting node.
1. When `left < right`, this tree was *turning*.
2. When `left == right`, this is a *divergent* node.
3. `left == 1 && right == 0` is used as a sentinel value for reaching maximum depth without
encountering any invalid trees (see [`REACHED_MAX_DEPTH`](@ref). All other `left > right`
values are disallowed.
"""
struct InvalidTree
    left::Int32
    right::Int32
end

InvalidTree(i::Integer) = InvalidTree(Base.unsafe_trunc(Int32, i), Base.unsafe_trunc(Int32, i))

is_divergent(invalid_tree::InvalidTree) = invalid_tree.left == invalid_tree.right

function Base.show(io::IO, invalid_tree::InvalidTree)
    msg = if is_divergent(invalid_tree)
        "divergence at position $(invalid_tree.left)"
    elseif invalid_tree == REACHED_MAX_DEPTH
        "reached maximum depth without divergence or turning"
    else
        @unpack left, right = invalid_tree
        "turning at positions $(left):$(right)"
    end
    print(io, msg)
end

"Sentinel value for reaching maximum depth."
const REACHED_MAX_DEPTH = InvalidTree(one(Int32), zero(Int32))

# on tree initialization, store

"""
    result, v = adjacent_tree(rng, trajectory, z, i, depth, is_forward)
Traverse the tree of given `depth` adjacent to point `z` in `trajectory`.
`is_forward` specifies the direction, `rng` is used for random numbers in
[`combine_proposals`](@ref). `i` is an integer position relative to the initial node (`0`).
The *first value* is either
1. an `InvalidTree`, indicating the first divergent node or turning subtree that was
encounteted and invalidated this tree.
2. a tuple of `(ζ, ω, τ, z′, i′), with
    - `ζ`: the proposal from the tree.
    - `ω`: the log weight of the subtree that corresponds to the proposal
    - `τ`: turn statistics
    - `z′`: the last node of the tree
    - `i′`: the position of the last node relative to the initial node.
The *second value* is always the visited node statistic.
"""
function adjacent_tree(rng, tree::Tree{P,T,L}, H, trajectory, z::PhasePoint{P,T,L}, i::Int32, depth::Int32, is_forward::Bool) where {P,T,L}
    i′ = i + (is_forward ? one(Int32) : -one(Int32) )
    # @show (1, bitstring(unsafe_load(reinterpret(Ptr{UInt32}, tree.root), 2)))
    # @show z.Q.q
    # @show depth, i′
    lb, ub = 5, 10
    # lb <= abs(i′) < ub && @show z
    if depth == zero(Int32) # moves from z into ζready
        z′ = move(tree, H, trajectory, z, is_forward)
        # lb <= abs(i′) < ub && @show logdensity(trajectory.H, z′), trajectory.π₀
        (ζ, ω, τ), v, invalid = leaf(tree, H, trajectory, z′, false)
        return (ζ, ω, τ, z′, i′), v, (invalid,InvalidTree(i′))
    else
        # “left” tree
        t₋, v₋, (invalid,it) = adjacent_tree(
            rng, tree, H, trajectory, z, i, depth - one(Int32), is_forward
        )
        # @show first(t₋)
        # @show t₋[4]
        invalid && return t₋, v₋, (invalid, it)
        ζ₋, ω₋, τ₋, z₋, i₋ = t₋

        # “right” tree — visited information from left is kept even if invalid
        t₊, v₊, (invalid,it) = adjacent_tree(
            rng, tree, H, trajectory, z₋, i₋, depth - one(Int32), is_forward
        )
        v = combine_visited_statistics(trajectory, v₋, v₊)
        invalid && return t₊, v, (invalid,it)
        ζ₊, ω₊, τ₊, z₊, i₊ = t₊

        # turning invalidates
        # try
        # @show (2, bitstring(unsafe_load(reinterpret(Ptr{UInt32}, tree.root), 2)))
        τ = combine_turn_statistics_in_direction(tree, trajectory, τ₋, τ₊, is_forward)
        # catch err
        #     @show
        #     rethrow(err)
        is_turning(trajectory, τ) && return t₊, v, (true, InvalidTree(i′, i₊))

        # valid subtree, combine proposals
        ζ, ω = combine_proposals_and_logweights(
            rng, tree, trajectory, ζ₋, ζ₊, ω₋, ω₊, is_forward, false
        )
        return (ζ, ω, τ, z₊, i₊), v, (false,REACHED_MAX_DEPTH)
    end
end

"""
$(SIGNATURES)
Sample a `trajectory` starting at `z`, up to `max_depth`. `directions` determines the tree
expansion directions.
Return the following values
- `ζ`: proposal from the tree
- `v`: visited node statistics
- `termination`: an `InvalidTree` (this includes the last doubling step turning, which is
  technically a valid tree) or `REACHED_MAX_DEPTH` when all subtrees were valid and no
  turning happens.
- `depth`: the depth of the tree that was sampled from. Doubling steps that lead to an
  invalid adjacent tree do not contribute to `depth`.
"""
function sample_trajectory(rng, tree::Tree, H, trajectory, z::PhasePoint{P,T,L}, max_depth::Integer, directions::Directions) where {P,T,L}
    #    @argcheck max_depth ≤ MAX_DIRECTIONS_DEPTH
    # original_flag = zᵢ.flag
    # protect_initial = true # Protect initial position by giving it a dummy flag
    # z = PhasePoint(zᵢ.Q, zᵢ.p, 0x00000000)
    # @show logdensity(trajectory.H, z), trajectory.π₀
    (ζ, ω, τ), v, invalid = leaf(tree, H, trajectory, z, true)
    z₋ = z₊ = z
    # z₋flag = z₊flag = 0x00000000
    depth = zero(Int32)
    termination = REACHED_MAX_DEPTH
    i₋ = i₊ = zero(Int32)
    # dealloc₋ = dealloc₊ = false
    while depth < max_depth
        is_forward, directions = next_direction(directions)
        # @show depth, is_forward
        if is_forward
            zᵢ, iᵢ = z₊, i₊
            alloc = z₋.flag
        else
            zᵢ, iᵢ = z₋, i₋
            alloc = z₊.flag
        end
        # ((ζ.flag === zᵢ) | (z₊.flag === z₋.flag)) || free_z!(tree, zᵢ.flag)
        VectorizationBase.store!(
            reinterpret(Ptr{UInt32}, tree.root), 0xffffffff ⊻ (alloc | ζ.flag)
        )
        # clear_all_but_z!( tree, alloc | ζ.flag )
        t′, v′, (invalid, it) = adjacent_tree(
            rng, tree, H, trajectory, zᵢ, iᵢ, depth, is_forward
        )

        v = combine_visited_statistics(trajectory, v, v′)

        # invalid adjacent tree: stop
        invalid && (termination = it; break)

        # extract information from adjacent tree
        ζ′, ω′, τ′, z′, i′ = t′
        # allocate!(tree, z′.flag)
        # ζ′g = z′.flag === ζ′.flag ? PhasePoint(ζ′.Q, ζ′.p, 0x00000000) : ζ′
        # update edges and combine proposals
        if is_forward
            z₊, i₊ = z′, i′
        else
            z₋, i₋ = z′, i′
        end

        # tree has doubled successfully
        ζ, ω = combine_proposals_and_logweights(
            rng, tree, trajectory, ζ, ζ′, ω, ω′, is_forward, true
        )
        depth += one(Int32)

        # when the combined tree is turning, stop
        τ = combine_turn_statistics_in_direction(tree, trajectory, τ, τ′, is_forward)
        is_turning(trajectory, τ) && (termination = InvalidTree(i₋, i₊); break)
    end
    # @show  original_flag, ζ.flag
    # clear_all_but_z!( tree, ζ.flag )# == 0x00000000 ? original_flag : ζ.flag )
    clear!( tree )
    ζ, v, termination, depth
end


text:
__precompile__() # this module is safe to precompile
module SolarRotation
using PyCall

const AstropyCoordinates = PyNULL()
const AstropyTime = PyNULL()
const Barycorrpy = PyNULL()
const SciPySpatialTransform = PyNULL()

function __init__()
	copy!(AstropyCoordinates , pyimport("astropy.coordinates") )
	copy!(AstropyTime , pyimport("astropy.time") )
	copy!(Barycorrpy , pyimport("barycorrpy") )
	# Uncomment below and remove pyimport from CalculateFWHMDifference_SolarRotation_from_loc_Ecliptic if start using it
	#copy!(SciPySpatialTransform , pyimport("scipy.spatial.transform") )
end

using Dates, LinearAlgebra
using DataFrames

"""
  `CalculateFWHMDifference_SolarRotation_from_loc_Ecliptic(jd; loc )`
Calculate the difference between the Observed Solar FWHM and Sidereal Solar FWHM
Inputs:
	jd: Julian Date (scalar)
	loc: Astropy Earth Location object. https://docs.astropy.org/en/stable/api/astropy.coordinates.EarthLocation.html

Output:
	Delta:  F_obs**2 - F_sid**2 [(km/s)^2]

Based on Colier Cameron et al. (2019)
Uses ses JPL ephemeris rather than assuming a Keplerian orbit
Performs calculations in ecliptic reference frame.

Python version by Shubham Kanodia March 2021 (with input from Jason Wright and Eric Ford)
Translated to Julia by Eric Ford March 9, 2021.
"""
function CalculateFWHMDifference_SolarRotation_from_loc_Ecliptic(jd::Real; loc::PyObject )
#function CalculateFWHMDifference_SolarRotationGen(loc::PyObject, jd::AVR ) where { T<:Real, AVR<:AbstractVector{T} } # JDUTC::PyObject)
	# Functions and constants from AstroPy or SciPy that we'll use multiple times
	get_body_barycentric_posvel = AstropyCoordinates.get_body_barycentric_posvel
	get_body_barycentric = AstropyCoordinates.get_body_barycentric
	speed_of_light = Barycorrpy.PhysicalConstants.c
	R_sun = 695700  # km # should we use better, but non-standard value, e.g, 695660?

	SciPySpatialTransform = pyimport("scipy.spatial.transform")
	SciPyRotation = SciPySpatialTransform.Rotation

	JDUTC = AstropyTime.Time(jd , format="jd", scale="utc")
	# Convert times to obtain TDB and TT
	JDTDB = JDUTC.tdb
	JDTT = JDUTC.tt

	################################
	######EARTH EPHEMERIS ############
	################################
	##### NUTATION, PRECESSION, ETC. #####

	# Observatory position wrt Geocenter
	r_pint, v_pint = Barycorrpy.PINT_erfautils.gcrs_posvel_from_itrf(loc, JDUTC, JDTT)
	r_eci = r_pint'  # [m]
	v_eci = v_pint'  # [m/s]

	##### EPHEMERIDES #####
	ephemeris = "de430"
	earth_geo = get_body_barycentric_posvel("earth", JDTDB, ephemeris=ephemeris) # [km]
	r_geo = earth_geo[1].xyz*1000.0 # [m]
	v_geo = earth_geo[2].xyz*1000.0/86400.0  # [m/s]

	PosVector_EarthSSB = r_eci + r_geo # [m]MethodError

	# Relativistic Addition of Velocities
	VelVector_EarthSSB = (v_eci + v_geo) / (1.0 + dot(v_eci,v_geo)/speed_of_light^2) # [m/s]

	################################
	######SOLAR EPHEMERIS ############
	################################
	solar_ephem = get_body_barycentric_posvel("sun", JDTDB, ephemeris=ephemeris)
	PosVector_SolSSB = solar_ephem[1].xyz*1000. #[m]
	VelVector_SolSSB = solar_ephem[2].xyz*1000.0/86400.0  # [m/s]

	################################
	####EQUATORIAL COORD VECTORS ####
	################################
	PosVector_EarthSol = PosVector_EarthSSB - PosVector_SolSSB
	PosMag_EarthSol = norm(PosVector_EarthSol)

	VelVector_EarthSol = (VelVector_EarthSSB - VelVector_SolSSB) / (1.0 + dot(VelVector_SolSSB,VelVector_EarthSSB)/speed_of_light^2)

	#OmegaVector = map(i->cross(PosVector_EarthSol[:,i], VelVector_EarthSol[:,i]) ./ dot(PosVector_EarthSol[:,i],PosVector_EarthSol[:,i]), 1:size(PosVector_EarthSol,2))
	OmegaVector = cross(vec(PosVector_EarthSol), vec(VelVector_EarthSol)) ./ dot(PosVector_EarthSol,PosVector_EarthSol)
	################################
	################################

	#DeltaCentury = map(i->(Dates.year(JDUTC.datetime) .- 2000)./100, 1:length(jd) )
	DeltaCentury = (Dates.year(JDUTC.datetime) - 2000)/100

	# https://www2.mps.mpg.de/homes/fraenz/systems/systems3art.pdf

	#Eqn 14
	SolarInclination = 7.25
	SolarLongitude = 75.76 + 1.397*DeltaCentury
	EclipticEpsilon = 23.44

	# Need to perform Extrinsic Euler Rotations to rotate from equatorial to ecliptic
	#REcliptic = R.from_euler("X",  -EclipticEpsilon, degrees=True).as_matrix()
	REcliptic = SciPyRotation.from_euler("X",  -EclipticEpsilon, degrees=true).as_dcm()
	# may need to change to .as_matrix() with SciPy 1.4?

	# Intrinsic rotation to go from solar axis to ecliptic
	#RObliquity = R.from_euler("xz",  [SolarInclination, SolarLongitude], degrees=true).as_matrix()
	RObliquity = SciPyRotation.from_euler("xz",  [SolarInclination, SolarLongitude], degrees=true).as_dcm()

	################################
	####ROTATED COORD VECTORS ####
	################################

	#RotatedPositionVector = map(i->REcliptic * PosVector_EarthSol[:,i], 1:size(PosVector_EarthSol,2) )
	#RotatedPositionHat = RotatedPositionVector ./ norm.(RotatedPositionVector)
	RotatedPositionVector = REcliptic * PosVector_EarthSol
	RotatedPositionHat = RotatedPositionVector / norm(RotatedPositionVector)

	#OmegaEarthVectorRotated = map(i->REcliptic * OmegaVector[i],1:size(OmegaVector,2))
	OmegaEarthVectorRotated = REcliptic * OmegaVector
	################################
	################################

	OmegaSolVector = [0.0, 0.0,  2.972 *1e-6]
	#OmegaSolHat = [0.0 , 0.0 , 1.0]

	# Rotated Solar rotation vector to ecliptic plane
	#OmegaSolVector = map(i->RObliquity * OmegaSolVector[:,i], 1:size(OmegaSolVector,2) )
	#OmegaSolHat = OmegaSolVector./norm.(OmegaSolVector)
	OmegaSolVector = RObliquity * OmegaSolVector
	OmegaSolHat = OmegaSolVector/norm(OmegaSolVector)

	sini = sqrt(1 - dot(OmegaSolHat,RotatedPositionHat)^2)
	Gamma = 1.04
	DeltaOmega = OmegaSolVector .- OmegaEarthVectorRotated

	Delta = ((Gamma* R_sun)^2) * (dot(DeltaOmega,DeltaOmega)*sini^2 - dot(OmegaSolVector,OmegaSolVector))
end

"""
   `CalculateFWHMDifference_SolarRotation_from_loc_Equatorial(jd; loc )`
Calculate the difference between the Observed Solar FWHM and Sidereal Solar FWHM
Inputs:
  jd: Julian Date (scalar)
  loc: Astropy Earth Location object. https://docs.astropy.org/en/stable/api/astropy.coordinates.EarthLocation.html

Output:
  Delta:  F_obs**2 - F_sid**2 [(km/s)^2]

Based on Colier Cameron et al. (2019)
Uses ses JPL ephemeris rather than assuming a Keplerian orbit.
Performs calculations in equatorial reference frame.

Python version by Shubham Kanodia (with input from Jason Wright and Eric Ford) March 2021
Translated to Julia by Eric Ford March 9, 2021.
"""
function CalculateFWHMDifference_SolarRotation_from_loc_Equatorial(jd::Real; loc::PyObject )
#function CalculateFWHMDifference_SolarRotationGen(loc::PyObject, jd::AVR ) where { T<:Real, AVR<:AbstractVector{T} } # JDUTC::PyObject)
	# Functions and constants from AstroPy or SciPy that we'll use multiple times
	get_body_barycentric_posvel = AstropyCoordinates.get_body_barycentric_posvel
	get_body_barycentric = AstropyCoordinates.get_body_barycentric
	speed_of_light = Barycorrpy.PhysicalConstants.c
	R_sun = 695700  # km # should we use better, but non-standard value, e.g, 695660?

	JDUTC = AstropyTime.Time(jd , format="jd", scale="utc")
	#return JDUTC
	# Convert times to obtain TDB and TT
	JDTDB = JDUTC.tdb
	JDTT = JDUTC.tt

	################################
	######EARTH EPHEMERIS ############
	################################
	##### NUTATION, PRECESSION, ETC. #####

	# Observatory position wrt Geocenter
	r_pint, v_pint = Barycorrpy.PINT_erfautils.gcrs_posvel_from_itrf(loc, JDUTC, JDTT)
	r_eci = r_pint'  # [m]
	v_eci = v_pint'  # [m/s]

	##### EPHEMERIDES #####
	ephemeris = "de430"
	earth_geo = get_body_barycentric_posvel("earth", JDTDB, ephemeris=ephemeris) # [km]
	r_geo = earth_geo[1].xyz*1000.0 # [m]
	v_geo = earth_geo[2].xyz*1000.0/86400.0  # [m/s]

	PosVector_EarthSSB = r_eci + r_geo # [m]

	# Relativistic Addition of Velocities
	VelVector_EarthSSB = (v_eci + v_geo) / (1.0 + dot(v_eci,v_geo)/speed_of_light^2) # [m/s]

	################################
	######SOLAR EPHEMERIS ############
	################################
	solar_ephem = get_body_barycentric_posvel("sun", JDTDB, ephemeris=ephemeris)
	PosVector_SolSSB = solar_ephem[1].xyz*1000. #[m]
	VelVector_SolSSB = solar_ephem[2].xyz*1000.0/86400.0  # [m/s]

	################################
	####EQUATORIAL COORD VECTORS ####
	################################
	PosVector_EarthSol = PosVector_EarthSSB - PosVector_SolSSB
	PosMag_EarthSol = norm(PosVector_EarthSol)
	PosHat_EarthSol = PosVector_EarthSol ./ PosMag_EarthSol
	VelVector_EarthSol = (VelVector_EarthSSB - VelVector_SolSSB) / (1.0 + dot(VelVector_SolSSB,VelVector_EarthSSB)/speed_of_light^2)

	#OmegaVector = map(i->cross(PosVector_EarthSol[:,i], VelVector_EarthSol[:,i]) ./ dot(PosVector_EarthSol[:,i],PosVector_EarthSol[:,i]), 1:size(PosVector_EarthSol,2))
	OmegaVector = cross(vec(PosVector_EarthSol), vec(VelVector_EarthSol)) ./ dot(PosVector_EarthSol,PosVector_EarthSol)

	################################
	################################
	# Rotation Axis of the Sun with respect to Equatorial System
	# https://www2.mps.mpg.de/homes/fraenz/systems/systems3art.pdf

	#Eqn 13
	alpha = 286.13 # Rotate along the z axis
	dec = 63.87 # Rotate along the x axis
	Theta = alpha * π/180
	Phi = dec * π/180
	EclipticEpsilon = 23.44
	SolarInclination = 7.25
	SolarLongitude = 75.76

	# Need to perform Extrinsic Euler Rotations
	# Transform to solar rotation axis
	REquatorial = [cos(Theta)*cos(Phi), sin(Theta)*cos(Phi), sin(Phi)]

	OmegaSolVector = [0.0, 0.0,  2.972 *1e-6]
	#OmegaSolHat = [0.0 , 0.0 , 1.0]

	OmegaSolHat = REquatorial
	OmegaSolVector = REquatorial * norm(OmegaSolVector)
	#OmegaSolHat = OmegaSolVector / norm(OmegaSolVector)

	sini = sqrt(1 - dot(OmegaSolHat,PosHat_EarthSol)^2)
	Gamma = 1.04
	DeltaOmega = OmegaSolVector .- OmegaVector

	Delta = ((Gamma* R_sun)^2) * (dot(DeltaOmega,DeltaOmega)*sini^2 - dot(OmegaSolVector,OmegaSolVector))

end

"""
	`CalculateFWHMDifference_SolarRotation_from_long_lat_alt(jd; long, lat, alt)`
Calculate the difference between the Observed Solar FWHM and Sidereal Solar FWHM
	Inputs:
		jd: Julian Date (scalar)
		long: longitude (degrees, west negative)
		lat: latitude (degrees)
		alt: altitude (m)
"""
function CalculateFWHMDifference_SolarRotation_from_long_lat_alt(jd::Real; long::Real, lat::Real, alt::Real )
	loc = AstropyCoordinates.EarthLocation.from_geodetic(long, lat, height=alt)
	delta = CalculateFWHMDifference_SolarRotation_from_loc_Equatorial(jd, loc=loc)
	#delta = CalculateFWHMDifference_SolarRotation_from_loc_Ecliptic(jd, loc=loc)
	return delta
end

"""
	`CalculateFWHMDifference_SolarRotation_from_obs(jd; obs)`
Calculate the difference between the Observed Solar FWHM and Sidereal Solar FWHM
Inputs:
	jd: Julian Date (scalar)
	obs: observatory name (symbol)
Currently has HARPS-N and NEID/WIYN coordinates.
"""
function CalculateFWHMDifference_SolarRotation_from_obs(jd::Real; obs::Symbol )
	if obs == :HARPSN
		longi = -17.88905
		lat = 28.754
		alt = 2387.2
	elseif  obs == :WIYN
		longi = -111.600562
		lat =  31.958092
		alt =  2091.0
	else
		@error("Don't have coordinates for obs = " * string(obs))
	end
	CalculateFWHMDifference_SolarRotation_from_long_lat_alt(jd, long=longi, lat=lat, alt=alt)
end

#=
function get_wiyn_loc()
	obsname = "WIYN"
	longi = -111.600562
	lat =  31.958092
	alt =  2091.0
	AstropyCoordinates.EarthLocation.from_geodetic(longi, lat, height=alt)
end

function get_harpsn_loc()
	obsname = "HARPS-N"
	longi = -17.88905
	lat = 28.754
	alt = 2387.2
	loc = AstropyCoordinates.EarthLocation.from_geodetic(longi, lat, height=alt)
end
=#

export CalculateFWHMDifference_SolarRotation_from_loc_Ecliptic, CalculateFWHMDifference_SolarRotation_from_loc_Equatorial
export CalculateFWHMDifference_SolarRotation_from_long_lat_alt, CalculateFWHMDifference_SolarRotation_from_obs

#export get_wiyn_loc, get_harpsn_loc

""" get_obs_loc(obs::Symbol)
 Returns a Dict with long & lat (degrees) and elevation (km)
 Warning: Currently only has info for :WIYN and :HARPSN.
"""
function get_obs_loc(obs::Symbol)
	valid_obs = [:WIYN, :HARPSN]
	@assert obs ∈ valid_obs
	if obs == :HARPSN
		loc = Dict("lon"=> -17.88905, "lat"=> 28.754, "elevation"=> 2.3872)
	elseif obs == :WIYN
		loc = Dict("lon"=> -111.600562, "lat"=> 31.958092, "elevation"=> 2.091)
	else
		@error("Don't have coordinates for obs = " * string(obs))
	end
	return loc
end

function get_lst(jd::Union{Real,Vector{<:Real}}; obs::Symbol)
   loc = get_obs_loc(obs)
   pyloc = AstropyCoordinates.EarthLocation.from_geodetic(loc["lon"], loc["lat"], height=loc["elevation"])
   TimeObs = AstropyTime.Time(jd, format="jd", scale="utc", location=pyloc)
   LST = TimeObs.sidereal_time("mean")
   if typeof(jd) <: Real
	   return first(LST)
   else
	   return LST
   end
end

function get_solar_ra_dec(jd::Union{Real,Vector{<:Real}}; obs::Symbol)
	jlloc = get_obs_loc(obs)
    pyloc = AstropyCoordinates.EarthLocation.from_geodetic(jlloc["lon"], jlloc["lat"], height=jlloc["elevation"]) #*1000)
    TimeObs = AstropyTime.Time(jd, format="jd", scale="utc", location=pyloc)
	SolarObj = AstropyCoordinates.get_sun(TimeObs)
	RA, Dec = SolarObj.ra, SolarObj.dec
	return (ra=vec(RA), dec=vec(Dec))
end

function get_solar_hour_angle(jd::Union{Real,Vector{<:Real}}; obs::Symbol)
	lst = get_lst(jd, obs=obs)
	ra = get_solar_ra_dec(jd, obs=obs).ra
	ha = lst .- ra./15
	ha[ha.<-12] .+= 24
	ha[ha.>12] .-= 24
	return ha
end

function get_earth_sun_dist(jd::Union{Real,Vector{<:Real}}; obs::Symbol)
	jlloc = get_obs_loc(obs)
    loc = AstropyCoordinates.EarthLocation.from_geodetic(jlloc["lon"], jlloc["lat"], height=jlloc["elevation"])
	ephemeris = "de430"
	JDUTC = AstropyTime.Time(jd , format="jd", scale="utc")
	# Convert times to obtain TDB and TT
	JDTDB = JDUTC.tdb
	JDTT = JDUTC.tt

	get_body_barycentric_posvel = AstropyCoordinates.get_body_barycentric_posvel
	earth_geo = get_body_barycentric_posvel("earth", JDTDB, ephemeris=ephemeris)[1].xyz * 1000.0 # [km->m]
	solar_ephem = get_body_barycentric_posvel("sun", JDTDB, ephemeris=ephemeris)[1].xyz * 1000.0 # [km->m]
	r_pint, v_pint = Barycorrpy.PINT_erfautils.gcrs_posvel_from_itrf(loc, JDUTC, JDTT) # [m]
	diff_vector = (earth_geo + r_pint') -solar_ephem #[m]
	dist = vec(sqrt.(sum(diff_vector.^2,dims=1)))
	AU = 149597870.7  # km
	dist = dist/(AU*1000)
	if type(jd) <:Real
		return first(dist)
	else
		return dist
	end
end

function get_solar_info(jd::Union{Real,Vector{<:Real}}; obs::Symbol)
	jlloc = get_obs_loc(obs)
    pyloc = AstropyCoordinates.EarthLocation.from_geodetic(jlloc["lon"], jlloc["lat"], height=jlloc["elevation"]) #*1000)
    TimeObs = AstropyTime.Time(jd, format="jd", scale="utc", location=pyloc)
	SolarObj = AstropyCoordinates.get_sun(TimeObs)
	frame_obsnight = AstropyCoordinates.AltAz(obstime=TimeObs, location=pyloc)
	SolarAltAz = SolarObj.transform_to(frame_obsnight)
	RA, Dec = SolarObj.ra, SolarObj.dec
	Alt, Az = SolarAltAz.alt, SolarAltAz.az
	Airmass = SolarAltAz.secz
	LST = TimeObs.sidereal_time("mean")

	ephemeris = "de430"
	JDUTC = AstropyTime.Time(jd , format="jd", scale="utc")
	# Convert times to obtain TDB and TT
	JDTDB = JDUTC.tdb
	JDTT = JDUTC.tt
	r_pint, v_pint = Barycorrpy.PINT_erfautils.gcrs_posvel_from_itrf(pyloc, JDUTC, JDTT) # [m]
	get_body_barycentric_posvel = AstropyCoordinates.get_body_barycentric_posvel
	earth_geo = get_body_barycentric_posvel("earth", JDTDB, ephemeris=ephemeris)[1].xyz * 1000.0 # [km->m]
	solar_ephem = get_body_barycentric_posvel("sun", JDTDB, ephemeris=ephemeris)[1].xyz * 1000.0 # [km->m]
	diff_vector = (earth_geo + r_pint') -solar_ephem #[m]
	dist = sqrt.(sum(diff_vector.^2,dims=1))
	AU = 149597870.7  # km
	dist = dist./(AU*1000)

	if typeof(jd) <: Real
		result = Dict(:ra=>first(RA), :dec=>first(Dec), :alt=>first(Alt), :az=>first(Az), :airmass=>first(Airmass), :lst=>first(LST), :sol_dist_au=>first(dist) )
	else
		result = DataFrame(:ra=>RA, :dec=>Dec, :alt=>Alt, :az=>Az, :airmass=>Airmass, :lst=>LST, :sol_dist_au=>dist )
	end
	return result
end

export get_obs_loc, get_lst, get_solar_hour_angle, get_solar_ra_dec, get_earth_sun_dist, get_solar_info

end # module SolarRotation


text:
using PGFPlots
using CSV
using Query
using DataFrames

# ld
filename = Pkg.dir("ContinuousPOMDPTreeSearchExperiments", "icaps_2018", "data", "ld_discretization_Tuesday_27_Feb_15_23.csv")
alldata = CSV.read(filename)

aggregated = by(alldata, [:solver, :d]) do df
    r = df[:reward]
    return DataFrame(mean=mean(r), sem=std(r)/sqrt(length(r)))
end

names = Dict("pomcpow"=>"POMCPOW\\textsuperscript{D}",
             "d_pomcp"=>"POMCP\\textsuperscript{D}",
             "d_despot"=>"DESPOT\\textsuperscript{D}")
plots = Plots.Plot[]
for k in ["d_pomcp", "d_despot", "pomcpow"]
    series = @from pt in aggregated begin
        @where pt.solver == k
        @select {pt.d, pt.mean, pt.sem}
        @collect DataFrame
    end
    push!(plots, Plots.Linear(series, x=:d, y=:mean, errorBars=ErrorBars(y=series[:sem]), legendentry=names[k]))
end
ds = sort(unique(aggregated[:d]))
push!(plots, Plots.Linear(ds,
                          fill!(similar(ds), 57.7),
                          errorBars=ErrorBars(y=fill!(similar(ds), 0.5)),
                          legendentry="POMCPOW\\phantom{\\textsuperscript{D}}"))
p = plot(Axis(plots,
              xlabel="Discretization bin size",
              ylabel="Mean accumulated reward",
              xmode="log",
              style="legend style={at={(1.06,0.3)},anchor=south east,font=\\small}"))
save("/tmp/ld_discretization.pdf", p)

# Subhunt

filename = Pkg.dir("ContinuousPOMDPTreeSearchExperiments", "icaps_2018", "data", "subhunt_discretization_Tuesday_27_Feb_11_38.csv")
alldata = CSV.read(filename)

aggregated = by(alldata, [:solver, :binsize]) do df
    r = df[:reward]
    return DataFrame(mean=mean(r), sem=std(r)/sqrt(length(r)))
end

# println(aggregated)

names = Dict("pomcpow"=>"POMCPOW\\textsuperscript{D}",
             "d_pomcp"=>"POMCP\\textsuperscript{D}",
             "d_despot"=>"DESPOT\\textsuperscript{D}")
plots = Plots.Plot[]
for k in ["d_pomcp", "d_despot", "pomcpow"]
    series = @from pt in aggregated begin
        @where pt.solver == k
        @select {pt.binsize, pt.mean, pt.sem}
        @collect DataFrame
    end
    push!(plots, Plots.Linear(series, x=:binsize, y=:mean, errorBars=ErrorBars(y=series[:sem]), legendentry=names[k]))
end
ds = sort(unique(aggregated[:binsize]))
push!(plots, Plots.Linear(ds,
                          fill!(similar(ds), 69.2),
                          errorBars=ErrorBars(y=fill!(similar(ds), 1.3)),
                          legendentry="POMCPOW\\phantom{\\textsuperscript{D}}"))
p = plot(Axis(plots,
              ymin=0,
              xmode="log",
              xlabel="Discretization bin size",
              ylabel="Mean accumulated reward",
              style="legend style={at={(0.01,0.45)},anchor=south west,font=\\small}"))
save("/tmp/subhunt_discretization.pdf", p)


text:
function LazyLoadedEop(arg0::DataProvidersManager)
    return LazyLoadedEop((DataProvidersManager,), arg0)
end

function add_default_eop1980_history_loaders(obj::LazyLoadedEop, arg0::JString, arg1::JString, arg2::JString, arg3::JString, arg4::JString, arg5::Supplier)
    return jcall(obj, "addDefaultEOP1980HistoryLoaders", void, (JString, JString, JString, JString, JString, Supplier), arg0, arg1, arg2, arg3, arg4, arg5)
end

function add_default_eop2000_history_loaders(obj::LazyLoadedEop, arg0::JString, arg1::JString, arg2::JString, arg3::JString, arg4::JString, arg5::Supplier)
    return jcall(obj, "addDefaultEOP2000HistoryLoaders", void, (JString, JString, JString, JString, JString, Supplier), arg0, arg1, arg2, arg3, arg4, arg5)
end

function add_eop_history_loader(obj::LazyLoadedEop, arg0::IERSConventions, arg1::EOPHistoryLoader)
    return jcall(obj, "addEOPHistoryLoader", void, (IERSConventions, EOPHistoryLoader), arg0, arg1)
end

function clear_eop_history_loaders(obj::LazyLoadedEop)
    return jcall(obj, "clearEOPHistoryLoaders", void, ())
end

function get_data_providers_manager(obj::LazyLoadedEop)
    return jcall(obj, "getDataProvidersManager", DataProvidersManager, ())
end

function get_eop_history(obj::LazyLoadedEop, arg0::IERSConventions, arg1::jboolean, arg2::TimeScales)
    return jcall(obj, "getEOPHistory", EOPHistory, (IERSConventions, jboolean, TimeScales), arg0, arg1, arg2)
end

function set_eop_continuity_threshold(obj::LazyLoadedEop, arg0::jdouble)
    return jcall(obj, "setEOPContinuityThreshold", void, (jdouble,), arg0)
end


text:
function plot_sum(rates,Color =nothing)
    prob = ODEProblem(f,u₀,tspan,rates)
    # prob = remake(prob, p=10 .^ rates)
    sol = solve(prob,Rosenbrock23())
    summ = reduce(vcat,sum(sol,dims=1))
    if Color == nothing
        plot!(sol.t,summ)
    else
        plot!(sol.t,summ,c=Color)
    end
end

function display_rates(rates)
    for rate in rates
        @show rate
    end
end


text:
# ---
# title: 1473. Paint House III
# id: problem1473
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/paint-house-iii/description/>
# hidden: true
# ---
#
# There is a row of `m` houses in a small city, each house must be painted with
# one of the `n` colors (labeled from 1 to `n`), some houses that has been
# painted last summer should not be painted again.
#
# A neighborhood is a maximal group of continuous houses that are painted with
# the same color. (For example: houses = [1,2,2,3,3,2,1,1] contains 5
# neighborhoods  [{1}, {2,2}, {3,3}, {2}, {1,1}]).
#
# Given an array `houses`, an `m * n` matrix `cost` and an integer `target`
# where:
#
#   * `houses[i]`: is the color of the house `i`, **0** if the house is not painted yet.
#   * `cost[i][j]`: is the cost of paint the house `i` with the color `j+1`.
#
# Return the minimum cost of painting all the remaining houses in such a way
# that there are exactly `target` neighborhoods, if not possible return **-1**.
#
#
#
# **Example 1:**
#
#
#
#     Input: houses = [0,0,0,0,0], cost = [[1,10],[10,1],[10,1],[1,10],[5,1]], m = 5, n = 2, target = 3
#     Output: 9
#     Explanation: Paint houses of this way [1,2,2,1,1]
#     This array contains target = 3 neighborhoods, [{1}, {2,2}, {1,1}].
#     Cost of paint all houses (1 + 1 + 1 + 1 + 5) = 9.
#
#
# **Example 2:**
#
#
#
#     Input: houses = [0,2,1,2,0], cost = [[1,10],[10,1],[10,1],[1,10],[5,1]], m = 5, n = 2, target = 3
#     Output: 11
#     Explanation: Some houses are already painted, Paint the houses of this way [2,2,1,2,2]
#     This array contains target = 3 neighborhoods, [{2,2}, {1}, {2,2}].
#     Cost of paint the first and last house (10 + 1) = 11.
#
#
# **Example 3:**
#
#
#
#     Input: houses = [0,0,0,0,0], cost = [[1,10],[10,1],[1,10],[10,1],[1,10]], m = 5, n = 2, target = 5
#     Output: 5
#
#
# **Example 4:**
#
#
#
#     Input: houses = [3,1,2,3], cost = [[1,1,1],[1,1,1],[1,1,1],[1,1,1]], m = 4, n = 3, target = 3
#     Output: -1
#     Explanation: Houses are already painted with a total of 4 neighborhoods [{3},{1},{2},{3}] different of target = 3.
#
#
#
#
# **Constraints:**
#
#   * `m == houses.length == cost.length`
#   * `n == cost[i].length`
#   * `1 <= m <= 100`
#   * `1 <= n <= 20`
#   * `1 <= target <= m`
#   * `0 <= houses[i] <= n`
#   * `1 <= cost[i][j] <= 10^4`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
include("dependencies_for_runtests.jl")

using Statistics
using Oceananigans.Architectures: arch_array
using Oceananigans.AbstractOperations: BinaryOperation
using Oceananigans.Fields: ReducedField, CenterField, ZFaceField, compute_at!, @compute
using Oceananigans.BoundaryConditions: fill_halo_regions!
using Oceananigans.Grids: halo_size

trilinear(x, y, z) = x + y + z

@testset "Fields computed by Reduction" begin
    @info "Testing Fields computed by reductions..."

    for arch in archs
        arch_str = string(typeof(arch))

        regular_grid = RectilinearGrid(arch, size = (2, 2, 2),
                                       x = (0, 2), y = (0, 2), z = (0, 2),
                                       topology = (Periodic, Periodic, Bounded))

        xy_regular_grid = RectilinearGrid(arch, size=(2, 2, 2),
                                          x=(0, 2), y=(0, 2), z=[0, 1, 2],
                                          topology = (Periodic, Periodic, Bounded))

        @testset "Averaged fields [$arch_str]" begin
            @info "  Testing averaged Fields [$arch_str]"

            for grid in (regular_grid, xy_regular_grid)

                Nx, Ny, Nz = size(grid)

                w = ZFaceField(grid)
                T = CenterField(grid)
                ζ = Field{Face, Face, Face}(grid)

                set!(T, trilinear)
                set!(w, trilinear)
                set!(ζ, trilinear)

                @compute Txyz = Field(Average(T, dims=(1, 2, 3)))

                # Note: halo regions must be *filled* prior to computing an average
                # if the average within halo regions is to be correct.
                fill_halo_regions!(T, arch)
                @compute Txy = Field(Average(T, dims=(1, 2)))

                fill_halo_regions!(T, arch)
                @compute Tx = Field(Average(T, dims=1))

                @compute wxyz = Field(Average(w, dims=(1, 2, 3)))
                @compute wxy = Field(Average(w, dims=(1, 2)))
                @compute wx = Field(Average(w, dims=1))

                @compute ζxyz = Field(Average(ζ, dims=(1, 2, 3)))
                @compute ζxy = Field(Average(ζ, dims=(1, 2)))
                @compute ζx = Field(Average(ζ, dims=1))

                for T′ in (Tx, Txy)
                    @test T′.operand.operand === T
                end

                for w′ in (wx, wxy)
                    @test w′.operand.operand === w
                end

                for ζ′ in (ζx, ζxy)
                    @test ζ′.operand.operand === ζ
                end

                for f in (wx, wxy, Tx, Txy, ζx, ζxy)
                    @test f.operand isa Reduction
                    @test f.operand.reduce! === mean!
                end

                @test Txyz.operand isa Reduction
                @test wxyz.operand isa Reduction
                @test ζxyz.operand isa Reduction

                # Different behavior for regular grid z vs not.
                if grid === regular_grid
                    @test Txyz.operand.reduce! === mean!
                    @test wxyz.operand.reduce! === mean!
                    @test Txyz.operand.operand === T
                    @test wxyz.operand.operand === w
                else
                    @test Txyz.operand.reduce! === sum!
                    @test wxyz.operand.reduce! === sum!
                    @test Txyz.operand.operand isa BinaryOperation
                    @test wxyz.operand.operand isa BinaryOperation
                end

                @test Tx.operand.dims === tuple(1)
                @test wx.operand.dims === tuple(1)
                @test Txy.operand.dims === (1, 2)
                @test wxy.operand.dims === (1, 2)
                @test Txyz.operand.dims === (1, 2, 3)
                @test wxyz.operand.dims === (1, 2, 3)

                @test Txyz[1, 1, 1] ≈ 3
                @test Array(interior(Txy))[1, 1, :] ≈ [2.5, 3.5]
                @test Array(interior(Tx))[1, :, :] ≈ [[2, 3] [3, 4]]
                @test wxyz[1, 1, 1] ≈ 3
                @test Array(interior(wxy))[1, 1, :] ≈ [2, 3, 4]
                @test Array(interior(wx))[1, :, :] ≈ [[1.5, 2.5] [2.5, 3.5] [3.5, 4.5]]

                @compute Txyz = Field(Average(T, condition=T.>3))
                @compute Txy = Field(Average(T, dims=(1, 2), condition=T.>3))
                @compute Tx = Field(Average(T, dims=1, condition=T.>2))
                @test Txyz[1,1,1] ≈ 3.75
                @test Array(interior(Txy))[1, 1, :] ≈ [3.5, 11.5/3]
                @test Array(interior(Tx))[1, :, :] ≈ [[2.5, 3] [3, 4]]

                @compute wxyz = Field(Average(w, condition=w.>3))
                @compute wxy = Field(Average(w, dims=(1, 2), condition=w.>2))
                @compute wx = Field(Average(w, dims=1, condition=w.>1))
                @test wxyz[1,1,1] ≈ 4.25
                @test Array(interior(wxy))[1, 1, :] ≈ [3, 10/3, 4]
                @test Array(interior(wx))[1, :, :] ≈ [[2, 2.5] [2.5, 3.5] [3.5, 4.5]]
            end

            # Test whether a race condition gets hit for averages over large fields
            big_grid = RectilinearGrid(arch,
                                       topology = (Periodic, Periodic, Bounded),
                                       size = (256, 256, 128),
                                       x = (0, 2),
                                       y = (0, 2),
                                       z = (0, 2))

            c = CenterField(big_grid)
            c .= 1

            C = Field(Average(c, dims=(1, 2)))

            # Test that the mean consistently returns 1 at every z for many evaluations
            # Warm up
            for i = 1:10
                sum!(C, C.operand.operand)
            end

            results = []
            for i = 1:10
                mean!(C, C.operand.operand)
                push!(results, all(interior(C) .== 1))
            end

            @test mean(results) == 1.0
        end

        @testset "Allocating reductions [$arch_str]" begin
            @info "  Testing allocating reductions"

            grid = RectilinearGrid(arch, size = (2, 2, 2),
                                   x = (0, 2), y = (0, 2), z = (0, 2),
                                   topology = (Periodic, Periodic, Bounded))

            w = ZFaceField(grid)
            T = CenterField(grid)
            set!(T, trilinear)
            set!(w, trilinear)
            fill_halo_regions!(T, arch)
            fill_halo_regions!(w, arch)

            @compute Txyz = Field(Average(T, dims=(1, 2, 3)))
            @compute Txy = Field(Average(T, dims=(1, 2)))
            @compute Tx = Field(Average(T, dims=1))

            @compute wxyz = Field(Average(w, dims=(1, 2, 3)))
            @compute wxy = Field(Average(w, dims=(1, 2)))
            @compute wx = Field(Average(w, dims=1))

            # Mean
            @test Txyz[1, 1, 1] == mean(T)
            @test interior(Txy) == interior(mean(T, dims=(1, 2)))
            @test interior(Tx) == interior(mean(T, dims=1))

            @test wxyz[1, 1, 1] == mean(w)
            @test interior(wxy) == interior(mean(w, dims=(1, 2)))
            @test interior(wx) == interior(mean(w, dims=1))

            # Maximum and minimum
            @test maximum(T) == maximum(interior(T))
            @test minimum(T) == minimum(interior(T))

            for dims in (1, (1, 2))
                @test interior(minimum(T; dims)) == minimum(interior(T); dims)
                @test interior(minimum(T; dims)) == minimum(interior(T); dims)
            end
        end

        @testset "Conditional computation of averaged Fields [$(typeof(arch))]" begin
            @info "  Testing conditional computation of averaged Fields [$(typeof(arch))]"
            for FT in float_types
                grid = RectilinearGrid(arch, FT, size=(2, 2, 2), extent=(1, 1, 1))
                c = CenterField(grid)

                for dims in (1, 2, 3, (1, 2), (2, 3), (1, 3), (1, 2, 3))
                    C = Field(Average(c, dims=dims))

                    @test !isnothing(C.status)

                    # Test conditional computation
                    set!(c, 1)
                    compute_at!(C, FT(1)) # will compute
                    @test all(interior(C) .== 1)
                    @test C.status.time == FT(1)

                    set!(c, 2)
                    compute_at!(C, FT(1)) # will not compute because status == 1
                    @test C.status.time == FT(1)
                    @test all(interior(C) .== 1)

                    compute_at!(C, FT(2))
                    @test C.status.time == FT(2)
                    @test all(interior(C) .== 2)
                end
            end
        end
    end
end


text:
type Description
  name::String
  fn::Function
  depth::Int
  children
  function Description(name::String,fn::Function,depth::Int)
    new(name,fn,depth,{})
  end
end

type Test
  name::String
  desc::Description
  fn::Function
  iit::Bool
  function Test(name::String, desc::Description, fn::Function, iit::Bool)
    new(name,desc,fn,iit)
  end
  function Test(name::String, desc::Description, fn::Function)
    new(name,desc,fn,false)
  end
end

type Error
  test::Test
  err
  trace
end

type DescriptionError
  desc::Description
  err
  trace
end

noob() = begin end
depth = 0
iitOn = false
iitTests= Test[]
descriptions = Description("", noob,0)
passes = Test[]
errors = Union(Error,DescriptionError)[]
currDesc = descriptions

function run(test::Test)
  iitOn && !test.iit && return
  fin = 0
  try
    t0 = time()
    test.fn()
    fin = time() - t0
  catch e
    err = Error(test,e,catch_backtrace())
    failedReport(err)
    return push!(errors, err)
  end
  passReport(test, toMilis(fin))
  push!(passes,test)
end

function run{T<:Array}(suite::T)
  for test in suite
    run(test)
  end
end

function run(desc::Description)
  global currDesc, depth
  currDesc = desc
  depth = desc.depth
  try
    desc.fn()
    report(desc)
    run(desc.children)
  catch e
    err = DescriptionError(desc,e,catch_backtrace())
    failedReport(err)
    push!(errors, err)
  end
end

function run()
  global descriptions, passes, errors, depth, iitOn, iitTests, currDesc
  t0 = time()
  run(descriptions)
  summaryReport(passes,errors, toMilis(time() - t0))

  # Reset the global variables
  tests = None
  depth = 0
  iitOn = false
  empty!(descriptions.children)
  empty!(passes)
  empty!(errors)
  currDesc = descriptions

  println('\n',"\033[33mAll test finished running\n", RESET)
  return length(errors)
end

function describe(fn::Function, name::String)
  global currDesc, depth
  push!(currDesc.children, Description(name,fn,depth + 1))
end

function it(fn::Function, name::String)
  push!(currDesc.children, Test(name,currDesc,fn))
end

function iit(fn::Function, name::String)
  global iitOn = true
  push!(currDesc.children, Test(name,currDesc,fn,true))
end

_describe(args...) = begin end
_it(args...) = begin end


text:
@show ARGS

@show readline()
@show readline()

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
#    Problem :
#    *********
#    A banded function with semi-bandwidth 10 and
#    negative curvature near the starting point
#
#    Origonal SIF Source: problem 8 in
#    L. Luksan, C. Matonoha and J. Vlcek
#    Modified CUTE problems for sparse unconstrained optimization,
#    Technical Report 1081,
#    Institute of Computer Science,
#    Academy of Science of the Czech Republic
#
#    CURLY10.SIF classification OUR2-AN-V-0
#
#    Number of variables, n, is variable such that n ≥ 2
#    The initial iterate x0 in the source is not used,
#    rather, x0 is defined as it is the CUTEst enviroment
#
# Daniel Henderson, 08/2021

f = (x) -> begin
	n = lastindex(x)
    q = similar(x)
	fx = 0.0
	for i in 1:n
        q[i] = sum(x[j] for j = i:min(i + 10, n))
		fx  += q[i]^4 - 20q[i]^2 - 0.1q[i]
	end
    return fx
end

g! = (g, x) -> begin
	n = lastindex(x)
    q = similar(x)
	for i in 1:n
        q[i] = sum(x[j] for j = i:min(i + 10, n))
		for j = i:min(i + 10, n)
        	g[j] += 4q[i]^3 - 40q[i] - 0.1
		end
	end
    return g
end

fg! = (g, x) -> begin
	fx = 0.0
	n = lastindex(x)
    q = similar(x)
	for i in 1:n
        q[i] = sum(x[j] for j = i:min(i + 10, n))
		fx  += q[i]^4 - 20q[i]^2 - 0.1q[i]
		for j = i:min(i + 10, n)
			g[j] += 4q[i]^3 - 40q[i] - 0.1
		end
	end
    return fx, g
end

init = (n::Int=10000) -> begin
	n < 2 && @warn("CURLY10: number of variables must be ≥ 2")
	n = max(n, 2)
	x0 = [0.0001*i/(n+ 1) for i in 1:n]
    return n, x0
end

TestSet["CURLY10"] = UncProgram("CURLY10", f, g!, fg!, init)


text:
using PyCall

@pyimport selenium.webdriver as wd

function gen_wrapper(class, objarg, doc_section; jl_wrap=(x,y)->y, curry=false)
    methods = doc_section[:find_elements_by_css_selector](".method")
    attr = doc_section[:find_elements_by_css_selector](".attribute")

    ms = map(f->gen_method(class, f, objarg; jl_wrap=jl_wrap, curry=curry), methods)
    as = map(f->gen_method(class, f, objarg, true; jl_wrap=jl_wrap, curry=curry), attr)
    names = vcat(map(first, ms), map(first, as))
    code = vcat(map(last, ms), map(last, as))
    "export " * join(names, ",") *"\n\n"* join(code, "\n\n")
end

find(x, q) = x[:find_element_by_css_selector](q)
findall(x, q) = x[:find_elements_by_css_selector](q)

cleanup_quotes(str) = replace(str, r"“|”|‘|’|'", "\"")

function strip_default(arg)
    idx = findfirst(arg, '=')
    maybe_splat(idx == 0 ? arg : arg[1:idx-1])
end

maybe_splat(x) =
    if startswith(x, "**")
        "; $(x[3:end])..."
    elseif startswith(x, "*")
        "$(x[2:end])..."
    else
        x
    end

function rewrite_default(arg)
    idx = findfirst(arg, '=')
    if idx == 0
        return maybe_splat(arg)
    else
        k = arg[1:idx-1]
        argval = rewrite_val(arg[idx+1:end])
        return "$k = $argval"
    end
end

function rewrite_val(v)
    if v == "None"
        "nothing"
    elseif v in ("False", "True")
        lowercase(v)
    else
        cleanup_quotes(v)
    end
end

function gen_method(class, el, objarg, attr=false; jl_wrap=(fn, x)->x, curry=false)
    fn = find(el, ".descname")[:text]
    args = map(x->x[:text], findall(el, "em"))
    doc = cleanup_quotes(find(el, "dd")[:text])
    expr = if attr
        """$objarg.o[:$fn]"""
    else
        """$objarg.o[:$fn]($(join(map(strip_default, args), ", ")))"""
    end

    expr = jl_wrap(fn, expr)
    method2 = if curry
        """

        $fn($(join(map(rewrite_default, args), ", "))) = $objarg -> $expr
        """
    else
        ""
    end

    fn => ("""
    \"\"\"
    $doc
    \"\"\"
    $fn($objarg::$class, $(join(map(rewrite_default, args), ", "))) = $expr
    """ * method2)

end

code = ""
#### Remote WebDriver

driver = wd.PhantomJS("phantomjs")

driver[:get]("http://selenium-python.readthedocs.io/api.html")

### Driver

function wrap_elem(fn, expr)
    if startswith(fn, "find_elements")
        "map(WebElement, $expr)"
    elseif startswith(fn, "find_element")
        "WebElement($expr)"
    else
        expr
    end
end

remote_section =
    driver[:find_element_by_id]("module-selenium.webdriver.remote.webdriver")

code *= gen_wrapper(:Driver, :driver, remote_section, jl_wrap=wrap_elem)

#### WebElement

webel_section =
    driver[:find_element_by_id]("module-selenium.webdriver.remote.webelement")

code *= gen_wrapper(:WebElement, :elem, webel_section; jl_wrap=wrap_elem)

#### Action chains

action_chain_section =
    driver[:find_element_by_id]("module-selenium.webdriver.common.action_chains")

function chain_wrapper(fn, x)
    if fn == "perform"
        x
    else
        "ActionChain($x)"
    end
end

code *= gen_wrapper(:ActionChain, :chain, action_chain_section; jl_wrap=chain_wrapper, curry=true)

### Special keys

keys_section =
    driver[:find_element_by_id]("module-selenium.webdriver.common.keys")

consts = map(findall(keys_section, ".attribute")) do k
    find(k, ".descname")[:text]
end

exports = "export " * join(consts, ",")

constdefs = join(["const $x = ks.Keys[:$x]" for x in consts], "\n")
code *= """
export Keys
module Keys
using PyCall

$exports

@pyimport selenium.webdriver.common.keys as ks
$constdefs
end
"""

println(code)


text:
@testset "SparseVoxelGrid tests" begin
    # Create points for a uniformly spaced 4x4x4 grid
    x = collect(linspace(0.0, 3.0, 3))
    cnt = 1
    points2d = zeros(2, 27)
    points3d = zeros(3, 27)
    for i in x, j in x, k in x
        points2d[1, cnt] = points3d[1, cnt] = i
        points2d[2, cnt] = points3d[2, cnt] = j
        points3d[3, cnt] = k
        cnt += 1
    end

    @testset "Rasterizer" begin
        # Test standard uniformly sized voxels
        d2d = rasterize_points(points2d, 1.0)
        [@test length(r) == 3 for r in values(d2d)]
        d3d = rasterize_points(points3d, 1.0)
        [@test length(r) == 3 for r in values(d3d)]

        @test length(keys(d2d)) == 9
        @test d2d[(0x00000000,0x00000000)] == UInt32[1, 2, 3]

        @test_throws DimensionMismatch rasterize_points([SVector{4, Float32}(rand(4)) for i= 1:4], 1.0)

        io = IOBuffer()
        show(io, d3d)
    end

    @testset "Voxelization" begin
        # Test standard uniformly sized voxels
        voxel_size = 1.0
        grid = SparseVoxelGrid(points3d, voxel_size)
        @test isempty(grid) == false
        @test haskey(grid, (0,0,0)) == true
        @test length(grid) == 27
        @test length(collect(grid)) == 27
        [@test length(collect(voxel)) == 1 for voxel in grid]
        @test voxel_center(grid, (1, 1, 1)) == SVector{3,Float64}(0.5, 0.5, 0.5)
        @test length(collect(SparseVoxelGrid(points3d, (2.0, 2.5, 4.0)))) == 4
        @test length(collect(SparseVoxelGrid(points3d, SVector(2.0, 2.5, 4.0)))) == 4

        vector_points = [SVector{3, Float64}(points3d[:,i]) for i = 1:size(points3d, 2)]
        grid = SparseVoxelGrid(vector_points, SVector(1, 1, 1))
        @test length(grid) == 27

        grid = SparseVoxelGrid(vector_points, 2.0)
        @test length(grid) == 8
    end

    @testset "Neighbouring voxel" begin
        radius = 1
        grid = SparseVoxelGrid(points3d, 2.0)
        for voxel in grid
            # Test using anonymous function method
            in_cuboid(grid, voxel, radius) do voxel
                @test haskey(grid, voxel.id)
            end
        end

        grid = SparseVoxelGrid(points3d, 1.5)
        voxel_list = [(1,1,1), (1,2,1), (2,2,2), (10,10,10)]
        @test length(collect(in_cuboid(grid, voxel_list[1], 2))) == 26
        @test length(collect(in_cuboid(grid, voxel_list[2], 2))) == 26
        @test length(collect(in_cuboid(grid, voxel_list[3], 2))) == 26
        @test length(collect(in_cuboid(grid, voxel_list[4], 2))) == 0
        @test length(collect(in_cuboid(grid, collect(grid)[1], 2))) == 26

        # Test show methods
        io = IOBuffer()
        voxel = collect(grid)[1]
        cuboid = in_cuboid(grid, (0,0,0), radius)
        show(io, grid)
        show(io, voxel)
        show(io, cuboid)
    end
end


text:
using Gadfly, DataArrays, RDatasets

plot(dataset("ggplot2", "diamonds"), x="Price", color="Cut", Geom.density)


text:
function d2tduftau(t,
		   fft_u,  fft_v,
		   fft_du, fft_dv,
		   A1, A2,
		   matr, conjmatr,
		   sigma, llambda, Ntaumm)

    sigma = 1

    u     = ifft(exp(1im * t * A1) * fft_u * matr)
    v     = ifft(exp(1im * t * A1) * fft_v * matr)
    dtu   = ifft(exp(1im * t * A1) * (1im * A1) * fft_u) * matr
    dtv   = ifft(exp(1im * t * A1) * (1im * A1) * fft_v) * matr
    d2tu  = ifft(exp(1im * t * A1) * (-A1 .^ 2) * fft_u) * matr
    d2tv  = ifft(exp(1im * t * A1) * (-A1 .^ 2) * fft_v) * matr

    du    = ifft(exp(1im * t * A1) * fft_du * matr)
    dv    = ifft(exp(1im * t * A1) * fft_dv * matr)
    dtdu  = ifft(exp(1im * t * A1) * (1im * A1) * fft_du) * matr
    dtdv  = ifft(exp(1im * t * A1) * (1im * A1) * fft_dv) * matr
    d2tdu = ifft(exp(1im * t * A1) * (-A1 .^ 2) * fft_du) * matr
    d2tdv = ifft(exp(1im * t * A1) * (-A1 .^ 2) * fft_dv) * matr

    z     = (u + conj(v)) / 2
    dtz   = (dtu + conj(dtv)) / 2
    d2tz  = (d2tu + conj(d2tv)) / 2
    dz    = (du + conj(dv)) / 2
    dtdz  = (dtdu + conj(dtdv)) / 2
    d2tdz = (d2tdu + conj(d2tdv)) / 2

    fz1 = (  2 * dz  * conj(z)  * d2tz
	   + 2 * z   * conj(dz) * d2tz
	   + 2 * z   * conj(z)  * d2tdz
	   + 4 * dtz * dtdz     * conj(z)
	   + 2 * dtz .^ 2 * conj(dz)
	   + 4 * dtdz * conj(dtz)  * z
	   + 4 * dtz  * conj(dtdz) * z
	   + 4 * dtz  * conj(dtz)  * dz
	   + 2 * z * dz * conj(d2tz)
	   +     z .^ 2 * conj(d2tdz))

    champu1 = -1im * llambda * A2 .* exp(-1im * t * A1) .* fft(conjmatr .* fz1)
    champv1 = -1im * llambda * A2 .* exp(-1im * t * A1) .* fft(conjmatr .* conj(fz1))

    fz1 = (  2 * dtz * conj(z) * dz
	   + 2 * z * conj(dtz) * dz
	   + 2 * abs(z) .^ 2 * dtdz
	   + 2 * z * dtz * conj(dz)
	   + z .^ 2 * conj(dtdz))

    champu2 = -1im * llambda .* A2 .* exp.(-1im * t * A1) .* (-1im * A1) .* fft(conjmatr .* fz1)
    champv2 = -1im * llambda .* A2 .* exp.(-1im * t * A1) .* (-1im * A1) .* fft(conjmatr .* conj(fz1))

    fz1 = 2 * abs.(z) .^ 2 * dz .+ z .^ 2 * conj(dz)
    champu3 = -1im * llambda * A2 .* exp.(-1im * t .* A1) .* (-A1 .^ 2) .* fft(conjmatr .* fz1)
    champv3 = -1im * llambda * A2 .* exp.(-1im * t .* A1) .* (-A1 .^ 2) .* fft(conjmatr .* conj(fz1))

    champu = champu1 .+ 2 * champu2 .+ champu3
    champv = champv1 .+ 2 * champv2 .+ champv3

    champu, champv

end


text:
struct Status
    name::String
    duration::Int
end

max_stars_repo_name: Flobii/discord-bot

max_issues_repo_name: Flobii/discord-bot

max_forks_repo_name: Flobii/discord-bot


text:
################################################################################
#
#	ABSTRACT TYPE
#
#   ReciprocalPoint{D}
#   --> D is the dimension of embedding space
#
#   FILE CONTAINS
#       - abstract type definition
#       - interface definition
#       - interface testing
#
################################################################################

################################################################################
#
#   ABSTRACT TYPE DEFINITION
#
################################################################################
abstract type AbstractReciprocalPoint{D} <: AbstractSite{String, D} end

# export the type
export AbstractReciprocalPoint

################################################################################
#
#	INTERFACING / ACCESSING RECIPROCAL POINTS
#	(functions have to be overwritten by concrete types)
#
################################################################################

# default constructor interface
# used for creation of new reciprocal points
function newReciprocalPoint(
            :: Type{R},
            point   :: Vector{<:Real},
            label   :: AbstractString
        ) :: R where {R<:AbstractReciprocalPoint{D} where D}

    # print an error because implementation for concrete type is missing
    error(  "not implemented function 'newReciprocalPoint' for concrete reciprocal point type " *
            string(R) * " in reciprocal space dimension " * string(D) )
end

# export the newSite interface
export newReciprocalPoint

# interface of AbstractSite
# so it DOES NOT HAVE TO BE OVERWRITTEN by concrete recipocal point
function newSite(
            :: Type{R},
            point   :: Vector{<:Real},
            label   :: AbstractString
        ) :: R where {D, R <: AbstractReciprocalPoint{D}}

    # return the respective function
    return newReciprocalPoint(R, point, label)
end

# export the newSite interface
export newSite

# get label
function label(
            p :: AbstractReciprocalPoint{D}
        ) :: String where {D}

    # print an error because implementation for concrete type is missing
    error("not implemented interface function 'label' for reciprocal point type " * string(typeof(p)))
end
# set label
function label!(
            p :: AbstractReciprocalPoint{D},
            l :: AbstractString
        ) where {D}

    # print an error because implementation for concrete type is missing
    error("not implemented interface function 'label!' for reciprocal point type " * string(typeof(p)))
end

# export the label interface
export label, label!

# get point
function point(
            p :: AbstractReciprocalPoint{D}
        ) :: Vector{Float64} where {D}

    # print an error because implementation for concrete type is missing
    error("not implemented interface function 'point' for reciprocal point type " * string(typeof(p)))
end
# set point
function point!(
            p :: AbstractReciprocalPoint{D},
            r :: Vector{<:Real}
        ) where {D}

    # print an error because implementation for concrete type is missing
    error("not implemented interface function 'point!' for reciprocal point type " * string(typeof(p)))
end

# export the point interface
export point, point!

# SIMILAR FUNCTION (can be overwritten but does not have to be overwritten)

# without new parameters
function similar(
            p :: R
        ) :: R where {D,R<:AbstractReciprocalPoint{D}}

    # return a new site object
    return newReciprocalPoint(R, deepcopy(point(p)), deepcopy(label(p)))
end
# with new parameters
function similar(
            p :: R,
            r :: Vector{<:Real},
            l :: AbstractString
        ) :: R where {D,R<:AbstractReciprocalPoint{D}}

    # create a new site object
    p_new = similar(p)
    # set parameters
    point!(p_new, r)
    label!(p_new, l)
    # return the new object
    return p_new
end

# export the similar interface
export similar

# LaTeX functionality

# get label
function labelLaTeX(
            p :: AbstractReciprocalPoint{D}
        ) :: String where {D}

    # return the usual label
    return label(p)
end
# set label
function labelLaTeX!(
            p :: AbstractReciprocalPoint{D},
            l :: AbstractString
        ) where {D}

    # set the usual label
    return label!(p,l)
end

# export the label interface
export labelLaTeX, labelLaTeX!

# new point with different LaTeX label
function newReciprocalPoint(
            :: Type{R},
            point      :: Vector{<:Real},
            label      :: AbstractString,
            labelLaTeX :: AbstractString
        ) :: R where {R<:AbstractReciprocalPoint{D} where D}

    # call the usual constructor by neglecting the LaTeX label
    return newReciprocalPoint(R, point, label)
end

# export the newSite interface
export newReciprocalPoint

################################################################################
#
#	OVERWRITING BASE.SHOW FOR RECIPROCAL POINTS
#
################################################################################
function show(io::IO, p::R) where {R<:AbstractReciprocalPoint{D} where D}
    print(io, R, " @", point(p), ": (\"", label(p), "\" / L\"", labelLaTeX(p), "\")")
end


text:
module TestIndexing
    using Base.Test
    using DataArrays
    using DataFrames

    #
    # DataVector and PooledDataVector Indexing
    #

    dv = DataArray([1, 2])
    pdv = PooledDataArray([1, 2])

    adv_indices = {1,
                   1.0,
                   1:2,
                   1:1:2,
                   1.0:2.0,
                   1.0:1.0:2.0,
                   [1, 2],
                   [1.0, 2.0],
                   [true, true],
                   trues(2),
                   @data([1, 2]),
                   @data([1.0, 2.0]),
                   @data([true, true])}

    single_adv_indices = {1,
                          1.0}

    multi_adv_indices = {1:2,
                         1:1:2,
                         1.0:2.0,
                         1.0:1.0:2.0,
                         [1, 2],
                         [1.0, 2.0],
                         [true, true],
                         trues(2),
                         @data([1, 2]),
                         @data([1.0, 2.0]),
                         @data([true, true])}

    # avd[Index]
    for adv_index in adv_indices
        dv[adv_index]
        pdv[adv_index]
    end

    # avd[SingleIndex] = Single Value
    for single_adv_index in single_adv_indices
        dv[single_adv_index] = NA
        dv[single_adv_index] = 1
        dv[single_adv_index] = 1.0
        pdv[single_adv_index] = NA
        pdv[single_adv_index] = 1
        pdv[single_adv_index] = 1.0
    end

    # avd[MultiIndex] = Multiple Values
    for multi_adv_index in multi_adv_indices
        dv[multi_adv_index] = [1, 2]
        dv[multi_adv_index] = [1.0, 2.0]
        dv[multi_adv_index] = 1:2
        dv[multi_adv_index] = 1:1:2
        dv[multi_adv_index] = 1.0:2.0
        dv[multi_adv_index] = 1.0:1.0:2.0
        pdv[multi_adv_index] = [1, 2]
        pdv[multi_adv_index] = [1.0, 2.0]
        pdv[multi_adv_index] = 1:2
        pdv[multi_adv_index] = 1:1:2
        pdv[multi_adv_index] = 1.0:2.0
        pdv[multi_adv_index] = 1.0:1.0:2.0
    end

    # avd[MultiIndex] = Single Value
    for multi_adv_index in multi_adv_indices
        dv[multi_adv_index] = 1
        dv[multi_adv_index] = 1.0
        pdv[multi_adv_index] = 1
        pdv[multi_adv_index] = 1.0
    end

    #
    # DataFrame indexing
    #

    df = DataFrame(A = 1:2, B = 3:4)

    row_indices = {1,
                   1.0,
                   1:2,
                   [1, 2],
                   [1.0, 2.0],
                   [true, false],
                   trues(2),
                   @data([1, 2]),
                   @data([1.0, 2.0]),
                   @data([true, false]),
                   :(A .== 1)}

    column_indices = {1,
                      1.0,
                      "A",
                      1:2,
                      [1, 2],
                      [1.0, 2.0],
                      [true, false],
                      trues(2),
                      ["A", "B"],
                      @data([1, 2]),
                      @data([1.0, 2.0]),
                      @data([true, false]),
                      @data(["A", "B"]),
                      :(names(_DF) .== "B")}

    #
    # getindex()
    #

    for column_index in column_indices
        df[column_index]
    end

    for row_index in row_indices
        for column_index in column_indices
            df[row_index, column_index]
        end
    end

    #
    # setindex!()
    #

    for column_index in column_indices
        df[column_index] = df[column_index]
    end

    for row_index in row_indices
        for column_index in column_indices
            df[row_index, column_index] = df[row_index, column_index]
        end
    end

    #
    # Broadcasting assignments
    #

    for column_index in column_indices
        df[column_index] = NA
        df[column_index] = 1
        df[column_index] = 1.0
        df[column_index] = "A"
        df[column_index] = DataArray([1 + 0im, 2 + 1im])
    end

    # Only assign into columns for which new value is type compatible
    for row_index in row_indices
        for column_index in column_indices
            df[row_index, column_index] = NA
            df[row_index, column_index] = 1
            df[row_index, column_index] = 1.0
            df[row_index, column_index] = "A"
            df[row_index, column_index] = DataArray([1 + 0im, 2 + 1im])
        end
    end
end


text:
## Iterator protocol for `Group`

Base.eltype(::Type{G}) where {G<:Group} =
    throw(InterfaceNotImplemented(:Iteration, "Base.eltype(::Type{$G})"))
Base.iterate(G::Group) =
    throw(InterfaceNotImplemented(:Iteration, "Base.iterate(::$(typeof(G)))"))
Base.iterate(G::Group, state) = throw(
    InterfaceNotImplemented(:Iteration, "Base.iterate(::$(typeof(G)), state)"),
)

#=
TODO: This is a convention; cannot define length unless Base.IteratorSize returns HasLength or HasShape!

You need to override the default
> `Base.IteratorSize(::Type{MyGroup}) = Base.HasLength()`
by
> `Base.IteratorSize(::Type{MyGroup}) = Base.IsInfinite()`
if your group is always infinite, or by
> `Base.IteratorSize(::Type{MyGroup}) = Base.SizeUnknown()`
if finiteness can not be easily established.

otherwise define

Base.length(G::MyGroup) = order(Int, G)

which is "best effort", cheap computation of length of the group iterator. For practical reasons the largest group you could iterate over in your lifetime is of order ~factorial(19) (assuming 10ns per element).
=#

## `Group` interface

### Obligatory methods for `Group`
@doc Markdown.doc"""
    one(G::Group)
Return the identity of the group.
"""
Base.one(G::Group) =
    throw(InterfaceNotImplemented(:Group, "Base.one(::$(typeof(G)))"))

@doc Markdown.doc"""
    order([BigInt, ]G::Group)
    order(I::Type{<:Integer}, g::Group)
Return the order of `g` as an instance of `I`.

Only arbitrary size integers are required to return mathematically correct answer.
"""
AbstractAlgebra.order(::Type{<:Integer}, G::Group) = throw(
    InterfaceNotImplemented(
        :Group,
        "GroupsCore.order(::Type{<:Integer}, ::$(typeof(G)))",
    ),
)
@doc Markdown.doc"""
    gens(G::Group)

Return a random-accessed collection of generators of `G`.

If a group does not come with a generating set (or it may be prohibitively expensive to compute), one needs to alter `GroupsCore.hasgens(::Group) = false`.
"""
AbstractAlgebra.gens(G::Group) =
    throw(InterfaceNotImplemented(:Group, "GroupsCore.gens(::$(typeof(G)))"))

function Base.rand(
    rng::Random.AbstractRNG,
    rs::Random.SamplerTrivial{G},
) where {G<:Group}
    throw(
        InterfaceNotImplemented(
            :Random,
            "Base.rand(::Random.AbstractRNG, ::Random.SamplerTrivial{$G}))",
        ),
    )
end

### Default implementations for `Group`

function Base.isfinite(G::Group)
    IS = Base.IteratorSize(G)
    IS isa Base.HasLength && return true
    IS isa Base.HasShape && return true
    # else : IS isa (Base.SizeUnknown, Base.IsInfinite, ...)
    return false
end

hasgens(G::Group) = true

AbstractAlgebra.order(G::Group) = order(BigInt, G)

@doc Markdown.doc"""
    elem_type(parent_type)
Given the type of a parent object return the type of its elements.
"""
AbstractAlgebra.elem_type(T::Type{<:Group}) = eltype(T)

function AbstractAlgebra.gens(G::Group, i::Integer)
    hasgens(G) && return gens(G)[i]
    # TODO: throw something more specific
    throw(
        "Group does not seem to have generators. Did you alter `hasgens(::$(typeof(G)))`?",
    )
end

function AbstractAlgebra.ngens(G::Group)
    hasgens(G) && return length(gens(G))
    # TODO: throw something more specific
    throw(
        "Group does not seem to have generators. Did you alter `hasgens(::$(typeof(G)))`?",
    )
end

pseudo_rand(G::Group, args...) = pseudo_rand(Random.default_rng(), G, args...)
pseudo_rand(rng::Random.AbstractRNG, G::Group, args...) = rand(rng, G, args...)


text:
using RTS

RTS.julia_main([""])