text:
#####
##### MsgPack conversions for Base types
#####

MsgPack.msgpack_type(::Type{Nanosecond}) = MsgPack.IntegerType()
MsgPack.from_msgpack(::Type{Nanosecond}, x::Integer) = Nanosecond(x)
MsgPack.to_msgpack(::MsgPack.IntegerType, x::Nanosecond) = x.value

MsgPack.msgpack_type(::Type{VersionNumber}) = MsgPack.StringType()
MsgPack.from_msgpack(::Type{VersionNumber}, x::String) = VersionNumber(x[2:end])
MsgPack.to_msgpack(::MsgPack.StringType, x::VersionNumber) = string('v', x)

MsgPack.msgpack_type(::Type{UUID}) = MsgPack.StringType()
MsgPack.from_msgpack(::Type{UUID}, x::String) = UUID(x)
MsgPack.to_msgpack(::MsgPack.StringType, x::UUID) = string(x)

MsgPack.msgpack_type(::Type{DataType}) = MsgPack.StringType()
function MsgPack.from_msgpack(::Type{DataType}, x::String)
    return julia_type_from_onda_sample_type(x)
end
function MsgPack.to_msgpack(::MsgPack.StringType, T::DataType)
    return onda_sample_type_from_julia_type(T)
end

#####
##### Julia DataType <--> Onda `sample_type` string
#####

function julia_type_from_onda_sample_type(t::AbstractString)
    t == "int8" && return Int8
    t == "int16" && return Int16
    t == "int32" && return Int32
    t == "int64" && return Int64
    t == "uint8" && return UInt8
    t == "uint16" && return UInt16
    t == "uint32" && return UInt32
    t == "uint64" && return UInt64
    return error("sample type ", t, " is not supported by Onda")
end

function onda_sample_type_from_julia_type(T::Type)
    T === Int8 && return "int8"
    T === Int16 && return "int16"
    T === Int32 && return "int32"
    T === Int64 && return "int64"
    T === UInt8 && return "uint8"
    T === UInt16 && return "uint16"
    T === UInt32 && return "uint32"
    T === UInt64 && return "uint64"
    return error("sample type ", T, " is not supported by Onda")
end

#####
##### annotations
#####

"""
    Annotation <: AbstractTimeSpan

A type representing an individual Onda annotation object. Instances contain
the following fields, following the Onda specification for annotation objects:

  - `value::String`
  - `start_nanosecond::Nanosecond`
  - `stop_nanosecond::Nanosecond`
"""
struct Annotation <: AbstractTimeSpan
    value::String
    start_nanosecond::Nanosecond
    stop_nanosecond::Nanosecond
    function Annotation(value::AbstractString, start::Nanosecond, stop::Nanosecond)
        _validate_timespan(start, stop)
        return new(value, start, stop)
    end
end

MsgPack.msgpack_type(::Type{Annotation}) = MsgPack.StructType()

function Annotation(value, span::AbstractTimeSpan)
    return Annotation(value, first(span), last(span))
end

Base.first(annotation::Annotation) = annotation.start_nanosecond

Base.last(annotation::Annotation) = annotation.stop_nanosecond

#####
##### signals
#####

"""
    Signal

A type representing an individual Onda signal object. Instances contain
the following fields, following the Onda specification for signal objects:

  - `channel_names::Vector{Symbol}`
  - `start_nanosecond::Nanosecond`
  - `stop_nanosecond::Nanosecond`
  - `sample_unit::Symbol`
  - `sample_resolution_in_unit::Float64`
  - `sample_offset_in_unit::Float64`
  - `sample_type::DataType`
  - `sample_rate::Float64`
  - `file_extension::Symbol`
  - `file_options::Union{Nothing,Dict{Symbol,Any}}`
"""
Base.@kwdef struct Signal
    channel_names::Vector{Symbol}
    start_nanosecond::Nanosecond
    stop_nanosecond::Nanosecond
    sample_unit::Symbol
    sample_resolution_in_unit::Float64
    sample_offset_in_unit::Float64
    sample_type::DataType
    sample_rate::Float64
    file_extension::Symbol
    file_options::Union{Nothing,Dict{Symbol,Any}}
    function Signal(channel_names, start_nanosecond, stop_nanosecond, sample_unit,
                    sample_resolution_in_unit, sample_offset_in_unit, sample_type,
                    sample_rate, file_extension, file_options)
        _validate_timespan(start_nanosecond, stop_nanosecond)
        return new(channel_names, start_nanosecond, stop_nanosecond, sample_unit,
                   sample_resolution_in_unit, sample_offset_in_unit, sample_type,
                   sample_rate, file_extension, file_options)
    end
end

function Base.:(==)(a::Signal, b::Signal)
    return all(name -> getfield(a, name) == getfield(b, name), fieldnames(Signal))
end

MsgPack.msgpack_type(::Type{Signal}) = MsgPack.StructType()

function is_valid(signal::Signal)
    return is_lower_snake_case_alphanumeric(string(signal.sample_unit)) &&
           all(n -> is_lower_snake_case_alphanumeric(string(n), ('-', '.')),
               signal.channel_names) &&
           onda_sample_type_from_julia_type(signal.sample_type) isa AbstractString
end

function file_option(signal::Signal, name, default)
    signal.file_options isa Dict && return get(signal.file_options, name, default)
    return default
end

"""
    signal_from_template(signal::Signal;
                         channel_names=signal.channel_names,
                         start_nanosecond=signal.start_nanosecond,
                         stop_nanosecond=signal.stop_nanosecond,
                         sample_unit=signal.sample_unit,
                         sample_resolution_in_unit=signal.sample_resolution_in_unit,
                         sample_offset_in_unit=signal.sample_offset_in_unit,
                         sample_type=signal.sample_type,
                         sample_rate=signal.sample_rate,
                         file_extension=signal.file_extension,
                         file_options=signal.file_options)

Return a `Signal` where each field is mapped to the corresponding keyword argument.
"""
function signal_from_template(signal::Signal; channel_names=signal.channel_names,
                              start_nanosecond=signal.start_nanosecond,
                              stop_nanosecond=signal.stop_nanosecond,
                              sample_unit=signal.sample_unit,
                              sample_resolution_in_unit=signal.sample_resolution_in_unit,
                              sample_offset_in_unit=signal.sample_offset_in_unit,
                              sample_type=signal.sample_type,
                              sample_rate=signal.sample_rate,
                              file_extension=signal.file_extension,
                              file_options=signal.file_options)
    return Signal(channel_names, start_nanosecond, stop_nanosecond, sample_unit,
                  sample_resolution_in_unit, sample_offset_in_unit, sample_type,
                  sample_rate, file_extension, file_options)
end

"""
    channel(signal::Signal, name::Symbol)

Return `i` where `signal.channel_names[i] == name`.
"""
function channel(signal::Signal, name::Symbol)
    return findfirst(isequal(name), signal.channel_names)
end

"""
    channel(signal::Signal, i::Integer)

Return `signal.channel_names[i]`.
"""
channel(signal::Signal, i::Integer) = signal.channel_names[i]

"""
    channel_count(signal::Signal)

Return `length(signal.channel_names)`.
"""
channel_count(signal::Signal) = length(signal.channel_names)

"""
    span(signal::Signal)

Return `TimeSpan(signal.start_nanosecond, signal.stop_nanosecond)`.
"""
span(signal::Signal) = TimeSpan(signal.start_nanosecond, signal.stop_nanosecond)

"""
    duration(signal::Signal)

Return `duration(span(signal))`.
"""
duration(signal::Signal) = duration(span(signal))

"""
    sample_count(signal::Signal)

Return the number of multichannel samples that fit within `duration(signal)`
given `signal.sample_rate`.
"""
function sample_count(signal::Signal)
    return index_from_time(signal.sample_rate, duration(signal)) - 1
end

"""
    sizeof_samples(signal::Signal)

Returns the expected size (in bytes) of the encoded `Samples` object corresponding
to the entirety of `signal`:

    sample_count(signal) * channel_count(signal) * sizeof(signal.sample_type)
"""
function sizeof_samples(signal::Signal)
    return sample_count(signal) * channel_count(signal) * sizeof(signal.sample_type)
end

#####
##### recordings
#####

"""
    Recording

A type representing an individual Onda recording object. Instances contain
the following fields, following the Onda specification for recording objects:

  - `signals::Dict{Symbol,Signal}`
  - `annotations::Set{Annotation}`
"""
struct Recording
    signals::Dict{Symbol,Signal}
    annotations::Set{Annotation}
end

function Base.:(==)(a::Recording, b::Recording)
    return all(name -> getfield(a, name) == getfield(b, name), fieldnames(Recording))
end

MsgPack.msgpack_type(::Type{Recording}) = MsgPack.StructType()

"""
    annotate!(recording::Recording, annotation::Annotation)

Returns `push!(recording.annotations, annotation)`.
"""
function annotate!(recording::Recording, annotation::Annotation)
    return push!(recording.annotations, annotation)
end

"""
    duration(recording::Recording)

Returns `maximum(s -> s.stop_nanosecond, values(recording.signals))`; throws an
`ArgumentError` if `recording.signals` is empty.
"""
function duration(recording::Recording)
    isempty(recording.signals) &&
        throw(ArgumentError("`duration(recording)` is not defined if `isempty(recording.signals)`"))
    return maximum(s -> s.stop_nanosecond, values(recording.signals))
end

"""
    set_span!(recording::Recording, name::Symbol, span::AbstractTimeSpan)

Replace `recording.signals[name]` with a copy that has the `start_nanosecond`
and `start_nanosecond` fields set to match the provided `span`. Returns the
newly constructed `Signal` instance.
"""
function set_span!(recording::Recording, name::Symbol, span::AbstractTimeSpan)
    signal = signal_from_template(recording.signals[name]; start_nanosecond=first(span),
                                  stop_nanosecond=last(span))
    recording.signals[name] = signal
    return signal
end

"""
    set_span!(recording::Recording, span::TimeSpan)

Return `Dict(name => set_span!(recording, name, span) for name in keys(recording.signals))`
"""
function set_span!(recording::Recording, span::AbstractTimeSpan)
    return Dict(name => set_span!(recording, name, span)
                for name in keys(recording.signals))
end

#####
##### reading/writing `recordings.msgpack.zst`
#####

struct Header
    onda_format_version::VersionNumber
    ordered_keys::Bool
end

MsgPack.msgpack_type(::Type{Header}) = MsgPack.StructType()

function read_recordings_file(path)
    file_path = joinpath(path, "recordings.msgpack.zst")
    bytes = zstd_decompress(read(file_path))
    io = IOBuffer(bytes)
    # `0x92` is the MessagePack byte prefix for 2-element array
    read(io, UInt8) == 0x92 || error("recordings.msgpack.zst has bad byte prefix")
    header = MsgPack.unpack(io, Header)
    if !is_supported_onda_format_version(header.onda_format_version)
        @warn("attempting to load `Dataset` with unsupported Onda version",
              supported = ONDA_FORMAT_VERSION, attempting = header.onda_format_version)
        @warn("consider upgrading old datasets via `Onda.upgrade_onda_format_from_v0_2_to_v0_3!`")
    end
    strict = header.ordered_keys ? (Recording,) : ()
    recordings = MsgPack.unpack(io, Dict{UUID,Recording}; strict=strict)
    return header, recordings
end

function write_recordings_file(path, header::Header, recordings::Dict{UUID,Recording})
    file_path = joinpath(path, "recordings.msgpack.zst")
    backup_file_path = joinpath(path, "_recordings.msgpack.zst.backup")
    isfile(file_path) && mv(file_path, backup_file_path)
    io = IOBuffer()
    MsgPack.pack(io, [header, recordings])
    bytes = zstd_compress(resize!(io.data, io.size))
    write(file_path, bytes)
    rm(backup_file_path; force=true)
    return nothing
end


text:
push!(LOAD_PATH,"../src/")
using Documenter, Diagonalizations

makedocs(
   sitename="Diagonalizations",
   authors="Marco Congedo",
   modules=[Diagonalizations],
   pages =
   [
      "index.md",
      "Diagonalizations" => "Diagonalizations.md",
      "Filters" => Any[
         "One dataset (m=1)" => Any[
            "PCA" => "pca.md",
            "Whitening" => "whitening.md",
            "CSP" => "csp.md",
            "CSTP" => "cstp.md",
            "AJD" => "ajd.md",
         ],
         "Two datasets (m=2)" => Any[
            "MCA" => "mca.md",
            "CCA" => "cca.md",
         ],
         "Several datasets (m>2)" => Any[
            "gMCA" => "gmca.md",
            "gCCA" => "gcca.md",
            "mAJD" => "majd.md",
         ],
      ],
      "Tools" => "tools.md",
      "Algorithms" => "algorithms.md",
   ]
)

deploydocs(
   # root
   # target = "build", # add this folder to .gitignore!
   repo = "github.com/Marco-Congedo/Diagonalizations.jl.git",
   # branch = "gh-pages",
   # osname = "linux",
   # deps = Deps.pip("pygments", "mkdocs"),
   # devbranch = "dev",
   # devurl = "dev",
   # versions = ["stable" => "v^", "v#.#", devurl => devurl],
)

ext: jl


text:
module MainTests

end # module


text:
"""
This can contain a maximum of `length(par)` parenting-relations
par is an array of `Int`, which is the index of the parent node.
"""
mutable struct DisjointSet
    par::Vector{Int}
    function DisjointSet(size)
        x=[i for i in 1:size]
        return new(x)
    end
end

"""
Find the ancestor of node `x`.
"""
function find(set::DisjointSet,x::Int)::Int
    if set.par[x]==x return x
    else return set.par[x]=find(set,set.par[x])
    end
end

function Base.merge!(set::DisjointSet,x::Int,y::Int)
    x=find(set,x)
    y=find(set,y)
    set.par[x]=y
end

max_stars_repo_name: KohRongSoon/Julia

max_issues_repo_name: KohRongSoon/Julia

max_forks_repo_name: KohRongSoon/Julia


text:
#!/usr/bin/env julia

using Gcd

for a = 1:5, b = 1:2:15
    println("gcd($a, $b) = $(my_gcd(a, b))")
end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: jl:light,ipynb
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.5.1
#   kernelspec:
#     display_name: Julia 1.4.2
#     language: julia
#     name: julia-1.4
# ---

# +
"""Tutorial: Overlap, Kinetic, and Dipole Integrals"""

__authors__   = ["D. Menendez", "Adam S. Abbott", "Boyi Zhang", "Justin M. Turney"]
__credit__    = ["D. Menendez", "Adam S. Abbott", "Boyi Zhang", "Justin M. Turney"]

__copyright__ = "(c) 2014-2020, The Psi4Julia Developers"
__license__   = "BSD-3-Clause"
__date__      = "2020-08-03"
# -

# # Overlap, Kinetic, and Dipole Integrals
#
# In this tutorial we will compute the overlap and kinetic energy integrals encountered in Hartree-Fock. We will also compute the dipole integrals to obtain the molecular dipole moment.
#
# ## Recurrence formula for one-electron integrals over Gaussian functions
# The direct calculation of every single integral in our arrays is not only a complicated procedure, but an expensive one too. Here, we will use the Obara-Saika recursion scheme for a simpler and more efficient implementation. For example, for the overlap integrals, given the value of $(s|s)$, one can recursively determine all other overlap integrals.
#
# ### Cartesian Gaussian functions
# Denote the origin of a 3-dimensional cartesian gaussian function by the coordinates $\mathbf{R} = (R_x, R_y, R_z)$.
# Let $\mathbf{r} = (x, y, z)$ be the coordinates of the electron, and $\alpha$ be the orbital exponent. We can now define an *unnormalized* Cartesian Gaussian function as:
#
# \begin{equation}
# \phi(\mathbf r; \alpha, \mathbf n, \mathbf R) = (x - R_x)^{n_x} (y - R_y)^{n_y} (z - R_z)^{n_z} \exp[-\alpha (\mathbf r - \mathbf R)^2]
# \end{equation}
#
# where $\alpha$ is the orbital exponent, and $\mathbf{n} = (n_x, n_y, n_z)$ is the angular momentum index vector. The sum $n_x + n_y + n_z = \lambda$ will hereafter be referred to as the **angular momentum**. We define a **shell** to be a set of functions (*components*) which share the same origin $\mathbf{R}$, and angular momentum $\lambda$.
#
# The shells with $\lambda$  equal to $0, 1, 2,...,$ are referred to as the $s, p, d, ...$ shell. Each shell has $(\lambda + 1) (\lambda + 2)/2$ components. The $s$ shell, with angular momentum $\lambda = 0$ has one component usually designated as $s$. The $p$ shell ($\lambda = 1$) has three components, designated as $p_x, p_y, p_z$. The $d$ shell ($\lambda = 2$) has six components, designated as $d_{xx}, d_{yy}, d_{zz}, d_{xy}, d_{xz}, d_{yz}$.
#
# In quantum chemistry, we typically represent a single component (an **atomic orbital**) by a linear combination of Gaussians ($c_1 \phi_1 + c_2 \phi_2 + c_3 \phi_3 ...$) hereon referred to as **primitives**. For example, the STO-3G basis set uses three primitives for each atomic orbital basis function. Each primitive is weighted by a coefficient $c$.
#
# Using our angular momentum index vector $\mathbf{n}$ we denote the single component of the $s$ shell to have angular momentum index $\mathbf{n} = \mathbf{0} = (0, 0, 0)$. Since the $p$ shell has three components, we may compactly express the angular momentum index vector as $\mathbf{1}_i$ where $i$ may be $x$, $y$, or $z$, and $\mathbf{1}_i = (\delta_{ix}, \delta_{iy}, \delta_{iz})$. For example, $p_x$ may be represented as $\mathbf{1}_x = (1, 0, 0)$. For the six components of the $d$ shell, we require a sum of two angular momentum index vectors $\mathbf{1}_i + \mathbf{1}_j$, where $(i,j = x,y,z)$. In this notation, the $d_{xy}$ component is $\mathbf{1}_x + \mathbf{1}_y = (1,0,0) + (0,1,0) = (1,1,0)$. To obtain higher order angular momentum components, we add the appropriate number of $\mathbf{1}_i$'s ($\mathbf{n}$'s).
#
#
# ### Two-center overlap integrals
# Two-center overlap integrals over unnormalized cartesian gaussian functions are of the form:
# \begin{equation}
# (\mathbf a|\mathbf b) = \int d\mathbf r\ \phi(\mathbf r; \alpha_a, \mathbf a, \mathbf A)\phi(\mathbf r; \alpha_b, \mathbf b, \mathbf B)
# \end{equation}
#
# Given $(\mathbf 0_A | \mathbf 0_B)$, we can use the Obara-Saika recursion relation to obtain overlap integrals between all basis functions. The overlap over $s$ functions is given by:
#
# \begin{equation}
# (\mathbf 0_A | \mathbf 0_B) = \left(\frac{\pi}{\alpha}\right)^{3/2} \exp[-\xi(\mathbf A-\mathbf B)^2]
# \end{equation}
#
# where $\alpha = \alpha_a + \alpha_b$ and $\zeta = \frac{\alpha_a\alpha_b}{\alpha}$.
# The recursion relations are given below. For a full derivation, see the appendix, or the
# [original paper](http://aip.scitation.org/doi/abs/10.1063/1.450106) by Obara and Saika. To increment the left side angular momentum:
#
# \begin{equation}
# (\mathbf a+\mathbf 1_i|\mathbf b) = (\mathbf{P - A})(\mathbf a|\mathbf b) + \frac{1}{2\alpha} N_i(\mathbf a)(\mathbf a-\mathbf 1_i|\mathbf b) + \frac{1}{2\alpha} N_i(\mathbf b)(\mathbf a|\mathbf b-\mathbf 1_i)
# \end{equation}
# and similarily, to increment the right side:
# \begin{equation}
# (\mathbf a|\mathbf b+\mathbf 1_i) = (\mathbf{P - B})(\mathbf a|\mathbf b) + \frac{1}{2\alpha} N_i(\mathbf a)(\mathbf a-\mathbf 1_i|\mathbf b) + \frac{1}{2\alpha} N_i(\mathbf b)(\mathbf a|\mathbf b-\mathbf 1_i)
# \end{equation}
#
#
# where \begin{equation}\mathbf{P} = \frac{\alpha_a \mathbf{A} + \alpha_b \mathbf{B}} {\alpha} \end{equation}
#
# and ${N}_i(\mathbf{a})$, ${N}_i(\mathbf{b})$ are just the angular momenta of $\mathbf{a}$ and $\mathbf{b}$.
#
# To fill in the first row, the seoncd term goes to zero, since we cannot have a negative angular momentum (cannot subtract $\mathbf 1_i$). To fill in the first column, the third term goes to zero for the same reason.
#
#

# # Implementation
# ### Write the recursion function
#
# The Obara-Saika recursion relationships depend on $\mathbf{P-A}$,  $\mathbf{P-B}$, $\alpha$, and the angular momentum values for $\mathbf a$ and $\mathbf b$, which we will denote as `PA`, `PB`, `alpha`, `AMa`, and `AMb`.
# Let's write a function that takes these parameters and returns three matrices containing the x, y and z components of our unnormalized overlap integrals. These same components also can be used to construct our kinetic energy and dipole integrals later, as we will see.
#
# For now, we will set $(\mathbf 0_A | \mathbf 0_B)$ to $1.0$ for simplicity. Later, we will set the value as defined above.

# +
using PyCall: pyimport
psi4 = pyimport("psi4")
np   = pyimport("numpy") # used only to cast to Psi4 arrays

function os_recursion(PA, PB, alpha, AMa, AMb)
    if length(PA) != 3 || length(PB) != 3
       error("PA and PB must be xyz coordinates.")
    end

    # Allocate space x, y, and z matrices
    # We add one because the equation for the kinetic energy
    # integrals require terms one beyond those in the overlap
    x = zeros(AMa + 1, AMb + 1)
    y = zeros(AMa + 1, AMb + 1)
    z = zeros(AMa + 1, AMb + 1)

    # Define 1/2alpha factor for convenience
    oo2a = 1 / (2alpha)

    # Set initial conditions (0a|0b) to 1.0 for each cartesian component
    x[1, 1] = y[1, 1] = z[1, 1] = 1.0

    # BEGIN RECURSION
    # Fill in the [0,1] position with PB
    if AMb > 0
        x[1, 2] = PB[1]
        y[1, 2] = PB[2]
        z[1, 2] = PB[3]
    end

    # Fill in the rest of row zero
    for b in 2:AMb
        x[1, b+1] = PB[1] * x[1, b] + b * oo2a * x[1, b - 1]
        y[1, b+1] = PB[2] * y[1, b] + b * oo2a * y[1, b - 1]
        z[1, b+1] = PB[3] * z[1, b] + b * oo2a * z[1, b - 1]
    end

    # Now, we have for each cartesian component
    # | 1.0  PB #  #|
    # |  0   0  0  0|
    # |  0   0  0  0|
    # |  0   0  0  0|

    # Upward recursion in a for all b's
    # Fill in the [1,0] position with PA
    if AMa > 0
        x[2, 1] = PA[1]
        y[2, 1] = PA[2]
        z[2, 1] = PA[3]

    # Now, we have for each cartesian component
    # | 1.0  PB #  #|
    # |  PA  0  0  0|
    # |  0   0  0  0|
    # |  0   0  0  0|

        # Fill in the rest of row one
        for b in 2:AMb
            x[2, b] = PA[1] * x[1, b] + b * oo2a * x[1, b - 1]
            y[2, b] = PA[2] * y[1, b] + b * oo2a * y[1, b - 1]
            z[2, b] = PA[3] * z[1, b] + b * oo2a * z[1, b - 1]
        end

        # Now, we have for each cartesian component
        # | 1.0  PB #  #|
        # |  PA  #  #  #|
        # |  0   0  0  0|
        # |  0   0  0  0|

        # Fill in the rest of column 0
        for a in 2:AMa
            x[a + 1, 1] = PA[1] * x[a, 1] + a * oo2a * x[a - 1, 1]
            y[a + 1, 1] = PA[2] * y[a, 1] + a * oo2a * y[a - 1, 1]
            z[a + 1, 1] = PA[3] * z[a, 1] + a * oo2a * z[a - 1, 1]

        # Now, we have for each cartesian component
        # | 1.0  PB #  #|
        # |  PA  #  #  #|
        # |  #   0  0  0|
        # |  #   0  0  0|

        # Fill in the rest of the a'th row
            for b in 2:AMb
                x[a + 1, b] = PA[1] * x[a, b] + a * oo2a * x[a - 1, b] + b * oo2a * x[a, b - 1]
                y[a + 1, b] = PA[2] * y[a, b] + a * oo2a * y[a - 1, b] + b * oo2a * y[a, b - 1]
                z[a + 1, b] = PA[3] * z[a, b] + a * oo2a * z[a - 1, b] + b * oo2a * z[a, b - 1]
            end
        end

        # Now, we have for each cartesian component
        # | 1.0  PB #  #|
        # |  PA  #  #  #|
        # |  #   #  #  #|
        # |  #   #  #  #|
    end

    # Return the results
    (x, y, z)
end
# -

# # Overlap Integrals
#
# Now that we have our recursion set up, we are ready to compute the integrals, starting with the overlap integrals. First we need a molecule and a basis set. We do a quick Hartree-Fock computation to obtain a wavefunction for easy access to the basis set. Later, we will also use the wavefunction's density matrix to compute the dipole moment.
#
# By default, Psi4 uses spherical harmonic Gaussians, sometimes referred to as "pure angular momentum", whereas we will be using Cartesian Gaussians. We set `puream` to be off so that we can compare our integrals to Psi4's.

# +
psi4.core.set_output_file("output.dat", false)

mol = psi4.geometry("""
                        O
                        H 1 1.1
                        H 1 1.1 2 104
                        symmetry c1
                        """)

psi4.set_options(Dict("basis"         => "sto-3g",
                      "scf_type"      => "pk",
                      "mp2_type"      => "conv",
                      "puream"        => 0,
                      "e_convergence" => 1e-8,
                      "d_convergence" => 1e-8))

scf_e, scf_wfn = psi4.energy("scf", return_wfn=true)

basis = scf_wfn.basisset()
# -

# The factored form of the two-center overlap integrals is
# \begin{equation}
# (\mathbf a | \mathbf b) = \space I_x(n_{ax},n_{bx},n_{cx})\ I_y(n_{ay},n_{by},n_{cy})\ I_z(n_{az},n_{bz},n_{cz})
# \end{equation}
#
# where $I_x$ is the x component of the unnormalized overlap matrix.
#
# To obtain the overlap matrix, we must loop over all shells, and for each shell, loop over the primitive gaussians.
#
# ~~~julia
# # make space to store the overlap integral matrix
# S = zeros(basis.nao(),basis.nao())
#
# # loop over the shells, basis.nshell is the number of shells
# for i in 1:basis.nshell(), j in 1:basis.nshell()
#     # basis.shell is a shell (1s, 2s, 2p, etc.)
#     # for water, there are 5 shells: (H: 1s, H: 1s, O: 1s, 2s, 2p)
#     ishell = basis.shell(i-1)
#     jshell = basis.shell(j-1)
#     # each shell has some number of primitives which make up each component of a shell
#     # sto-3g has 3 primitives for every component of every shell.
#     nprimi = ishell.nprimitive
#     nprimj = jshell.nprimitive
# ~~~
#
# At this point we access the information for each primitive basis set $\mathbf{a}$ and $\mathbf{b}$ in $(\mathbf{a}|\mathbf{b})$. The `basis.shell()` (which we have set to `ishell` and `jshell`) is a `psi4.core.GaussianShell` object, which contains all the information we need about the primitives within a shell, such as the exponent `(shell().exp())`, the coefficient `(shell().coef())`, the angular momentum `(shell().am)` and the atom number it is centered over `(shell().ncenter)`, which we can use to get the coordinates of the primitive center ($\mathbf{A}$ or $\mathbf{B}$) (*whew!*). We store all of these parameters as variables in our for-loops for easy access. Continuing our sample script:
#
# ~~~julia
#     # loop over the primitives within a shell
#     for a in 1:nprimi, b in 1:nprimj
#         expa = ishell.exp(a-1) # primitive exponents
#         expb = jshell.exp(b-1)
#         coefa = ishell.coef(a-1)  # primitive coefficients
#         coefb = jshell.coef(b-1)
#         AMa = ishell.am  # angular momenta
#         AMb = jshell.am
#         # defining centers for each basis function:
#         # mol.x() returns the x coordinate of the atom given by ishell.ncenter
#         # we use this to define a coordinate vector for our centers
#         A = [mol.x(ishell.ncenter), mol.y(ishell.ncenter), mol.z(ishell.ncenter)]
#         B = [mol.x(jshell.ncenter), mol.y(jshell.ncenter), mol.z(jshell.ncenter)]
# ~~~
#
# We are now ready to define `PA` and `PB` which are arguments for our recursion function.
#
#
# ~~~julia
#         alpha = expa + expb
#         zeta = (expa * expb) / alpha
#         P = (expa * A + expb * B) / alpha
#         PA = P - A
#         PB = P - B
#         AB = A - B
# ~~~
#
# We also define the value of the $(s|s)$ overlap `start`, originally set to 1.0 in the recursion function.
#
# \begin{equation}
# s = (\mathbf 0_A | \mathbf 0_B) = \left(\frac{\pi}{\alpha}\right)^{3/2} \exp[-\zeta(\mathbf A-\mathbf B)^2]
# \end{equation}
#
# ~~~julia
#        start = (π / alpha)^(3 / 2) * exp(-zeta * (AB[1]^2 + AB[2]^2 + AB[3]^2))
# ~~~
#
# Form the x-component, y-component, and z-component overlap matrices with our recursion function. Recall that we set $\mathbf{(0|0)} = 1.0$ and have ignored using the coefficients thus far. We will account for this later when we build each element of our overlap matrix.
#
# ~~~julia
#        # call the recursion
#        x, y, z = os_recursion(PA, PB, alpha, AMa+1, AMb+1)
# ~~~
#
#
# We are now ready to assign the elements of our overlap matrix. Currently, we are only looping over the shells, but not the individual components of each shell (e.g. we are not distinguishing between $p_x$, $p_y$, and $p_z$). Given the total angular momentum of **a** (`AMa`) and **b** (`AMb`), we loop over all possible allocations of this total angular momentum to each angular momentum index in $\mathbf{n} = (n_x, n_y, n_z)$ which we denote as `l`, `m`, and `n`.
#
# ~~~julia
#         # Basis function index where the shell begins
#         i_idx = ishell.function_index
#         j_idx = jshell.function_index
#
#         # We use counters to keep track of which component (e.g., p_x, p_y, p_z)
#         # within the shell we are on
#         counta = 1
#
#         for p in 1:AMa
#             la = AMa - p                    # Let l take on all values, and p be the leftover a.m.
#             for q in 1:p
#                 ma = p - q                  # distribute all leftover a.m. to m and n
#                 na = q
#                 countb = 1
#                 for r in 1:AMb
#                     lb = AMb - r            # Let l take on all values, and r the leftover a.m.
#                     for s in 1:r
#                         mb = r - s          # distribute all leftover a.m. to m and n
#                         nb = s
#
#                         # set the value in the full overlap matrix
#                         S[i_idx + counta, j_idx + countb] += start    *
#                                                              coefa    *
#                                                              coefb    *
#                                                              x[la+1,lb+1] *
#                                                              y[ma+1,mb+1] *
#                                                              z[na+1,nb+1]
#
#                         countb += 1
#                     end
#                 end
#                 counta += 1
#             end
#         end
#     end
# end
# ~~~

# ## Computing the Full Overlap Integral Matrix
#
# Putting all the above code snippets together we form our full overlap matrix:

# +
# make space to store the overlap integral matrix
S = zeros(basis.nao(),basis.nao())

# loop over the shells, basis.nshell is the number of shells
for i in 1:basis.nshell(), j in 1:basis.nshell()
    # basis.shell is a shell (1s, 2s, 2p, etc.)
    # for water, there are 5 shells: (H: 1s, H: 1s, O: 1s, 2s, 2p)
    ishell = basis.shell(i-1)
    jshell = basis.shell(j-1)
    # each shell has some number of primitives which make up each component of a shell
    # sto-3g has 3 primitives for every component of every shell.
    nprimi = ishell.nprimitive
    nprimj = jshell.nprimitive
    # loop over the primitives within a shell
    for a in 1:nprimi, b in 1:nprimj
        expa = ishell.exp(a-1) # exponents
        expb = jshell.exp(b-1)
        coefa = ishell.coef(a-1)  # coefficients
        coefb = jshell.coef(b-1)
        AMa = ishell.am  # angular momenta
        AMb = jshell.am
        # defining centers for each basis function
        # mol.x() returns the x coordinate of the atom given by ishell.ncenter
        # we use this to define a coordinate vector for our centers
        A = [mol.x(ishell.ncenter), mol.y(ishell.ncenter), mol.z(ishell.ncenter)]
        B = [mol.x(jshell.ncenter), mol.y(jshell.ncenter), mol.z(jshell.ncenter)]

        alpha = expa + expb
        zeta = (expa * expb) / alpha
        P = (expa * A + expb * B) / alpha
        PA = P - A
        PB = P - B
        AB = A - B
        start = (π / alpha)^(3 / 2) * exp(-zeta * (AB[1]^2 + AB[2]^2 + AB[3]^2))

        # call the recursion
        x, y, z = os_recursion(PA, PB, alpha, AMa+1, AMb+1)

        # Basis function index where the shell begins
        i_idx = ishell.function_index
        j_idx = jshell.function_index

        # We use counters to keep track of which component (e.g., p_x, p_y, p_z)
        # within the shell we are on
        counta = 1

        for p in 1:AMa
            la = AMa - p                    # Let l take on all values, and p be the leftover a.m.
            for q in 1:p
                ma = p - q                  # distribute all leftover a.m. to m and n
                na = q
                countb = 1
                for r in 1:AMb
                    lb = AMb - r            # Let l take on all values, and r the leftover a.m.
                    for s in 1:r
                        mb = r - s          # distribute all leftover a.m. to m and n
                        nb = s

                        # set the value in the full overlap matrix
                        S[i_idx + counta, j_idx + countb] += start    *
                                                             coefa    *
                                                             coefb    *
                                                             x[la+1,lb+1] *
                                                             y[ma+1,mb+1] *
                                                             z[na+1,nb+1]

                        countb += 1
                    end
                end
                counta += 1
            end
        end
    end
end
# -

# Check our overlap against Psi4:

mints = psi4.core.MintsHelper(basis)
Spsi4 = np.asarray(mints.ao_overlap())
@assert np.allclose(S, Spsi4, 6)

#
# # Kinetic energy integrals
# Before proceeding to an analysis of the kinetic energy integral, it will prove convenient to establish a short-hand notation for integrals related to the overlap integral:
# \begin{equation}
# (0 \mid 0) = \int G(\alpha_a, \mathbf A, l_a, m_a, n_a) G(\alpha_b, \mathbf B, l_b, m_b, n_b)\ d\tau
# \end{equation}
# The symbol $(+1 \mid 0)_{\mathrm{x}}$ will denote an integral of the form given by the above equation, except that the quantum number $l_a$ has been incremented by 1. Similarly, $(0 \mid +1)_{\mathrm{x}}$ increments $l_b$ by 1. Quantum numbers $m$ and $n$ are incremented in the same way with subscripts $y$ and $z$.
#
# The kinetic energy integral is defined as
#
# \begin{equation}
# T = - \frac{1}{2}  \int G(\alpha_a, \mathbf A, l_a, m_a, n_a) \left(\frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2} + \frac{\partial^2}{\partial z^2}\right) G(\alpha_b, \mathbf B, l_b, m_b, n_b)\ d\tau
# \end{equation}
#
# The equation for this integral is "unsymmetric" because the quantum numbers of the Gaussian-type function centered on $\mathbf B$ are altered by the Laplacian while those of the Gaussian-type function centered on $\mathbf A$ are not. The symmetric $x$, $y$, and $z$ components are given by:
#
#
# \begin{equation}
# T_x = \frac{1}{2}\Big\{l_a l_b (-1|-1)_{\mathrm{x}} + 4\alpha_a \alpha_b (+1|+1)_{\mathrm{x}} - 2\alpha_a l_b (+1|-1)_{\mathrm{x}} - 2\alpha_b l_a (-1|+1)_{\mathrm{x}}\Big\}
# \end{equation}
#
# \begin{equation}
# T_y = \frac{1}{2}\Big\{l_a l_b (-1|-1)_{\mathrm{y}} + 4\alpha_a \alpha_b (+1|+1)_{\mathrm{y}} - 2\alpha_a l_b (+1|-1)_{\mathrm{y}} - 2\alpha_b l_a (-1|+1)_{\mathrm{y}}\Big\}
# \end{equation}
#
# \begin{equation}
# T_z = \frac{1}{2}\Big\{l_a l_b (-1|-1)_{\mathrm{z}} + 4\alpha_a \alpha_b (+1|+1)_{\mathrm{z}} - 2\alpha_a l_b (+1|-1)_{\mathrm{z}} - 2\alpha_b l_a (-1|+1)_{\mathrm{z}}\Big\}
# \end{equation}
#
# and the full kinetic energy integral is:
#
# \begin{equation}
# T =  (T_x + T_y + T_z)
# \end{equation}
#
#
# Within our same for-loop structure as for the overlap integrals, we can compute the kinetic energy integrals with:
#
# ~~~julia
# Tx = (1 / 2) * (la * lb * x[lam, lbm] + 4 * expa * expb * x[la + 2, lb + 2] -
#        2 * expa * lb * x[la + 2, lbm] - 2 * expb * la * x[lam, lb + 2])   *
#        y[ma+1, mb+1] * z[na+1, nb+1]
#
# Ty = (1 / 2) * (ma * mb * y[mam, mbm] + 4 * expa * expb * y[ma + 2, mb + 2] -
#        2 * expa * mb * y[ma + 2, mbm] - 2 * expb * ma * y[mam, mb + 2])   *
#        x[la+1, lb+1] * z[na+1, nb+1]
#
# Tz = (1 / 2) * (na * nb * z[nam, nbm] + 4 * expa * expb * z[na + 2, nb + 2] -
#        2 * expa * nb * z[na + 2, nbm] - 2 * expb * na * z[nam, nb + 2])   *
#        x[la+1, lb+1] * y[ma+1, mb+1]
#
# # incorporate the value of (0|0) and coefficients
#
# T[i_idx + counta, j_idx + countb] += start * coefa * coefb * (Tx + Ty + Tz)
# ~~~
#
# # Dipole moment integrals
# We discuss these non-energy integrals here because they are frequently used and closely related to overlap and kinetic energy integrals. The dipole moment is defined with respect to a point in space $\mathbf C$, which is almost always taken to be the center of mass. Fortunately for our purposes, Psi4 automatically moves the molecule to be centered on the center of mass and thus $\mathbf C$ is $\mathbf 0$. The dipole moment integral is written for the $x$ direction as
#
# \begin{equation}
# d_{\mathrm{x}} =  \int G(\alpha_a, \mathbf A, l_a, m_a, n_a) G(\alpha_b, \mathbf B, l_b, m_b, n_b) \mathrm{x_\mathbf{C}}\ d\tau
# \end{equation}
#
# and similarly for the operators $\mathrm{y_\mathbf{C}}$ and $\mathrm{z_\mathbf{C}}$. Here, $x_\mathbf{C}$ is
#
# \begin{equation}
# \mathrm{x_\mathbf{C}} = \mathrm{x} - \mathrm{C_x}
# \end{equation}
#
#
# A convenient procedure is to redefine $\mathrm{x_\mathbf{C}}$ in terms of  $\mathrm{x_\mathbf{A}}$ or $\mathrm{x_\mathbf{B}}$ (we will use $\mathrm{x_\mathbf{A}}$):
#
# \begin{equation}
# \mathrm{x_\mathbf{C}}= (\mathrm{x - A_x}) + (\mathrm{A_x - C_x}) = \mathrm{x_\mathbf{A}} + (\mathrm{A_x - C_x})
# \end{equation}
#
# and we find
#
# \begin{equation}
# d_\mathrm{x} = (+1|0)_{\mathrm{x}} + (\mathrm{A_x - C_x}) (0|0)
# \end{equation}
#
# and the expressions are similar for $d_{\mathrm{y}}$ and $d_{\mathrm{z}}$. Putting this into code:
#
#
# ~~~julia
# dx = (x[la + 2, lb+1] + A[1] * x[la+1, lb+1]) * y[ma+1, mb+1] * z[na+1, nb+1]
# dy = (y[ma + 2, mb+1] + A[2] * y[ma+1, mb+1]) * x[la+1, lb+1] * z[na+1, nb+1]
# dz = (z[na + 2, nb+1] + A[3] * z[na+1, nb+1]) * x[la+1, lb+1] * y[ma+1, mb+1]
#
# # incorporate the (0|0) value and the coefficients:
#
# Dx[i_idx + counta, j_idx + countb] += start * coefa * coefb * dx
# Dy[i_idx + counta, j_idx + countb] += start * coefa * coefb * dy
# Dz[i_idx + counta, j_idx + countb] += start * coefa * coefb * dz
# ~~~
#
#

# # Full Implementation
#
# Putting everything together, we get the overlap, kinetic, and dipole integrals:

# +
# make space to store the overlap, kinetic, and dipole integral matrices
S = zeros(basis.nao(),basis.nao())
T = zeros(basis.nao(),basis.nao())
Dx = zeros(basis.nao(),basis.nao())
Dy = zeros(basis.nao(),basis.nao())
Dz = zeros(basis.nao(),basis.nao())

# loop over the shells, basis.nshell is the number of shells
for i in 1:basis.nshell(), j in 1:basis.nshell()
    # basis.shell is a shell (1s, 2s, 2p, etc.)
    # for water, there are 5 shells: (H: 1s, H: 1s, O: 1s, 2s, 2p)
    ishell = basis.shell(i-1)
    jshell = basis.shell(j-1)
    # each shell has some number of primitives which make up each component of a shell
    # sto-3g has 3 primitives for every component of every shell.
    nprimi = ishell.nprimitive
    nprimj = jshell.nprimitive
    # loop over the primitives within a shell
    for a in 1:nprimi, b in 1:nprimj
        expa = ishell.exp(a-1) # exponents
        expb = jshell.exp(b-1)
        coefa = ishell.coef(a-1)  # coefficients
        coefb = jshell.coef(b-1)
        AMa = ishell.am  # angular momenta
        AMb = jshell.am
        # defining centers for each basis function
        # mol.x() returns the x coordinate of the atom given by ishell.ncenter
        # we use this to define a coordinate vector for our centers
        A = [mol.x(ishell.ncenter), mol.y(ishell.ncenter), mol.z(ishell.ncenter)]
        B = [mol.x(jshell.ncenter), mol.y(jshell.ncenter), mol.z(jshell.ncenter)]
        alpha = expa + expb
        zeta = (expa * expb) / alpha
        P = (expa * A + expb * B) / alpha
        PA = P - A
        PB = P - B
        AB = A - B
        start = (π / alpha)^(3 / 2) * exp(-zeta * (AB[1]^2 + AB[2]^2 + AB[3]^2))
        # call the recursion
        x, y, z = os_recursion(PA, PB, alpha, AMa+1, AMb+1)

        # Basis function index where the shell begins
        i_idx = ishell.function_index
        j_idx = jshell.function_index

        # We use counters to keep track of which component (e.g., p_x, p_y, p_z)
        # within the shell we are on
        counta = 1

        for p in 1:AMa
            la = AMa - p                    # Let l take on all values, and p be the leftover a.m.
            for q in 1:p
                ma = p - q                  # distribute all leftover a.m. to m and n
                na = q
                countb = 1
                for r in 1:AMb
                    lb = AMb - r            # Let l take on all values, and r the leftover a.m.
                    for s in 1:r
                        mb = r - s          # distribute all leftover a.m. to m and n
                        nb = s

                        S[i_idx + counta, j_idx + countb] += start    *
                                                             coefa    *
                                                             coefb    *
                                                             x[la+1,lb+1] *
                                                             y[ma+1,mb+1] *
                                                             z[na+1,nb+1]

                        lam = la < 1 ? lastindex(x,1) : la
                        lbm = lb < 1 ? lastindex(x,1) : lb
                        nam = na < 1 ? lastindex(x,1) : na
                        nbm = nb < 1 ? lastindex(x,1) : nb
                        mam = ma < 1 ? lastindex(x,1) : ma
                        mbm = mb < 1 ? lastindex(x,1) : mb
                        Tx = (1 / 2) * (la * lb * x[lam, lbm] + 4 * expa * expb * x[la + 2, lb + 2] -
                               2 * expa * lb * x[la + 2, lbm] - 2 * expb * la * x[lam, lb + 2])   *
                               y[ma+1, mb+1] * z[na+1, nb+1]

                        Ty = (1 / 2) * (ma * mb * y[mam, mbm] + 4 * expa * expb * y[ma + 2, mb + 2] -
                               2 * expa * mb * y[ma + 2, mbm] - 2 * expb * ma * y[mam, mb + 2])   *
                               x[la+1, lb+1] * z[na+1, nb+1]

                        Tz = (1 / 2) * (na * nb * z[nam, nbm] + 4 * expa * expb * z[na + 2, nb + 2] -
                               2 * expa * nb * z[na + 2, nbm] - 2 * expb * na * z[nam, nb + 2])   *
                               x[la+1, lb+1] * y[ma+1, mb+1]

                        T[i_idx + counta, j_idx + countb] += start * coefa * coefb * (Tx + Ty + Tz)

                        dx = (x[la + 2, lb+1] + A[1] * x[la+1, lb+1]) * y[ma+1, mb+1] * z[na+1, nb+1]
                        dy = (y[ma + 2, mb+1] + A[2] * y[ma+1, mb+1]) * x[la+1, lb+1] * z[na+1, nb+1]
                        dz = (z[na + 2, nb+1] + A[3] * z[na+1, nb+1]) * x[la+1, lb+1] * y[ma+1, mb+1]

                        Dx[i_idx + counta, j_idx + countb] += start * coefa * coefb * dx
                        Dy[i_idx + counta, j_idx + countb] += start * coefa * coefb * dy
                        Dz[i_idx + counta, j_idx + countb] += start * coefa * coefb * dz

                        countb += 1
                    end
                end
                counta += 1
            end
        end
    end
end
# -

# Check the kinetic energy integrals:

Tpsi4 = np.asarray(mints.ao_kinetic())
@assert np.allclose(T, Tpsi4, 6)

# # Computing the Dipole Moment

# The x-component of the dipole moment is given by
#
# \begin{align}
# \mu_x =& - \mu_{\mathrm{elec}} + \mu_{\mathrm{nuclear}}\\
# \mu_x =& - \sum_{\mu \nu}^{\mathrm{AO}} D_{\mu \nu} d_{\mu \nu}^{x} + \sum_A^N Z_A X_A
# \end{align}
#
# You can obtain the nuclear contribution to the dipole $\mu_{\mathrm{nuclear}}$ from Psi4 using
#
#     psi4.core.molecule.nuclear_dipole()
#
# The total dipole moment is given by
#
# \begin{equation}
# \mu = \sqrt{\mu_x^2 + \mu_y^2 + \mu_z^2}
# \end{equation}

# +
#Code to compute dipole moment

# Get density matrices
Da = np.asarray(scf_wfn.Da())
Db = np.asarray(scf_wfn.Db())
D = Da + Db

# Get nuclear dipole
nuc_dipole = mol.nuclear_dipole()

# Compute dipole moment components
mux = -sum(D .* Dx) + get(nuc_dipole,1)
muy = -sum(D .* Dy) + get(nuc_dipole,2)
muz = -sum(D .* Dz) + get(nuc_dipole,3)

# Compute the dipole moment
mu = sqrt(mux^2 + muy^2 + muz^2)

# Get Psi4's dipole moment components
mux_psi4 = psi4.core.variable("SCF DIPOLE X")
muy_psi4 = psi4.core.variable("SCF DIPOLE Y")
muz_psi4 = psi4.core.variable("SCF DIPOLE Z")

# Compute Psi4 dipole moment
mu_psi4 = sqrt(mux_psi4^2 + muy_psi4^2 + muz_psi4^2)

# Psi4 prints in Debye, ours is in a.u.
# 2.54174623 Debye in 1 a.u.
mu_psi4 *= 1/2.54174623

# Compare Psi4 dipole moment to ours

@assert np.allclose(mu, mu_psi4, 4)

println("Psi4 calculated molecular dipole magnitude: ", mu_psi4)
println("Here calculated molecular dipole magnitude: ", mu)
# -

# # References
#
# 1. Recursion scheme:
#     > [[Obara:1986](http://aip.scitation.org/doi/abs/10.1063/1.450106)] S. Obara and A. Saika, *J. Chem. Phys* **84**, 3963 (1986)
#
# 2. T. Helgaker, P. Jorgensen, and J. Olsen, *Molecular Electronic Structure Theory*, John Wiley & Sons Inc, 2000.


text:
import Base: position
using FFTW

"""
    PositionBasis(xmin, xmax, Npoints)
    PositionBasis(b::MomentumBasis)

Basis for a particle in real space.

For simplicity periodic boundaries are assumed which means that
the rightmost point defined by `xmax` is not included in the basis
but is defined to be the same as `xmin`.

When a [`MomentumBasis`](@ref) is given as argument the exact values
of ``x_{min}`` and ``x_{max}`` are due to the periodic boundary conditions
more or less arbitrary and are chosen to be
``-\\pi/dp`` and ``\\pi/dp`` with ``dp=(p_{max}-p_{min})/N``.
"""
struct PositionBasis{T,X1,X2} <: Basis
    shape::Vector{T}
    xmin::Float64
    xmax::Float64
    N::T
    function PositionBasis{X1,X2}(xmin::Real, xmax::Real, N::T) where {X1,X2,T<:Int}
        @assert isa(X1, Real) && isa(X2, Real)
        new{T,X1,X2}([N], xmin, xmax, N)
    end
end
PositionBasis(xmin::Real, xmax::Real, N::Int) = PositionBasis{xmin,xmax}(xmin,xmax,N)

"""
    MomentumBasis(pmin, pmax, Npoints)
    MomentumBasis(b::PositionBasis)

Basis for a particle in momentum space.

For simplicity periodic boundaries are assumed which means that
`pmax` is not included in the basis but is defined to be the same as `pmin`.

When a [`PositionBasis`](@ref) is given as argument the exact values
of ``p_{min}`` and ``p_{max}`` are due to the periodic boundary conditions
more or less arbitrary and are chosen to be
``-\\pi/dx`` and ``\\pi/dx`` with ``dx=(x_{max}-x_{min})/N``.
"""
struct MomentumBasis{P1,P2} <: Basis
    shape::Vector{Int}
    pmin::Float64
    pmax::Float64
    N::Int
    function MomentumBasis{P1,P2}(pmin::Real, pmax::Real, N::Int) where {P1,P2}
        @assert isa(P1, Real) && isa(P2, Real)
        new([N], pmin, pmax, N)
    end
end
MomentumBasis(pmin::Real, pmax::Real, N::Int) = MomentumBasis{pmin,pmax}(pmin, pmax, N)

PositionBasis(b::MomentumBasis) = (dp = (b.pmax - b.pmin)/b.N; PositionBasis(-pi/dp, pi/dp, b.N))
MomentumBasis(b::PositionBasis) = (dx = (b.xmax - b.xmin)/b.N; MomentumBasis(-pi/dx, pi/dx, b.N))

==(b1::PositionBasis, b2::PositionBasis) = b1.xmin==b2.xmin && b1.xmax==b2.xmax && b1.N==b2.N
==(b1::MomentumBasis, b2::MomentumBasis) = b1.pmin==b2.pmin && b1.pmax==b2.pmax && b1.N==b2.N

"""
    gaussianstate(b::PositionBasis, x0, p0, sigma)
    gaussianstate(b::MomentumBasis, x0, p0, sigma)

Create a Gaussian state around `x0` and` p0` with width `sigma`.

In real space the gaussian state is defined as

```math
\\Psi(x) = \\frac{1}{\\pi^{1/4}\\sqrt{\\sigma}}
            e^{i p_0 (x-\\frac{x_0}{2}) - \\frac{(x-x_0)^2}{2 \\sigma^2}}
```

and is connected to the momentum space definition

```math
\\Psi(p) = \\frac{\\sqrt{\\sigma}}{\\pi^{1/4}}
            e^{-i x_0 (p-\\frac{p_0}{2}) - \\frac{1}{2}(p-p_0)^2 \\sigma^2}
```

via a Fourier-transformation

```math
\\Psi(p) = \\frac{1}{\\sqrt{2\\pi}}
            \\int_{-\\infty}^{\\infty} e^{-ipx}\\Psi(x) \\mathrm{d}x
```

The state has the properties

* ``⟨p⟩ = p_0``
* ``⟨x⟩ = x_0``
* ``\\mathrm{Var}(x) = \\frac{σ^2}{2}``
* ``\\mathrm{Var}(p) = \\frac{1}{2 σ^2}``

Due to the numerically necessary discretization additional scaling
factors ``\\sqrt{Δx}`` and ``\\sqrt{Δp}`` are used so that
``\\langle x_i|Ψ\\rangle = \\sqrt{Δ x} Ψ(x_i)`` and ``\\langle p_i|Ψ\\rangle = \\sqrt{Δ p} Ψ(p_i)`` so
that the resulting Ket state is normalized.
"""
function gaussianstate(b::PositionBasis, x0::Real, p0::Real, sigma::Real)
    psi = Ket(b)
    dx = spacing(b)
    alpha = 1.0/(pi^(1/4)*sqrt(sigma))*sqrt(dx)
    x = b.xmin
    for i=1:b.N
        psi.data[i] = alpha*exp(1im*p0*(x-x0/2) - (x-x0)^2/(2*sigma^2))
        x += dx
    end
    return psi
end

function gaussianstate(b::MomentumBasis, x0::Real, p0::Real, sigma::Real)
    psi = Ket(b)
    dp = spacing(b)
    alpha = sqrt(sigma)/pi^(1/4)*sqrt(dp)
    p = b.pmin
    for i=1:b.N
        psi.data[i] = alpha*exp(-1im*x0*(p-p0/2) - (p-p0)^2/2*sigma^2)
        p += dp
    end
    return psi
end

"""
    spacing(b::PositionBasis)

Difference between two adjacent points of the real space basis.
"""
spacing(b::PositionBasis) = (b.xmax - b.xmin)/b.N
"""
    spacing(b::MomentumBasis)

Momentum difference between two adjacent points of the momentum basis.
"""
spacing(b::MomentumBasis) = (b.pmax - b.pmin)/b.N

"""
    samplepoints(b::PositionBasis)

x values of the real space basis.
"""
samplepoints(b::PositionBasis) = (dx = spacing(b); Float64[b.xmin + i*dx for i=0:b.N-1])
"""
    samplepoints(b::MomentumBasis)

p values of the momentum basis.
"""
samplepoints(b::MomentumBasis) = (dp = spacing(b); Float64[b.pmin + i*dp for i=0:b.N-1])

"""
    position(b::PositionBasis)

Position operator in real space.
"""
position(b::PositionBasis) = SparseOperator(b, sparse(Diagonal(complex(samplepoints(b)))))

"""
    position(b:MomentumBasis)

Position operator in momentum space.
"""
function position(b::MomentumBasis)
    b_pos = PositionBasis(b)
    transform(b, b_pos)*dense(position(b_pos))*transform(b_pos, b)
end

"""
    momentum(b:MomentumBasis)

Momentum operator in momentum space.
"""
momentum(b::MomentumBasis) = SparseOperator(b, sparse(Diagonal(complex(samplepoints(b)))))

"""
    momentum(b::PositionBasis)

Momentum operator in real space.
"""
function momentum(b::PositionBasis)
    b_mom = MomentumBasis(b)
    transform(b, b_mom)*dense(momentum(b_mom))*transform(b_mom, b)
end

"""
    potentialoperator(b::PositionBasis, V(x))

Operator representing a potential ``V(x)`` in real space.
"""
function potentialoperator(b::PositionBasis, V::Function)
    x = samplepoints(b)
    diagonaloperator(b, V.(x))
end

"""
    potentialoperator(b::MomentumBasis, V(x))

Operator representing a potential ``V(x)`` in momentum space.
"""
function potentialoperator(b::MomentumBasis, V::Function)
    b_pos = PositionBasis(b)
    transform(b, b_pos)*dense(potentialoperator(b_pos, V))*transform(b_pos, b)
end

"""
    potentialoperator(b::CompositeBasis, V(x, y, z, ...))

Operator representing a potential ``V`` in more than one dimension.

# Arguments
* `b`: Composite basis consisting purely either of `PositionBasis` or
    `MomentumBasis`. Note, that calling this with a composite basis in
    momentum space might consume a large amount of memory.
* `V`: Function describing the potential. ATTENTION: The number of arguments
    accepted by `V` must match the spatial dimension. Furthermore, the order
    of the arguments has to match that of the order of the tensor product of
    bases (e.g. if `b=bx⊗by⊗bz`, then `V(x,y,z)`).
"""
function potentialoperator(b::CompositeBasis, V::Function)
    if isa(b.bases[1], PositionBasis)
        potentialoperator_position(b, V)
    elseif isa(b.bases[1], MomentumBasis)
        potentialoperator_momentum(b, V)
    else
        throw(IncompatibleBases())
    end
end
function potentialoperator_position(b::CompositeBasis, V::Function)
    for base=b.bases
        @assert isa(base, PositionBasis)
    end

    points = [samplepoints(b1) for b1=b.bases]
    dims = length.(points)
    n = length(b.bases)
    data = Array{ComplexF64}(undef, dims...)
    @inbounds for i=1:length(data)
        index = Tuple(CartesianIndices(data)[i])
        args = (points[j][index[j]] for j=1:n)
        data[i] = V(args...)
    end

    diagonaloperator(b, data[:])
end
function potentialoperator_momentum(b::CompositeBasis, V::Function)
    bases_pos = []
    for base=b.bases
        @assert isa(base, MomentumBasis)
        push!(bases_pos, PositionBasis(base))
    end
    b_pos = tensor(bases_pos...)
    transform(b, b_pos)*dense(potentialoperator_position(b_pos, V))*transform(b_pos, b)
end

"""
    FFTOperator

Abstract type for all implementations of FFT operators.
"""
abstract type FFTOperator{BL<:Basis, BR<:Basis, T} <: AbstractOperator{BL,BR} end

Base.eltype(x::FFTOperator) = promote_type(eltype(x.mul_before), eltype(x.mul_after))

"""
    FFTOperators

Operator performing a fast fourier transformation when multiplied with a state
that is a Ket or an Operator.
"""
mutable struct FFTOperators{BL<:Basis,BR<:Basis,T<:Array,P1,P2,P3,P4} <: FFTOperator{BL, BR, T}
    basis_l::BL
    basis_r::BR
    fft_l!::P1
    fft_r!::P2
    fft_l2!::P3
    fft_r2!::P4
    mul_before::T
    mul_after::T
    function FFTOperators(b1::BL, b2::BR,
        fft_l!::P1,
        fft_r!::P2,
        fft_l2!::P3,
        fft_r2!::P4,
        mul_before::T,
        mul_after::T) where {BL<:Basis,BR<:Basis,T,P1,P2,P3,P4}
        new{BL,BR,T,P1,P2,P3,P4}(b1, b2, fft_l!, fft_r!, fft_l2!, fft_r2!, mul_before, mul_after)
    end
end

"""
    FFTKets

Operator that can only perform fast fourier transformations on Kets.
This is much more memory efficient when only working with Kets.
"""
mutable struct FFTKets{BL<:Basis,BR<:Basis,T<:Array,P1,P2} <: FFTOperator{BL, BR, T}
    basis_l::BL
    basis_r::BR
    fft_l!::P1
    fft_r!::P2
    mul_before::T
    mul_after::T
    function FFTKets(b1::BL, b2::BR,
        fft_l!::P1,
        fft_r!::P2,
        mul_before::T,
        mul_after::T) where {BL<:Basis,BR<:Basis, T, P1, P2}
        new{BL, BR, T, P1, P2}(b1, b2, fft_l!, fft_r!, mul_before, mul_after)
    end
end

"""
    transform(b1::MomentumBasis, b2::PositionBasis)
    transform(b1::PositionBasis, b2::MomentumBasis)

Transformation operator between position basis and momentum basis.
"""
function transform(basis_l::MomentumBasis, basis_r::PositionBasis; ket_only::Bool=false)
    Lx = (basis_r.xmax - basis_r.xmin)
    dp = spacing(basis_l)
    dx = spacing(basis_r)
    if basis_l.N != basis_r.N || abs(2*pi/dp - Lx)/Lx > 1e-12
        throw(IncompatibleBases())
    end
    mul_before = exp.(-1im*basis_l.pmin*(samplepoints(basis_r) .- basis_r.xmin))
    mul_after = exp.(-1im*basis_r.xmin*samplepoints(basis_l))/sqrt(basis_r.N)
    x = Vector{ComplexF64}(undef, length(basis_r))
    if ket_only
        FFTKets(basis_l, basis_r, plan_bfft!(x), plan_fft!(x), mul_before, mul_after)
    else
        A = Matrix{ComplexF64}(undef, length(basis_r), length(basis_r))
        FFTOperators(basis_l, basis_r, plan_bfft!(x), plan_fft!(x), plan_bfft!(A, 2), plan_fft!(A, 1), mul_before, mul_after)
    end
end

"""
    transform(b1::CompositeBasis, b2::CompositeBasis)

Transformation operator between two composite bases. Each of the bases
has to contain bases of type PositionBasis and the other one a corresponding
MomentumBasis.
"""
function transform(basis_l::PositionBasis, basis_r::MomentumBasis; ket_only::Bool=false)
    Lx = (basis_l.xmax - basis_l.xmin)
    dp = spacing(basis_r)
    dx = spacing(basis_l)
    if basis_l.N != basis_r.N || abs(2*pi/dp - Lx)/Lx > 1e-12
        throw(IncompatibleBases())
    end
    mul_before = exp.(1im*basis_l.xmin*(samplepoints(basis_r) .- basis_r.pmin))
    mul_after = exp.(1im*basis_r.pmin*samplepoints(basis_l))/sqrt(basis_r.N)
    x = Vector{ComplexF64}(undef, length(basis_r))
    if ket_only
        FFTKets(basis_l, basis_r, plan_fft!(x), plan_bfft!(x), mul_before, mul_after)
    else
        A = Matrix{ComplexF64}(undef, length(basis_r), length(basis_r))
        FFTOperators(basis_l, basis_r, plan_fft!(x), plan_bfft!(x), plan_fft!(A, 2), plan_bfft!(A, 1), mul_before, mul_after)
    end
end

function transform(basis_l::CompositeBasis, basis_r::CompositeBasis; ket_only::Bool=false, index::Vector{Int}=Int[])
    @assert length(basis_l.bases) == length(basis_r.bases)
    if length(index) == 0
        check_pos = [isa.(basis_l.bases, PositionBasis)...]
        check_mom = [isa.(basis_l.bases, MomentumBasis)...]
        if any(check_pos) && !any(check_mom)
            index = [1:length(basis_l.bases);][check_pos]
        elseif any(check_mom) && !any(check_pos)
            index = [1:length(basis_l.bases);][check_mom]
        else
            throw(IncompatibleBases())
        end
    end
    if all(isa.(basis_l.bases[index], PositionBasis))
        @assert all(isa.(basis_r.bases[index], MomentumBasis))
        transform_xp(basis_l, basis_r, index; ket_only=ket_only)
    elseif all(isa.(basis_l.bases[index], MomentumBasis))
        @assert all(isa.(basis_r.bases[index], PositionBasis))
        transform_px(basis_l, basis_r, index; ket_only=ket_only)
    else
        throw(IncompatibleBases())
    end
end

function transform_xp(basis_l::CompositeBasis, basis_r::CompositeBasis, index::Vector{Int}; ket_only::Bool=false)
    n = length(basis_l.bases)
    Lx = [(b.xmax - b.xmin) for b=basis_l.bases[index]]
    dp = [spacing(b) for b=basis_r.bases[index]]
    dx = [spacing(b) for b=basis_l.bases[index]]
    N = [length(b) for b=basis_l.bases]
    for i=1:n
        if N[i] != length(basis_r.bases[i])
            throw(IncompatibleBases())
        end
    end
    for i=1:length(index)
        if abs(2*pi/dp[i] - Lx[i])/Lx[i] > 1e-12
            throw(IncompatibleBases())
        end
    end

    if index[1] == 1
        mul_before = exp.(1im*basis_l.bases[1].xmin*(samplepoints(basis_r.bases[1]) .- basis_r.bases[1].pmin))
        mul_after = exp.(1im*basis_r.bases[1].pmin*samplepoints(basis_l.bases[1]))/sqrt(basis_r.bases[1].N)
    else
        mul_before = ones(N[1])
        mul_after = ones(N[1])
    end
    for i=2:n
        if any(i .== index)
            mul_before = kron(exp.(1im*basis_l.bases[i].xmin*(samplepoints(basis_r.bases[i]) .- basis_r.bases[i].pmin)), mul_before)
            mul_after = kron(exp.(1im*basis_r.bases[i].pmin*samplepoints(basis_l.bases[i]))/sqrt(basis_r.bases[i].N), mul_after)
        else
            mul_before = kron(ones(N[i]), mul_before)
            mul_after = kron(ones(N[i]), mul_after)
        end
    end
    mul_before = reshape(mul_before, (N...,))
    mul_after = reshape(mul_after, (N...,))

    x = Array{ComplexF64}(undef, N...)
    if ket_only
        FFTKets(basis_l, basis_r, plan_fft!(x, index), plan_bfft!(x, index), mul_before, mul_after)
    else
        A = Array{ComplexF64}(undef, [N; N]...)
        FFTOperators(basis_l, basis_r, plan_fft!(x, index), plan_bfft!(x, index), plan_fft!(A, [n + 1:2n;][index]), plan_bfft!(A, [1:n;][index]), mul_before, mul_after)
    end
end

function transform_px(basis_l::CompositeBasis, basis_r::CompositeBasis, index::Vector{Int}; ket_only::Bool=false)
    n = length(basis_l.bases)
    Lx = [(b.xmax - b.xmin) for b=basis_r.bases[index]]
    dp = [spacing(b) for b=basis_l.bases[index]]
    dx = [spacing(b) for b=basis_r.bases[index]]
    N = [length(b) for b=basis_l.bases]
    for i=1:n
        if N[i] != length(basis_r.bases[i])
            throw(IncompatibleBases())
        end
    end
    for i=1:length(index)
        if abs(2*pi/dp[i] - Lx[i])/Lx[i] > 1e-12
            throw(IncompatibleBases())
        end
    end

    if index[1] == 1
        mul_before = exp.(-1im*basis_l.bases[1].pmin*(samplepoints(basis_r.bases[1]) .- basis_r.bases[1].xmin))
        mul_after = exp.(-1im*basis_r.bases[1].xmin*samplepoints(basis_l.bases[1]))/sqrt(N[1])
    else
        mul_before = ones(N[1])
        mul_after = ones(N[1])
    end
    for i=2:n
        if i in index
            mul_before = kron(exp.(-1im*basis_l.bases[i].pmin*(samplepoints(basis_r.bases[i]) .- basis_r.bases[i].xmin)), mul_before)
            mul_after = kron(exp.(-1im*basis_r.bases[i].xmin*samplepoints(basis_l.bases[i]))/sqrt(N[i]), mul_after)
        else
            mul_before = kron(ones(N[i]), mul_before)
            mul_after = kron(ones(N[i]), mul_after)
        end
    end
    mul_before = reshape(mul_before, (N...,))
    mul_after = reshape(mul_after, (N...,))

    x = Array{ComplexF64}(undef, N...)
    if ket_only
        FFTKets(basis_l, basis_r, plan_bfft!(x, index), plan_fft!(x, index), mul_before, mul_after)
    else
        A = Array{ComplexF64}(undef, [N; N]...)
        FFTOperators(basis_l, basis_r, plan_bfft!(x, index), plan_fft!(x, index), plan_bfft!(A, [n + 1:2n;][index]), plan_fft!(A, [1:n;][index]), mul_before, mul_after)
    end
end

DenseOperator(op::FFTOperator) = op*identityoperator(DenseOpType, op.basis_r)

dagger(op::FFTOperators) = transform(op.basis_r, op.basis_l)
dagger(op::FFTKets) = transform(op.basis_r, op.basis_l; ket_only=true)

tensor(A::FFTOperators, B::FFTOperators) = transform(tensor(A.basis_l, B.basis_l), tensor(A.basis_r, B.basis_r))
tensor(A::FFTKets, B::FFTKets) = transform(tensor(A.basis_l, B.basis_l), tensor(A.basis_r, B.basis_r); ket_only=true)

function mul!(result::Ket{B1},M::FFTOperator{B1,B2},b::Ket{B2},alpha_,beta_) where {B1<:Basis,B2<:Basis}
    alpha = convert(ComplexF64, alpha_)
    beta = convert(ComplexF64, beta_)
    N::Int = length(M.basis_r)
    if beta==Complex(0.)
        @inbounds for i=1:N
            result.data[i] = M.mul_before[i] * b.data[i]
        end
        M.fft_r! * reshape(result.data, size(M.mul_before))
        @inbounds for i=1:N
            result.data[i] *= M.mul_after[i] * alpha
        end
    else
        psi_ = Ket(M.basis_l, copy(b.data))
        @inbounds for i=1:N
            psi_.data[i] *= M.mul_before[i]
        end
        M.fft_r! * reshape(psi_.data, size(M.mul_before))
        @inbounds for i=1:N
            result.data[i] = beta*result.data[i] + alpha * psi_.data[i] * M.mul_after[i]
        end
    end
    result
end

function mul!(result::Bra{B2},b::Bra{B1},M::FFTOperator{B1,B2},alpha_,beta_) where {B1<:Basis,B2<:Basis}
    alpha = convert(ComplexF64, alpha_)
    beta = convert(ComplexF64, beta_)
    N::Int = length(M.basis_l)
    if beta==Complex(0.)
        @inbounds for i=1:N
            result.data[i] = conj(M.mul_after[i]) * conj(b.data[i])
        end
        M.fft_l! * reshape(result.data, size(M.mul_after))
        @inbounds for i=1:N
            result.data[i] = conj(result.data[i]) * M.mul_before[i] * alpha
        end
    else
        psi_ = Bra(M.basis_r, conj(b.data))
        @inbounds for i=1:N
            psi_.data[i] *= conj(M.mul_after[i])
        end
        M.fft_l! * reshape(psi_.data, size(M.mul_after))
        @inbounds for i=1:N
            result.data[i] = beta*result.data[i] + alpha * conj(psi_.data[i]) * M.mul_before[i]
        end
    end
    result
end

function mul!(result::Operator{B1,B3,T},A::Operator{B1,B2},B::FFTOperators{B2,B3},alpha_,beta_) where {B1<:Basis,B2<:Basis,B3<:Basis,T}
    alpha = convert(ComplexF64, alpha_)
    beta = convert(ComplexF64, beta_)
    if beta != Complex(0.)
        data = similar(result.data, size(result.data, 1), size(result.data, 2))
    else
        data = result.data
    end
    copyto!(data, A.data)
    @inbounds for j=1:length(B.mul_after), i=1:length(B.mul_after)
        data[i, j] *= B.mul_after[j]
    end
    conj!(data)
    n = size(B.mul_after)
    B.fft_l2! * reshape(data, n..., n...)
    conj!(data)
    N = prod(n)
    @inbounds for j=1:N, i=1:N
        data[i, j] *= B.mul_before[j]
    end
    if alpha != Complex(1.)
        lmul!(alpha, data)
    end
    if beta != Complex(0.)
        rmul!(result.data, beta)
        result.data += data
    end
    result
end

function mul!(result::Operator{B1,B3,T},A::FFTOperators{B1,B2},B::Operator{B2,B3},alpha_,beta_) where {B1<:Basis,B2<:Basis,B3<:Basis,T}
    alpha = convert(ComplexF64, alpha_)
    beta = convert(ComplexF64, beta_)
    if beta != Complex(0.)
        data = similar(result.data, size(result.data, 1), size(result.data, 2))
    else
        data = result.data
    end
    copyto!(data, B.data)
    @inbounds for j=1:length(A.mul_before), i=1:length(A.mul_before)
        data[i, j] *= A.mul_before[i]
    end
    n = size(A.mul_before)
    A.fft_r2! * reshape(data, n...,n...)
    N = prod(n)
    @inbounds for j=1:N, i=1:N
        data[i, j] *= A.mul_after[i]
    end
    if alpha != Complex(1.)
        lmul!(alpha, data)
    end
    if beta != Complex(0.)
        rmul!(result.data, beta)
        result.data += data
    end
    result
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

@testset "padding (lpad and rpad)" begin
    @test lpad("foo", 2) == "foo"
    @test rpad("foo", 2) == "foo"
    @test lpad("foo", 3) == "foo"
    @test rpad("foo", 3) == "foo"
    @test lpad("foo", 4) == " foo"
    @test rpad("foo", 4) == "foo "
    @test lpad("foo", 5) == "  foo"
    @test rpad("foo", 5) == "foo  "
    @test lpad("foo", 2, "123") == "foo"
    @test rpad("foo", 2, "123") == "foo"
    @test lpad("foo", 3, "123") == "foo"
    @test rpad("foo", 3, "123") == "foo"
    @test lpad("foo", 4, "123") == "1foo"
    @test rpad("foo", 4, "123") == "foo1"
    @test lpad("foo", 5, "123") == "12foo"
    @test rpad("foo", 5, "123") == "foo12"
    @test lpad("foo", 6, "123") == "123foo"
    @test rpad("foo", 6, "123") == "foo123"
    @test lpad("foo", 7, "123") == "1231foo"
    @test rpad("foo", 7, "123") == "foo1231"
    @test lpad("foo", 8, "123") == "12312foo"
    @test rpad("foo", 8, "123") == "foo12312"
    @test lpad("foo", 9, "123") == "123123foo"
    @test rpad("foo", 9, "123") == "foo123123"
    @test lpad("αβ", 2, "¹₂³") == "αβ"
    @test rpad("αβ", 2, "¹₂³") == "αβ"
    @test lpad("αβ", 3, "¹₂³") == "¹αβ"
    @test rpad("αβ", 3, "¹₂³") == "αβ¹"
    @test lpad("αβ", 4, "¹₂³") == "¹₂αβ"
    @test rpad("αβ", 4, "¹₂³") == "αβ¹₂"
    @test lpad("αβ", 5, "¹₂³") == "¹₂³αβ"
    @test rpad("αβ", 5, "¹₂³") == "αβ¹₂³"
    @test lpad("αβ", 6, "¹₂³") == "¹₂³¹αβ"
    @test rpad("αβ", 6, "¹₂³") == "αβ¹₂³¹"
    @test lpad("αβ", 7, "¹₂³") == "¹₂³¹₂αβ"
    @test rpad("αβ", 7, "¹₂³") == "αβ¹₂³¹₂"
    @test lpad("αβ", 8, "¹₂³") == "¹₂³¹₂³αβ"
    @test rpad("αβ", 8, "¹₂³") == "αβ¹₂³¹₂³"
    @test lpad("αβ", 9, "¹₂³") == "¹₂³¹₂³¹αβ"
    @test rpad("αβ", 9, "¹₂³") == "αβ¹₂³¹₂³¹"
end

# string manipulation
@testset "lstrip/rstrip/strip" begin
    @test strip("") == ""
    @test strip(" ") == ""
    @test strip("  ") == ""
    @test strip("   ") == ""
    @test strip("\t  hi   \n") == "hi"
    @test strip("foobarfoo", ['f','o']) == "bar"
    @test strip("foobarfoo", ('f','o')) == "bar"

    for s in ("", " ", " abc", "abc ", "  abc  "),
        f in (lstrip, rstrip, strip)

        fs = f(s)
        for T = (String, GenericString)
            local t, b
            t = convert(T,s)
            ft = f(t)
            @test s == t
            @test fs == ft
            @test typeof(ft) == SubString{T}

            b = convert(SubString{T}, t)
            fb = f(b)
            @test s == b
            @test fs == fb
            @test typeof(fb) == SubString{T}
        end
    end
end

@testset "rsplit/split" begin
    @test split("foo,bar,baz", 'x') == ["foo,bar,baz"]
    @test split("foo,bar,baz", ',') == ["foo","bar","baz"]
    @test split("foo,bar,baz", ",") == ["foo","bar","baz"]
    @test split("foo,bar,baz", r",") == ["foo","bar","baz"]
    @test split("foo,bar,baz", ','; limit=0) == ["foo","bar","baz"]
    @test split("foo,bar,baz", ','; limit=1) == ["foo,bar,baz"]
    @test split("foo,bar,baz", ','; limit=2) == ["foo","bar,baz"]
    @test split("foo,bar,baz", ','; limit=3) == ["foo","bar","baz"]
    @test split("foo,bar", "o,b") == ["fo","ar"]

    @test split("", ',') == [""]
    @test split(",", ',') == ["",""]
    @test split(",,", ',') == ["","",""]
    @test split("", ','  ; keepempty=false) == []
    @test split(",", ',' ; keepempty=false) == []
    @test split(",,", ','; keepempty=false) == []

    @test split("a b c") == ["a","b","c"]
    @test split("a  b \t c\n") == ["a","b","c"]

    @test split("a b c"; limit=2) == ["a","b c"]
    @test split("a  b \t c\n"; limit=3) == ["a","b","\t c\n"]
    @test split("a b c"; keepempty=true) == ["a","b","c"]
    @test split("a  b \t c\n"; keepempty=true) == ["a","","b","","","c",""]

    @test rsplit("foo,bar,baz", 'x') == ["foo,bar,baz"]
    @test rsplit("foo,bar,baz", ',') == ["foo","bar","baz"]
    @test rsplit("foo,bar,baz", ",") == ["foo","bar","baz"]
    @test rsplit("foo,bar,baz", ','; limit=0) == ["foo","bar","baz"]
    @test rsplit("foo,bar,baz", ','; limit=1) == ["foo,bar,baz"]
    @test rsplit("foo,bar,baz", ','; limit=2) == ["foo,bar","baz"]
    @test rsplit("foo,bar,baz", ','; limit=3) == ["foo","bar","baz"]
    @test rsplit("foo,bar", "o,b") == ["fo","ar"]

    @test rsplit("", ',') == [""]
    @test rsplit(",", ',') == ["",""]
    @test rsplit(",,", ',') == ["","",""]
    @test rsplit(",,", ','; limit=2) == [",",""]
    @test rsplit("", ','  ; keepempty=false) == []
    @test rsplit(",", ',' ; keepempty=false) == []
    @test rsplit(",,", ','; keepempty=false) == []

    @test rsplit("a b c") == ["a","b","c"]
    @test rsplit("a  b \t c\n") == ["a","b","c"]

    @test rsplit("a b c"; limit=2) == ["a b", "c"]
    @test rsplit("a  b \t c\n"; limit=3) == ["a ","b","c"]
    @test rsplit("a b c"; keepempty=true) == ["a","b","c"]
    @test rsplit("a  b \t c\n"; keepempty=true) == ["a","","b","","","c",""]

    let str = "a.:.ba..:..cba.:.:.dcba.:."
    @test split(str, ".:.") == ["a","ba.",".cba",":.dcba",""]
    @test split(str, ".:."; keepempty=false) == ["a","ba.",".cba",":.dcba"]
    @test split(str, ".:.") == ["a","ba.",".cba",":.dcba",""]
    @test split(str, r"\.(:\.)+") == ["a","ba.",".cba","dcba",""]
    @test split(str, r"\.(:\.)+"; keepempty=false) == ["a","ba.",".cba","dcba"]
    @test split(str, r"\.+:\.+") == ["a","ba","cba",":.dcba",""]
    @test split(str, r"\.+:\.+"; keepempty=false) == ["a","ba","cba",":.dcba"]

    @test rsplit(str, ".:.") == ["a","ba.",".cba.:","dcba",""]
    @test rsplit(str, ".:."; keepempty=false) == ["a","ba.",".cba.:","dcba"]
    @test rsplit(str, ".:."; limit=2) == ["a.:.ba..:..cba.:.:.dcba", ""]
    @test rsplit(str, ".:."; limit=3) == ["a.:.ba..:..cba.:", "dcba", ""]
    @test rsplit(str, ".:."; limit=4) == ["a.:.ba.", ".cba.:", "dcba", ""]
    @test rsplit(str, ".:."; limit=5) == ["a", "ba.", ".cba.:", "dcba", ""]
    @test rsplit(str, ".:."; limit=6) == ["a", "ba.", ".cba.:", "dcba", ""]
    end

    # zero-width splits
    @test split("", "") == rsplit("", "") == [""]
    @test split("abc", "") == rsplit("abc", "") == ["a","b","c"]
    @test rsplit("abc", "", limit=2) == ["ab","c"]
    @test split("abc", "", limit=2) == ["a","bc"]

    @test split("", r"") == [""]
    @test split("abc", r"") == ["a","b","c"]
    @test split("abcd", r"b?") == ["a","c","d"]
    @test split("abcd", r"b*") == ["a","c","d"]
    @test split("abcd", r"b+") == ["a","cd"]
    @test split("abcd", r"b?c?") == ["a","d"]
    @test split("abcd", r"[bc]?") == ["a","","d"]
    @test split("abcd", r"a*") == ["","b","c","d"]
    @test split("abcd", r"a+") == ["","bcd"]
    @test split("abcd", r"d*") == ["a","b","c",""]
    @test split("abcd", r"d+") == ["abc",""]
    @test split("abcd", r"[ad]?") == ["","b","c",""]

    # multi-byte unicode characters (issue #26225)
    @test split("α β γ", " ") == rsplit("α β γ", " ") ==
          split("α β γ", isspace) == rsplit("α β γ", isspace) == ["α","β","γ"]
    @test split("ö.", ".") == rsplit("ö.", ".") == ["ö",""]
    @test split("α β γ", "β") == rsplit("α β γ", "β") == ["α "," γ"]
end

@testset "replace" begin
    @test replace("\u2202", '*' => '\0') == "\u2202"

    @test replace("foobar", 'o' => '0') == "f00bar"
    @test replace("foobar", 'o' => '0', count=1) == "f0obar"
    @test replace("foobar", 'o' => "") == "fbar"
    @test replace("foobar", 'o' => "", count=1) == "fobar"
    @test replace("foobar", 'f' => 'F') == "Foobar"
    @test replace("foobar", 'r' => 'R') == "foobaR"

    @test replace("foofoofoo", "foo" => "bar") == "barbarbar"
    @test replace("foobarfoo", "foo" => "baz") == "bazbarbaz"
    @test replace("barfoofoo", "foo" => "baz") == "barbazbaz"

    @test replace("", "" => "") == ""
    @test replace("", "" => "x") == "x"
    @test replace("", "x" => "y") == ""

    @test replace("abcd", "" => "^") == "^a^b^c^d^"
    @test replace("abcd", "b" => "^") == "a^cd"
    @test replace("abcd", r"b?" => "^") == "^a^c^d^"
    @test replace("abcd", r"b+" => "^") == "a^cd"
    @test replace("abcd", r"b?c?" => "^") == "^a^d^"
    @test replace("abcd", r"[bc]?" => "^") == "^a^^d^"

    @test replace("foobarfoo", r"(fo|ba)" => "xx") == "xxoxxrxxo"
    @test replace("foobarfoo", r"(foo|ba)" => "bar") == "barbarrbar"

    @test replace("foobar", 'o' => 'ø') == "føøbar"
    @test replace("foobar", 'o' => 'ø', count=1) == "føobar"
    @test replace("føøbar", 'ø' => 'o') == "foobar"
    @test replace("føøbar", 'ø' => 'o', count=1) == "foøbar"
    @test replace("føøbar", 'ø' => 'ö') == "fööbar"
    @test replace("føøbar", 'ø' => 'ö', count=1) == "föøbar"
    @test replace("føøbar", 'ø' => "") == "fbar"
    @test replace("føøbar", 'ø' => "", count=1) == "føbar"
    @test replace("føøbar", 'f' => 'F') == "Føøbar"
    @test replace("ḟøøbar", 'ḟ' => 'F') == "Føøbar"
    @test replace("føøbar", 'f' => 'Ḟ') == "Ḟøøbar"
    @test replace("ḟøøbar", 'ḟ' => 'Ḟ') == "Ḟøøbar"
    @test replace("føøbar", 'r' => 'R') == "føøbaR"
    @test replace("føøbaṙ", 'ṙ' => 'R') == "føøbaR"
    @test replace("føøbar", 'r' => 'Ṙ') == "føøbaṘ"
    @test replace("føøbaṙ", 'ṙ' => 'Ṙ') == "føøbaṘ"

    @test replace("ḟøøḟøøḟøø", "ḟøø" => "bar") == "barbarbar"
    @test replace("ḟøøbarḟøø", "ḟøø" => "baz") == "bazbarbaz"
    @test replace("barḟøøḟøø", "ḟøø" => "baz") == "barbazbaz"

    @test replace("foofoofoo", "foo" => "ƀäṙ") == "ƀäṙƀäṙƀäṙ"
    @test replace("fooƀäṙfoo", "foo" => "baz") == "bazƀäṙbaz"
    @test replace("ƀäṙfoofoo", "foo" => "baz") == "ƀäṙbazbaz"

    @test replace("foofoofoo", "foo" => "bar") == "barbarbar"
    @test replace("foobarfoo", "foo" => "ƀäż") == "ƀäżbarƀäż"
    @test replace("barfoofoo", "foo" => "ƀäż") == "barƀäżƀäż"

    @test replace("ḟøøḟøøḟøø", "ḟøø" => "ƀäṙ") == "ƀäṙƀäṙƀäṙ"
    @test replace("ḟøøƀäṙḟøø", "ḟøø" => "baz") == "bazƀäṙbaz"
    @test replace("ƀäṙḟøøḟøø", "ḟøø" => "baz") == "ƀäṙbazbaz"

    @test replace("ḟøøḟøøḟøø", "ḟøø" => "bar") == "barbarbar"
    @test replace("ḟøøbarḟøø", "ḟøø" => "ƀäż") == "ƀäżbarƀäż"
    @test replace("barḟøøḟøø", "ḟøø" => "ƀäż") == "barƀäżƀäż"

    @test replace("ḟøøḟøøḟøø", "ḟøø" => "ƀäṙ") == "ƀäṙƀäṙƀäṙ"
    @test replace("ḟøøƀäṙḟøø", "ḟøø" => "ƀäż") == "ƀäżƀäṙƀäż"
    @test replace("ƀäṙḟøøḟøø", "ḟøø" => "ƀäż") == "ƀäṙƀäżƀäż"

    @test replace("", "" => "ẍ") == "ẍ"
    @test replace("", "ẍ" => "ÿ") == ""

    @test replace("äƀçđ", "" => "π") == "πäπƀπçπđπ"
    @test replace("äƀçđ", "ƀ" => "π") == "äπçđ"
    @test replace("äƀçđ", r"ƀ?" => "π") == "πäπçπđπ"
    @test replace("äƀçđ", r"ƀ+" => "π") == "äπçđ"
    @test replace("äƀçđ", r"ƀ?ç?" => "π") == "πäπđπ"
    @test replace("äƀçđ", r"[ƀç]?" => "π") == "πäππđπ"

    @test replace("foobarfoo", r"(fo|ba)" => "ẍẍ") == "ẍẍoẍẍrẍẍo"

    @test replace("ḟøøbarḟøø", r"(ḟø|ba)" => "xx") == "xxøxxrxxø"
    @test replace("ḟøøbarḟøø", r"(ḟøø|ba)" => "bar") == "barbarrbar"

    @test replace("fooƀäṙfoo", r"(fo|ƀä)" => "xx") == "xxoxxṙxxo"
    @test replace("fooƀäṙfoo", r"(foo|ƀä)" => "ƀäṙ") == "ƀäṙƀäṙṙƀäṙ"

    @test replace("ḟøøƀäṙḟøø", r"(ḟø|ƀä)" => "xx") == "xxøxxṙxxø"
    @test replace("ḟøøƀäṙḟøø", r"(ḟøø|ƀä)" => "ƀäṙ") == "ƀäṙƀäṙṙƀäṙ"

    @test replace("foo", "oo" => uppercase) == "fOO"

    # Issue 13332
    @test replace("abc", 'b' => 2.1) == "a2.1c"

    # test replace with a count for String and GenericString
    # check that replace is a no-op if count==0
    for s in ["aaa", Test.GenericString("aaa")]
        # @test replace("aaa", 'a' => 'z', count=0) == "aaa" # enable when undeprecated
        @test replace(s, 'a' => 'z', count=1) == "zaa"
        @test replace(s, 'a' => 'z', count=2) == "zza"
        @test replace(s, 'a' => 'z', count=3) == "zzz"
        @test replace(s, 'a' => 'z', count=4) == "zzz"
        @test replace(s, 'a' => 'z', count=typemax(Int)) == "zzz"
        @test replace(s, 'a' => 'z')    == "zzz"
    end

    # Issue 25741
    @test replace("abc", ['a', 'd'] => 'A') == "Abc"

    # for Char pattern call Char replacement function
    @test replace("a", "a" => typeof) == "SubString{String}"
    @test replace("a", r"a" => typeof) == "SubString{String}"
    @test replace("a", 'a' => typeof) == "Char"
    @test replace("a", in("a") => typeof) == "Char"
    @test replace("a", ['a'] => typeof) == "Char"

end

@testset "chomp/chop" begin
    @test chomp("foo\n") == "foo"
    @test chomp("fo∀\n") == "fo∀"
    @test chomp("foo\r\n") == "foo"
    @test chomp("fo∀\r\n") == "fo∀"
    @test chomp("fo∀") == "fo∀"
    @test chop("fooε") == "foo"
    @test chop("foεo") == "foε"
    @test chop("∃∃∃∃") == "∃∃∃"
    @test chop("∀ϵ∃Δ", head=0, tail=0) == "∀ϵ∃Δ"
    @test chop("∀ϵ∃Δ", head=0, tail=1) == "∀ϵ∃"
    @test chop("∀ϵ∃Δ", head=0, tail=2) == "∀ϵ"
    @test chop("∀ϵ∃Δ", head=0, tail=3) == "∀"
    @test chop("∀ϵ∃Δ", head=0, tail=4) == ""
    @test chop("∀ϵ∃Δ", head=0, tail=5) == ""
    @test chop("∀ϵ∃Δ", head=1, tail=0) == "ϵ∃Δ"
    @test chop("∀ϵ∃Δ", head=2, tail=0) == "∃Δ"
    @test chop("∀ϵ∃Δ", head=3, tail=0) == "Δ"
    @test chop("∀ϵ∃Δ", head=4, tail=0) == ""
    @test chop("∀ϵ∃Δ", head=5, tail=0) == ""
    @test chop("∀ϵ∃Δ", head=1, tail=1) == "ϵ∃"
    @test chop("∀ϵ∃Δ", head=2, tail=2) == ""
    @test chop("∀ϵ∃Δ", head=3, tail=3) == ""
    @test_throws ArgumentError chop("∀ϵ∃Δ", head=-3, tail=3)
    @test_throws ArgumentError chop("∀ϵ∃Δ", head=3, tail=-3)
    @test_throws ArgumentError chop("∀ϵ∃Δ", head=-3, tail=-3)

    @test isa(chomp("foo"), SubString)
    @test isa(chop("foo"), SubString)
end

@testset "bytes2hex and hex2bytes" begin
    hex_str = "d7a8fbb307d7809469ca9abcb0082e4f8d5651e46d3cdb762d02d0bf37c9e592"
    bin_val = hex2bytes(hex_str)

    @test div(length(hex_str), 2) == length(bin_val)
    @test hex_str == bytes2hex(bin_val)

    bin_val = hex2bytes("07bf")
    @test bin_val[1] == 7
    @test bin_val[2] == 191
    @test typeof(bin_val) == Array{UInt8, 1}
    @test length(bin_val) == 2

    # all valid hex chars
    @test "0123456789abcdefabcdef" == bytes2hex(hex2bytes("0123456789abcdefABCDEF"))

    # odd size
    @test_throws ArgumentError hex2bytes("0123456789abcdefABCDEF0")

    #non-hex characters
    @test_throws ArgumentError hex2bytes("0123456789abcdefABCDEFGH")

    @testset "Issue 23161" begin
        arr = b"0123456789abcdefABCDEF"
        arr1 = Vector{UInt8}(undef, length(arr) >> 1)
        @test hex2bytes!(arr1, arr) === arr1 # check in-place
        @test "0123456789abcdefabcdef" == bytes2hex(arr1)
        @test hex2bytes("0123456789abcdefABCDEF") == hex2bytes(arr)
        @test_throws ArgumentError hex2bytes!(arr1, b"") # incorrect arr1 length
        @test hex2bytes(b"") == UInt8[]
        @test hex2bytes(view(b"012345",1:6)) == UInt8[0x01,0x23,0x45]
        @test begin
            s = view(b"012345ab",1:6)
            d = view(zeros(UInt8, 10),1:3)
            hex2bytes!(d,s) == UInt8[0x01,0x23,0x45]
        end
        # odd size
        @test_throws ArgumentError hex2bytes(b"0123456789abcdefABCDEF0")

        #non-hex characters
        @test_throws ArgumentError hex2bytes(b"0123456789abcdefABCDEFGH")
    end
end

# b"" should be immutable
let testb() = b"0123"
    b = testb()
    @test eltype(b) === UInt8
    @test b isa AbstractVector
    @test_throws ErrorException b[4] = '4'
    @test testb() == UInt8['0','1','2','3']
end

max_stars_repo_name: djsegal/julia-fork

max_issues_repo_name: djsegal/julia-fork

max_forks_repo_name: djsegal/julia-fork


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
using BinaryProvider

name = "Xorg"
version = v"7.5"
# url = "https://mirrors.edge.kernel.org/pub/linux/utils/util-linux/v2.34/util-linux-2.34-rc1.tar.gz"
# path = download(url)
# hash = open(path) do io
#     bytes2hex(BinaryProvider.sha256(io))
# end

# libs = [
#     "bigreqsproto-1.1.2.tar.bz2",
#     "compositeproto-0.4.2.tar.bz2",
#     "damageproto-1.2.1.tar.bz2",
#     "dmxproto-2.3.1.tar.bz2",
#     "dri2proto-2.8.tar.bz2",
#     "dri3proto-1.0.tar.bz2",
#     "fixesproto-5.0.tar.bz2",
#     "fontsproto-2.1.3.tar.bz2",
#     "glproto-1.4.17.tar.bz2",
#     "inputproto-2.3.2.tar.bz2",
#     "kbproto-1.0.7.tar.bz2",
#     "presentproto-1.1.tar.bz2",
#     "randrproto-1.5.0.tar.bz2",
#     "recordproto-1.14.2.tar.bz2",
#     "renderproto-0.11.1.tar.bz2",
#     "resourceproto-1.2.0.tar.bz2",
#     "scrnsaverproto-1.2.2.tar.bz2",
#     "videoproto-2.3.3.tar.bz2",
#     "xcmiscproto-1.2.2.tar.bz2",
#     "xextproto-7.3.0.tar.bz2",
#     "xf86bigfontproto-1.2.0.tar.bz2",
#     "xf86dgaproto-2.1.tar.bz2",
#     "xf86driproto-2.1.1.tar.bz2",
#     "xf86vidmodeproto-2.3.1.tar.bz2",
#     "xineramaproto-1.2.1.tar.bz2",
#     "xproto-7.0.31.tar.bz2",
# ]
# open("deps.jl", "w") do fio
#     println(fio, "sources = [")
#     for lib in libs
#         url = "https://www.x.org/pub/individual/proto/" * lib
#         println(url)
#         path = download(url)
#         hash = open(path) do io
#             bytes2hex(BinaryProvider.sha256(io))
#         end
#         println(fio, repr(url), " => ", repr(hash), ",")
#     end
#     println(fio, "]")
# end
#> Collection of sources required to build Fontconfig
#

sources = [
    "http://xcb.freedesktop.org/dist/libpthread-stubs-0.3.tar.bz2" => "35b6d54e3cc6f3ba28061da81af64b9a92b7b757319098172488a660e3d87299",
    "https://ftp.gnu.org/gnu/pth/pth-2.0.7.tar.gz" => "72353660c5a2caafd601b20e12e75d865fd88f6cf1a088b306a3963f0bc77232",
    "http://xorg.freedesktop.org/releases/individual/util/util-macros-1.19.0.tar.bz2" => "2835b11829ee634e19fa56517b4cfc52ef39acea0cd82e15f68096e27cbed0ba",
    "https://ftp.gnu.org/gnu/gettext/gettext-0.19.tar.gz" => "a9f9c7764bc68782aca57687a1f16c228bc321069c538325923513b841bd847c",
    # xorg protocol headers
    "https://www.x.org/pub/individual/proto/bigreqsproto-1.1.2.tar.bz2" => "462116ab44e41d8121bfde947321950370b285a5316612b8fce8334d50751b1e",
    "https://www.x.org/pub/individual/proto/compositeproto-0.4.2.tar.bz2" => "049359f0be0b2b984a8149c966dd04e8c58e6eade2a4a309cf1126635ccd0cfc",
    "https://www.x.org/pub/individual/proto/damageproto-1.2.1.tar.bz2" => "5c7c112e9b9ea8a9d5b019e5f17d481ae20f766cb7a4648360e7c1b46fc9fc5b",
    "https://www.x.org/pub/individual/proto/dmxproto-2.3.1.tar.bz2" => "e72051e6a3e06b236d19eed56368117b745ca1e1a27bdc50fd51aa375bea6509",
    "https://www.x.org/pub/individual/proto/dri2proto-2.8.tar.bz2" => "f9b55476def44fc7c459b2537d17dbc731e36ed5d416af7ca0b1e2e676f8aa04",
    "https://www.x.org/pub/individual/proto/dri3proto-1.0.tar.bz2" => "01be49d70200518b9a6b297131f6cc71f4ea2de17436896af153226a774fc074",
    "https://www.x.org/pub/individual/proto/fixesproto-5.0.tar.bz2" => "ba2f3f31246bdd3f2a0acf8bd3b09ba99cab965c7fb2c2c92b7dc72870e424ce",
    "https://www.x.org/pub/individual/proto/fontsproto-2.1.3.tar.bz2" => "259046b0dd9130825c4a4c479ba3591d6d0f17a33f54e294b56478729a6e5ab8",
    "https://www.x.org/pub/individual/proto/glproto-1.4.17.tar.bz2" => "adaa94bded310a2bfcbb9deb4d751d965fcfe6fb3a2f6d242e2df2d6589dbe40",
    "https://www.x.org/pub/individual/proto/inputproto-2.3.2.tar.bz2" => "893a6af55733262058a27b38eeb1edc733669f01d404e8581b167f03c03ef31d",
    "https://www.x.org/pub/individual/proto/kbproto-1.0.7.tar.bz2" => "f882210b76376e3fa006b11dbd890e56ec0942bc56e65d1249ff4af86f90b857",
    "https://www.x.org/pub/individual/proto/presentproto-1.1.tar.bz2" => "f69b23a8869f78a5898aaf53938b829c8165e597cda34f06024d43ee1e6d26b9",
    "https://www.x.org/pub/individual/proto/randrproto-1.5.0.tar.bz2" => "4c675533e79cd730997d232c8894b6692174dce58d3e207021b8f860be498468",
    "https://www.x.org/pub/individual/proto/recordproto-1.14.2.tar.bz2" => "a777548d2e92aa259f1528de3c4a36d15e07a4650d0976573a8e2ff5437e7370",
    "https://www.x.org/pub/individual/proto/renderproto-0.11.1.tar.bz2" => "06735a5b92b20759204e4751ecd6064a2ad8a6246bb65b3078b862a00def2537",
    "https://www.x.org/pub/individual/proto/resourceproto-1.2.0.tar.bz2" => "3c66003a6bdeb0f70932a9ed3cf57cc554234154378d301e0c5cfa189d8f6818",
    "https://www.x.org/pub/individual/proto/scrnsaverproto-1.2.2.tar.bz2" => "8bb70a8da164930cceaeb4c74180291660533ad3cc45377b30a795d1b85bcd65",
    "https://www.x.org/pub/individual/proto/videoproto-2.3.3.tar.bz2" => "c7803889fd08e6fcaf7b68cc394fb038b2325d1f315e571a6954577e07cca702",
    "https://www.x.org/pub/individual/proto/xcmiscproto-1.2.2.tar.bz2" => "b13236869372256c36db79ae39d54214172677fb79e9cdc555dceec80bd9d2df",
    "https://www.x.org/pub/individual/proto/xextproto-7.3.0.tar.bz2" => "f3f4b23ac8db9c3a9e0d8edb591713f3d70ef9c3b175970dd8823dfc92aa5bb0",
    "https://www.x.org/pub/individual/proto/xf86bigfontproto-1.2.0.tar.bz2" => "ba9220e2c4475f5ed2ddaa7287426b30089e4d29bd58d35fad57ba5ea43e1648",
    "https://www.x.org/pub/individual/proto/xf86dgaproto-2.1.tar.bz2" => "ac5ef65108e1f2146286e53080975683dae49fc94680042e04bd1e2010e99050",
    "https://www.x.org/pub/individual/proto/xf86driproto-2.1.1.tar.bz2" => "9c4b8d7221cb6dc4309269ccc008a22753698ae9245a398a59df35f1404d661f",
    "https://www.x.org/pub/individual/proto/xf86vidmodeproto-2.3.1.tar.bz2" => "45d9499aa7b73203fd6b3505b0259624afed5c16b941bd04fcf123e5de698770",
    "https://www.x.org/pub/individual/proto/xineramaproto-1.2.1.tar.bz2" => "977574bb3dc192ecd9c55f59f991ec1dff340be3e31392c95deff423da52485b",
    "https://www.x.org/pub/individual/proto/xproto-7.0.31.tar.bz2" => "c6f9747da0bd3a95f86b17fb8dd5e717c8f3ab7f0ece3ba1b247899ec1ef7747",

    "http://downloads.sourceforge.net/expat/expat-2.1.0.tar.gz" => "823705472f816df21c8f6aa026dd162b280806838bb55b3432b0fb1fcca7eb86",
    "https://mirrors.edge.kernel.org/pub/linux/utils/util-linux/v2.34/util-linux-2.34-rc1.tar.gz" => "f97d7725509f53c690c3c65c4432f97857ca398651fb06fded61fa17c19bef9f",

    "http://xcb.freedesktop.org/dist/xcb-proto-1.11.tar.bz2" => "b4aceee6502a0ce45fc39b33c541a2df4715d00b72e660ebe8c5bb444771e32e",
    "http://xorg.freedesktop.org/releases/individual/lib/libXau-1.0.8.tar.bz2" => "fdd477320aeb5cdd67272838722d6b7d544887dfe7de46e1e7cc0c27c2bea4f2",
    "http://xcb.freedesktop.org/dist/libxcb-1.11.tar.bz2" => "03635d70045b9ede90778e67516135828a57de87ac508f987024f43c03620ff7",
    "http://xorg.freedesktop.org/releases/individual/lib/xtrans-1.3.4.tar.bz2" => "054d4ee3efd52508c753e9f7bc655ef185a29bd2850dd9e2fc2ccc33544f583a",
    "http://xorg.freedesktop.org/releases/individual/lib/libX11-1.6.2.tar.bz2" => "2aa027e837231d2eeea90f3a4afe19948a6eb4c8b2bec0241eba7dbc8106bd16",
    "http://xorg.freedesktop.org/releases/individual/lib/libXext-1.3.3.tar.bz2" => "b518d4d332231f313371fdefac59e3776f4f0823bcb23cf7c7305bfb57b16e35",
    "http://xorg.freedesktop.org/releases/individual/lib/libFS-1.0.6.tar.bz2" => "9007fcfbd37bfa6ede4000ec4b5c2dc6dc0f52789358361b43f68e9ae4d6aed7",
    "http://xorg.freedesktop.org/releases/individual/lib/libICE-1.0.9.tar.bz2" => "8f7032f2c1c64352b5423f6b48a8ebdc339cc63064af34d66a6c9aa79759e202",
    "http://xorg.freedesktop.org/releases/individual/lib/libSM-1.2.2.tar.bz2" => "0baca8c9f5d934450a70896c4ad38d06475521255ca63b717a6510fdb6e287bd",
    "http://xorg.freedesktop.org/releases/individual/lib/libXScrnSaver-1.2.2.tar.bz2" => "8ff1efa7341c7f34bcf9b17c89648d6325ddaae22e3904e091794e0b4426ce1d",
    "http://xorg.freedesktop.org/releases/individual/lib/libXt-1.1.4.tar.bz2" => "843a97a988f5654872682a4120486d987d853a71651515472f55519ffae2dd57",
    "http://xorg.freedesktop.org/releases/individual/lib/libXmu-1.1.2.tar.bz2" => "756edc7c383254eef8b4e1b733c3bf1dc061b523c9f9833ac7058378b8349d0b",
    "http://xorg.freedesktop.org/releases/individual/lib/libXpm-3.5.11.tar.bz2" => "c5bdafa51d1ae30086fac01ab83be8d47fe117b238d3437f8e965434090e041c",
    "http://xorg.freedesktop.org/releases/individual/lib/libXaw-1.0.12.tar.bz2" => "96fc314874fce9979556321d1d6ee00b5baf32fb333b7278853b4983bc3cdbf6",
    "http://xorg.freedesktop.org/releases/individual/lib/libXfixes-5.0.1.tar.bz2" => "63bec085084fa3caaee5180490dd871f1eb2020ba9e9b39a30f93693ffc34767",
    "http://xorg.freedesktop.org/releases/individual/lib/libXcomposite-0.4.4.tar.bz2" => "ede250cd207d8bee4a338265c3007d7a68d5aca791b6ac41af18e9a2aeb34178",
    "http://xorg.freedesktop.org/releases/individual/lib/libXrender-0.9.8.tar.bz2" => "1d14b02f0060aec5d90dfdcf16a996f17002e515292906ed26e3dcbba0f4fc62",
    "http://xorg.freedesktop.org/releases/individual/lib/libXcursor-1.1.14.tar.bz2" => "9bc6acb21ca14da51bda5bc912c8955bc6e5e433f0ab00c5e8bef842596c33df",
    "http://xorg.freedesktop.org/releases/individual/lib/libXdamage-1.1.4.tar.bz2" => "7c3fe7c657e83547f4822bfde30a90d84524efb56365448768409b77f05355ad",
    "http://xorg.freedesktop.org/releases/individual/lib/libfontenc-1.1.2.tar.bz2" => "a9a4efed3359b2e80161bb66b65038fac145137fa134e71335264cbc23b02f62",
    "http://xorg.freedesktop.org/releases/individual/lib/libXfont-1.5.0.tar.bz2" => "3a3c52c4adf9352b2160f07ff0596af17ab14f91d6509564e606678a1261c25f",
    "http://xorg.freedesktop.org/releases/individual/lib/libXft-2.3.2.tar.bz2" => "f5a3c824761df351ca91827ac221090943ef28b248573486050de89f4bfcdc4c",
    "http://xorg.freedesktop.org/releases/individual/lib/libXi-1.7.4.tar.bz2" => "2cffc2686618dc1803725636cd92b36342c512dc60a7a35cba34bf7192a42244",
    "http://xorg.freedesktop.org/releases/individual/lib/libXinerama-1.1.3.tar.bz2" => "7a45699f1773095a3f821e491cbd5e10c887c5a5fce5d8d3fced15c2ff7698e2",
    "http://xorg.freedesktop.org/releases/individual/lib/libXrandr-1.4.2.tar.bz2" => "caa7b31ac769be51a532343c65376f1d4df3f307afaed58e34fb5e82e8b825ad",
    "http://xorg.freedesktop.org/releases/individual/lib/libXres-1.0.7.tar.bz2" => "26899054aa87f81b17becc68e8645b240f140464cf90c42616ebb263ec5fa0e5",
    "http://xorg.freedesktop.org/releases/individual/lib/libXtst-1.2.2.tar.bz2" => "ef0a7ffd577e5f1a25b1663b375679529663a1880151beaa73e9186c8309f6d9",
    "http://xorg.freedesktop.org/releases/individual/lib/libXv-1.0.10.tar.bz2" => "55fe92f8686ce8612e2c1bfaf58c057715534419da700bda8d517b1d97914525",
    "http://xorg.freedesktop.org/releases/individual/lib/libXvMC-1.0.8.tar.bz2" => "5e1a401efa433f959d41e17932b8c218c56b931348f494b8fa4656d7d798b204",
    "http://xorg.freedesktop.org/releases/individual/lib/libXxf86dga-1.1.4.tar.bz2" => "8eecd4b6c1df9a3704c04733c2f4fa93ef469b55028af5510b25818e2456c77e",
    "http://xorg.freedesktop.org/releases/individual/lib/libXxf86vm-1.1.3.tar.bz2" => "da5e86c32ee2069b9e6d820e4c2e4242d4877cb155a2b2fbf2675a1480ec37b8",
    "http://xorg.freedesktop.org/releases/individual/lib/libdmx-1.1.3.tar.bz2" => "c97da36d2e56a2d7b6e4f896241785acc95e97eb9557465fd66ba2a155a7b201",
    "http://xorg.freedesktop.org/releases/individual/lib/libpciaccess-0.13.2.tar.bz2" => "abab8c2b050afb89bc29280e9c6b50ec5867174960d88bfb81a01036ec20de19",
    "http://xorg.freedesktop.org/releases/individual/lib/libxkbfile-1.0.8.tar.bz2" => "8aa94e19c537c43558f30906650cea6e15fa012591445d9f927658c3b32a8f3a",
    "http://xorg.freedesktop.org/releases/individual/lib/libxshmfence-1.1.tar.bz2" => "dbc2db2925ca9f216fd1e9c63d0974db9f4d49aaf5877ffb606d2d8d7e58cebe",
]

# # Bash recipe for building across all platforms
script = sprint() do io
    println(io, raw"cd $WORKSPACE/srcdir")
    println(io, raw"""
    pathappend () {
        local PATHVARIABLE=${2:-PATH}
        export $PATHVARIABLE="${!PATHVARIABLE:+${!PATHVARIABLE}:}$1"
    }
    export XORG_PREFIX=$prefix
    export XORG_CONFIG="--prefix=$prefix --host=$target --sysconfdir=/etc --localstatedir=/var --disable-static LANG=C"

    pathappend $XORG_PREFIX/bin             PATH
    pathappend $XORG_PREFIX/lib/pkgconfig   PKG_CONFIG_PATH
    pathappend $XORG_PREFIX/share/pkgconfig PKG_CONFIG_PATH

    pathappend $XORG_PREFIX/lib             LIBRARY_PATH
    pathappend $XORG_PREFIX/include         C_INCLUDE_PATH
    pathappend $XORG_PREFIX/include         CPLUS_INCLUDE_PATH
    LINGUAS="de nl en"

    ACLOCAL='aclocal -I $XORG_PREFIX/share/aclocal'

    export PATH PKG_CONFIG_PATH ACLOCAL LIBRARY_PATH C_INCLUDE_PATH CPLUS_INCLUDE_PATH
    """)
    for (url, hash) in sources
        name = replace(basename(url), r"(.tar.bz2)|(.tar.gz)" => "")
        println(io, "cd $name")
        if match(r"pth-[0-9]*", name) !== nothing
            println(io, raw"sed -i 's#$(LOBJS): Makefile#$(LOBJS): pth_p.h Makefile#' Makefile.in")
        end
        if name == "libxb"
            println(io, raw"""
            sed "s/pthread-stubs//" -i configure
            ./configure $XORG_CONFIG --enable-xinput --docdir='${datadir}'/doc/libxcb-1.11
            """)
        elseif match(r"libXfont-[0-9]*", name) !== nothing
            println(io, raw"./configure $XORG_CONFIG --disable-devel-docs")
        elseif match(r"libXt-[0-9]*", name) !== nothing
            println(io, raw"./configure $XORG_CONFIG --with-appdefaultdir=/etc/X11/app-defaults")
        else
            println(io, raw" ./configure $XORG_CONFIG")
        end
        println(io, raw"""
        make
        make install
        cd ..
        """)
    end
end

# println(script)
# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:x86_64, libc=:glibc)
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libxorg", :libxorg),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "https://github.com/JuliaGraphics/FreeTypeBuilder/releases/download/v2.9.1-1/build_FreeType2.v2.9.1.jl",
    "https://github.com/bicycle1885/ZlibBuilder/releases/download/v1.0.2/build_Zlib.v1.2.11.jl",
    "https://github.com/bjarthur/FontconfigBuilder/releases/download/v2.13.1/build_Fontconfig.v2.13.1.jl",
]

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

max_stars_repo_name: JuliaGraphics/XorgBuilder

max_issues_repo_name: JuliaGraphics/XorgBuilder

max_forks_repo_name: JuliaGraphics/XorgBuilder


text:
module MomentExpansions
    using DataStructures
    import SymPy
    import SymEngine
    import Base.zero, Base.one, Base.show
    using Reexport
    @reexport using DiffEqBase
    @reexport using DiffEqBiological

    include("types.jl")
    include("utils.jl")
    include("meanmoments.jl")
    include("centralmoments.jl")

    """
		approx(network::AbstractReactionNetwork, order::Int; algebra::Algebra)

	Given an AbstractReactionNetwork, and an expansion order, returns the standard moment expansion.
    To switch between using SymPy or SymEngine as the backend, use the kwargs `algebra=sympy()` or `algebra=symengine()`.
	"""
    function approx(network::DiffEqBase.AbstractReactionNetwork; order::Int=2, algebra::Algebra=symengine())
        momentstruct = gen_moment_syms(network, order, algebra=algebra)
        dmean = dμ_dt(network, momentstruct, algebra=algebra)
        dcen = dcen_dt(network, momentstruct, dmean, algebra=algebra)
        ode_expr, ode_func = get_ode(network, momentstruct, dmean, dcen)
        MomentExpansion(network, momentstruct, dmean, dcen, ode_expr, ode_func)
    end

    function DiffEqBase.ODEProblem(momentexpansion::MomentExpansion, u0::Union{AbstractArray, Number}, args...; kwargs...)
        DiffEqBase.ODEProblem(momentexpansion.ode, u0, args...; kwargs...)
    end

    export approx, sympy, symengine
end


text:
using Immerse, Colors
using Test

include("figs.jl")
include("properties.jl")
include("gadfly_panzoom.jl")
include("hittesting.jl")
#include("faces.jl")


text:
using JuliaFormatter
using JSON

function main()
    # Data read from `test/test.jl`
    text = JSON.json(Dict(
        "method" => "format",
        "params" => Dict("text" => String(read(abspath(joinpath(@__DIR__, "..", "tests", "test.jl")))))
    ))
    data = JSON.parse(String(text))
    text = data["params"]["text"]
    output = text
    indent = typemax(Int64)
    format_options = Dict{Symbol,Any}()
    for line in split(text, '\n')
        if length(line) > 0
            indent = min(length(line) - length(lstrip(line)), indent)
        end
    end
    output = format_text(text; format_options...)
    data["status"] = "success"
    data["params"]["text"] =
        [rstrip(lpad(l, length(l) + indent)) for l in split(output, "\n")]
    println(stdout, JSON.json(data))
end

main()


text:
# Autogenerated wrapper script for DecFP_jll for x86_64-apple-darwin
export libbid

JLLWrappers.@generate_wrapper_header("DecFP")
JLLWrappers.@declare_library_product(libbid, "@rpath/libbid.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libbid,
        "lib/libbid.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()