text:
export Horn, Detector, InstrumentDB, BandshapeInfo, SpectrumInfo, NoiseTemperatureInfo
export InstrumentDB, defaultdbfolder, parsefpdict, parsedetdict
export sensitivity_tant, t_to_trj, trj_to_t, deltat_to_deltatrj, deltatrj_to_deltat
export detector, bandpass, bandshape, spectrum, fknee_hz, tnoise

import YAML
import Stripeline
import Base: show
import RecipesBase

using Printf

@doc raw"""
Information about a STRIP horn

This structure holds a number of parameters relative to each feed horn in the
STRIP focal plane.

You should initialize `Horn` objects via the `InstrumentDB` constructor, which
loads their definition from a STRIP instrument database in YAML format.

Field              | Type           | Meaning
:----------------- |:-------------- |:------------------------------------------------------------
`name`             | String         | Name of the horn, e.g., `I0`
`id`               | Int            | Unique number of the horn, starting from 1
`polid`            | Int            | Unique ID of the polarimeter associated with the horn
`polarizerid`      | Int            | Unique ID of the polarizer+OMT associated with the horn
`moduleid`         | Int            | Number of the horn within the module, from 0 to 6
`color`            | String         | Name of the color associated with the module
`orientation`      | Array{Float64} | 3D vector containing the orientation of the horn in the sky
`fwhm_x_deg`       | Float64        | FWHM of the beam along the X axis, in degrees
`fwhm_y_deg`       | Float64        | FWHM of the beam along the Y axis, in degrees
`main_spillover`   | Float64        | Main reflector spillover
`sub_spillover`    | Float64        | Sub-reflector spillover
`xpd_db`           | Float64        | Cross-polarization, in dB
`directivity_dbi`  | Float64        | Directivity, in dBi
`ellipticity`      | Float64        | Ellipticity
"""
struct Horn
    name::String
    id::Int
    polid::Int
    polarizerid::Int
    moduleid::Int
    color::String
    orientation::Array{Float64,1}
    fwhm_x_deg::Float64
    fwhm_y_deg::Float64
    main_spillover::Float64
    sub_spillover::Float64
    xpd_db::Float64
    directivity_dbi::Float64
    ellipticity::Float64
end

function Base.show(io::IO, horn::Horn)
    if get(io, :compact, false)
        print(io, "Horn $(horn.name), module $(horn.moduleid) ($(horn.color))")
    else
        @printf(io, """
            Horn %s, module %d (%s):
                Orientation: [%.4f, %.4f, %.4f]
                FWHM (X/Y): %.4f × %.4f °
                Spillover: %f (main), %f (sub)
                Cross-polarization: %.2f dB
                Directivity: %.2f dBi
                Ellipticity: %.4f""",
            horn.name, horn.moduleid, horn.color,
            horn.orientation[1], horn.orientation[2], horn.orientation[3],
            horn.fwhm_x_deg, horn.fwhm_y_deg,
            horn.main_spillover, horn.sub_spillover,
            horn.xpd_db,
            horn.directivity_dbi,
            horn.ellipticity)
    end
end

@doc raw"""
    BandshapeInfo

Information about the spectral band response of a polarimeter.

Field                      | Type             | Meaning
:------------------------- |:---------------- |:-------------------------------------------------------
`center_frequency_hz`      | Float64          | Estimate for the center frequency, in Hz
`center_frequency_err_hz`  | Float64          | Estimated error on the center frequency, in Hz
`bandwidth_hz`             | Float64          | Estimated bandwidth, in Hz
`bandwidth_err_hz`         | Float64          | Estimated error on the bandwidth, in Hz
`lowest_frequency_hz`      | Float64          | Lowest frequency of the bandshape in `response`, in Hz
`highest_frequency_hz`     | Float64          | Highest frequency of the bandshape in `response`, in Hz
`num_of_frequencies`       | Int              | Number of samples in `response`
`bandshape`                | Array{Float64,1} | Profile of the bandshape (pure numbers)
`bandshape_error`          | Array{Float64,1} | Estimated error on the profile of the bandshape
`test_id`                  | Array{Int,1}     | ID of the unit-level test used to characterize the bandshape
`analysis_id`              | Int              | ID of the unit-level analysis used to characterize the bandshape

You can plot a `BandshapeInfo` object by importing `Plots` and using `plot`:

```julia
db = InstrumentDB()
plot(bandpass(db, "I0"), show_error = true)
```

The following keywords are recognized in the call to `plot`:

- `show_error` (default: `true`): include an error bar.
- `show_centerfreq` (default: `false`): include a vertical bar showing
  the position of the center frequency

"""
struct BandshapeInfo
    center_frequency_hz::Float64
    center_frequency_err_hz::Float64
    bandwidth_hz::Float64
    bandwidth_err_hz::Float64
    lowest_frequency_hz::Float64
    highest_frequency_hz::Float64
    num_of_frequencies::Int
    bandshape::Array{Float64,1}
    bandshape_error::Array{Float64,1}
    test_id::Array{Int,1}
    analysis_id::Int
end

function Base.show(io::IO, band::BandshapeInfo)
    if get(io, :compact, false)
        @printf(io, "BandshapeInfo(ν0=%.2f GHz, Δν=%.2f GHz)",
                band.center_frequency_hz * 1e-9,
                band.bandwidth_hz * 1e-9)
    else
        @printf(io, """
            Bandshape:
                Center frequency: %.2f ± %.2f GHz
                Bandwidth: %.2f ± %.2f GHz
                Frequency range: [%.2f, %.2f] GHz (%d points)
                Test ID: [%s]
                Analysis ID: %d""",
            band.center_frequency_hz * 1e-9, band.center_frequency_err_hz * 1e-9,
            band.bandwidth_hz * 1e-9, band.bandwidth_err_hz * 1e-9,
            band.lowest_frequency_hz * 1e-9, band.highest_frequency_hz * 1e-9, band.num_of_frequencies,
            join(["$x" for x in band.test_id], ", "),
            band.analysis_id)
    end
end

RecipesBase.@recipe function plot(band::BandshapeInfo; show_error = true, show_centerfreq = false)
    seriestype --> :path
    xguide --> "Frequency [GHz]"

    if show_centerfreq
        RecipesBase.@series begin
            seriestype --> :path
            color --> :gray
            label --> ""
            let centerfreq = band.center_frequency_hz * 1e-9
                [centerfreq, centerfreq], [0, 1]
            end
        end
    end

    if show_error
        ribbon --> (band.bandshape_error, band.bandshape_error)
        fillalpha --> 0.3
    end

    ν, profile = bandshape(band)
    ν .* 1e-9, profile
end

@doc raw"""
    BandshapeInfo()

Initialize a BandshapeInfo object with all values set to zero.
"""
BandshapeInfo() = BandshapeInfo(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0, Float64[], 0, 0)

function bandshape(bandinfo::BandshapeInfo)
    ν = range(
        bandinfo.lowest_frequency_hz,
        stop = bandinfo.highest_frequency_hz,
        length = bandinfo.num_of_frequencies,
    )

    @assert length(ν) == length(bandinfo.bandshape)
    (ν, bandinfo.bandshape, bandinfo.bandshape_error)
end

@doc raw"""
    bandshape(bandinfo::BandshapeInfo) -> Tuple{Array{Float64, 1}, Array{Float64, 1}}
    bandshape(db::InstrumentDB, polid::Integer) -> Tuple{Array{Float64, 1}, Array{Float64, 1}}
    bandshape(db::InstrumentDB, horn_name::AbstractString) -> Tuple{Array{Float64, 1}, Array{Float64, 1}}

Return a pair `(ν_hz, B, Berr)` containing the shape of the bandpass
in `bandinfo` (first form), or the bandpass taken from the instrument
database (second and third form). The two elements of the tuple
`(ν_hz, B)` are two arrays of the same length containing the
frequencies (in Hz) and the bandpass response at the same frequency
(pure number), and they are suitable to be plotted, like in the
following example:

```julia
db = InstrumentDB()
x, y, err = bandshape(db, "G2")
plot(x, y, ribbon=(err, err))   # Plot the bandpass and the error bar
```

However, it is easier just to use `plot` on a [`BandshapeInfo`](@ref)
object.

"""
bandshape

@doc raw"""
    SpectrumInfo

Information about the noise spectrum of the output of a polarimeter.

Field                | Type     | Meaning
:------------------- |:-------- |:-------------------------------------------------------
`slope_i`            | Float64  | The slope ($\alpha$) of the 1/f component of the noise in the I signal
`slope_i_err`        | Float64  | Error associated with the value of `slope_i`
`slope_q`            | Float64  | Same as `slope_i`, but for the Q signal
`slope_q_err`        | Float64  | Error associated with the value of `slope_q`
`slope_u`            | Float64  | Same as `slope_i`, but for the U signal
`slope_u_err`        | Float64  | Error associated with the value of `slope_u`
`fknee_i_hz`         | Float64  | Knee frequency of the I signal, in Hz
`fknee_i_err_hz`     | Float64  | Error associated with the value of `fknee_i_hz`
`fknee_q_hz`         | Float64  | Knee frequency of the Q signal, in Hz
`fknee_q_err_hz`     | Float64  | Error associated with the value of `fknee_q_hz`
`fknee_u_hz`         | Float64  | Knee frequency of the U signal, in Hz
`fknee_u_err_hz`     | Float64  | Error associated with the value of `fknee_u_hz`
`wn_i_k2_hz`         | Float64  | White noise level for the I signal, in K^2 Hz
`wn_i_err_k2_hz`     | Float64  | Error associated with the value of `wn_i_k2_hz`
`wn_q_k2_hz`         | Float64  | White noise level for the Q signal, in K^2 Hz
`wn_q_err_k2_hz`     | Float64  | Error associated with the value of `wn_q_k2_hz`
`wn_u_k2_hz`         | Float64  | White noise level for the U signal, in K^2 Hz
`wn_u_err_k2_hz`     | Float64  | Error associated with the value of `wn_u_k2_hz`
`load_temperature_k` | Float64  | System brightness temperature used during the tests (in K)
`test_id`            | Int      | ID of the unit-level test used to characterize the bandshape
`analysis_id`        | Int      | ID of the unit-level analysis used to characterize the bandshape

You can quickly plot the theoretical shape of the noise power spectrum
using `plot` on a `SpectrumInfo` object.

"""
struct SpectrumInfo
    slope_i::Float64
    slope_q::Float64
    slope_u::Float64
    slope_i_err::Float64
    slope_q_err::Float64
    slope_u_err::Float64
    fknee_i_hz::Float64
    fknee_q_hz::Float64
    fknee_u_hz::Float64
    fknee_i_err_hz::Float64
    fknee_q_err_hz::Float64
    fknee_u_err_hz::Float64
    wn_i_k2_hz::Float64
    wn_q_k2_hz::Float64
    wn_u_k2_hz::Float64
    wn_i_err_k2_hz::Float64
    wn_q_err_k2_hz::Float64
    wn_u_err_k2_hz::Float64
    load_temperature_k::Float64
    test_id::Int
    analysis_id::Int
end

function Base.show(io::IO, spec::SpectrumInfo)
    if get(io, :compact, false)
        @printf(io, "SpectrumInfo(α=[%.4f, %.4f, %.4f], fknee=[%.3f, %.3f, %.3f], wn=[%.3f, %.3f, %.3f])",
                spec.slope_i, spec.slope_q, spec.slope_u,
                spec.fknee_i_hz, spec.fknee_q_hz, spec.fknee_u_hz,
                spec.wn_i_k2_hz, spec.wn_q_k2_hz, spec.wn_u_k2_hz)
    else
        @printf(io, """
            Noise spectrum:
                Slope: I = %.4f ± %.4f, Q = %.4f ± %.4f, U = %.4f ± %.4f
                Knee frequency: I = %.1f ± %.1f mHz, Q = %.1f ± %.1f mHz, U = %.1f ± %.1f mHz
                White noise: Q = %.1f ± %.1f mK^2 Hz, U = %.1f ± %.1f mK^2 Hz
                System brightness temperature: %.1f K
                Test ID: %d
                Analysis ID: %d""",
            spec.slope_i, spec.slope_i_err,
            spec.slope_q, spec.slope_q_err,
            spec.slope_u, spec.slope_u_err,
            spec.fknee_i_hz * 1e3, spec.fknee_i_err_hz * 1e3,
            spec.fknee_q_hz * 1e3, spec.fknee_q_err_hz * 1e3,
            spec.fknee_u_hz * 1e3, spec.fknee_u_err_hz * 1e3,
            spec.wn_q_k2_hz * 1e6, spec.wn_q_err_k2_hz * 1e6,
            spec.wn_u_k2_hz * 1e6, spec.wn_u_err_k2_hz * 1e6,
            spec.load_temperature_k,
            spec.test_id,
            spec.analysis_id)
    end
end

RecipesBase.@recipe function plot(spec::SpectrumInfo)
    model = (freq, α, fknee, wn) -> wn * (1 + fknee/freq)^α
    fknee = Float64[]
    spec.fknee_q_hz > 0 && push!(fknee, spec.fknee_q_hz)
    spec.fknee_u_hz > 0 && push!(fknee, spec.fknee_u_hz)

    min_freq, max_freq = if length(fknee) ≥ 1
        1e-4 * minimum(fknee), 1e+4 * maximum(fknee)
    else
        0.01, 1.0
    end
    freqs = 10 .^ range(log10(min_freq), log10(max_freq), length = 30)

    scale --> :log10

    RecipesBase.@series begin
        label --> "Q"
        (freqs, model.(freqs, Ref(spec.slope_q), Ref(spec.fknee_q_hz), Ref(spec.wn_q_k2_hz)))
    end

    label --> "U"
    xlabel --> "Frequency [Hz]"
    ylabel --> "Power [K^2/Hz]"
    (freqs, model.(freqs, Ref(spec.slope_u), Ref(spec.fknee_u_hz), Ref(spec.wn_u_k2_hz)))
end

@doc raw"""
    SpectrumInfo()

Initialize a SpectrumInfo object with all values set to zero.
"""
SpectrumInfo() = SpectrumInfo(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0, 0)

@doc raw"""
    NoiseTemperatureInfo

Information about the noise temperature of a polarimeter. This structure is used
for the field `tnoise` of the `Detector` struct.

Field          | Type             | Meaning
:------------- |:---------------- |:------------------------------------------------------------------------------
`tnoise_k`     | Float64          | Noise temperature computed from `tnoise_values_k`, in K
`tnoise_err_k` | Float64          | Error associated with `tnoise_k`, computed from `tnoise_values_k`
`test_ids`     | Array{Int,1}     | List of unit-level test IDs used to estimate the noise temperature
`analysis_ids` | Array{Int,1}     | List of unit-level analysis report IDs used to estimate the noise temperature
`values_k`     | Array{Float64,1} | List of noise temperatures estimated from the tests
"""
struct NoiseTemperatureInfo
    tnoise_k::Float64
    tnoise_err_k::Float64
    test_ids::Array{Int,1}
    analysis_ids::Array{Int,1}
    values_k::Array{Float64,1}
end

function Base.show(io::IO, tnoise::NoiseTemperatureInfo)
    if get(io, :compact, false)
        @printf(io, "NoiseTemperatureInfo(tnoise=%.1f)",
                tnoise.tnoise_k)
    else
        @printf(io, """
            Noise temperature:
                Tnoise: %.1f ± %.1f K
                Estimates: [%s] K
                Test IDs: [%s]
                Analysis IDs: [%s]""",
            tnoise.tnoise_k, tnoise.tnoise_err_k,
            join([@sprintf("%.1f", x) for x in tnoise.values_k], ", "),
            join(["$x" for x in tnoise.test_ids], ", "),
            join(["$x" for x in tnoise.analysis_ids], ", ")
            )
    end
end

@doc raw"""
    NoiseTemperatureInfo()

Initialize a NoiseTemperatureInfo object with all values set to zero.
"""
NoiseTemperatureInfo() = NoiseTemperatureInfo(0.0, 0.0, Int[], Int[], Float64[])

@doc raw"""
Information about a STRIP detector

This structure holds information about a STRIP polarimeter.

You should initialize `Detector` objects via the `InstrumentDB` constructor,
which loads their definition from a local STRIP instrument database.

Field       | Type                 | Meaning
:---------- |:-------------------- |:-----------------------------------------------------------------
`id`        | Int                  | Integer ID of the polarimeter, e.g., `2` for `STRIP02`
`name`      | String               | Full name of the polarimeter, e.g., `STRIP02`
`band`      | String               | Band: it can either be `Q` or `W`
`bandshape` | BandshapeInfo        | Information about the bandpass response
`spectrum`  | SpectrumInfo         | Information about the noise spectrum (white noise and 1/f noise)
`tnoise`    | NoiseTemperatureInfo | Information about the noise temperature
"""
struct Detector
    id::Int
    name::String
    band::String

    bandshape::BandshapeInfo
    spectrum::SpectrumInfo
    tnoise::NoiseTemperatureInfo
end

function Base.show(io::IO, det::Detector)
    if get(io, :compat, false)
        print(io, "Detector($(det.name), $(det.band) band)")
    else
        @printf(io, """
            Detector %s (%s band):
                Center frequency: %.2f ± %.2f GHz
                Bandwidth: %.2f ± %.2f GHz
                Noise temperature: %.1f ± %.2f K
                Knee frequency: %.1f ± %.1f mHz (Q), %.1f ± %.1f mHz (U)""",
            det.name, det.band,
            det.bandshape.center_frequency_hz * 1e-9, det.bandshape.center_frequency_err_hz * 1e-9,
            det.bandshape.bandwidth_hz * 1e-9, det.bandshape.bandwidth_err_hz * 1e-9,
            det.tnoise.tnoise_k, det.tnoise.tnoise_err_k,
            det.spectrum.fknee_q_hz * 1e3, det.spectrum.fknee_q_err_hz * 1e3,
            det.spectrum.fknee_u_hz * 1e3, det.spectrum.fknee_u_err_hz * 1e3)
    end
end

@doc raw"""
STRIP instrument database

The "database" contains information about feed horns and polarimeters:

- The field `focalplane` is a dictionary (mapping) associating the string
  identifying a horn (e.g., `I0`) with a [`Horn`](@ref) structure;

- The field `detectors` is a dictionary associating the ID of the
  polarimeter (e.g., 2 stands for `STRIP02`) with a [`Detector`](@ref) structure.

You should usually create an object of this kind using the default constructor,
which parses a set of YAML files containing the real parameters of the
instrument.

# Examples

```jldoctest
julia> db = InstrumentDB();

julia> print("Number of horns in the database: $(length(keys(db.focalplane)))")
Number of horns in the database: 55

julia> print("Number of polarimeters in the database: $(length(keys(db.detectors)))")
Number of polarimeters in the database: 66
```

# Visualization

You can produce a table describing the contents of the instrument
database using `show` and passing `text/markdown` as MIME type:

````julia
db = InstrumentDB()
show(stdout, MIME("text/markdown"), db)
````

The table can be converted to other formats (HTML, LaTeX, Microsoft
Word, …) using commonly-available tools, e.g.,
[Pandoc](https://pandoc.org/).

"""
struct InstrumentDB
    focalplane::Dict{String,Horn}
    detectors::Dict{Int,Detector}
end

function Base.show(io::IO, db::InstrumentDB)
    @printf(io, "InstrumentDB(%d horns, %d detectors)",
            length(keys(db.focalplane)),
            length(keys(db.detectors)))
end

function Base.show(io::IO, ::MIME"text/markdown", db::InstrumentDB)
    print(io, """# Horns

| Horn | Band | Polarimeter | FWHM (x) [deg] | FWHM (y) [deg] |
|------|------|-------------|----------------|----------------|
""")

    # Present the horns in a well-defined order (alphabetic order
    # would put W-band before the Y module).
    for mod in ["I", "V", "B", "G", "Y", "O", "R", "W"]
        for horn in 0:6
            cur_key = "$mod$horn"
            cur_key in keys(db.focalplane) || continue
            cur_horn = db.focalplane[cur_key]
            @printf(
                "| %s | %s | %s | %.2f | %.2f |\n",
                cur_key,
                db.detectors[cur_horn.polid].name,
                db.detectors[cur_horn.polid].band,
                cur_horn.fwhm_x_deg,
                cur_horn.fwhm_y_deg,
            )
        end
    end

    print(io, """

# Polarimeters

| Name | Band | Noise temperature [K] | Center Frequency [GHz] | Bandpass [GHz] |
|------|------|-----------------------|------------------------|----------------|
""")

    det_keys = sort(keys(db.detectors) |> collect)
    for cur_key in det_keys
        cur_horn = db.detectors[cur_key]
        @printf(
            "| %s | %s | %.1f | %.2f | %.2f |\n",
            cur_horn.name,
            cur_horn.band,
            cur_horn.tnoise.tnoise_k,
            cur_horn.bandshape.center_frequency_hz * 1e-9,
            cur_horn.bandshape.bandwidth_hz * 1e-9,
        )
    end

end

@doc raw"""
    parsefpdict(fpdict)

Return a dictionary associating an horn name (e.g., `I0`) to a `Horn`
object containing information about some horn in the STRIP focal
plane. The information are parsed from `fpdict`, which should be a
dictionary loaded from a YAML file. The default YAML file to be used
is located in the folder returned by [`defaultdbfolder`](@ref) and is
usually named `strip_focal_plane.yaml`

"""
function parsefpdict(fpdict::Dict{Any,Any})
    focalplane = Dict{String,Horn}()
    for (key, value) in fpdict["horns"]
        refhorn = if haskey(fpdict["pairs"], key)
            refhornname = fpdict["pairs"][key]
            fpdict["horns"][refhornname]
        else
            value
        end

        focalplane[key] = Horn(key,
                               value["id"],
                               value["polarimeter_id"],
                               value["polarizer_id"],
                               value["module_id"],
                               value["color"],
                               value["orientation"],
                               refhorn["fwhm_x_deg"],
                               refhorn["fwhm_y_deg"],
                               refhorn["main_spillover"],
                               refhorn["sub_spillover"],
                               refhorn["xpd_db"],
                               refhorn["directivity_dbi"],
                               refhorn["ellipticity"])
    end

    focalplane
end

function parsebandshape(banddict::Dict{Any,Any})
    BandshapeInfo(get(banddict, "center_frequency_hz", 0.0),
        get(banddict, "center_frequency_err_hz", 0.0),
        get(banddict, "bandwidth_hz", 0.0),
        get(banddict, "bandwidth_err_hz", 0.0),
        get(banddict, "lowest_frequency_hz", 0.0),
        get(banddict, "highest_frequency_hz", 0.0),
        get(banddict, "num_of_frequencies", 0),
        get(banddict, "bandshape", Float64[]),
        get(banddict, "bandshape_error", Float64[]),
        get(banddict, "test_id", Int[]),
        get(banddict, "analysis_id", 0))
end

function parsespectrum(specdict::Dict{Any,Any})
    SpectrumInfo(get(specdict, "I_slope", 0.0),
        get(specdict, "Q_slope", 0.0),
        get(specdict, "U_slope", 0.0),
        get(specdict, "I_slope_err", 0.0),
        get(specdict, "Q_slope_err", 0.0),
        get(specdict, "U_slope_err", 0.0),
        get(specdict, "I_fknee_hz", 0.0),
        get(specdict, "Q_fknee_hz", 0.0),
        get(specdict, "U_fknee_hz", 0.0),
        get(specdict, "I_fknee_err_hz", 0.0),
        get(specdict, "Q_fknee_err_hz", 0.0),
        get(specdict, "U_fknee_err_hz", 0.0),
        get(specdict, "I_wn_level_k2_hz", 0.0),
        get(specdict, "Q_wn_level_k2_hz", 0.0),
        get(specdict, "U_wn_level_k2_hz", 0.0),
        get(specdict, "I_wn_level_err_k2_hz", 0.0),
        get(specdict, "Q_wn_level_err_k2_hz", 0.0),
        get(specdict, "U_wn_level_err_k2_hz", 0.0),
        get(specdict, "load_average_temperature_k", 20.0),
        get(specdict, "test_id", 0),
        get(specdict, "analysis_id", 0))
end

function parsetnoise(tnoisedict::Dict{Any,Any})
    NoiseTemperatureInfo(get(tnoisedict, "tnoise_k", 0.0),
        get(tnoisedict, "tnoise_err_k", 0.0),
        get(tnoisedict, "tnoise_test_ids", Int[]),
        get(tnoisedict, "analysis_ids", Int[]),
        get(tnoisedict, "values_k", Float64[]))
end

@doc raw"""
    parsedetdict(detdict)

Return a dictionary associating an integer number to a `Detector`
object containing information about the STRIP detector with the
corresponding number. The information are parsed from `detdict`, which
should be a dictionary loaded from a YAML file. The default YAML file
to be used is located in the folder returned by
[`defaultdbfolder`](@ref) and is usually named `strip_detectors.yaml`

"""
function parsedetdict(detdict)
    detectors = Dict{Int,Detector}()
    for curdet in detdict
        id = curdet["id"]

        bandpass = haskey(curdet, "bandpass") ? parsebandshape(curdet["bandpass"]) : BandshapeInfo()
        spectrum = haskey(curdet, "spectrum") ? parsespectrum(curdet["spectrum"]) : SpectrumInfo()
        tnoise = haskey(curdet, "tnoise") ? parsetnoise(curdet["tnoise"]) : NoiseTemperatureInfo()

        detectors[id] = Detector(id,
                                 get(curdet, "name", @sprintf("STRIP%02d", id)),
                                 get(curdet, "band", id ≤ 70 ? "Q" : "W"),
                                 bandpass,
                                 spectrum,
                                 tnoise)
    end

    detectors
end

@doc raw"""
    InstrumentDB(dbpath::AbstractString)

Load the STRIP instrument database from the specified path. Return an
instance of a InstrumentDB object.

"""
function InstrumentDB(dbpath::AbstractString)
    focalplanedict = open(joinpath(dbpath, "strip_focal_plane.yaml")) do f
        YAML.load(f)
    end

    detectordict = open(joinpath(dbpath, "strip_detectors.yaml")) do f
        YAML.load(f)
    end

    InstrumentDB(parsefpdict(focalplanedict), parsedetdict(detectordict))
end

@doc raw"""
    defaultdbfolder()

Return a string containing the (local) full path to the YAML files
containing the reference instrument DB.

"""
defaultdbfolder() = joinpath(dirname(pathof(Stripeline)), "..", "instrumentdb")

@doc raw"""
    InstrumentDB()

Load the STRIP instrument database from the directory returned by
[`defaultdbfolder`](@ref). Return an instance of a InstrumentDB object.
"""
InstrumentDB() = InstrumentDB(defaultdbfolder())

################################################################################

detector(db::InstrumentDB, polid::Integer) = db.detectors[polid]
detector(db::InstrumentDB, horn_name::AbstractString) = detector(db, db.focalplane[horn_name].polid)

@doc raw"""
    detector(db::InstrumentDB, polid::Integer) -> Detector
    detector(db::InstrumentDB, horn_name::AbstractString) -> Detector

Return a [`Detector`](@ref) structure, taken from the instrument
database. If the form with `polid` is used, `polid` is the progressive
number of the polarimeter; e.g., for `STRIP02`, `polid == 2`. In the
second form, you pass the string identifying the horn on the focal
plane, e.g., `I0`, `W3`, etc.

```julia
db = InstrumentDB()
pol1 = detector(db, 16)   # Get information about STRIP16
pol2 = detector(db, "V4") # Get information about the detector connected to horn V4
```

"""
detector

# These functions are already documented above
bandshape(db::InstrumentDB, polid::Integer) = bandshape(bandpass(db, polid))
bandshape(db::InstrumentDB, horn_name::AbstractString) = bandshape(bandpass(db, horn_name))

bandpass(db::InstrumentDB, polid::Integer) = detector(db, polid).bandshape
bandpass(db::InstrumentDB, horn_name::AbstractString) = bandpass(db, db.focalplane[horn_name].polid)

@doc raw"""
    bandpass(db::InstrumentDB, polid::Integer) -> BandshapeInfo
    bandpass(db::InstrumentDB, horn_name::AbstractString) -> BandshapeInfo

Return a pair `(ν_hz, B)` containing the bandpass `B` for the horn
with the specified ID (`polid`) or associated to some horn
(`horn_name`). To understand how `polid` and `horn_name` work, see the
documentation for [`detector`](@ref).

The two elements of the tuple `(ν_hz, B)` are two arrays of the same
length containing the frequencies (in Hz) and the bandpass response at
the same frequency (pure number).

```julia
db = InstrumentDB()
x, y = bandpass(db, "G2")
plot(x, y)   # Plot the bandpass
```

See also [`bandshape`](@ref).

"""
bandpass

spectrum(db::InstrumentDB, polid::Integer) = detector(db, polid).spectrum
spectrum(db::InstrumentDB, horn_name::AbstractString) = spectrum(db, db.focalplane[horn_name].polid)

@doc raw"""
    spectrum(db::InstrumentDB, polid::Integer) -> SpectrumInfo
    spectrum(db::InstrumentDB, horn_name::AbstractString) -> SpectrumInfo

Return a [`SpectrumInfo`](@ref) object, taken from the instrument
database. The meaning of the parameters `polid` and `horn_name` is
explained in the documentation for [`detector`](@ref).

"""
spectrum

tnoise(db::InstrumentDB, polid::Integer) = detector(db, polid).tnoise
tnoise(db::InstrumentDB, horn_name::AbstractString) = tnoise(db, db.focalplane[horn_name].polid)

@doc raw"""
    tnoise(db::InstrumentDB, polid::Integer) -> NoiseTemperatureInfo
    tnoise(db::InstrumentDB, horn_name::AbstractString) -> NoiseTemperatureInfo

Return a [`NoiseTemperatureInfo`](@ref) object, taken from the
instrument database. The meaning of the parameters `polid` and
`horn_name` is explained in the documentation for [`detector`](@ref).

"""
tnoise

function fknee_hz(db::InstrumentDB, polid::Integer; tsys_k = missing)
    specinfo = spectrum(db, polid)

    tsys_k === missing && return (specinfo.fknee_i_hz, specinfo.fknee_q_hz, specinfo.fknee_u_hz)

    # Correct the value of fknee depending on the ratio between the system
    # temperature now and the temperature used during the characterization of the
    # polarimeter

    (α_i, α_q, α_u) = (specinfo.slope_i, specinfo.slope_q, specinfo.slope_u)
    load_ratio = specinfo.load_temperature_k / tsys_k

    (specinfo.fknee_i_hz * load_ratio^(1 / α_i),
        specinfo.fknee_q_hz * load_ratio^(1 / α_q),
        specinfo.fknee_u_hz * load_ratio^(1 / α_u),)
end

fknee_hz(
    db::InstrumentDB,
    horn_name::AbstractString;
    tsys_k = missing,
) = fknee_hz(db, db.focalplane[horn_name].polid; tsys_k = tsys_k)

@doc raw"""
    fknee_hz(db::InstrumentDB, polid::Integer; tsys_k = missing) -> Tuple{Float64, Float64, Float64}
    fknee_hz(db::InstrumentDB, horn_name::AbstractString; tsys_k = missing) -> Tuple{Float64, Float64, Float64}

Return the knee frequency for the selected detector, taken from the
instrument database. The meaning of the parameters `polid` and
`horn_name` is explained in the documentation for [`detector`](@ref).

If `tsys_k` is specified, the system temperature is rescaled to the
desired temperature of the load feeding the polarimeter, so that the
1/f component of the noise remains unchanged but the white noise
plateau raises/lowers by an appropriate amount. Otherwise, the
function returns the raw frequency taken from the instrument database.

"""
fknee_hz

################################################################################

@doc raw"""
    sensitivity_tant(db::InstrumentDB, load_tant; modules = Set([0, 1, 2, 3, 4, 5, 6]))

Calculate the white-noise sensitivity of an array of detectors,
measured in K⋅√s, given some antenna temperature for the load. The
result takes in account only those horns belonging to the modules
listed in the keyword `modules` (the W-band horns belong to module
`-1`). By default, only the Q-band modules are considered.

The result assumes the radiometer equation: ``σ⋅√τ =
\frac{T_{sys}}{√2β}``, where `T_{sys}` is the system temperature, `β`
is the bandwidth, and `τ` is the acquisition time. The factor 2 comes
from the way Strip polarimeters operate. The system temperature is
assumed to be the noise temperature of each detector, plus the term
`load_tant`, which should take into account all the other sources of
power entering the system (e.g., telescope, atmosphere, etc.). The
term `load_tant` should be expressed as an antenna temperature.

"""
function sensitivity_tant(db::InstrumentDB, load_tant; modules = Set([0, 1, 2, 3, 4, 5, 6]))
    result = 0.0
    num = 0

    for (name, horn) in db.focalplane
        (horn.moduleid ∈ modules) || continue

        polarimeter = db.detectors[horn.polid]
        tsys = load_tant + polarimeter.tnoise.tnoise_k
        result += 2 * polarimeter.bandshape.bandwidth_hz / (tsys^2)
        num += 1
    end

    (1.0 / sqrt(result), num)
end

const kb = 1.3806503e-23  # Boltzmann constant
const hplanck = 6.62607015e-34 # Planck constant

@doc raw"""
    t_to_trj(temperature_k, nu_hz)

Convert a thermodynamic temperature (in K) into a Rayleigh-Jeans temperature, given some
specified frequency `nu_hz` (in Hz).

See also [`trj_to_t`](@ref) for the inverse transformation.
"""
t_to_trj(temperature_k, nu_hz) = hplanck * nu_hz / kb / (exp(hplanck * nu_hz / (kb * temperature_k)) - 1)

function bisect(fn, range)
    # Plain old bisection method, with a fixed precision

    min_x, max_x = range
    threshold = 1e-9

    while true
        mid_x = (min_x + max_x) / 2
        abs(max_x - min_x) < threshold && return mid_x

        min_y, mid_y, max_y = (fn(x) for x in (min_x, mid_x, max_x))
        if sign(min_y * mid_y) < 0
            max_x = mid_x
        else
            min_x = mid_x
        end
    end
end

@doc raw"""
    trj_to_t(temperature_k, nu_hz)

Convert a Rayleigh-Jeans temperature (in K) into a thermodynamic temperature, given some
specified frequency `nu_hz` (in Hz).

See also [`t_to_trj`](@ref) for the inverse transformation.

"""
function trj_to_t(temperature_k, nu_hz)
    bisect(x->t_to_trj(x, nu_hz) - temperature_k,
        (temperature_k / 10.0, temperature_k * 10.0),
    )
end

@doc raw"""
    deltat_to_deltatrj(temperature_k, deltat_k, nu_hz)

Convert a small temperature fluctuation `deltat_k` around temperature
`temperature_k` from thermodynamic temperature to Rayleigh-Jeans (RJ)
temperature. This function can be used to convert sensitivities expressed as
thermodynamic temperatures in RJ sensitivities.

See also `deltatrj_to_deltat` for the inverse function.
"""
function deltat_to_deltatrj(temperature_k, deltat_k, nu_hz)
    # The formula we are using here is ∂T_RJ/∂T × δT, but we are
    # using the analytical formula for the derivative

    exponential = exp(hplanck * nu_hz / (kb * temperature_k))
    (hplanck * nu_hz / (temperature_k * kb))^2 * exponential / (exponential - 1)^2 * deltat_k
end

@doc raw"""
    deltat_to_deltatrj(temperature_k, deltat_k, nu_hz)

Convert a small temperature fluctuation `deltat_k` around temperature
`temperature_k` from Rayleigh-Jeans (RJ) temperature to thermodynamic
temperature. This function can be used to convert sensitivities expressed as
RJ temperatures in thermodynamic sensitivities.

See also `deltatrj_to_deltat` for the inverse function.
"""
function deltatrj_to_deltat(temperature_k, deltat_k, nu_hz)
    tcmb_k = trj_to_t(temperature_k, nu_hz)

    exponential = exp(hplanck * nu_hz / (kb * tcmb_k))
    (tcmb_k * kb / (hplanck * nu_hz))^2 * (exponential - 1)^2 / exponential * deltat_k
end


text:
module HarmonicOrthogonalPolynomials
using FastTransforms, LinearAlgebra, ClassicalOrthogonalPolynomials, ContinuumArrays, DomainSets,
        BlockArrays, BlockBandedMatrices, InfiniteArrays, StaticArrays, QuasiArrays, Base, SpecialFunctions
import Base: OneTo, oneto, axes, getindex, convert, to_indices, _maybetail, tail, eltype, *, ==, ^, copy, -, abs, resize!
import BlockArrays: block, blockindex, unblock, BlockSlice
import DomainSets: indomain, Sphere
import LinearAlgebra: norm, factorize
import QuasiArrays: to_quasi_index, SubQuasiArray, *
import ContinuumArrays: TransformFactorization, @simplify, ProjectionFactorization
import ClassicalOrthogonalPolynomials: checkpoints, _sum, cardinality, increasingtruncations
import BlockBandedMatrices: BlockRange1, _BandedBlockBandedMatrix
import FastTransforms: Plan, interlace
import QuasiArrays: LazyQuasiMatrix, LazyQuasiArrayStyle
import InfiniteArrays: InfStepRange, RangeCumsum

export SphericalHarmonic, UnitSphere, SphericalCoordinate, RadialCoordinate, Block, associatedlegendre, RealSphericalHarmonic, sphericalharmonicy, Laplacian, AbsLaplacianPower, abs, -, ^, AngularMomentum
cardinality(::Sphere) = ℵ₁

include("multivariateops.jl")
include("spheretrav.jl")
include("coordinates.jl")

increasingtruncations(::BlockedUnitRange{<:RangeCumsum{Int,<:AbstractRange}}) = broadcast(n -> Block.(oneto(n)), 4:2:∞)

checkpoints(::UnitSphere{T}) where T = [SphericalCoordinate{T}(0.1,0.2), SphericalCoordinate{T}(0.3,0.4)]

abstract type AbstractSphericalHarmonic{T} <: MultivariateOrthogonalPolynomial{3,T} end
struct RealSphericalHarmonic{T} <: AbstractSphericalHarmonic{T} end
struct SphericalHarmonic{T} <: AbstractSphericalHarmonic{T} end
SphericalHarmonic() = SphericalHarmonic{ComplexF64}()
RealSphericalHarmonic() = RealSphericalHarmonic{Float64}()
copy(a::AbstractSphericalHarmonic) = a

axes(S::AbstractSphericalHarmonic{T}) where T = (Inclusion{SphericalCoordinate{real(T)}}(UnitSphere{real(T)}()), blockedrange(1:2:∞))

associatedlegendre(m) = ((-1)^m*prod(1:2:(2m-1)))*(UltrasphericalWeight((m+1)/2).*Ultraspherical(m+1/2))
lgamma(n) = logabsgamma(n)[1]

function sphericalharmonicy(ℓ, m, θ, φ)
    m̃ = abs(m)
    exp((lgamma(ℓ+m̃+1)+lgamma(ℓ-m̃+1)-2lgamma(ℓ+1))/2)*sqrt((2ℓ+1)/(4π)) * exp(im*m*φ) * sin(θ/2)^m̃ * cos(θ/2)^m̃ * jacobip(ℓ-m̃,m̃,m̃,cos(θ))
end

function getindex(S::SphericalHarmonic{T}, x::SphericalCoordinate, K::BlockIndex{1}) where T
    ℓ = Int(block(K))
    k = blockindex(K)
    m = k-ℓ
    convert(T, sphericalharmonicy(ℓ-1, m, x.θ, x.φ))::T
end

==(::SphericalHarmonic{T},::SphericalHarmonic{T}) where T = true
==(::RealSphericalHarmonic{T},::RealSphericalHarmonic{T}) where T = true

# function getindex(S::RealSphericalHarmonic{T}, x::ZSphericalCoordinate, K::BlockIndex{1}) where T
#     # sorts entries by ...-2,-1,0,1,2... scheme
#     ℓ = Int(block(K))
#     k = blockindex(K)
#     m = k-ℓ
#     m̃ = abs(m)
#     indepm = (-1)^m̃*exp((lgamma(ℓ-m̃)-lgamma(ℓ+m̃))/2)*sqrt((2ℓ-1)/(2π))*associatedlegendre(m̃)[x.z,ℓ-m̃]
#     m>0 && return cos(m*x.φ)*indepm
#     m==0 && return cos(m*x.φ)/sqrt(2)*indepm
#     m<0 && return sin(m̃*x.φ)*indepm
# end

function getindex(S::RealSphericalHarmonic{T}, x::SphericalCoordinate, K::BlockIndex{1}) where T
    # starts with m=0, then alternates between sin and cos terms (beginning with sin).
    ℓ = Int(block(K))
    m = blockindex(K)-1
    z = cos(x.θ)
    if iszero(m)
        return sqrt((2ℓ-1)/(4*π))*associatedlegendre(0)[z,ℓ]
    elseif isodd(m)
        m = (m+1)÷2
        return sin(m*x.φ)*(-1)^m*exp((lgamma(ℓ-m)-lgamma(ℓ+m))/2)*sqrt((2ℓ-1)/(2*π))*associatedlegendre(m)[z,ℓ-m]
    else
        m = m÷2
        return cos(m*x.φ)*(-1)^m*exp((lgamma(ℓ-m)-lgamma(ℓ+m))/2)*sqrt((2ℓ-1)/(2*π))*associatedlegendre(m)[z,ℓ-m]
    end
end

getindex(S::AbstractSphericalHarmonic, x::StaticVector{3}, K::BlockIndex{1}) = S[SphericalCoordinate(x), K]
getindex(S::AbstractSphericalHarmonic, x::StaticVector{3}, K::Block{1}) = S[x, axes(S,2)[K]]
getindex(S::AbstractSphericalHarmonic, x::StaticVector{3}, KR::BlockOneTo) = mortar([S[x, K] for K in KR])
getindex(S::AbstractSphericalHarmonic, x::StaticVector{3}, k::Int) = S[x, findblockindex(axes(S,2), k)]
getindex(S::AbstractSphericalHarmonic, x::StaticVector{3}, kr::AbstractUnitRange{Int}) = [S[x, k] for k in kr]

# @simplify *(Ac::QuasiAdjoint{<:Any,<:SphericalHarmonic}, B::SphericalHarmonic) =

##
# Expansion
##

const FiniteSphericalHarmonic{T} = SubQuasiArray{T,2,SphericalHarmonic{T},<:Tuple{<:Inclusion,<:BlockSlice{BlockRange1{OneTo{Int}}}}}
const FiniteRealSphericalHarmonic{T} = SubQuasiArray{T,2,RealSphericalHarmonic{T},<:Tuple{<:Inclusion,<:BlockSlice{BlockRange1{OneTo{Int}}}}}
copy(a::FiniteRealSphericalHarmonic) = a
copy(a::FiniteSphericalHarmonic) = a

function grid(S::FiniteSphericalHarmonic)
    T = real(eltype(S))
    N = blocksize(S,2)
    # The colatitudinal grid (mod $\pi$):
    θ = ((1:N) .- one(T)/2)/N
    # The longitudinal grid (mod $\pi$):
    M = 2*N-1
    φ = (0:M-1)*2/convert(T, M)
    SphericalCoordinate.(π*θ, π*φ')
end
function grid(S::FiniteRealSphericalHarmonic)
    T = real(eltype(S))
    N = blocksize(S,2)
    # The colatitudinal grid (mod $\pi$):
    θ = ((1:N) .- one(T)/2)/N
    # The longitudinal grid (mod $\pi$):
    M = 2*N-1
    φ = (0:M-1)*2/convert(T, M)
    SphericalCoordinate.(π*θ, π*φ')
end

struct SphericalHarmonicTransform{T} <: Plan{T}
    sph2fourier::FastTransforms.FTPlan{T,2,FastTransforms.SPINSPHERE}
    analysis::FastTransforms.FTPlan{T,2,FastTransforms.SPINSPHEREANALYSIS}
end
struct RealSphericalHarmonicTransform{T} <: Plan{T}
    sph2fourier::FastTransforms.FTPlan{T,2,FastTransforms.SPHERE}
    analysis::FastTransforms.FTPlan{T,2,FastTransforms.SPHEREANALYSIS}
end

SphericalHarmonicTransform{T}(N::Int) where T<:Complex = SphericalHarmonicTransform{T}(plan_spinsph2fourier(T, N, 0), plan_spinsph_analysis(T, N, 2N-1, 0))
RealSphericalHarmonicTransform{T}(N::Int) where T<:Real = RealSphericalHarmonicTransform{T}(plan_sph2fourier(T, N), plan_sph_analysis(T, N, 2N-1))

*(P::SphericalHarmonicTransform{T}, f::Matrix{T}) where T = SphereTrav(P.sph2fourier \ (P.analysis * f))
*(P::RealSphericalHarmonicTransform{T}, f::Matrix{T}) where T = RealSphereTrav(P.sph2fourier \ (P.analysis * f))

factorize(S::FiniteSphericalHarmonic{T}) where T =
    TransformFactorization(grid(S), SphericalHarmonicTransform{T}(blocksize(S,2)))
factorize(S::FiniteRealSphericalHarmonic{T}) where T =
    TransformFactorization(grid(S), RealSphericalHarmonicTransform{T}(blocksize(S,2)))

function _sum(A::AbstractSphericalHarmonic{T}, dims) where T
    @assert dims == 1
    PseudoBlockArray(Hcat(sqrt(4convert(T, π)), Zeros{T}(1,∞)), (Base.OneTo(1),axes(A,2)))
end

include("laplace.jl")
include("angularmomentum.jl")

end # module

ext: jl


text:
# converters from selected RSEXPREC to Hash
# They are used to translate SEXPREC attributes into Hash

function Base.convert(::Type{Hash}, pl::RPairList)
    res = Hash()
    for i in 1:length(pl.items)
        setindex!(res, pl.items[i], pl.tags[i])
    end
    res
end

##############################################################################
##
## Conversion of intermediate R objects into DataArray and DataFrame objects
##
##############################################################################

namask(rl::RLogicalVector) = BitArray(rl.data .== R_NA_INT32)
namask(ri::RIntegerVector) = BitArray(ri.data .== R_NA_INT32)
namask(rn::RNumericVector) = BitArray(map(isna_float64, reinterpret(UInt64, rn.data)))
# if re or im is NA, the whole complex number is NA
# FIXME avoid temporary Vector{Bool}
namask(rc::RComplexVector) = BitArray([isna_float64(v.re) || isna_float64(v.im) for v in reinterpret(Complex{UInt64}, rc.data)])
namask(rv::RNullableVector) = rv.na

DataArrays.data(rv::RVEC) = DataArray(rv.data, namask(rv))

function DataArrays.data(ri::RIntegerVector)
    if !isfactor(ri) return DataArray(ri.data, namask(ri)) end
    # convert factor into PooledDataArray
    pool = getattr(ri, "levels", emptystrvec)
    sz = length(pool)
    REFTYPE = sz <= typemax(UInt8)  ? UInt8 :
              sz <= typemax(UInt16) ? UInt16 :
              sz <= typemax(UInt32) ? UInt32 :
                                      UInt64
    dd = ri.data
    dd[namask(ri)] = 0
    refs = convert(Vector{REFTYPE}, dd)
    return PooledDataArray(DataArrays.RefArray(refs), pool)
end

function sexp2julia(rex::RSEXPREC)
    warn("Conversion of $(typeof(rex)) to Julia is not implemented")
    return nothing
end

function sexp2julia(rv::RVEC)
    # FIXME dimnames
    # FIXME forceDataArrays option to always convert to DataArray
    nas = namask(rv)
    hasna = any(nas)
    if hasnames(rv)
        # if data has no NA, convert to simple Vector
        return DictoVec(hasna ? DataArray(rv.data, nas) : rv.data, names(rv))
    else
        hasdims = hasdim(rv)
        if !hasdims && length(rv.data)==1
            # scalar
            # FIXME handle NAs
            # if hasna
            return rv.data[1]
        elseif !hasdims
            # vectors
            return hasna ? DataArray(rv.data, nas) : rv.data
        else
            # matrices and so on
            dims = tuple(convert(Vector{Int64}, getattr(rv, "dim"))...)
            return hasna ? DataArray(reshape(rv.data, dims), reshape(nas, dims)) :
                         reshape(rv.data, dims)
        end
    end
end

function sexp2julia(rl::RList)
    if isdataframe(rl)
        # FIXME remove Any type assertion workaround
        DataFrame(Any[data(col) for col in rl.data], map(identifier, names(rl)))
    elseif hasnames(rl)
        DictoVec(Any[sexp2julia(item) for item in rl.data], names(rl))
    else
        # FIXME return DictoVec if forceDictoVec is on
        map(sexp2julia, rl.data)
    end
end


text:
using QBase
using Documenter

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

makedocs(;
    modules=[QBase],
    authors="Brian Doolittle <brian.d.doolittle@gmail.com> and contributors",
    repo="https://github.com/ChitambarLab/QBase.jl/blob/{commit}{path}#L{line}",
    sitename="QBase.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://ChitambarLab.github.io/QBase.jl",
        assets=String["assets/custom.css"],
    ),
    pages=[
        "Home" => "index.md",
        "Tutorial" => "tutorial.md",
        "Bras and Kets" => "brakets.md",
        "Operators" => "operators.md",
        "States" => "states.md",
        "Evolution" => "evolution.md",
        "Measurements" => "measurements.md",
        "Probabilities" => "probabilities.md",
        "Information Theory" => "information_theory.md",
        "Math Utilities" => "math_utilities.md",
    ],
)

deploydocs(;
    repo="github.com/ChitambarLab/QBase.jl",
    devbranch = "main",
)


text:
module Categories

export Category, Object, Morphism, Composition, Ob, Hom, category, id, label, domain, codomain, compose

struct Object{O}
    object::O
end

obtype(::Object{O}) where {O} = O

abstract type AbstractMorphism{O,M} end

obtype(::AbstractMorphism{O}) where {O} = O

homtype(::AbstractMorphism{O,M}) where {O,M} = M

struct Id{O,M} <: AbstractMorphism{O,M} end
(id::Id{O})(x::O) where {O} = x

struct Morphism{O,M} <: AbstractMorphism{O,M}
    morphism::M
    domain::Object{O}
    codomain::Object{O}
end

morphism(m::Morphism) = m.morphism

domain(m::Morphism) = m.domain

codomain(m::Morphism) = m.codomain

struct Composition{O,M,N} <: AbstractMorphism{O,M}
    morphisms::NTuple{N,Morphism{O,M}}
end

domain(c::Composition) = first(c.morphisms).domain

codomain(c::Composition) = last(c.morphisms).codomain

struct Hom{O,M}
    domain::Object{O}
    codomain::Object{O}
end

domain(hom::Hom) = hom.domain

codomain(hom::Hom) = hom.codomain

struct Category{O,M}
    ob::Type{Object{O}}
    hom::Type{Morphism{O,M}}
    id::Id{O,M}

    Category(::Type{O},::Type{M}) where {O,M} = new{O,M}(Object{O},Morphism{O,M},Id{O,M}())
end

Ob(C::Category) = C.ob

Hom(C::Category)= C.hom

Hom(::Category{O,M},X::Object{O},Y::Object{O}) where {O,M} = Hom{O,M}(X,Y)

id(C::Category) = C.id

Base.in(::Object{O1}, ::Type{Object{O2}}) where {O1,O2} = O1 <: O2

Base.in(::AbstractMorphism{O1,M1}, ::Type{Morphism{O2,M2}}) where {O1,O2,M1,M2} = (O1 <: O2) && (M1 <: M2)

function Base.in(m::AbstractMorphism{O1,M1}, hom::Hom{O2,M2}) where {O1,O2,M1,M2}
    (O1 <: O2) &&
    (M1 <: M2) &&
    (domain(m) === domain(hom)) &&
    (codomain(m) === codomain(hom))
end

function compose(m1::Morphism{O,M},m2::Morphism{O,M}) where {O,M}
    m1.codomain === m2.domain || error("Domain mismatch")
    return Composition((m1,m2))
end

function compose(m::Morphism{O,M},c::Composition{O,M}) where {O,M}
    m.codomain === domain(c) || error("Domain mismatch")
    Composition((m,c.morphisms...))
end

function compose(c::Composition{O,M},m::Morphism{O,M}) where {O,M}
    codomain(c) === m.domain || error("Domain mismatch")
    Composition((c.morphisms...,m))
end

compose(m::AbstractMorphism{O,M},::Id{O,M}) where {O,M} = m

compose(::Id{O,M},m::AbstractMorphism{O,M}) where {O,M} = m

Base.:∘(m1::AbstractMorphism,m2::AbstractMorphism) = compose(m1,m2)

function (comp::Composition{O})(x::O) where {O}
    domain(comp).object === x || error("Domain mismatch")
    for mor in comp.morphisms
        x = mor(x)
    end
    return x
end

Base.show(io::IO, obj::Object) = print(io, obj.object)

Base.show(io::IO, mor::Morphism) = print(io, mor.morphism, ": ", mor.domain," → ", mor.codomain)

function Base.show(io::IO, comp::Composition)
    print(
        io,
        join(morphism.(comp.morphisms)," ∘ "),
        ": ",
        domain(first(comp.morphisms)),
        " → ",
        codomain(last(comp.morphisms))
    )
end

end # module


text:
using Base.Test
using ForwardNets


text:
export LinearQuadratic, LinearQuadratic2BVP
export LinearQuadraticQuasiMetricSpace, DoubleIntegrator
export steer, steer_pairwise

### QuasiMetric Typedef/Evaluation
immutable LinearQuadratic2BVP{T<:AbstractFloat,
                              FGinv<:Function,
                              FexpAt<:Function,
                              Fcdrift<:Function,
                              Fcost<:Function,
                              Fdcost<:Function,
                              Fddcost<:Function,
                              Fu<:Function,
                              Fx<:Function}
    A::Matrix{T}
    B::Matrix{T}
    c::Vector{T}
    R::Matrix{T}
    Ginv::FGinv
    expAt::FexpAt
    cdrift::Fcdrift
    cost::Fcost
    dcost::Fdcost
    ddcost::Fddcost
    u::Fu
    x::Fx
end
type LinearQuadratic{T<:AbstractFloat,L<:LinearQuadratic2BVP} <: QuasiMetric
    bvp::L
    cmax::T
end
LinearQuadratic{T}(A::Matrix{T}, B::Matrix{T}, c::Vector{T}, R::Matrix{T}, cmax::T = T(1)) =
    LinearQuadratic(LinearQuadratic2BVP(A, B, c, R), cmax)
setup_steering(d::LinearQuadratic, r) = (d.cmax = r)
steer(d::LinearQuadratic, v::AbstractVector, w::AbstractVector) = steer(d.bvp, v, w, d.cmax)

evaluate(d::LinearQuadratic, v::AbstractVector, w::AbstractVector) = steer(d, v, w)[1]
changeprecision{T<:AbstractFloat}(::Type{T}, d::LinearQuadratic) =
    LinearQuadratic(LinearQuadratic2BVP((map(T,x) for x in (d.bvp.A, d.bvp.B, d.bvp.c, d.bvp.R))...), T(d.cmax))

### QuasiMetric Space Instantiation
function LinearQuadraticQuasiMetricSpace(lo::AbstractVector, hi::AbstractVector,
                                         A::Matrix, B::Matrix, c::Vector, R::Matrix, C::Matrix)
    BoundedStateSpace(SVector(lo), SVector(hi), LinearQuadratic(A, B, c, R), OutputMatrix(C))
end
function DoubleIntegrator(d::Int, lo = zeros(d), hi = ones(d); vmax = 1.5, r = 1.)
    A = [zeros(d,d) eye(d); zeros(d,2d)]
    B = [zeros(d,d); eye(d)]
    c = zeros(2d)
    R = r*eye(d)
    C = [eye(d) zeros(d,d)]
    LinearQuadraticQuasiMetricSpace(SVector([lo; -vmax*ones(d)]), SVector([hi; vmax*ones(d)]), A, B, c, R, C)
end
function WebbJvdB13quad10d()
    g = 9.8
    m = .5
    # l =
    A32 = [0 g; -g 0; 0 0]
    A = [zeros(3,3) eye(3)      zeros(3,2)  zeros(3,2);
         zeros(3,3) zeros(3,3)  A32         zeros(3,2);
         zeros(2,3) zeros(2,3)  zeros(2,2)  eye(2);
         zeros(2,10)]
    B = [zeros(5, 3);
         ]
    c = zeros(10)
end

function helper_data_structures{S}(V::Vector{S}, M::LinearQuadratic, batchsize = 1001)
    N = length(V)
    DUpairs = [steer_pairwise(M.bvp, V, V[rng], M.cmax) for rng in [i:min(i+batchsize-1, N) for i in 1:batchsize:N]]
    Dmat = hcat([D for (D,U) in DUpairs]...)
    Umat = hcat([U for (D,U) in DUpairs]...)
    DSF = BruteDistanceDS(Dmat')
    DSB = BruteDistanceDS(Dmat)
    US = EmptyControlCache()   # TODO: non-empty ControlDS
    DSF, US, DSB, US
end

### Steering
propagate(d::LinearQuadratic, v::AbstractVector, u::DurationAndTargetControl) = typeof(v)(u.x)
propagate(d::LinearQuadratic, v::AbstractVector, u::DurationAndTargetControl, s::AbstractFloat) =
    s <= 0 ? v : s >= u.t ? typeof(v)(u.x) : typeof(v)(d.bvp.x(v, u.x, u.t, s))

steering_control(d::LinearQuadratic, v::AbstractVector, w::AbstractVector) = DurationAndTargetControl(steer(d,v,w)[2],w)
function collision_waypoints(d::LinearQuadratic, v::AbstractVector, w::AbstractVector)
    t = steer(d, v, w)[2]
    [d.bvp.x(v, w, t, s) for s in linspace(typeof(t)(0), t, 5)]
end

### LQ2BVP Nuts and Bolts
import SymPy

## Symbolic stuff
function expAt(A::Matrix, t::SymPy.Sym)
    n = size(A,1)
    maximum(abs(A^n)) > 0 && error("TODO: implement more cases than nilpotent A! (e.g. diagonalizable)")
    sum([A^i*(t^i/factorial(i)) for i in 0:n-1])
end
function Sym2Function(s::SymPy.Sym, args::Union{Symbol, Expr} = :t, replace_rules = ())
    sliteral = string(SymPy.simplify(SymPy.expand(s)))
    for rr in replace_rules
        sliteral = replace(sliteral, rr...)
    end
    @eval $args -> $(parse(sliteral))
end
function Sym2Function(s::Vector{SymPy.Sym}, args::Union{Symbol, Expr} = :t, replace_rules = (); sv = false)
    vliteral = (sv ? "SVector((" : "[") *
               join(map(x -> string(SymPy.simplify(SymPy.expand(x))), s), ", ") *
               (sv ? ",))" : "]")
    for rr in replace_rules
        vliteral = replace(vliteral, rr...)
    end
    @eval $args -> $(parse(vliteral))
end
function Sym2Function(s::Matrix{SymPy.Sym}, args::Union{Symbol, Expr} = :t, replace_rules = ())
    mliteral = "[" *
               join(mapslices(rw -> join(rw, " "), map(x -> string(SymPy.simplify(SymPy.expand(x))), s), 2), "; ") *
               "]"
    for rr in replace_rules
        mliteral = replace(mliteral, rr...)
    end
    @eval $args -> $(parse(mliteral))
end

## 2BVP Utility Functions
function LinearQuadratic2BVP{T}(A::Matrix{T}, B::Matrix{T}, c::Vector{T}, R::Matrix{T})
    n = size(A, 1)
    t = SymPy.symbols("t", real=true)
    s = SymPy.symbols("s", real=true)
    xS = [SymPy.symbols(join(["x$i" for i in 1:n], ", "), real = true)...]
    yS = [SymPy.symbols(join(["y$i" for i in 1:n], ", "), real = true)...]
    replace_rules = ((r"[xy]\d+", s -> s[1:1] * "[" * s[2:end] * "]"),
                     (r"\d+\.\d*", x -> "$T("*x*")"))

    expAtS = expAt(A, t)
    expAsS = expAt(A, s)
    GS = SymPy.integrate(expAtS*B*inv(R)*B'*expAtS', t)
    GinvS = inv(GS)
    cdriftS = SymPy.integrate(expAtS, t)*c
    xbarS = expAtS*xS + cdriftS
    costS = t + ((yS - xbarS)'*GinvS*(yS - xbarS))[1]
    dcostS = SymPy.diff(costS, t)
    ddcostS = SymPy.diff(costS, t, 2)
    uS = inv(R)*B'*expAt(A, t - s)'*GinvS*(yS-xbarS)
    xS = expAsS*xS + SymPy.integrate(expAsS, s)*c +
         SymPy.integrate(expAsS*B*inv(R)*B'*expAsS', s)*expAt(A, t - s)'*GinvS*(yS-xbarS)

    LinearQuadratic2BVP(A, B, c, R,
                        Sym2Function(GinvS, :(t::$T), replace_rules),
                        Sym2Function(expAtS, :(t::$T), replace_rules),
                        Sym2Function(cdriftS, :(t::$T), replace_rules),
                        Sym2Function(costS, :(x::AbstractVector, y::AbstractVector, t::$T), replace_rules),
                        Sym2Function(dcostS, :(x::AbstractVector, y::AbstractVector, t::$T), replace_rules),
                        Sym2Function(ddcostS, :(x::AbstractVector, y::AbstractVector, t::$T), replace_rules),
                        Sym2Function(uS, :(x::AbstractVector, y::AbstractVector, t::$T, s::$T), replace_rules, sv=true),
                        Sym2Function(xS, :(x::AbstractVector, y::AbstractVector, t::$T, s::$T), replace_rules, sv=true))
end

## Time-Optimal 2BVP
function topt_bisection{T}(dc, x0::AbstractVector, x1::AbstractVector, tm::T; tol = T(1e-3))
    # Bisection
    b = tm
    dc(x0, x1, b) < 0 && return tm
    a = tm / 100
    while dc(x0, x1, a) > 0; a /= 2; end
    m = T(0)
    cdval = T(1)
    while abs(cdval) > tol && abs(a - b) > tol
        m = (a+b)/2
        cdval = dc(x0, x1, m)
        cdval > 0 ? b = m : a = m
    end
    m
end
function topt_newton{T}(dc, ddc, x0::AbstractVector, x1::AbstractVector, tm::T; tol = T(1e-6))
    # Bisection / Newton's method combo
    b = tm
    dc(x0, x1, b) < 0 && return tm
    a = tm / 100
    while dc(x0, x1, a) > 0; a /= 2; end
    t = tm / 2
    cdval = dc(x0, x1, t)
    while abs(cdval) > tol && abs(a - b) > tol
        t = t - cdval / ddc(x0, x1, t)
        (t < a || t > b) && (t = (a+b)/2)
        cdval = dc(x0, x1, t)
        cdval > 0 ? b = t : a = t
    end
    t
end
function steer(L::LinearQuadratic2BVP, x0::AbstractVector, x1::AbstractVector, r::AbstractFloat)
    x0 == x1 && return zero(r), zero(r)
    t = isa(r, Float64) ? topt_newton(L.dcost, L.ddcost, x0, x1, r) : topt_bisection(L.dcost, x0, x1, r)
    L.cost(x0, x1, t), t
end
function steer_pairwise{T,S<:AbstractVector}(L::LinearQuadratic2BVP{T}, V::Vector{S}, W::Vector{S}, r::AbstractFloat)
    M = length(V)
    N = length(W)
    Vmat = statevec2mat(V)
    Wmat = statevec2mat(W)
    Vbarmat = L.expAt(r)*Vmat .+ L.cdrift(r)
    Ginv = L.Ginv(r)
    BRB = L.B*inv(L.R)*L.B'

    cd = pairwise(SqMahalanobis(Ginv*BRB*Ginv), Vbarmat, Wmat)
    LHT = Ginv*(L.A*Wmat .+ L.c)

    T1 = Distances.dot_percol(Wmat, LHT)
    BLAS.gemm!('T', 'N', T(-2), Vbarmat, LHT, T(1), cd)
    broadcast!(.+, cd, cd, 2*T1')
    broadcast!(.-, cd, 1, cd)

    cands = cd .> 0
    if is(V, W)
        for i in 1:N
            cands[i,i] = false
        end
    end
    IS, JS = findn(cands)
    VTS = [steer(L, V[i], W[j], r) for (i,j) in zip(IS, JS)]
    II = find(T[v for (v,s) in VTS] .<= r)
    DS = sparse(IS[II], JS[II], T[v for (v,s) in VTS[II]], M, N)
    US = sparse(IS[II], JS[II], [DurationAndTargetControl(VTS[i][2], W[JS[i]]) for i in II], M, N)
    DS, US
end


text:
include("test_minhash.jl")
include("test_simhash.jl")


text:
# Variant of vtk_grid with 2-D array "points".
# Size: (dim, num_points)
function vtk_grid{T<:AbstractFloat}(
        filename_noext::AbstractString,
        points::AbstractMatrix{T}, cells::Vector{MeshCell};
        compress::Bool=true, append::Bool=true)

    xvtk = XMLDocument()

    dim, Npts = size(points)

    # Reshape to 3D (if its not already)
    if dim == 1
        _points = zeros(3,Npts)
        _points[1,:] = points
    elseif dim == 2
        _points = zeros(3,Npts)
        _points[1:2,:] = points
    elseif dim == 3
        _points = points
    else
        throw(ArgumentError("points must be of size (dim,Npts) where dim = 1, 2 or 3 and Npts the number of points."))
    end

    Ncls = length(cells)

    vtk = DatasetFile(xvtk, filename_noext*".vtu", "UnstructuredGrid",
                      Npts, Ncls, compress, append)

    # VTKFile node
    xroot = vtk_xml_write_header(vtk)

    # UnstructuredGrid node
    xGrid = new_child(xroot, vtk.gridType_str)

    # Piece node
    xPiece = new_child(xGrid, "Piece")
    set_attribute(xPiece, "NumberOfPoints", vtk.Npts)
    set_attribute(xPiece, "NumberOfCells",  vtk.Ncls)

    # Points node
    xPoints = new_child(xPiece, "Points")

    # DataArray node
    data_to_xml(vtk, xPoints, _points, "Points", 3)

    # Cells node (below the Piece node)
    xCells = new_child(xPiece, "Cells")

    # Create data arrays.
    offsets = Array(Int32, Ncls)
    types = Array(UInt8, Ncls)

    Nconn = 0     # length of the connectivity array
    if Ncls >= 1  # it IS possible to have no cells
        offsets[1] = length(cells[1].connectivity)
    end

    for (n, c) in enumerate(cells)
        Npts_cell = length(c.connectivity)
        Nconn += Npts_cell
        types[n] = c.ctype.vtk_id
        if n >= 2
            offsets[n] = offsets[n-1] + Npts_cell
        end
    end

    # Create connectivity array.
    conn = Array(Int32, Nconn)
    const ONE = one(Int32)
    n = 1
    for c in cells
        for i in c.connectivity
            # We transform to zero-based indexing, required by VTK.
            conn[n] = i - ONE
            n += 1
        end
    end

    # Add arrays to the XML file (DataArray nodes).
    data_to_xml(vtk, xCells, conn,    "connectivity")
    data_to_xml(vtk, xCells, offsets, "offsets"     )
    data_to_xml(vtk, xCells, types,   "types"       )

    return vtk::DatasetFile
end

# Variant of vtk_grid with 1-D arrays x, y, z.
# Size of each array: (num_points)
function vtk_grid{T<:AbstractFloat}(
        filename_noext::AbstractString,
        x::AbstractVector{T}, y::AbstractVector{T}, z::AbstractVector{T},
        cells::Vector{MeshCell};
        compress::Bool=true, append::Bool=true)
    @assert length(x) == length(y) == length(z)
    Npts = length(x)
    points = Array(T, 3, Npts)
    for n = 1:Npts
        points[1, n] = x[n]
        points[2, n] = y[n]
        points[3, n] = z[n]
    end
    return vtk_grid(filename_noext, points, cells;
                    compress=compress, append=append)::DatasetFile
end

# 2D version
vtk_grid{T<:AbstractFloat}(
        filename_noext::AbstractString,
        x::AbstractVector{T}, y::AbstractVector{T},
        cells::Vector{MeshCell};
        compress::Bool=true, append::Bool=true) =
    vtk_grid(filename_noext, x, y, zero(x), cells, compress=compress, append=append)

# 1D version
vtk_grid{T<:AbstractFloat}(
        filename_noext::AbstractString,
        x::AbstractVector{T},
        cells::Vector{MeshCell};
        compress::Bool=true, append::Bool=true) =
    vtk_grid(filename_noext, x, zero(x), zero(x), cells, compress=compress, append=append)

# Variant with 4-D Array
function vtk_grid{T<:AbstractFloat}(
        filename_noext::AbstractString,
        points::AbstractArray{T,4}, cells::Vector{MeshCell};
        compress::Bool=true, append::Bool=true)

    dim, Ni, Nj, Nk = size(points)

    return vtk_grid(filename_noext, reshape(points,(dim,Ni*Nj*Nk)), cells, compress=compress, append=append)
end


text:
module MathLink

using Printf

export @W_str, @W_cmd, weval

include("../deps/deps.jl")
include("types.jl")
include("consts.jl")
include("init.jl")
include("link.jl")
include("wstp.jl")
include("extras.jl")
include("eval.jl")

end


text:
#-------------vortex.jl--------------------------------------------------------#
#
# Purpose: Find a field on the outside of a superfluid vortex
#
#------------------------------------------------------------------------------#

# Field type will be used for velocity field around vortex
type Field
   xDim::Int64
   yDim::Int64
   zDim::Int64

   dx::Float64
   dy::Float64
   dz::Float64

   # x, y, and z are all of size xDim*yDim*zDim
   x::Array{Float64, 1}
   y::Array{Float64, 1}
   z::Array{Float64, 1}

end

# position type will be used to find the location of vortices
type Position
    x::Float64
    y::Float64
    z::Float64
end

# function to initialize our field
function init_field(xDim::Int64, yDim::Int64, zDim::Int64,
                    dx::Float64, dy::Float64, dz::Float64)
    gsize = xDim * yDim * zDim
    empty_x = fill(0.0, gsize)
    empty_y = fill(0.0, gsize)
    empty_z = fill(0.0, gsize)
    return Field(xDim, yDim, zDim, dx, dy, dz,
                 empty_x, empty_y, empty_z)
end

# Function to find distance between two positions
function distance(loc1::Position, loc2::Position)
    x = loc1.x - loc2.x
    y = loc1.y - loc2.y
    z = loc1.z - loc2.z

    return sqrt(x*x + y*y + z*z)
end

# Function for the dot product between two positions
function dot(loc1::Position, loc2::Position)
    return loc1.x*loc2.x + loc1.y*loc2.y + loc1.z*loc2.z
end

# Function to find the magnitude of a position
function mag(loc::Position)
    return sqrt(loc.x*loc.x + loc.y*loc.y + loc.z*loc.z)
end

# Function to find the angle between two positions
function anglexy(loc1::Position, loc2::Position, radius)
    x = loc2.x - loc1.x
    y = loc2.y - loc1.y
    #return acos(dot(loc1,loc2)/(mag(loc1)*mag(loc2)))

    angle = 0
    if (sign(x) <= 0 && sign(y) > 0)
        angle = atan(-x/y) + 0.5*pi
    elseif (sign(x) < 0 && sign(y) <= 0)
        angle = atan(y/x) + pi
    elseif (sign(x) >= 0 && sign(y) < 0)
        angle = atan(-x/y) + 1.5 * pi
    else
        angle = atan(y/x)
    end
    return angle

end

# Function to find velocity around vortex in 2D
function find_field2d(location::Vector{Position}, xDim::Int64, yDim::Int64,
                      dx::Float64, dy::Float64)

    # initialize the field
    vel = init_field(xDim, yDim, 1, dx, dy, 0.0)

    # index is for internal keeping of the indices
    index = 0

    for i = 0:xDim - 1
        for j = 1:yDim
            index = j + i * yDim

            # Going through all vortex locations to find velocity field
            for k = 1:length(location)
                # First find the distance and angle to the vortex center
                point = Position(Float64((i)) * dx, Float64(j - 1)*dy, 0)
                radius = distance(point, location[k])
                theta = anglexy(point, location[k], radius)
                #println(theta)
                vel.x[index] = (1/(2pi * radius)) * cos(theta)
                vel.y[index] = (1/(2pi * radius)) * sin(theta)
            end
        end
    end

    return vel
end

# function to output data in a gnuplot fashion
function output2d(vel::Field)

    # opening files
    xfile = open("x.dat", "w")
    yfile = open("y.dat", "w")

    for i = 0:vel.xDim - 1
        for j = 1:vel.yDim
            index = j + i * vel.yDim
            x = Float64((i+1)) * vel.dx
            y = Float64(j) * vel.dy
            write(xfile, "$x\t$y\t$(vel.x[index])\n")
            write(yfile, "$x\t$y\t$(vel.y[index])\n")
        end
    end

    close(xfile)
    close(yfile)
end

'''
# main function
function main()
    xDim = 64
    yDim = 64
    zDim = 1

    dx = 1 / xDim
    dy = 1 / yDim
    dz = 0.0

    vortices = [Position(0.5,0.5,0)]

    vel = find_field2d(vortices, xDim, yDim, dx, dy)

    output2d(vel)
end

main()
'''

max_issues_repo_name: shiffman/simuleios

max_forks_repo_name: shiffman/simuleios


text:
function cleanup_queues()
    for q in sqs_list_queues(aws, "ocaws-jl-test-queue-")
        sqs_delete_queue(q)
    end
end

@testset "FIFO Queue" begin
    create_fifo_queue(; queue_name="ocaws-jl-test-queue-" * lowercase(Dates.format(now(Dates.UTC), "yyyymmddTHHMMSSZ")) * ".fifo") = sqs_create_queue(aws, queue_name; Dict(:FifoQueue => "true")...)

    @testset "Create FIFO Queue" begin
        queue_name = "ocaws-jl-test-queue-" * lowercase(Dates.format(now(Dates.UTC), "yyyymmddTHHMMSSZ")) * ".fifo"
        created_queue = create_fifo_queue(queue_name=queue_name)
        fetched_queue = sqs_get_queue(aws, queue_name)

        @test created_queue[:resource] == fetched_queue[:resource]
    end

    @testset "Send Message" begin
        queue = create_fifo_queue()

        options = Dict(
            :MessageGroupId => "1aaa11aa-11aa-11a1-aa1a-111aa1a1a111",
            :MessageDeduplicationId => "1aaa11aa-11aa-11a1-aa1a-111aa1a1a111"
        )
        sqs_send_message(queue, "Hello!", options...)
        message = sqs_receive_message(queue)

        @test message[:message] == "Hello!"
    end

    @testset "Queue Attributes" begin
        queue = create_fifo_queue()
        info = sqs_get_queue_attributes(queue)

        @test info["FifoQueue"] == "true"
    end
end

@testset "Queue" begin
    create_queue(; queue_name="ocaws-jl-test-queue-" * lowercase(Dates.format(now(Dates.UTC), "yyyymmddTHHMMSSZ"))) = sqs_create_queue(aws, queue_name)

    @testset "Create Queue" begin
        queue_name = "ocaws-jl-test-queue-" * lowercase(Dates.format(now(Dates.UTC), "yyyymmddTHHMMSSZ"))
        created_queue = create_queue(queue_name=queue_name)
        fetched_queue = sqs_get_queue(aws, queue_name)

        @test created_queue[:resource] == fetched_queue[:resource]
    end

    @testset "Send Message" begin
        queue = create_queue()
        sqs_send_message(queue, "Hello!")
        message = sqs_receive_message(queue)

        @test message[:message] == "Hello!"
    end

    @testset "Send Batch Messages" begin
        queue = create_queue()
        num_messages = 10

        msgs = repeat(["test message"], outer=num_messages)
        sqs_send_message_batch(queue, msgs)

        message_count = sqs_count(queue)
        @test message_count == num_messages

        while (m=sqs_receive_message(queue)) != nothing
            @test m[:message] == "test message"
        end
    end

    @testset "Queue Attributes" begin
        queue = create_queue()
        info = sqs_get_queue_attributes(queue)

        @test info["ApproximateNumberOfMessages"] == "0"
    end

    @testset "Delete Queue" begin
        queue_name = "ocaws-jl-test-queue-" * lowercase(Dates.format(now(Dates.UTC), "yyyymmddTHHMMSSZ"))
        queue = create_queue(queue_name=queue_name)
        sqs_delete_queue(queue)
        deleted_queue = sqs_get_queue(aws, queue_name)

        @test deleted_queue == nothing
    end

@testset "Message Visibility" begin
        queue = create_queue()
        message = "Hello!"

        @testset "Default" begin
            sqs_send_message(queue, message)
            response = sqs_receive_message(queue)

            @test response[:message] == message

            response = sqs_receive_message(queue)

            @test response == nothing
        end

        @testset "0 seconds" begin
            sqs_send_message(queue, message)
            response = sqs_receive_message(queue)

            @test response[:message] == message

            sqs_change_message_visibility(queue, response, 0)
            response = sqs_receive_message(queue)

            @test response[:message] == message
        end

        @testset "3 seconds" begin
            sqs_send_message(queue, message)
            response = sqs_receive_message(queue)

            @test response[:message] == message

            sqs_change_message_visibility(queue, response, 3)
            response = sqs_receive_message(queue)

            @test response == nothing

            sleep(3)
            response = sqs_receive_message(queue)

            @test response[:message] == message
        end
    end
end

cleanup_queues()


text:
function bar_reorder(args...; node::Int64=Int[], n::Int64=Int[],m::Int64=Int[],
    e1::Array{Int64, 1}=Int[], e2::Array{Int64, 1}=Int[],
    neworder::Array{Int64, 1}=Int[], iii::Int64=Int[],
    L::Array{Int64, 1}=Int[], pos::Array{Int64, 1}=Int[])

    for i  = 1:length(pos) #1
        for u = 1:pos[i]
            j = pos[i] - 1
            while j >= 0
                neworder[iii] = L[i+m*j] + 1
                iii=iii-1
                j = j - 1  # ????
            end
            k = e2[L[i+m*(u-1)]+1]
            if k <= n
                break
            else
                break
            end
        end
    end
    return neworder
end

#  function bar_reorder(args...; node::Int64=Int[], n::Int64=Int[],m::Int64=Int[],
#  	e1::Array{Int64, 1}=Int[], e2::Array{Int64, 1}=Int[],
#  	neworder::Array{Int64, 1}=Int[], iii::Int64=Int[],
#      L::Array{Int64, 1}=Int[], pos::Array{Int64, 1}=Int[], neworder2::Array{Int64, 1}=Int[],
#      j::Array{Int64, 1}=Int[],k::Array{Int64, 1}=Int[],i::Array{Int64, 1}=Int[])
#
#      i  = node - n #1
#      j = pos[i] - 1

#      while j >= 0
#      	neworder[iii] = L[i+m*j] + 1
#      	iii=iii-1
#      	j = j - 1  # ????
#      end
#      neworder
#
#      for j = 1:pos[i]
#      	k = e2[L[i+m*(j-1)]+1]
#      	if k > n
#      	   #bar_reorder(node=k,n=n,m=m,e1=e1,e2=e2,neworder=neworder,iii=iii,L=L,pos=pos,neworder2=neworder2)
#      	end
#      end
#      return neworder2
#  end


text:
using AlgebraPDF
using Test

@testset "Construction of pdf with parametric function" begin
    d0 = Normalized(FunctionWithParameters((x;p)->x^2+p.a; p=(a=2.0,)), (-2, 4))
    func(d0, 3) == 11
    pars(d0) == (a=2.0,)
    d1 = updatepars(d0, (a=3.0,))
    func(d1, 3) == 12
end

@testset "short cut: Normalized(function), fixedshapepdf" begin
    d = Normalized((e;p)->e^2+p.a; lims=(-1,2), p=(a=1.0,))
    @test typeof(pars(d)) <: NamedTuple
    @test length(freepars(d)) == 1
    @test length(fixedpars(d)) == 0
    #
    d1 = fixedshapepdf(x->exp(-(4x)^2), (-1, 2))
    @test length(freepars(d1)) == 0
end

@testset "fix parameters example" begin
    d0 = Normalized((e;p)->e^2+p.a; p=Ext(a=1.0), lims=(-1,2))
    @test d0(1.0) ≈ 2.0/(9/3+3*pars(d0).a)
    #
    d1 = fixpar(d0, :a, 2.9)
    @test d1(1.0) ≈ 3.9/(9/3+3*pars(d1).a)
    #
    d1′ = fixpars(d0, (:a,))
    @test freepars(d1′) == freepars(fixpars(d0, [:a]))
    @test freepars(d1′) == freepars(fixpars(d0, (a=1.0,)))
end

@testset "noparsnormf " begin
    d = Normalized((e;p)->e^2+p.a; p=Ext(a=1.0), lims=(-1,2))
    f = noparsnormf(d)
    d(1.1) == f(1.1)
end

@testset "integral " begin
    d1 = Normalized((e;p)->e^2+p.a; p=(a=1.0,), lims=(-1,2))
    @test integral(d1, (0,2)) < 1
    d2 = Normalized((e;p)->-e^2+p.b; p=(b=4.0,), lims=(-1,2))
    @test integral(d2, lims(d2)) ≈ 1
    s = d1+d2
    @test integral(s) == 2
    @test 1.8 < integral(s, (-0.9,1.9)) < 2
    s2 = updatepars(s,(α1=2,α2=4))
    #
    @test integral(s2[1]) == 2
    @test integral(s2[2]) == 4
    #
    @test 1.7 < integral(s2[1], (-0.9,1.9)) < 2
end

# implementation with NAMES of parameters build into the funciton call
struct nBW1{P,L} <: AbstractPDF{1}
    p::P
    lims::L
end
import AlgebraPDF:func
func(bw::nBW1, x::NumberOrTuple; p=pars(bw)) = abs2(p.m*p.Γ/(p.m^2-x^2-1im*p.m*p.Γ))

@testset "User-def type with fixed name" begin
    bw = nBW1((m=3.1,Γ=0.1), (0,4))
    @test pars(bw).m == 3.1
    @test pars(bw).Γ == 0.1
    @test bw(3.1) != 1.0 # because of normalization
    @test lims(bw) == (0,4)
end

# implementation with ORDER of parameters build into the funciton call
struct nBW2{P,L} <: AbstractPDF{1}
    p::P
    lims::L
end
import AlgebraPDF:func
function func(bw::nBW2, x::NumberOrTuple; p=pars(bw))
    m,Γ = (getproperty(p,s) for s in keys(bw.p))
    m*Γ/(m^2-x^2-1im*m*Γ)
end

@testset "User-def type with fixed order" begin
    bw_i = nBW2((m_i=3.1,Γ_i=0.1), (0.0,3.0))
    @test pars(bw_i).m_i == 3.1
    @test pars(bw_i).Γ_i == 0.1
    @test lims(bw_i) == (0.0,3.0)
    bw_j = nBW2((m_j=4.1,Γ_j=0.2), (1,4))
    @test pars(bw_j).m_j == 4.1
    @test pars(bw_j).Γ_j == 0.2
    @test lims(bw_j) == (1,4)
end

# same as the first option but with a macro
@makepdftype nBW3(y;p) = abs2(p.m*p.Γ/(p.m^2-y^2-1im*p.m*p.Γ))

@testset "User-def type with fixed name" begin
    bw = nBW3((m=3.1,Γ=0.1), (0,4))
    @test pars(bw).m == 3.1
    @test pars(bw).Γ == 0.1
    @test bw(3.1) != 1.0 # because of normalization
    @test lims(bw) == (0,4)
end


text:
#=  ADMM Algorithmic Regularization Paths for Sparse Statistical Machine Learning
#   Yue Hu, Eric C. Chi and Genevera I. Allen
#   arXiv:1504.06637v1
# 3.1 Sparse Linear Regression
  Goal:  fit for B in Y = X*B+eps
=#

function soft_threshold(x::Number, gamma::Float64)
  # Soft threshold of x to be writen to z
  gamma_by_2 = gamma/2
    if x > gamma_by_2
      return x - gamma_by_2
    elseif x < -gamma_by_2
      return x + gamma_by_2
    else
      return 0.0
    end
end

function soft_threshold!(x::Matrix, gamma::Float64, z::Matrix)
  @assert size(x) == size(z)
  for i in 1:size(x,1)
    for j in 1:size(x,2)
      z[i,j] = soft_threshold(x[i,j],gamma)
      #println("i=",i," j=",j," x=",x[i,j]," z=",z[i,j], " gamma=",gamma)
    end
  end
  return z
end

function fit_arp_slr_2d{T<:Number}(X::Array{T,2}, Y::Array{T,2}; stepsize::Float64 =1.05, gamma_init::Float64 =0.01, max_itterations::Integer=200, min_active::Integer=0, verbose::Bool=false)
  @assert size(X,1) == size(Y,1)
  n = size(X,1)
  p = size(X,2)
  q = size(Y,2)
  gamma_list = zeros(T,max_itterations)
  chisq_list = zeros(T,max_itterations)
  Z_list = zeros(max_itterations,p,q)
  B = zeros(p,q)  # Preallocate working matrices
  Z = zeros(p,q)
  U = zeros(p,q)
  H = inv(X'*X/n+diagm(ones(p)))
  HXtY = H*X'*Y
  local gamma = gamma_init
  last_itteration = max_itterations
  #println("# fit_arp_slr_2d")
  for k in 1:max_itterations
    gamma = gamma * stepsize
    B = HXtY+H*(Z-U)
    soft_threshold!(B+U,gamma,Z)
    gamma_list[k]= gamma
    Z_list[k,:,:] = Z
    U = U + B - Z
    num_active = length(find(Z))
    if verbose
      chisq = sum(abs2,Y-X*Z)
      print("# ", k, " ", gamma,": ")
      print("Z norms= ", num_active, " ", sum(abs,Z), " ", sum(abs2,Z), " chi^2= ", chisq, " ")
      #print("B=",B," Z=",Z," U=",U," ")
      println("")
    end
    if sum(abs2,Z) == 0. || num_active < min_active
      last_itteration = k
      break
    end
  end
  resize!(gamma_list,last_itteration)
  resize!(chisq_list,last_itteration)
  # resize to last_itteration?
  return (gamma_list,Z_list[1:last_itteration,:,:])
end

# Simple linear regression to update the non-zero coefficients
function refit_arp_slr_2d{T<:Number}(X::Array{T,2}, Y::Array{T,2}, B::Array{T,2})
  @assert size(X,1) == size(Y,1)
  p = size(X,2)
  Z = copy(B)
  for col in 1:size(B,2)         # Separate mutliple parameter regression problems into 1 at a time
    idx_active = find(B[:,col])  # find non-zero entries
    Xact = X[:,idx_active]
    Bact = (Xact'*Xact) \ (Xact'*Y[:,col])
    Z[idx_active,col] = vec(Bact)
  end
  pred = X*Z
  chisq = sum(abs2,Y-pred)
  return ( Z,chisq )
end

function refit_arp_slr_2d{T<:Number}(X::Array{T,2}, Y::Array{T,2}, B::Array{T,2}, Xcv::Array{T,2}, Ycv::Array{T,2} )
  @assert size(X,1) == size(Y,1)
  @assert size(X) == size(Xcv)
  @assert size(Y) == size(Ycv)
  (Z, chisq_refit) = refit_arp_slr_2d(X,Y,B)
  chisq_cv = sum(abs2,Ycv-Xcv*Z)
  return (Z,chisq_refit, chisq_cv )
end

#=
using Optim

function refit_arp_slr_2d_nonlin{T<:Number}(X::Array{T,2}, Y::Array{T,2}, B::Array{T,2})
  @assert size(X,1) == size(Y,1)
  p = size(X,2)

  # Function to minimize
  score(Z::Array{T,2}) = sum(abs2,Y-X*Z)
  idx_active = find(B)
  function score(param::Array{Float64,1})
     Z = copy(B)
     Z[idx_active] = param
     score(Z)
  end
  # This is much slower than linear regression.  Left code in case the problem becomes weakly non-linear in future
  opt_result = optimize(score,B[idx_active],LBFGS())
  Z = copy(B)
  Z[idx_active] = Optim.minimizer(opt_result)
  return ( Z,Optim.minimum(opt_result) )
end
=#

max_stars_repo_name: eford/RVSpectraKitLearn

max_issues_repo_name: eford/RVSpectraKitLearn

max_forks_repo_name: eford/RVSpectraKitLearn


text:
function FieldBracketingNthOrderBrentSolver(arg0::RealFieldElement, arg1::RealFieldElement, arg2::RealFieldElement, arg3::jint)
    return FieldBracketingNthOrderBrentSolver((RealFieldElement, RealFieldElement, RealFieldElement, jint), arg0, arg1, arg2, arg3)
end

function get_absolute_accuracy(obj::FieldBracketingNthOrderBrentSolver)
    return jcall(obj, "getAbsoluteAccuracy", RealFieldElement, ())
end

function get_evaluations(obj::FieldBracketingNthOrderBrentSolver)
    return jcall(obj, "getEvaluations", jint, ())
end

function get_function_value_accuracy(obj::FieldBracketingNthOrderBrentSolver)
    return jcall(obj, "getFunctionValueAccuracy", RealFieldElement, ())
end

function get_max_evaluations(obj::FieldBracketingNthOrderBrentSolver)
    return jcall(obj, "getMaxEvaluations", jint, ())
end

function get_maximal_order(obj::FieldBracketingNthOrderBrentSolver)
    return jcall(obj, "getMaximalOrder", jint, ())
end

function get_relative_accuracy(obj::FieldBracketingNthOrderBrentSolver)
    return jcall(obj, "getRelativeAccuracy", RealFieldElement, ())
end

function solve(obj::FieldBracketingNthOrderBrentSolver, arg0::jint, arg1::RealFieldUnivariateFunction, arg2::RealFieldElement, arg3::RealFieldElement, arg4::AllowedSolution)
    return jcall(obj, "solve", RealFieldElement, (jint, RealFieldUnivariateFunction, RealFieldElement, RealFieldElement, AllowedSolution), arg0, arg1, arg2, arg3, arg4)
end

function solve(obj::FieldBracketingNthOrderBrentSolver, arg0::jint, arg1::RealFieldUnivariateFunction, arg2::RealFieldElement, arg3::RealFieldElement, arg4::RealFieldElement, arg5::AllowedSolution)
    return jcall(obj, "solve", RealFieldElement, (jint, RealFieldUnivariateFunction, RealFieldElement, RealFieldElement, RealFieldElement, AllowedSolution), arg0, arg1, arg2, arg3, arg4, arg5)
end

function solve_interval(obj::FieldBracketingNthOrderBrentSolver, arg0::jint, arg1::RealFieldUnivariateFunction, arg2::RealFieldElement, arg3::RealFieldElement, arg4::RealFieldElement)
    return jcall(obj, "solveInterval", BracketedRealFieldUnivariateSolver_Interval, (jint, RealFieldUnivariateFunction, RealFieldElement, RealFieldElement, RealFieldElement), arg0, arg1, arg2, arg3, arg4)
end


text:
function NAND(A, B)
    nand = 0
    # Your code goes here!
    return nand
end

max_stars_repo_name: basimr/lovelace

max_issues_repo_name: basimr/project-lovelace

max_forks_repo_name: basimr/project-lovelace


text:
module TestPaths
  using Test
  using EcologicalNetworks

  # Generate some data

  A = UnipartiteProbabilisticNetwork([0.0 0.1; 0.2 0.8])

  P = number_of_paths(A, n=2)

  @test P[1,1] ≈ 0.02
  @test P[1,2] ≈ 0.08
  @test P[2,2] ≈ 0.66

  m = rand([0,1], (5,5))
  U = UnipartiteQuantitativeNetwork(m)
  @test number_of_paths(U) == number_of_paths(UnipartiteNetwork(m.>0))
  @test shortest_path(U) == shortest_path(UnipartiteNetwork(m.>0))

  N = unipartitemotifs()[:S1]
  b = bellman_ford(N)
  @test length(b) == 3
  @test length(filter(x -> x.to == :c, b)) == 2

  @test sort(dijkstra(N, :a)) == sort(bellman_ford(N, :a))
  @test sort(dijkstra(N)) == sort(bellman_ford(N))
  d = dijkstra(N, :a)
  @test length(d) == 2
  @test length(filter(x -> x.from == :a, d)) == 2

end


text:
struct MonteCarloProblem{T,T2,T3,T4,T5} <: AbstractMonteCarloProblem
  prob::T
  prob_func::T2
  output_func::T3
  reduction::T4
  u_init::T5
end

MonteCarloProblem(prob::DEProblem;
    output_func = (sol,i)-> (sol,false),
    prob_func= (prob,i,repeat)->prob,
    reduction = (u,data,I)->(append!(u,data),false),
    u_init = []) =
    MonteCarloProblem(prob,prob_func,output_func,reduction,u_init)

MonteCarloProblem(;prob,
    output_func = (sol,i)-> (sol,false),
    prob_func= (prob,i,repeat)->prob,
    reduction = (u,data,I)->(append!(u,data),false),
    u_init = [], p = nothing) =
    MonteCarloProblem(prob,prob_func,output_func,reduction,u_init)


text:
# types from SymEngine to Julia
## CSetBasic
mutable struct CSetBasic
    ptr::Ptr{Cvoid}
end

function CSetBasic()
    z = CSetBasic(ccall((:setbasic_new, libsymengine), Ptr{Cvoid}, ()))
    _finalizer(CSetBasic_free, z)
    z
end

function CSetBasic_free(x::CSetBasic)
    if x.ptr != C_NULL
        ccall((:setbasic_free, libsymengine), Nothing, (Ptr{Cvoid},), x.ptr)
        x.ptr = C_NULL
    end
end

function Base.length(s::CSetBasic)
    ccall((:setbasic_size, libsymengine), UInt, (Ptr{Cvoid},), s.ptr)
end

function Base.getindex(s::CSetBasic, n::UInt)
    result = Basic()
    ccall((:setbasic_get, libsymengine), Nothing, (Ptr{Cvoid}, UInt, Ref{Basic}), s.ptr, n, result)
    result
end

function Base.convert(::Type{Vector}, x::CSetBasic)
    n = Base.length(x)
    [x[i-1] for i in 1:n]
end
Base.convert(::Type{Set}, x::CSetBasic) = Set(convert(Vector, x))

## VecBasic Need this for get_args...

mutable struct CVecBasic
    ptr::Ptr{Cvoid}
end

function CVecBasic()
    z = CVecBasic(ccall((:vecbasic_new, libsymengine), Ptr{Cvoid}, ()))
    _finalizer(CVecBasic_free, z)
    z
end

function CVecBasic_free(x::CVecBasic)
    if x.ptr != C_NULL
        ccall((:vecbasic_free, libsymengine), Nothing, (Ptr{Cvoid},), x.ptr)
        x.ptr = C_NULL
    end
end

function Base.length(s::CVecBasic)
    ccall((:vecbasic_size, libsymengine), UInt, (Ptr{Cvoid},), s.ptr)
end

function Base.getindex(s::CVecBasic, n)
    result = Basic()
    ccall((:vecbasic_get, libsymengine), Nothing, (Ptr{Cvoid}, UInt, Ref{Basic}), s.ptr, UInt(n), result)
    result
end

function Base.convert(::Type{Vector}, x::CVecBasic)
    n = Base.length(x)
    [x[i-1] for i in 1:n]
end

start(s::CVecBasic) = 0
done(s::CVecBasic, i) = (i == Base.length(s))
next(s::CVecBasic, i) = s[i], i+1
if VERSION < v"0.7.0-rc1"
    Base.start(s::CVecBasic) = start(s)
    Base.done(s::CVecBasic, i) = done(s, i)
    Base.next(s::CVecBasic, i) = next(s, i)
else
    function Base.iterate(s::CVecBasic, i=start(s))
        done(s, i) && return nothing
        next(s, i)
    end
end

## CMapBasicBasic
mutable struct CMapBasicBasic
    ptr::Ptr{Cvoid}
end

function CMapBasicBasic()
    z = CMapBasicBasic(ccall((:mapbasicbasic_new, libsymengine), Ptr{Cvoid}, ()))
    _finalizer(CMapBasicBasic_free, z)
    z
end

function CMapBasicBasic(dict::Dict)
    c = CMapBasicBasic()
    for (key, value) in dict
        c[Basic(key)] = Basic(value)
    end
    return c
end

function CMapBasicBasic_free(x::CMapBasicBasic)
    if x.ptr != C_NULL
        ccall((:mapbasicbasic_free, libsymengine), Nothing, (Ptr{Cvoid},), x.ptr)
        x.ptr = C_NULL
    end
end

function Base.length(s::CMapBasicBasic)
    ccall((:mapbasicbasic_size, libsymengine), UInt, (Ptr{Cvoid},), s.ptr)
end

function Base.getindex(s::CMapBasicBasic, k::Basic)
    result = Basic()
    ret = ccall((:mapbasicbasic_get, libsymengine), Cint, (Ptr{Cvoid}, Ref{Basic}, Ref{Basic}), s.ptr, k, result)
    if ret == 0
        throw(KeyError("Key not found"))
    end
    result
end

function Base.setindex!(s::CMapBasicBasic, k::Basic, v::Basic)
    ccall((:mapbasicbasic_insert, libsymengine), Nothing, (Ptr{Cvoid}, Ref{Basic}, Ref{Basic}), s.ptr, k, v)
end

Base.convert(::Type{CMapBasicBasic}, x::Dict{Any, Any}) = CMapBasicBasic(x)

## Dense matrix

mutable struct CDenseMatrix <: DenseArray{Basic, 2}
    ptr::Ptr{Cvoid}
end

Base.promote_rule(::Type{CDenseMatrix}, ::Type{Matrix{T}} ) where {T <: Basic} = CDenseMatrix

function CDenseMatrix_free(x::CDenseMatrix)
    if x.ptr != C_NULL
        ccall((:dense_matrix_free, libsymengine), Nothing, (Ptr{Cvoid},), x.ptr)
        x.ptr = C_NULL
    end
end

function CDenseMatrix()
    z = CDenseMatrix(ccall((:dense_matrix_new, libsymengine), Ptr{Cvoid}, ()))
    _finalizer(CDenseMatrix_free, z)
    z
end

function CDenseMatrix(m::Int, n::Int)
    z = CDenseMatrix(ccall((:dense_matrix_new_rows_cols, libsymengine), Ptr{Cvoid}, (Int, Int), m, n))
    _finalizer(CDenseMatrix_free, z)
    z
end

function CDenseMatrix(x::Array{T, 2}) where T
    r,c = size(x)
    M = CDenseMatrix(r, c)
    for j in 1:c
        for i in 1:r ## fill column by column
            M[i,j] = x[i,j]
        end
    end
    M
end

function Base.convert(::Type{Matrix}, x::CDenseMatrix)
    m,n = Base.size(x)
    [x[i,j] for i in 1:m, j in 1:n]
end

Base.convert(::Type{CDenseMatrix}, x::Array{T, 2}) where {T} = CDenseMatrix(x)
Base.convert(::Type{CDenseMatrix}, x::Array{T, 1}) where {T} = convert(CDenseMatrix, reshape(x, length(x), 1))

function toString(b::CDenseMatrix)
    a = ccall((:dense_matrix_str, libsymengine), Cstring, (Ptr{Cvoid}, ), b.ptr)
    string = unsafe_string(a)
    ccall((:basic_str_free, libsymengine), Nothing, (Cstring, ), a)
    string = replace(string, "**", "^") # de pythonify
    return string
end

function Base.show(io::IO, m::CDenseMatrix)
    r, c = size(m)
    println(io, "CDenseMatrix: $r x $c")
    println(io, toString(m))
end


text:
# Autogenerated wrapper script for PETSc_jll for x86_64-unknown-freebsd-libgfortran3
export libpetsc, libpetsc_Float32_Complex_Int32, libpetsc_Float32_Complex_Int64, libpetsc_Float32_Real_Int32, libpetsc_Float32_Real_Int64, libpetsc_Float64_Complex_Int32, libpetsc_Float64_Complex_Int64, libpetsc_Float64_Real_Int32, libpetsc_Float64_Real_Int64

using OpenBLAS32_jll
using MPICH_jll
using MicrosoftMPI_jll
using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("PETSc")
JLLWrappers.@declare_library_product(libpetsc, "libpetsc.so.3.15")
JLLWrappers.@declare_library_product(libpetsc_Float32_Complex_Int32, "libpetsc.so.3.15")
JLLWrappers.@declare_library_product(libpetsc_Float32_Complex_Int64, "libpetsc.so.3.15")
JLLWrappers.@declare_library_product(libpetsc_Float32_Real_Int32, "libpetsc.so.3.15")
JLLWrappers.@declare_library_product(libpetsc_Float32_Real_Int64, "libpetsc.so.3.15")
JLLWrappers.@declare_library_product(libpetsc_Float64_Complex_Int32, "libpetsc.so.3.15")
JLLWrappers.@declare_library_product(libpetsc_Float64_Complex_Int64, "libpetsc.so.3.15")
JLLWrappers.@declare_library_product(libpetsc_Float64_Real_Int32, "libpetsc.so.3.15")
JLLWrappers.@declare_library_product(libpetsc_Float64_Real_Int64, "libpetsc.so.3.15")
function __init__()
    JLLWrappers.@generate_init_header(OpenBLAS32_jll, MPICH_jll, MicrosoftMPI_jll, CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libpetsc,
        "lib/libpetsc_double_real_Int32.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libpetsc_Float32_Complex_Int32,
        "lib/libpetsc_single_complex_Int32.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libpetsc_Float32_Complex_Int64,
        "lib/libpetsc_single_complex_Int64.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libpetsc_Float32_Real_Int32,
        "lib/libpetsc_single_real_Int32.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libpetsc_Float32_Real_Int64,
        "lib/libpetsc_single_real_Int64.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libpetsc_Float64_Complex_Int32,
        "lib/libpetsc_double_complex_Int32.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libpetsc_Float64_Complex_Int64,
        "lib/libpetsc_double_complex_Int64.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libpetsc_Float64_Real_Int32,
        "lib/libpetsc_double_real_Int32.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libpetsc_Float64_Real_Int64,
        "lib/libpetsc_double_real_Int64.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
function test_mapped_dicts()
    # b = Fourier(20)
    # x1 = 1.0+im
    # x2 = 2.0+3im
    # m1 = interval_map(0.0, 1.0, x1, x2)
    # m2 = embedding_map(Float64, Complex{Float64})
    # m = m1 ∘ m2
    # c = mapped_dict(b, m)
    #
    # @test !in_support(c, 10, 0.5)
    # @test in_support(c, 10, x1)
    # @test in_support(c, 10, (x1+x2)/2)
    #
    # u = approximate(c, exp)
    # z = 1/3*x1+2/3*x2
    # @test abs(u(z)-exp(z)) < 0.2

    d = UnitCircle()
    m = parameterization(d)
    b = Fourier(20)
    c = mapped_dict(b, m)
    f1(x,y) = exp(x+im*y)
    u = approximate(c, f1)
    @test in_support(c, 1, SVector(1.0, 0.0))
    @test abs(u(SVector(cos(10),sin(10)))-f1(cos(10),sin(10))) < 1e-6
end


text:
import Statistics: quantile, mean
using GR

function draw_empty(xmin, xmax, ymin, ymax; xmargin=0.0, ymargin=0.05, fontscale=5)
    xdiff = (xmax-xmin)*xmargin; ydiff = (ymax-ymin)*ymargin
    # Empty plot with right axes:
    setcharheight(fontscale*0.027)
    plot(; xlim=[xmin-xdiff,xmax+xdiff],
           ylim=[ymin-ydiff,ymax+ydiff], ylabel="Cost")
    nothing
end

function setcolors_grayscale(Lx)
    # Color 1 = black
    setcolorrep(1, 0.0, 0.0, 0.0)
    # From 2=white to darker, in cycles
    for k in 1:Lx
        gr = 1.0 - rem(2*(k-1)/Lx, 1.0)*0.8
        setcolorrep(k+1, gr, gr, gr)
    end
end

function boxbars(x, dx, Q_)
    #verrorbars(x, Q_[:,3], Q_[:,1], Q_[:,5])
    for k in 1:length(x)
        x_ = x[k]
        setlinecolorind(1); setmarkercolorind(1)
        setfillcolorind(k+1)
        polyline([x_,x_], [Q_[k,1], Q_[k,5]])
        fillrect(x_-dx,x_+dx,Q_[k,2],Q_[k,4])
        drawrect(x_-dx,x_+dx,Q_[k,2],Q_[k,4])
        polyline([x_-dx,x_+dx], Q_[k,3]*ones(2))
        polyline([x_-dx,x_+dx], Q_[k,1]*ones(2))
        polyline([x_-dx,x_+dx], Q_[k,5]*ones(2))
    end
end

function plot_all_estimates(xs, E; group_width=0.4, box_width=0.8,
    max_y=Inf)
    n_xs = length(xs)
    n_ns = size(E)[2]
    dx = (xs[2]-xs[1])*group_width/n_xs
    setfillintstyle(GR.INTSTYLE_SOLID)
    setmarkertype(GR.MARKERTYPE_SOLID_CIRCLE)
    setmarkersize(1.0)
    xmin = xs[1]-dx*n_ns/2; xmax = xs[end]+dx*(n_ns+2)/2
    ymin = Inf; ymax = -Inf
    Q_ = Vector(undef,n_xs)
    M_ = Vector(undef, n_xs)
    for k in 1:n_xs
        Q_[k] = mapslices(x -> quantile(x, [0.05,0.25,0.5,0.75,0.95]), E[k,:,:], dims=2)
        ymin = min(ymin, minimum(Q_[k]))
        ymax = max(ymax, maximum(Q_[k]))
    end
    ymax = min(max_y, ymax)
    draw_empty(xmin, xmax, ymin, ymax)
    setcolors_grayscale(n_ns)
    for k in 1:n_xs
        setlinecolorind(1); setmarkercolorind(1)
        #setfillcolorind(rem(n_xs-k,3)+2)
        xs_ = xs[k] .+ (collect(1:n_ns) .- n_ns/2)*dx
        boxbars(xs_, dx/2*box_width, Q_[k])
        # Add mean:
        M_ = vec(mapslices(x->mean(x), E[k,:,:], dims=2))
        setmarkertype(GR.MARKERTYPE_DIAGONAL_CROSS)
        polymarker(xs_, M_)
    end
end

using JLD2
@load "paper/output/lgDemoCompareOptimalCost_results.jld2"

figure(figsize=(15,5))
beginprint("paper/output/lgDemoCompareOptimalCost.pdf")
plot_all_estimates(ts, [BS_cost AS_cost]; max_y=12e3)
endprint()


text:
#NOTE: distribute does not work yet on this version

using MPIClusterManagers
using Distributed
using Plasmo

graph = OptiGraph()

#Add nodes to a GraphModel
@optinode(graph,n1)
@optinode(graph,n2)
@optinode(graph,n3)
@optinode(graph,n4)

@variable(n1,0 <= x <= 2)
@variable(n1,0 <= y <= 3)
@variable(n1, z >= 0)
@constraint(n1,x+y+z >= 4)
@objective(n1,Min,y)

@variable(n2,x >= 0)
@NLconstraint(n2,ref,exp(x) >= 2)
@variable(n2,z >= 0)
@constraint(n2,z + x >= 4)
@objective(n2,Min,x)

@variable(n3,x[1:5] >= 0)
@NLconstraint(n3,ref,exp(x[3]) >= 5)
@constraint(n3,sum(x[i] for i = 1:5) == 10)
@objective(n3,Min,x[1] + x[2] + x[3])

@variable(n4,x[1:5])
@constraint(n4,sum(x[i] for i = 1:5) >= 10)
@NLconstraint(n4,ref,exp(x[2]) >= 4)
@objective(n4,Min,x[2])

#Link constraints take the same expressions as the JuMP @constraint macro
@linkconstraint(graph,n1[:x] == n2[:x])
@linkconstraint(graph,n2[:x] == n3[:x][3])
@linkconstraint(graph,n3[:x][1] == n4[:x][1])

manager=MPIManager(np=2)
addprocs(manager)

@everywhere begin
    using Pkg
    using Revise
    Pkg.activate((@__DIR__)*"/..")
    using Plasmo
    using PipsSolver
end

#Distribute the modelgraph to the Julia workers
julia_workers = collect(values(manager.mpi2j))
remote_graphs = distribute(graph,julia_workers,remote_name = :pipsgraph)  #create the variable pipsgraph on each worker

r1 = remote_graphs[1] #reference to the model graph on worker 1
r2 = remote_graphs[2] #reference to the model graph on worker 2

#Solve with PIPS-NLP
@mpi_do manager begin
    using MPI
    pipsnlp_solve(pipsgraph)  #this works because the vairable :pipsgraph was defined on each worker (mpirank)
end

println("Fetching solution")
rank_zero = manager.mpi2j[0] #get the julia process representing rank 0
solution = fetch(@spawnat(rank_zero, getfield(Main, :pipsgraph))) #transfer the solution from rank 0 to our local Julia process

for node in getnodes(solution)
    println(node.model.ext[:colVal])
end


text:
using Documenter, ModelConstructors

makedocs(modules = [ModelConstructors],
         clean = false,
         format = Documenter.HTML(),
         sitename = "ModelConstructors.jl",
         authors = "FRBNY-DSGE",
         linkcheck = false,
         strict = false,
         pages = Any[
                     "Home"                                 => "index.md",
                     "Model Design"                         => "model_design.md",
                     "Creating Models"                      => "example_model.md",
                     "Implementation Details"               =>  "implementation_details.md",
                     "Contributing to ModelConstructors.jl" => "contributing.md",
                     "License"                              => "license.md"
         ],
         doctest = false # for now
)

deploydocs(
    repo = "github.com/FRBNY-DSGE/ModelConstructors.jl.git",
    target = "build",
    deps = nothing,
    devbranch = "master",
    branch = "gh-pages",
    # versions = "v#",
    # julia = "0.7",
    # osname = "osx",
    make = nothing
)


text:
using Distributed

using CSV
using DataFrames
@everywhere using FileIO
using JLD2
using ProgressMeter
using Statistics

@assert length(ARGS) >= 1 "no experiment path given"

prefix_path = ARGS[1]
num_trajectories = length(ARGS)==2 ? parse(Int, ARGS[2]) : 1000

df = CSV.read(joinpath(prefix_path, "parameters_map.csv")) |> DataFrame
output_dirs = map( x-> joinpath(prefix_path, x, "output"), df.path)

results = @showprogress pmap(output_dirs) do path
  load(joinpath(path, "summary.jld2"), "last_infections", "num_infections")
end

last_infections = getindex.(results, 1)
num_infections = getindex.(results, 2)

df.mean_last_infection_time = last_infections .|> mean
df.mean_num_infected = num_infections .|> mean

save(joinpath(prefix_path,"summary.jld2"),
  "summary", df,
  "last_infections", last_infections,
  "num_infections", num_infections,
)

save(joinpath(prefix_path,"summary.csv"), df)


text:
# Helper for parity function. Default criteria to calculate parity
_parity_func(ϵ) = x -> isnan(x) ? NaN : (1-ϵ) <= x <= 1/(1-ϵ)

"""
    parity(df, ϵ=nothing; func=nothing)

Takes the dataframe `df` returned from `disparity` function and adds columns for parity values corresponding to each disparity column.
It then calculates the parity values for measures that were passed for disparity calculation.

Parity is a boolean value indicating whether a fairness constraint is satisfied by disparity values.

The default fairness criteria for a disparity value `x` and fairness threshold `ϵ` for a group is:

    (1-ϵ) <= x <= 1/(1-ϵ)

Here ϵ is the fairness threshold which is required if default fairness constraint is used.

## Keywords

* `func=nothing` : single argument function specifying custom Fairness constraint for disparity instead of the default criteria
"""
function parity(df::DataFrame, ϵ=nothing; func=nothing)
    if func==nothing
        ϵ!=nothing || throw(ArgumentError("For using default fairness criteria, ϵ i.e. fairness threshold should be passed."))
        func = _parity_func(ϵ)
    end
    cols = names(df)
    M = string.(cols)
    for i in 1:length(cols)
        if !endswith(M[i], "_disparity")
            continue
        end
        m = replace(M[i], "disparity" => "parity")
        m = Symbol(m)
        df[:, m] = func.(df[:, cols[i]])
    end
    return df
end