text:
@testset "underscore" begin
    data = [1,2,3]
    @test @_(data |> sum(_)) == 6
    a = 10
    @test @_(data |> _ .+ a) == data .+ a
    @test @_(data |> _.^2 .+ _ |> sum) == sum(data.^2 .+ data)

    two_outputs(x) = (x, x.+1)
    @test @_(data |> two_outputs |> _[1] |> sum) == sum(data)
    @test @_(data |> two_outputs |> _[2] |> sum) == sum(data.+1)
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

"""
    struct MCMCMultiCycleBurnin <: MCMCBurninAlgorithm

A multi-cycle MCMC burn-in algorithm.

Constructors:

* ```$(FUNCTIONNAME)(; fields...)```

Fields:

$(TYPEDFIELDS)
"""
@with_kw struct MCMCMultiCycleBurnin <: MCMCBurninAlgorithm
    nsteps_per_cycle::Int64 = 10000
    max_ncycles::Int = 30
    nsteps_final::Int64 = div(nsteps_per_cycle, 10)
end

export MCMCMultiCycleBurnin

function mcmc_burnin!(
    outputs::Union{AbstractVector{<:DensitySampleVector},Nothing},
    tuners::AbstractVector{<:AbstractMCMCTunerInstance},
    chains::AbstractVector{<:MCMCIterator},
    burnin_alg::MCMCMultiCycleBurnin,
    convergence_test::MCMCConvergenceTest,
    strict_mode::Bool,
    nonzero_weights::Bool,
    callback::Function
)
    nchains = length(chains)

    @info "Begin tuning of $nchains MCMC chain(s)."

    cycles = zero(Int)
    successful = false
    while !successful && cycles < burnin_alg.max_ncycles
        cycles += 1

        new_outputs = DensitySampleVector.(chains)

        next_cycle!.(chains)

        tuning_reinit!.(tuners, chains, burnin_alg.nsteps_per_cycle)

        mcmc_iterate!(
            new_outputs, chains, tuners,
            max_nsteps = burnin_alg.nsteps_per_cycle,
            nonzero_weights = nonzero_weights,
            callback = callback
        )

        tuning_update!.(tuners, chains, new_outputs)
        isnothing(outputs) || append!.(outputs, new_outputs)

        ct_result = check_convergence!(convergence_test, chains, new_outputs)

        ntuned = count(c -> c.info.tuned, chains)
        nconverged = count(c -> c.info.converged, chains)
        successful = (ntuned == nconverged == nchains)

        callback(Val(:mcmc_burnin), tuners, chains)

        @info "MCMC Tuning cycle $cycles finished, $nchains chains, $ntuned tuned, $nconverged converged."
    end

    tuning_finalize!.(tuners, chains)

    if successful
        @info "MCMC tuning of $nchains chains successful after $cycles cycle(s)."
    else
        msg = "MCMC tuning of $nchains chains aborted after $cycles cycle(s)."
        if strict_mode
            throw(ErrorException(msg))
        else
            @warn msg
        end
    end

    if burnin_alg.nsteps_final > 0
        @info "Running post-tuning stabilization steps for $nchains MCMC chain(s)."

        next_cycle!.(chains)

        mcmc_iterate!(
            outputs, chains,
            max_nsteps = burnin_alg.nsteps_final,
            nonzero_weights = nonzero_weights,
            callback = callback
        )
    end

    successful
end

ext: jl


text:
module PulsarSearch
using StatsBase
include("stats.jl")

"""
`z_n_binned(profile, n) --> zsq`

``Z^2_n`` statistic for pulse profiles from binned events

See Bachetti+2021, arXiv:2012.11397

Parameters
----------
`profile` : array of floats
    The folded pulse profile (containing the number of
    photons falling in each pulse bin)

`n` : int
    Number of harmonics, including the fundamental

Returns
-------
`zsq` : float
    The value of the statistic

# Examples
```jldoctest
julia> using PulsarSearch

julia> z_n_binned(zeros(10), 2)
0.0

julia> z_n_binned(ones(10), 2) < 1e-30
true

julia> z_n_binned([10., 0., 0., 0., 0.], 2)
40.0
```
"""
function z_n_binned(profile::AbstractVector{}, n::Integer)::Float64
    total = sum(profile)
    N = length(profile)
    phase = range(0, stop = N - 1)  * (2pi/ N)

    if iszero(total)
        return 0.0
    end

    z = zero(Float64)
    for k in range(1, stop = n)
        s = 0.0
        c = 0.0
        for i in eachindex(profile)
            sk, ck = sincos(k * phase[i])
            s += profile[i] * sk
            c += profile[i] * ck
        end
        s = s^2
        c = c^2
        z += c + s
    end
    return z * 2 / total
end

"""
`z_n(phases, n) --> zsq`

``Z^2_n`` statistics, à la Buccheri+83, A&A, 128, 245, eq. 2.

Parameters
----------
`phase` : array of floats
    The phases of the events

`n` : int, default 2
    Number of harmonics, including the fundamental

Returns
-------
`zsq` : float
    The ``Z^2_n`` statistic

# Examples
```jldoctest
julia> using PulsarSearch

julia> z_n([10., 0., 0., 0., 0.], 2)
20.0
julia> z_n(ones(10), 2)
40.0
julia> z_n(Array([0.5]), 2)
0.0
```
"""
function z_n(phases::AbstractVector{T}, n::Integer)::T where {T<:AbstractFloat}
    N = length(phases)
    if N < 2
        return zero(T)
    end
    z = zero(T)
    twopiphase = 2 * pi * phases
    for k in range(1, stop = n)
        s = zero(T)
        c = zero(T)
        for ph in twopiphase
            sk, ck = sincos(k * ph)
            s += sk
            c += ck
        end
        s = s^2
        c = c^2
        z += c + s
    end
    return z * 2 / N
end

"""
`z_n_search(times, n, fmin, fmax [,oversample]) --> freqs, zsq_stat`

Calculate the ``Z^2_n`` statistics at trial frequencies in photon data.

Parameters
----------
`times` : array-like
    the event arrival times

`n` : int
    the number of harmonics in ``Z^2_n``

Other Parameters
----------------
`fmin` : float
    minimum pulse frequency to search

`fmax` : float
    maximum pulse frequency to search

`oversample` : float
    Oversampling factor with respect to the usual 1/T/n rule

Returns
-------
`fgrid` : array-like
    frequency grid of the epoch folding periodogram

`zsq_stat` : array-like
    the Z^2_n statistics corresponding to each frequency bin.
"""
function z_n_search(
    times::AbstractVector{T},
    n::Integer,
    fmin::Number,
    fmax::Number;
    oversample::Number = 2,
) where {T <: AbstractFloat}
    t0 = first(times)
    t1 = last(times)
    df = 1 / (t1 - t0) / oversample
    freqs = fmin:df:fmax
    N = length(freqs)
    stats = Vector{T}(undef, N)
    for i in eachindex(freqs)
        phases = times * freqs[i]
        phases .-= floor.(phases)
        stats[i] = z_n(phases, n)
    end
    return freqs, stats
end

"""
`z_n_search_hist(times, n, fmin, fmax [,oversample]) --> freqs, zsq_stat`

Calculate the ``Z^2_n`` statistics at trial frequencies in photon data.
Pre-bins the data using a histogram.
At the moment, it is _not_ faster. It will be after I add the 2-d hist
+ shift-and-add. But it allowed to test that `z_n_binned works` with
binned data.

Parameters
----------
`times` : array-like
    the event arrival times

`n` : int
    the number of harmonics in ``Z^2_n``

Other Parameters
----------------
`fmin` : float
    minimum pulse frequency to search

`fmax` : float
    maximum pulse frequency to search

`oversample` : float
    Oversampling factor with respect to the usual 1/T/n rule

Returns
-------
`fgrid` : array-like
    frequency grid of the epoch folding periodogram

`zsq_stat` : array-like
    the Z^2_n statistics corresponding to each frequency bin.
"""
function z_n_search_hist(
    times::AbstractVector{T},
    n::Integer,
    fmin::Number,
    fmax::Number;
    oversample::Number = 2,
    nbin::Integer = 16
) where {T <: AbstractFloat}
    t0 = first(times)
    t1 = last(times)
    df = 1 / (t1 - t0) / oversample
    freqs = fmin:df:fmax
    N = length(freqs)
    stats = Vector{T}(undef, N)
    for i in eachindex(freqs)
        phases = times * freqs[i]
        phases .-= floor.(phases)
        hist = fit(Histogram, phases, 0.0:(1/nbin):1.0)
        stats[i] = z_n_binned(hist.weights, n)
    end
    return freqs, stats
end

export z_n
export z_n_binned
export z_n_search
export z_n_search_hist
# Write your package code here.

end


text:
const NUCLEIC_IUPAC	= alphabet(nt"ACTGBDHKMNRSVWXY-.",false)

const AMINO_IUPAC	= alphabet(aa"ABCDEFGHIJKLMNOPQRSTUVWXYZ*-.",false)

const AMINO_20_UPPERCASE= alphabet(aa"ACDEFGHIKLMNPQRSTVWY",true)

const DNA_COMPLEMENT = zeros(Nucleotide,256)
      DNA_COMPLEMENT[Int('.')] = '.'
      DNA_COMPLEMENT[Int('-')] = '-'
      DNA_COMPLEMENT[Int('A')] = 'T'
      DNA_COMPLEMENT[Int('B')] = 'V'
      DNA_COMPLEMENT[Int('C')] = 'G'
      DNA_COMPLEMENT[Int('D')] = 'H'
      DNA_COMPLEMENT[Int('G')] = 'C'
      DNA_COMPLEMENT[Int('H')] = 'D'
      DNA_COMPLEMENT[Int('K')] = 'M'
      DNA_COMPLEMENT[Int('M')] = 'K'
      DNA_COMPLEMENT[Int('T')] = 'A'
      DNA_COMPLEMENT[Int('V')] = 'B'
      DNA_COMPLEMENT[Int('S')] = 'S'
      DNA_COMPLEMENT[Int('W')] = 'W'
      DNA_COMPLEMENT[Int('R')] = 'Y'
      DNA_COMPLEMENT[Int('Y')] = 'R'
      DNA_COMPLEMENT[Int('N')] = 'N'
      DNA_COMPLEMENT[Int('X')] = 'X'
      DNA_COMPLEMENT[Int('a')] = 't'
      DNA_COMPLEMENT[Int('b')] = 'v'
      DNA_COMPLEMENT[Int('c')] = 'g'
      DNA_COMPLEMENT[Int('d')] = 'h'
      DNA_COMPLEMENT[Int('g')] = 'c'
      DNA_COMPLEMENT[Int('h')] = 'd'
      DNA_COMPLEMENT[Int('k')] = 'm'
      DNA_COMPLEMENT[Int('m')] = 'k'
      DNA_COMPLEMENT[Int('t')] = 'a'
      DNA_COMPLEMENT[Int('v')] = 'b'
      DNA_COMPLEMENT[Int('s')] = 's'
      DNA_COMPLEMENT[Int('w')] = 'w'
      DNA_COMPLEMENT[Int('r')] = 'y'
      DNA_COMPLEMENT[Int('y')] = 'r'
      DNA_COMPLEMENT[Int('x')] = 'x'
      DNA_COMPLEMENT[Int('n')] = 'n'

const RNA_COMPLEMENT = zeros(Nucleotide,256)
      RNA_COMPLEMENT[Int('.')] = '.'
      RNA_COMPLEMENT[Int('-')] = '-'
      RNA_COMPLEMENT[Int('A')] = 'U'
      RNA_COMPLEMENT[Int('B')] = 'V'
      RNA_COMPLEMENT[Int('C')] = 'G'
      RNA_COMPLEMENT[Int('D')] = 'H'
      RNA_COMPLEMENT[Int('G')] = 'C'
      RNA_COMPLEMENT[Int('H')] = 'D'
      RNA_COMPLEMENT[Int('K')] = 'M'
      RNA_COMPLEMENT[Int('M')] = 'K'
      RNA_COMPLEMENT[Int('U')] = 'A'
      RNA_COMPLEMENT[Int('V')] = 'B'
      RNA_COMPLEMENT[Int('S')] = 'S'
      RNA_COMPLEMENT[Int('W')] = 'W'
      RNA_COMPLEMENT[Int('R')] = 'Y'
      RNA_COMPLEMENT[Int('Y')] = 'R'
      RNA_COMPLEMENT[Int('N')] = 'N'
      RNA_COMPLEMENT[Int('X')] = 'X'
      RNA_COMPLEMENT[Int('a')] = 'u'
      RNA_COMPLEMENT[Int('b')] = 'v'
      RNA_COMPLEMENT[Int('c')] = 'g'
      RNA_COMPLEMENT[Int('d')] = 'h'
      RNA_COMPLEMENT[Int('g')] = 'c'
      RNA_COMPLEMENT[Int('h')] = 'd'
      RNA_COMPLEMENT[Int('k')] = 'm'
      RNA_COMPLEMENT[Int('m')] = 'k'
      RNA_COMPLEMENT[Int('u')] = 'a'
      RNA_COMPLEMENT[Int('v')] = 'b'
      RNA_COMPLEMENT[Int('s')] = 's'
      RNA_COMPLEMENT[Int('w')] = 'w'
      RNA_COMPLEMENT[Int('r')] = 'y'
      RNA_COMPLEMENT[Int('y')] = 'r'
      RNA_COMPLEMENT[Int('x')] = 'x'
      RNA_COMPLEMENT[Int('n')] = 'n'

const AMINO_3LETTERS_TO_1 = Dict{ASCIIString,AminoAcid}(
		  "ALA"	=>	aminoacid('A'),
		  "ASX"	=>	aminoacid('B'),
		  "CYS"	=>	aminoacid('C'),
 		  "ASP"	=>  aminoacid('D'),
		  "GLU"	=>	aminoacid('E'),
		  "PHE"	=>	aminoacid('F'),
		  "GLY"	=>	aminoacid('G'),
		  "HYS"	=>	aminoacid('H'),
		  "ILE"	=>	aminoacid('I'),
		  "XLE"	=>	aminoacid('J'),
		  "LYS"	=>	aminoacid('K'),
		  "LEU"	=>	aminoacid('L'),
		  "MET"	=>	aminoacid('M'),
		  "ASN"	=>	aminoacid('N'),
		  "PYL"	=>	aminoacid('O'),
		  "PRO"	=>	aminoacid('P'),
		  "GLN"	=>	aminoacid('Q'),
		  "ARG"	=>	aminoacid('R'),
		  "SER"	=>	aminoacid('S'),
		  "THR"	=>	aminoacid('T'),
		  "SEC"	=>	aminoacid('U'),
		  "VAL"	=>	aminoacid('V'),
		  "TPR"	=>	aminoacid('W'),
		  "XXX"	=>	aminoacid('X'),
		  "TYR"	=>	aminoacid('Y'),
		  "GLX"	=>	aminoacid('Z') )

const AMINO_1LETTER_TO_3 = Dict{AminoAcid,ASCIIString}(
		  aminoacid('A')	=>	"ALA",
		  aminoacid('B')	=>	"ASX",
		  aminoacid('C')	=>	"CYS",
		  aminoacid('D')	=>	"ASP",
		  aminoacid('E')	=>	"GLU",
		  aminoacid('F')	=>	"PHE",
		  aminoacid('G')	=>	"GLY",
		  aminoacid('H')	=>	"HYS",
		  aminoacid('I')	=>	"ILE",
		  aminoacid('J')	=>	"XLE",
		  aminoacid('K')	=>	"LYS",
		  aminoacid('L')	=>	"LEU",
		  aminoacid('M')	=>	"MET",
		  aminoacid('N')	=>	"ASN",
		  aminoacid('O')	=>	"PYL",
		  aminoacid('P')	=>	"PRO",
		  aminoacid('Q')	=>	"GLN",
		  aminoacid('R')	=>	"ARG",
		  aminoacid('S')	=>	"SER",
		  aminoacid('T')	=>	"THR",
		  aminoacid('U')	=>	"SEC",
		  aminoacid('V')	=>	"VAL",
		  aminoacid('W')	=>	"TPR",
		  aminoacid('X')	=>	"XXX",
		  aminoacid('Y')	=>	"TYR",
		  aminoacid('Z')	=>	"GLX" )

const NUCLEIC_IUPAC_AMBIGUOUS_DICT = Dict{Nucleotide,Array{Nucleotide,1}}(
			nucleotide('.')	=>	nt".-",
			nucleotide('-')	=>	nt".-",
			nucleotide('A')	=>	nt"A",
			nucleotide('C')	=>	nt"C",
			nucleotide('G')	=>	nt"G",
			nucleotide('T')	=>	nt"TU",
			nucleotide('U')	=>	nt"TU",
			nucleotide('R')	=>	nt"GA",
			nucleotide('Y')	=>	nt"TCU",
			nucleotide('M')	=>	nt"AC",
			nucleotide('K')	=>	nt"GTU",
			nucleotide('S')	=>	nt"GC",
			nucleotide('W')	=>	nt"ATU",
			nucleotide('H')	=>	nt"ACTU",
			nucleotide('B')	=>	nt"GCTU",
			nucleotide('V')	=>	nt"ACG",
			nucleotide('D')	=>	nt"AGTU",
			nucleotide('N')	=>	nt"ACTUG")

const AMINO_IUPAC_EXTENDED_DICT = Dict{AminoAcid,Array{AminoAcid,1}}(
		  aminoacid('-')	=>	aminoacid(['-'; '.']),
		  aminoacid('.')	=>	aminoacid(['-'; '.']),
		  aminoacid('*')	=>	aminoacid(['*']),
		  aminoacid('A')	=>	aminoacid(['A']),
		  aminoacid('B')	=>	aminoacid(['N'; 'D']),
		  aminoacid('C')	=>	aminoacid(['C']),
		  aminoacid('D')	=>	aminoacid(['D']),
		  aminoacid('E')	=>	aminoacid(['E']),
		  aminoacid('F')	=>	aminoacid(['F']),
		  aminoacid('G')	=>	aminoacid(['G']),
		  aminoacid('H')	=>	aminoacid(['H']),
		  aminoacid('I')	=>	aminoacid(['I']),
		  aminoacid('J')	=>	aminoacid(['I'; 'L']),
		  aminoacid('K')	=>	aminoacid(['K']),
		  aminoacid('L')	=>	aminoacid(['L']),
		  aminoacid('M')	=>	aminoacid(['M']),
		  aminoacid('N')	=>	aminoacid(['N']),
		  aminoacid('O')	=>	aminoacid(['O']),
		  aminoacid('P')	=>	aminoacid(['P']),
		  aminoacid('Q')	=>	aminoacid(['Q']),
		  aminoacid('R')	=>	aminoacid(['R']),
		  aminoacid('S')	=>	aminoacid(['S']),
		  aminoacid('T')	=>	aminoacid(['T']),
		  aminoacid('U')	=>	aminoacid(['U']),
		  aminoacid('V')	=>	aminoacid(['V']),
		  aminoacid('W')	=>	aminoacid(['W']),
		  aminoacid('X')	=>	aminoacid(['A'; 'C'; 'D'; 'E'; 'F'; 'G'; 'H'; 'I'; 'K';'L'; 'M'; 'N'; 'P'; 'Q'; 'R'; 'S'; 'T'; 'V'; 'W'; 'Y']),
		  aminoacid('Y')	=>	aminoacid(['Y']),
		  aminoacid('Z')	=>	aminoacid(['Q'; 'E']))

## For Use in Regex ##
const _AMBIGUOUS_NUCLEIC_IUPAC = Dict{Nucleotide,Array{UInt8,1}}(
			nucleotide('T')	=>	UInt8['T';'U'],
			nucleotide('U')	=>	UInt8['T';'U'],
			nucleotide('R')	=>	UInt8['R';'G';'A'],
			nucleotide('Y')	=>	UInt8['Y';'T';'C';'U'],
			nucleotide('M')	=>	UInt8['M';'A';'C'],
			nucleotide('K')	=>	UInt8['K';'G';'T';'U'],
			nucleotide('S')	=>	UInt8['S';'G';'C'],
			nucleotide('W')	=>	UInt8['W';'A';'T';'U'],
			nucleotide('H')	=>	UInt8['W';'M';'Y';'H';'A';'C';'T';'U'],
			nucleotide('B')	=>	UInt8['S';'K';'B';'G';'C';'T';'U'],
			nucleotide('V')	=>	UInt8['S';'M';'R';'V';'A';'C';'G'],
			nucleotide('D')	=>	UInt8['W';'K';'R';'D';'A';'G';'T';'U'],
			nucleotide('N')	=>	UInt8['X';'W';'S';'K';'M';'Y';'R';'N';'A';'C';'T';'U';'G'],
			nucleotide('X')	=>	UInt8['X';'W';'S';'K';'M';'Y';'R';'N';'A';'C';'T';'U';'G'])

## For Use in Regex ##
const _AMBIGUOUS_AMINO_IUPAC = Dict{AminoAcid,Array{UInt8,1}}(
		  aminoacid('B')	=>	UInt8['B'; 'N'; 'D'],
		  aminoacid('X')	=>	UInt8['X'; 'B'; 'Z'; 'A'; 'C'; 'D'; 'E'; 'F'; 'G'; 'H'; 'I'; 'K';'L'; 'M'; 'N'; 'P'; 'Q'; 'R'; 'S'; 'T'; 'V'; 'W'; 'Y'; 'J' ],
		  aminoacid('Z')	=>	UInt8['Z'; 'Q'; 'E'],
		  aminoacid('J')	=>	UInt8['J'; 'L'; 'I'])


text:
module FilteringWrapper

using JavaCall

const Class = @jimport java.lang.Class
const ElevationFilter = @jimport org.orekit.estimation.measurements.filtering.ElevationFilter
const GroundStation = @jimport org.orekit.estimation.measurements.GroundStation
const JString = @jimport java.lang.String
const MeasurementFilter = @jimport org.orekit.estimation.measurements.filtering.MeasurementFilter
const Object = @jimport java.lang.Object
const ObservedMeasurement = @jimport org.orekit.estimation.measurements.ObservedMeasurement
const ResidualFilter = @jimport org.orekit.estimation.measurements.filtering.ResidualFilter
const SpacecraftState = @jimport org.orekit.propagation.SpacecraftState

export ElevationFilter
export MeasurementFilter
export ResidualFilter
export equals
export filter
export get_class
export hash_code
export notify
export notify_all
export to_string
export wait

include("elevation_filter.jl")
include("measurement_filter.jl")
include("residual_filter.jl")

end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Functions related to the quaternions.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export dquat, norm, vect

################################################################################
#                                 Constructors
################################################################################

"""
    Quaternion(q0::T0, q1::T1, q2::T2, q3::T3) where {T0, T1, T2, T3}

Create the following quaternion:

    q0 + q1.i + q2.j + q3.k

in which:

- `q0` is the real part of the quaternion.
- `q1` is the X component of the quaternion vectorial part.
- `q2` is the Y component of the quaternion vectorial part.
- `q3` is the Z component of the quaternion vectorial part.

!!! note
    The quaternion type is obtained by promoting `T0`, `T1`, `T2`, and `T3`.

# Examples

```jldoctest
julia> Quaternion(1, 0, 0, 0)
Quaternion{Int64}:
  + 1 + 0⋅i + 0⋅j + 0⋅k

julia> Quaternion(1, 0, 0, 0.0)
Quaternion{Float64}:
  + 1.0 + 0.0⋅i + 0.0⋅j + 0.0⋅k
```

---

    Quaternion(v::AbstractVector)

If the vector `v` has 3 components, then create a quaternion in which the real
part is `0` and the vectorial or imaginary part has the same components of the
vector `v`. In other words:

    q = 0 + v[1].i + v[2].j + v[3].k

Otherwise, if the vector `v` has 4 components, then create a quaternion in which
the elements match those of the input vector:

    q = v[1] + v[2].i + v[3].j + v[4].k

!!! note
    If the length of `v` is not 3 or 4, then an error is thrown.

# Examples

```jldoctest
julia> Quaternion([0, cosd(45), sind(45)])
Quaternion{Float64}:
  + 0.0 + 0.0⋅i + 0.707107⋅j + 0.707107⋅k

julia> Quaternion([cosd(45), 0, sind(45), 0])
Quaternion{Float64}:
  + 0.707107 + 0.0⋅i + 0.707107⋅j + 0.707107⋅k
```

---

    Quaternion(r::Number, v::AbstractVector)

Create a quaternion with real part `r` and vectorial or imaginary part `v`:

    r + v[1].i + v[2].j + v[3].k

!!! note
    The quaternion type is obtained by promoting the type of `r` and the
    elements of `v`.

# Examples

```jldoctest
julia> Quaternion(cosd(45), [0, sind(45), 0])
Quaternion{Float64}:
  + 0.707107 + 0.0⋅i + 0.707107⋅j + 0.707107⋅k
```

---

    Quaternion(u::UniformScaling{T}) where T
    Quaternion{T}(u::UniformScaling) where T
    Quaternion(u::UniformScaling, Q::Quaternion{T}) where T

Create the quaternion `u.λ + 0.i + 0.j + 0.k`.

If a quaternion is passed as in the third signature, then the new quaternion
will have the same type.

# Examples

```jldoctest
julia> Quaternion(I)
Quaternion{Bool}:
  + true + false⋅i + false⋅j + false⋅k

julia> Quaternion(1.0I)
Quaternion{Float64}:
  + 1.0 + 0.0⋅i + 0.0⋅j + 0.0⋅k

julia> q = Quaternion{Float32}(I)
Quaternion{Float32}:
  + 1.0 + 0.0⋅i + 0.0⋅j + 0.0⋅k

julia> Quaternion(I, q)
Quaternion{Float32}:
  + 1.0 + 0.0⋅i + 0.0⋅j + 0.0⋅k
```
"""
function Quaternion(q0::T0, q1::T1, q2::T2, q3::T3) where {T0, T1, T2, T3}
    T = promote_type(T0, T1, T2, T3)
    return Quaternion{T}(q0, q1, q2, q3)
end

function Quaternion(v::AbstractVector)
    # The vector must have 3 or 4 components.
    if length(v) != 3 && length(v) != 4
        throw(ArgumentError("The input vector must have 3 or 4 components."))
    end

    if length(v) == 3
        return Quaternion(0, v[1], v[2], v[3])
    else
        return Quaternion(v[1], v[2], v[3], v[4])
    end
end

Quaternion(r::Number, v::AbstractVector) = Quaternion(r, v[1], v[2], v[3])
Quaternion(u::UniformScaling{T}) where T = Quaternion{T}(T(u.λ), T(0), T(0), T(0))
Quaternion{T}(u::UniformScaling) where T = Quaternion{T}(T(u.λ), T(0), T(0), T(0))
Quaternion(::UniformScaling, ::Quaternion{T}) where T = Quaternion{T}(I)

################################################################################
#                                  Operations
################################################################################

# Operation: +
# ==============================================================================

"""
    +(qa::Quaternion, qb::Quaternion)

Compute `qa + qb`.

If one of the operands is a `UniformScaling`:

    +(u::UniformScaling, q::Quaternion)
    +(q::Quaternion, u::UniformScaling)

then it is considered as the quaternion `u.λ + 0 ⋅ i + 0 ⋅ j + 0 ⋅ k`.

# Examples

```jldoctest
julia> q1 = Quaternion(1, 0, 0, 0)
Quaternion{Int64}:
  + 1 + 0⋅i + 0⋅j + 0⋅k

julia> q2 = Quaternion(0, cosd(45), 0, sind(45))
Quaternion{Float64}:
  + 0.0 + 0.707107⋅i + 0.0⋅j + 0.0⋅k

julia> q1 + q2
Quaternion{Float64}:
  + 1.0 + 0.707107⋅i + 0.0⋅j + 0.0⋅k

julia> q1 + I
Quaternion{Int64}:
  + 2 + 0⋅i + 0⋅j + 0⋅k
```
"""
@inline function +(qa::Quaternion, qb::Quaternion)
    return Quaternion(qa.q0 + qb.q0, qa.q1 + qb.q1, qa.q2 + qb.q2, qa.q3 + qb.q3)
end

@inline +(u::UniformScaling, q::Quaternion) = Quaternion(u.λ + q.q0, q.q1, q.q2, q.q3)
@inline +(q::Quaternion, u::UniformScaling) = u+q

# Operation: -
# ==============================================================================

"""
    -(q::Quaternion)

Return the quaterion `-q`.

# Examples

```jldoctest
julia> q = Quaternion(1, 0, 0, 0)
Quaternion{Int64}:
  + 1 + 0⋅i + 0⋅j + 0⋅k

julia> -q
Quaternion{Int64}:
  - 1 + 0⋅i + 0⋅j + 0⋅k
```
"""
@inline -(q::Quaternion) = -1*q

"""
    -(qa::Quaternion, qb::Quaternion)

Compute `qa - qb`.

If one of the operands is a `UniformScaling`:

    -(u::UniformScaling, q::Quaternion)
    -(q::Quaternion, u::UniformScaling)

then it is considered as the quaternion `u.λ + 0 ⋅ i + 0 ⋅ j + 0 ⋅ k`.

# Examples

```jldoctest
julia> q1 = Quaternion(1, 0, 0, 0)
Quaternion{Int64}:
  + 1 + 0⋅i + 0⋅j + 0⋅k

julia> q2 = Quaternion(0, cosd(45), 0, sind(45))
Quaternion{Float64}:
  + 0.0 + 0.707107⋅i + 0.0⋅j + 0.0⋅k

julia> q1 - q2
Quaternion{Float64}:
  + 1.0 - 0.707107⋅i + 0.0⋅j + 0.0⋅k

julia> q1 - I
Quaternion{Int64}:
  + 0 + 0⋅i + 0⋅j + 0⋅k
```
"""
@inline function -(qa::Quaternion, qb::Quaternion)
    return Quaternion(qa.q0 - qb.q0, qa.q1 - qb.q1, qa.q2 - qb.q2, qa.q3 - qb.q3)
end

@inline -(u::UniformScaling, q::Quaternion) = Quaternion(u.λ - q.q0, -q.q1, -q.q2, -q.q3)
@inline -(q::Quaternion, u::UniformScaling) = (-u) + q

# Operation: *
# ==============================================================================

"""
    *(λ::Number, q::Quaternion)
    *(q::Quaternion, λ::Number)

Compute `λ * q` or `q * λ`, in which `λ` is a scalar.

# Examples

```jldoctest
julia> q = Quaternion(1, 0, 0, 0)
Quaternion{Int64}:
  + 1 + 0⋅i + 0⋅j + 0⋅k

julia> 2 * q
Quaternion{Int64}:
  + 2 + 0⋅i + 0⋅j + 0⋅k
```
"""
@inline *(λ::Number, q::Quaternion) = Quaternion(λ * q.q0, λ * q.q1, λ * q.q2, λ * q.q3)
@inline *(q::Quaternion, λ::Number) = Quaternion(λ * q.q0, λ * q.q1, λ * q.q2, λ * q.q3)

"""
    *(q1::Quaternion, q2::Quaternion)

Compute the quaternion multiplication `q1 * q2` (Hamilton product).

If one of the operands is a `UniformScaling`:

    *(u::UniformScaling, q::Quaternion)
    *(q::Quaternion, u::UniformScaling)

then it is considered as the quaternion `u.λ + 0 ⋅ i + 0 ⋅ j + 0 ⋅ k`.

# Examples

```jldoctest
julia> q1 = Quaternion(cosd(30), 0, sind(30), 0)
Quaternion{Float64}:
  + 0.866025 + 0.0⋅i + 0.5⋅j + 0.5⋅k

julia> q2 = Quaternion(cosd(60), 0, sind(60), 0)
Quaternion{Float64}:
  + 0.5 + 0.0⋅i + 0.866025⋅j + 0.866025⋅k

julia> q1 * q2
Quaternion{Float64}:
  + 0.0 + 0.0⋅i + 1.0⋅j + 1.0⋅k

julia> I * q1
Quaternion{Float64}:
  + 0.866025 + 0.0⋅i + 0.5⋅j + 0.5⋅k
```
"""
@inline function *(q1::Quaternion, q2::Quaternion)
    return Quaternion(
        q1.q0 * q2.q0 - q1.q1 * q2.q1 - q1.q2 * q2.q2 - q1.q3 * q2.q3,
        q1.q0 * q2.q1 + q1.q1 * q2.q0 + q1.q2 * q2.q3 - q1.q3 * q2.q2,
        q1.q0 * q2.q2 - q1.q1 * q2.q3 + q1.q2 * q2.q0 + q1.q3 * q2.q1,
        q1.q0 * q2.q3 + q1.q1 * q2.q2 - q1.q2 * q2.q1 + q1.q3 * q2.q0
    )
end

@inline *(u::UniformScaling, q::Quaternion) = Quaternion(u)*q
@inline *(q::Quaternion, u::UniformScaling) = q*Quaternion(u)

"""
    *(v::AbstractVector, q::Quaternion)
    *(q::Quaternion, v::AbstractVector)

Compute the multiplication `qv * q` or `q * qv` in which `qv` is a quaternion
with real part `0` and vectorial/imaginary part `v` (Hamilton product).

# Examples

```jldoctest
julia> q = Quaternion(1, 0, 0, 0)
Quaternion{Int64}:
  + 1 + 0⋅i + 0⋅j + 0⋅k

julia> v = [0, cosd(60), sind(60)]
3-element Vector{Float64}:
 0.0
 0.5
 0.8660254037844386

julia> q * v
Quaternion{Float64}:
  + 0.0 + 0.0⋅i + 0.5⋅j + 0.5⋅k
```
"""
@inline function *(v::AbstractVector, q::Quaternion)
    return Quaternion(
        -v[1] * q.q1 - v[2] * q.q2 - v[3] * q.q3,
        +v[1] * q.q0 + v[2] * q.q3 - v[3] * q.q2,
        -v[1] * q.q3 + v[2] * q.q0 + v[3] * q.q1,
        +v[1] * q.q2 - v[2] * q.q1 + v[3] * q.q0
    )
end

@inline function *(q::Quaternion, v::AbstractVector)
    return Quaternion(
        - q.q1 * v[1] - q.q2 * v[2] - q.q3 * v[3],
        + q.q0 * v[1] + q.q2 * v[3] - q.q3 * v[2],
        + q.q0 * v[2] - q.q1 * v[3] + q.q3 * v[1],
        + q.q0 * v[3] + q.q1 * v[2] - q.q2 * v[1]
    )
end

# Operation: /
# ==============================================================================

"""
    /(λ::Number, q::Quaternion)
    /(q::Quaternion, λ::Number)

Compute the division `λ / q` or `q / λ`, in which `λ` is a scalar.

# Examples

```jldoctest
julia> q = Quaternion(2, 0, 0, 0)
Quaternion{Int64}:
  + 2 + 0⋅i + 0⋅j + 0⋅k

julia> q / 2
Quaternion{Float64}:
  + 1.0 + 0.0⋅i + 0.0⋅j + 0.0⋅k
```
"""
@inline function /(λ::Number, q::Quaternion)
    # Compute `λ*(1/q)`.
    norm_q² = q.q0 * q.q0 + q.q1 * q.q1 + q.q2 * q.q2 + q.q3 * q.q3

    return Quaternion(
        +λ * q.q0 / norm_q²,
        -λ * q.q1 / norm_q²,
        -λ * q.q2 / norm_q²,
        -λ * q.q3 / norm_q²
    )
end

@inline /(q::Quaternion, λ::Number) = q * (1 / λ)

"""
    /(q1::Quaternion, q2::Quaternion)

Compute `q1 * inv(q2)` (Hamilton product).

If one of the operands is a `UniformScaling`:

    /(u::UniformScaling, q::Quaternion)
    /(q::Quaternion, u::UniformScaling)

then it is considered as the quaternion `u.λ + 0 ⋅ i + 0 ⋅ j + 0 ⋅ k`.

# Examples

```jldoctest
julia> q1 = Quaternion(cosd(75), 0, sind(75), 0)
Quaternion{Float64}:
  + 0.258819 + 0.0⋅i + 0.965926⋅j + 0.965926⋅k

julia> q2 = Quaternion(cosd(30), 0, sind(30), 0)
Quaternion{Float64}:
  + 0.866025 + 0.0⋅i + 0.5⋅j + 0.5⋅k

julia> q1 / q2
Quaternion{Float64}:
  + 0.707107 + 0.0⋅i + 0.707107⋅j + 0.707107⋅k

julia> q1 / (2 * I)
Quaternion{Float64}:
  + 0.12941 + 0.0⋅i + 0.482963⋅j + 0.482963⋅k
```
"""
@inline /(q1::Quaternion, q2::Quaternion) = q1 * inv(q2)
@inline /(u::UniformScaling, q::Quaternion) = Quaternion(u) / q
@inline /(q::Quaternion, u::UniformScaling) = q / Quaternion(u)

# Operation: \
# ==============================================================================

"""
    \\(q1::Quaternion, q2::Quaternion)

Compute `inv(q1) * q2`.

If one of the operands is a `UniformScaling`:

    \\(u::UniformScaling, q::Quaternion)
    \\(q::Quaternion, u::UniformScaling)

then it is considered as the quaternion `u.λ + 0 ⋅ i + 0 ⋅ j + 0 ⋅ k`.

# Examples

```jldoctest
julia> q1 = Quaternion(cosd(75), 0, sind(75), 0)
Quaternion{Float64}:
  + 0.258819 + 0.0⋅i + 0.965926⋅j + 0.965926⋅k

julia> q2 = Quaternion(cosd(30), 0, sind(30), 0)
Quaternion{Float64}:
  + 0.866025 + 0.0⋅i + 0.5⋅j + 0.5⋅k

julia> q2 \\ q1
Quaternion{Float64}:
  + 0.707107 + 0.0⋅i + 0.707107⋅j + 0.707107⋅k
```
"""
@inline \(q1::Quaternion, q2::Quaternion) = inv(q1) * q2
@inline \(u::UniformScaling, q::Quaternion) = inv(Quaternion(u)) * q
@inline \(q::Quaternion, u::UniformScaling) = inv(q) * Quaternion(u)

"""
    \\(v::AbstractVector, q::Quaternion)
    \\(q::Quaternion, v::AbstractVector)

Compute the division `qv \\ q` or `q \\ qv` in which `qv` is a quaternion with
real part `0` and vectorial/imaginary part `v` (Hamilton product).

# Examples

```jldoctest
julia> q = Quaternion(1, 0, 0, 0)
Quaternion{Int64}:
  + 1 + 0⋅i + 0⋅j + 0⋅k

julia> v = [0, cosd(60), sind(60)]
3-element Vector{Float64}:
 0.0
 0.5
 0.8660254037844386

julia> v \\ q
Quaternion{Float64}:
  + 0.0 + 0.0⋅i - 0.5⋅j - 0.5⋅k
```
"""
@inline \(q::Quaternion, v::AbstractVector) = inv(q)*v
@inline \(v::AbstractVector, q::Quaternion) = inv(Quaternion(v))*q

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

"""
    conj(q::Quaternion)

Compute the complex conjugate of the quaternion `q`:

    q0 - q1.i - q2.j - q3.k

See also: [`inv`](@ref)

# Examples

```jldoctest
julia> q = Quaternion(1, cosd(75), 0, sind(75))
Quaternion{Float64}:
  + 1.0 + 0.258819⋅i + 0.0⋅j + 0.0⋅k

julia> conj(q)
Quaternion{Float64}:
  + 1.0 - 0.258819⋅i - 0.0⋅j - 0.0⋅k
```
"""
@inline conj(q::Quaternion) = Quaternion(q.q0, -q.q1, -q.q2, -q.q3)

"""
    copy(q::Quaternion{T}) where T

Create a copy of the quaternion `q`.

"""
@inline copy(q::Quaternion{T}) where T = Quaternion{T}(q.q0, q.q1, q.q2, q.q3)
# TODO: Do we really need a copy functions since Quaternion is not mutable?
# Maybe it is necessary for DifferentialEquations.jl

"""
    imag(q::Quaternion)

Return the vectorial or imaginary part of the quaternion `q` represented by a
3 × 1 vector of type `SVector{3}`.

See also: [`real`](@ref), [`vect`](@ref)

# Examples

```jldoctest
julia> q = Quaternion(cosd(75), 0, sind(75), 0)
Quaternion{Float64}:
  + 0.258819 + 0.0⋅i + 0.965926⋅j + 0.965926⋅k

julia> imag(q)
3-element StaticArrays.SVector{3, Float64} with indices SOneTo(3):
 0.0
 0.9659258262890683
 0.0
```
"""
@inline imag(q::Quaternion) = SVector{3}(q.q1, q.q2, q.q3)

"""
    inv(q::Quaternion)

Compute the inverse of the quaternion `q`:

    conj(q)
    -------
      |q|²

See also: [`conj`](@ref)

# Examples

```jldoctest
julia> q = Quaternion(1, 0, cosd(75), sind(75))
Quaternion{Float64}:
  + 1.0 + 0.0⋅i + 0.258819⋅j + 0.258819⋅k

julia> inv(q)
Quaternion{Float64}:
  + 0.5 - 0.0⋅i - 0.12941⋅j - 0.12941⋅k
```
"""
@inline function inv(q::Quaternion)
    # Compute the inverse of the quaternion.
    norm_q² = q.q0 * q.q0 + q.q1 * q.q1 + q.q2 * q.q2 + q.q3 * q.q3

    return Quaternion(
        +q.q0 / norm_q²,
        -q.q1 / norm_q²,
        -q.q2 / norm_q²,
        -q.q3 / norm_q²
    )
end

"""
    norm(q::Quaternion)

Compute the Euclidean norm of the quaternion `q`:

    √(q0² + q1² + q2² + q3²)

# Examples

```jldoctest
julia> q = Quaternion(cosd(75), 0, sind(75), 0)
Quaternion{Float64}:
  + 0.258819 + 0.0⋅i + 0.965926⋅j + 0.965926⋅k

julia> norm(q)
1.0
```
"""
@inline norm(q::Quaternion) = √(q.q0 * q.q0 + q.q1 * q.q1 + q.q2 * q.q2 + q.q3 * q.q3)

"""
    real(q::Quaternion)

Return the real part of the quaternion `q`: `q0`.

See also: [`imag`](@ref), [`vect`](@ref)

# Examples

```jldoctest
julia> q = Quaternion(cosd(75), 0, sind(75), 0)
Quaternion{Float64}:
  + 0.258819 + 0.0⋅i + 0.965926⋅j + 0.965926⋅k

julia> real(q)
0.25881904510252074
```
"""
@inline real(q::Quaternion) = q.q0

"""
    vect(q::Quaternion)

Return the vectorial or imaginary part of the quaternion `q` represented by a
3 × 1 vector of type `SVector{3}`.

See also: [`imag`](@ref), [`real`](@ref)

# Examples

```jldoctest
julia> q = Quaternion(cosd(75), 0, sind(75), 0)
Quaternion{Float64}:
  + 0.258819 + 0.0⋅i + 0.965926⋅j + 0.965926⋅k

julia> vect(q)
3-element StaticArrays.SVector{3, Float64} with indices SOneTo(3):
 0.0
 0.9659258262890683
 0.0
```
"""
@inline vect(q::Quaternion) = SVector{3}(q.q1, q.q2, q.q3)

"""
    zeros(::Type{Quaternion{T}}) where T

Create the null quaternion: `0 + 0 ⋅ i + 0 ⋅ j + 0 ⋅ k`.

The type of the new quaternion will be `T`. If the type `T` is omitted, then it
defaults to `Float64`.

# Example

```jldoctest
julia> zeros(Quaternion{Float32})
Quaternion{Float32}:
  + 0.0 + 0.0⋅i + 0.0⋅j + 0.0⋅k

julia> zeros(Quaternion)
Quaternion{Float64}:
  + 0.0 + 0.0⋅i + 0.0⋅j + 0.0⋅k
```

---

    zeros(q::Quaternion{T}) where T

Create the null quaternion with the same type `T` of another quaternion `q`.

# Examples

```jldoctest
julia> q = Quaternion{Float32}(1, 0, 0, 0)
Quaternion{Float32}:
  + 1.0 + 0.0⋅i + 0.0⋅j + 0.0⋅k

julia> zeros(q)
Quaternion{Float32}:
  + 0.0 + 0.0⋅i + 0.0⋅j + 0.0⋅k
```
"""
@inline zeros(::Type{Quaternion{T}}) where T = Quaternion{T}(T(0), T(0), T(0), T(0))
@inline zeros(::Type{Quaternion}) = Quaternion{Float64}(0, 0, 0, 0)
@inline zeros(q::Quaternion{T}) where T = zeros(Quaternion{T})

# The following functions make sure that a quaternion is an iterable object.
# This allows broadcasting without allocations.
Base.IndexStyle(::Type{<:Quaternion}) = IndexLinear()
@inline Base.size(::Quaternion) = (4,)

@inline function Base.getindex(q::Quaternion, i::Int)
    if i == 1
        return q.q0
    elseif i == 2
        return q.q1
    elseif i == 3
        return q.q2
    elseif i == 4
        return q.q3
    else
        throw(BoundsError(q, i))
    end
end

################################################################################
#                                      IO
################################################################################

function show(io::IO, q::Quaternion{T}) where T
    # Check if the user wants compact printing, defaulting to `true`.
    compact_printing = get(io, :compact, true)::Bool

    # Get the absolute values using `print`.
    aq0 = sprint(print, abs(q.q0), context = :compact => compact_printing)
    aq1 = sprint(print, abs(q.q1), context = :compact => compact_printing)
    aq2 = sprint(print, abs(q.q2), context = :compact => compact_printing)
    aq3 = sprint(print, abs(q.q3), context = :compact => compact_printing)

    # Get the signs.
    sq0 = signbit(q.q0) ? "-" : "+"
    sq1 = signbit(q.q1) ? "-" : "+"
    sq2 = signbit(q.q2) ? "-" : "+"
    sq3 = signbit(q.q3) ? "-" : "+"

    print(io, "Quaternion{$(T)}: ")
    print(io,
        sq0, " ", aq0, " ",
        sq1, " ", aq1, "⋅i ",
        sq2, " ", aq2, "⋅j ",
        sq3, " ", aq3, "⋅k"
    )

    return nothing
end

function show(io::IO, mime::MIME"text/plain", q::Quaternion{T}) where T
    # Check if the user wants colors.
    color = get(io, :color, false)::Bool

    # Check if the user wants compact printing, defaulting to `true`.
    compact_printing = get(io, :compact, true)::Bool

    # Assemble the context.
    context = IOContext(io, :color => color, :compact => compact_printing)

    b = color ? string(_b) : ""
    d = color ? string(_d) : ""

    # Get the absolute values using `print`.
    aq0 = sprint(print, abs(q.q0), context = context)
    aq1 = sprint(print, abs(q.q1), context = context)
    aq2 = sprint(print, abs(q.q2), context = context)
    aq3 = sprint(print, abs(q.q3), context = context)

    # Get the signs.
    sq0 = signbit(q.q0) ? "-" : "+"
    sq1 = signbit(q.q1) ? "-" : "+"
    sq2 = signbit(q.q2) ? "-" : "+"
    sq3 = signbit(q.q3) ? "-" : "+"

    println(io, "Quaternion{$(T)}:")
    print(io,
        "  ",
        sq0, " ", aq0, " ",
        sq1, " ", aq1, "⋅", b, "i", d, " ",
        sq2, " ", aq2, "⋅", b, "j", d, " ",
        sq2, " ", aq2, "⋅", b, "k", d
    )

    return nothing
end

################################################################################
#                                  Kinematics
################################################################################

"""
    dquat(qba::Quaternion, wba_b::AbstractVector)

Compute the time-derivative of the quaternion `qba` that rotates a reference
frame `a` into alignment to the reference frame `b` in which the angular
velocity of `b` with respect to `a`, and represented in `b`, is `wba_b`.

# Examples

```jldoctest
julia> q = Quaternion(1.0I);

julia> dquat(q,[1;0;0])
Quaternion{Float64}:
  - 0.0 + 0.5⋅i + 0.0⋅j + 0.0⋅k
```
"""
function dquat(qba::Quaternion, wba_b::AbstractVector)
    # Auxiliary variable.
    w = wba_b

    # Check the dimensions.
    if length(wba_b) != 3
        throw(ArgumentError("The angular velocity vector must have three components."))
    end

    # Return the time-derivative.
    #         1         x
    #   dq = --- (wba_b) . q
    #         2

    return Quaternion(
                       -w[1]/2*qba.q1 -w[2]/2*qba.q2 -w[3]/2*qba.q3,
        +w[1]/2*qba.q0                +w[3]/2*qba.q2 -w[2]/2*qba.q3,
        +w[2]/2*qba.q0 -w[3]/2*qba.q1                +w[1]/2*qba.q3,
        +w[3]/2*qba.q0 +w[2]/2*qba.q1 -w[1]/2*qba.q2
    )
end


text:
using Base.Filesystem
using CSV
using DataFrames

"Save objective values of a population to a file"
function log_population(population::Vector{Solution}, gen, log_dir)
    if ~isdir(log_dir)
        mkpath(log_dir)
    end

    popsize = length(population)
    nobjs = population[1].problem.nobjs

    df = DataFrame()
    df[!, "iteration"] = [gen for p = 1:popsize]
    for i=1:nobjs
        colname = "obj_$i"
        df[!, colname] = [population[p].objectives[i] for p = 1:popsize]
    end
    df[!, "born"] = [population[p].generation for p = 1:popsize]
    df[!, "rank"] = [population[p].rank for p = 1:popsize]

    output_file = Base.Filesystem.joinpath(log_dir, "log_$gen.csv")
    CSV.write(output_file, df)
    @show(gen)

end

max_stars_repo_name: phawthorne/echidna

max_issues_repo_name: phawthorne/echidna

max_forks_repo_name: phawthorne/echidna


text:
"""
    inbounds(x, max, overflow)

Check grid boundaries for a single coordinate and max value or a tuple
of coorinates and max values.

Returns a tuple containing the coordinate(s) followed by a boolean `true`
if the cell is in bounds, `false` if not.

Overflow of type [`RemoveOverflow`](@ref) returns the coordinate and `false` to skip
coordinates that overflow outside of the grid.
[`WrapOverflow`](@ref) returns a tuple with the current position or it's
wrapped equivalent, and `true` as it is allways in-bounds.
"""
@inline inbounds(xs::Tuple, maxs::Tuple, overflow) = begin
    a, inbounds_a = inbounds(xs[1], maxs[1], overflow)
    b, inbounds_b = inbounds(xs[2], maxs[2], overflow)
    (a, b), inbounds_a & inbounds_b
end
@inline inbounds(x::Number, max::Number, overflow::RemoveOverflow) =
    x, isinbounds(x, max, overflow)
@inline inbounds(x::Number, max::Number, overflow::WrapOverflow) =
    if x < oneunit(x)
        max + rem(x, max), true
    elseif x > max
        rem(x, max), true
    else
        x, true
    end

@inline isinbounds(xs::Tuple, maxs::Tuple, overflow) = all(isinbounds.(xs, maxs, Ref(overflow)))
@inline isinbounds(x::Number, max::Number, overflow::RemoveOverflow) =
    x > zero(x) && x <= max

"""
    sizefromradius(radius)

Get the size of a neighborhood dimension from its radius,
which is always 2r + 1.
"""
hoodsize(radius::Integer) = 2radius + 1

"""
Return a tuple of the base types of the rules in the ruleset
"""
ruletypes(ruleset::Ruleset) = ruletypes(typeof(Ruleset.rules))
ruletypes(t::Type) = t.name.wrapper
ruletypes(ts::Type{<:Tuple}) = (ruletypes.(ts.parameters)...,)

@inline ismasked(data::AbstractSimData, i...) = ismasked(mask(data), i...)
@inline ismasked(mask::Nothing, i...) = false
@inline ismasked(mask::AbstractArray, i...) = @inbounds return !mask[i...]


text:
using DataFrames: names!
using StatsBase: mean, std, sem
import StatsBase: sem
import Base.size
import Base.names

struct ChainDataFrame
    name::String
    df::DataFrame
end

ChainDataFrame(df::DataFrame) = ChainDataFrame("", df)

Base.size(c::ChainDataFrame) = size(c.df)
Base.names(c::ChainDataFrame) = names(c.df)
function Base.show(io::IO, c::ChainDataFrame)
    println(io, c.name)
    show(io, c.df, summary = false, allrows=true, allcols=true)
end

Base.getindex(c::ChainDataFrame, args...) = getindex(c.df, args...)
Base.getindex(c::ChainDataFrame, s::Union{Symbol, Vector{Symbol}}) = c.df[s]
Base.isequal(cs1::Vector{ChainDataFrame}, cs2::Vector{ChainDataFrame}) = isequal.(cs1, cs2)
Base.isequal(c1::ChainDataFrame, c2::ChainDataFrame) = isequal(c1, c2)

function Base.show(io::IO, cs::Vector{ChainDataFrame})
    println(io, summary(cs))
    for i in cs
        println(io)
        println(io, i.name)
        show(io, i.df, allrows=true, summary=false)
        println(io)
    end
end

# Allows overriding of `display`
function Base.show(io::IO, ::MIME"text/plain", cs::Vector{ChainDataFrame})
    show(io, cs)
end

function Base.getindex(c::ChainDataFrame, s::Union{Symbol, Vector{Symbol}}, m)
    s = s isa AbstractArray ? s : [s]
    s_ind = indexin(s, c.df[:, :parameters])
    return c.df[s_ind, m]
end

function Base.getindex(c::ChainDataFrame,
        s1::Vector{Symbol},
        s2::Union{Symbol, Vector{Symbol}})
    return c.df[map(x -> x in s1, c.df.parameters), s2]
end

function Base.getindex(c::ChainDataFrame,
        s1::Symbol,
        s2::Union{Symbol, Vector{Symbol}})
    return c.df[c.df.parameters .== s1, s2]
end

Base.lastindex(c::ChainDataFrame, i::Integer) = lastindex(c.df, i)

function Base.convert(::Type{T}, cs::Array{ChainDataFrame,1}) where T<:Array
    arrs = [convert(T, cs[j].df[:, 2:end]) for j in 1:length(cs)]
    return cat(arrs..., dims = 3)
end
"""

# Summarize a Chains object formatted as a DataFrame

Summarize method for an MCMCChains.Chains object.

### Method
```julia
  summarize(
    chn::MCMCChains.AbstractChains,
    funs...;
    sections::Vector{Symbol}=[:parameters],
    func_names=[],
    etype=:bm
  )
```

### Required arguments
```julia
* `chn` : Chains object to convert to a DataFrame-formatted summary
```

### Optional arguments
```julia
* `funs...` : zero or more vector functions, e.g. mean, std, etc.
* `sections = [:parameters]` : Sections from the Chains object to be included
* `etype = :bm` : Default for df_mcse
```

### Examples
```julia
* `summarize(chns)` : Complete chain summary
* `summarize(chns[[:parm1, :parm2]])` : Chain summary of selected parameters
* `summarize(chns, sections=[:parameters])`  : Chain summary of :parameters section
* `summarize(chns, sections=[:parameters, :internals])` : Chain summary for multiple sections
```

"""
function summarize(chn::Chains, funs...;
        sections::Union{Symbol, Vector{Symbol}}=Symbol[:parameters],
        func_names=[],
        append_chains::Bool=true,
        showall::Bool=false,
        name::String="",
        additional_df=nothing)
    # Check that we actually have :parameters.
    showall = showall ? true : !in(:parameters, keys(chn.name_map))

    # If we weren't given any functions, fall back on summary stats.
    if length(funs) == 0
        return summarystats(chn,
            sections=sections,
            showall=showall)
    end

    # Generate a dataframe to work on.
    df = DataFrame(chn, sections, showall=showall, append_chains=append_chains)

    # If no function names were given, make a new list.
    func_names = length(func_names) == 0 ?
        handle_funs(funs) : Symbol.(func_names)

    # Do all the math, make columns.
    columns = if append_chains
        vcat([names(df)],
             [[f(col) for col = eachcol(df, false)] for f in funs])
    else
        [vcat([names(df[1])],
              [[f(col) for col = eachcol(i, false)] for f in funs]) for i in df]
    end

    # Make a vector of column names.
    colnames = vcat([:parameters], func_names)

    # Build the dataframes.
    ret_df = if append_chains
        DataFrame(columns, colnames)
    else
        [DataFrame(i, colnames) for i in columns]
    end

    if additional_df != nothing
        if append_chains
            ret_df = join(ret_df, additional_df, on=:parameters)
        else
            ret_df = [join(r, additional_df, on=:parameters) for r in ret_df]
        end
    end

    ret_df = if append_chains
        ChainDataFrame(name, ret_df)
    else
        [ChainDataFrame(name, r) for r in ret_df]
    end

    return ret_df
end

function handle_funs(fns)
    tmp =  [string(f) for f in fns]
    Symbol.([split(tmp[i], ".")[end] for i in 1:length(tmp)])
end


text:
@benchgroup "traversals" begin
  @benchgroup "graphs" begin
    for (name, g) in GRAPHS
        @bench "$(name): bfs_tree" LightGraphs.bfs_tree($g, 1)
        @bench "$(name): dfs_tree" LightGraphs.dfs_tree($g, 1)
    end
  end # graphs
  @benchgroup "digraphs" begin
    for (name, g) in DIGRAPHS
        @bench "$(name): bfs_tree" LightGraphs.bfs_tree($g, 1)
        @bench "$(name): dfs_tree" LightGraphs.dfs_tree($g, 1)
    end
  end # digraphs
end # traversals


text:
using EntanglementDist
using SQLite
using Gadfly
using DataFrames

import Base.convert
convert(x::Type{Nullable{Float64}}, y::Int64) = convert(Nullable{Float64}, convert(Float64, y))

# The plots generated by this script are stored in:
# - `./pdfplots` in PDF format
# - `./plots` in PNG format

### settings, change these for different plots ###
# This can be any parameter, it depends on how it has been saved in the database
ps = collect(0.4:0.1:1)

# The collection of states that are to be plotted
states = AbstractString["pptRelax: Double Werner"]

# The database where the data is stored
db = SQLite.DB("test.sqlite")

# Title of the plot as a function of the parameter
plottitle(p) = "Double Werner (p = $p)"

# Due to a bug in SQLite it might be the case that
# an error is thrown because a certain column in the database
# includes both floats and integers. Untill this is resolved
# by the SQLite.jl team, one might resort to adding a very
# small number (e.g. 10^-8) to the relevant columns. Typically
# this happens to delta_min if `0`, `1` are present.
#
# Example query:
# table = SQLite.query(db, "UPDATE `RainsProb` SET `delta_max` = `delta_max` + 0.00000001")
for p in ps
  for state in states
    # calculate the performance of the protocols
    if state == "pptRelax: Double Werner" || state == "Triple Werner"
      stateFn = wernerState
    elseif state == "Double Ronald"
      stateFn = sState
    elseif state == "Double Ronald 2"
      stateFn = rState
    elseif state == "Ronald 2 averaged phase" || state == "ronald2StateCorrPhase"
      stateFn = rStateCorrPhase
    end

    # get data from DB
    q = string("SELECT * FROM `pptRelax` WHERE
      (`state` = '", state, "') AND
      `p` = '", p, "'")
      println(q)
    table = SQLite.query(db, q)
    # get rid of the nullables introduced by SQLite,
    # as they cannot be interpreted by the plotting
    # software
    for c in 1:ncol(table)
      table[c] = table[c].values
    end
    PPTdf = copy(table)

    q = string("SELECT * FROM `pptRelax` WHERE
      (`state` = '1ext ", state, "') AND
      `p` = '", p, "'")
      println(q)
    table = SQLite.query(db, q)
    # get rid of the nullables introduced by SQLite
    for c in 1:ncol(table)
      table[c] = table[c].values
    end
    kextdf = copy(table)

    q = string("SELECT * FROM `pptRelax` WHERE
      (`state` = '1ext ", state, "') AND
      `p` = '", p, "'")
      println(q)
    table = SQLite.query(db, q)
    # get rid of the nullables introduced by SQLite
    for c in 1:ncol(table)
      table[c] = table[c].values
    end
    onlySymdf = copy(table)

    PPTLayer = layer(PPTdf,
      x=:p_succ,
      y=:fidelity,
      Geom.line,
      # Geom.point,
      Theme(default_color=colorant"blue"))

    layers = typeof(PPTLayer)[]
    names = AbstractString[]
    colors = AbstractString[]
    if nrow(PPTdf) > 0 # only include this layer if it is nonempty
      push!(layers, PPTLayer)
      push!(names, "PPT")
      push!(colors, "blue")
    end

    if nrow(kextdf) > 0 && false # only include this layer if it is nonempty
      kextLayer = layer(kextdf,
        x=:p_succ,
        y=:fidelity,
        Geom.point,
        Theme(default_color=colorant"red"))
      push!(layers, kextLayer)
      push!(names, "1-ext")
      push!(colors, "red")
    end

    if nrow(onlySymdf) > 0 # only include this layer if it is nonempty
      onlySymLayer = layer(onlySymdf,
      x=:p_succ,
      y=:fidelity,
      Geom.point,
      Theme(default_color=colorant"green"))
      push!(layers, onlySymLayer)
      push!(names, "1-ext")
      push!(colors, "green")
    end

    # Add for some input states the EPL protocol
      F, p_succ = EPLParam(sortAB(copies(stateFn(p),2),2,2))
      eplLayer = layer(y = [F], x = [p_succ], Geom.point, Theme(default_color=colorant"purple"))
      unshift!(layers, eplLayer)
      unshift!(names, "EPL protocol")
      unshift!(colors, "purple")

    # Include schemes in plots
    F, p_succ = DEJMPSParam(stateFn(p))
    deutschLayer = layer(y = [F], x = [p_succ], Geom.point, Theme(default_color=colorant"cyan"))
    F, p_succ = BBPSSWParam(stateFn(p))
    bennettLayer = layer(y = [F], x = [p_succ], Geom.point, Theme(default_color=colorant"orange"))

    unshift!(layers, deutschLayer)
    unshift!(layers, bennettLayer)
    unshift!(names, "DEJMPS protocol")
    unshift!(names, "BBPSSW protocol")
    unshift!(colors, "cyan")
    unshift!(colors, "orange")

    pl = plot(
      Theme(panel_fill=colorant"white"),
      layers...,
      Stat.xticks(ticks=collect(0:0.1:1)),
      Stat.yticks(ticks=collect(0.5:0.1:1)),
      Guide.xlabel("Probability of success"),
      Guide.ylabel("Fidelity"),
      Guide.title(plottitle(p)),
      Guide.manual_color_key("Legend",
        names,
        colors)
        )

    filename = replace(replace(string(state,p), " ", ""), ".", "")

    draw(PNG(string("plots/", filename ,".png"), 600px, 300px), pl)
    draw(PDF(string("pdfplots/", filename ,".pdf"), 4inch, 3inch), pl)
  end
end


text:
@testset "cv_iht normal" begin
	# Since my code seems to work, putting in some output as they can be verified by comparing with simulation

    #simulat data with k true predictors, from distribution d and with link l.
	n = 1000
	p = 10000
	k = 10
	d = Normal
	l = canonicallink(d())

	#set random seed
	Random.seed!(2019)

	#construct snpmatrix, covariate files, and true model b
	x = simulate_random_snparray(n, p, undef)
	xbm = SnpBitMatrix{Float64}(x, model=ADDITIVE_MODEL, center=true, scale=true);
	z = ones(n, 1) # the intercept

	# simulate response, true model b, and the correct non-0 positions of b
	y, true_b, correct_position = simulate_random_response(x, xbm, k, d, l)

	#specify path and folds
	path = collect(1:20)
	num_folds = 3
	folds = rand(1:num_folds, size(x, 1))

	# run threaded IHT
	mses = cv_iht(d(), l, x, z, y, 1, path, num_folds, folds=folds, init=false, use_maf=false, debias=true, parallel=true)

	@test argmin(mses) == 11
	@test all(mses ≈ [1927.0765190526674;1443.8787350418297;1080.0413162424918;
		898.2611345915562;700.2502074095313;507.3948750413776;391.9679112461845;
		381.68778672828626;368.53087109541184;359.9090161930197;344.3626513420314;
		353.22416803442127;357.75754456737286;363.5298542018518;379.4998479342846;
		365.9493360314551;387.4562377561079;386.80694552347995;398.9880370887734;
		403.1125856821553])
end

@testset "cv_iht bernoulli" begin
	# Since my code seems to work, putting in some output as they can be verified by comparing with simulation

    #simulat data with k true predictors, from distribution d and with link l.
	n = 1000
	p = 10000
	k = 10
	d = Bernoulli
	l = canonicallink(d())

	#set random seed
	Random.seed!(2019)

	#construct snpmatrix, covariate files, and true model b
	x = simulate_random_snparray(n, p, undef)
	xbm = SnpBitMatrix{Float64}(x, model=ADDITIVE_MODEL, center=true, scale=true);
	z = ones(n, 1) # the intercept

	# simulate response, true model b, and the correct non-0 positions of b
	y, true_b, correct_position = simulate_random_response(x, xbm, k, d, l)

	#specify path and folds
	path = collect(1:20)
	num_folds = 3
	folds = rand(1:num_folds, size(x, 1))

	# run threaded IHT
	mses = cv_iht(d(), l, x, z, y, 1, path, num_folds, folds=folds, init=false, use_maf=false, debias=true, parallel=true)
	@test argmin(mses) == 6

	#not sure why mses[end] can change..... is this a bug?
	@test all(mses[1:end-1] .≈ [433.1645990291242; 396.4218354530673; 354.03582998371684; 337.2781542152179; 321.40856422243473;
			 300.13602523499037; 300.7545574376133; 304.0569785536426; 313.5244382980213; 315.21199531008796;
			 328.7337680232285; 355.3073847457196; 360.42023935617004; 354.7111928546249; 374.7397814724943;
			 380.5160395899203; 404.44221489146645; 404.77237379170714;412.59479205229724])
end

@testset "cv_iht Poisson" begin
	# Since my code seems to work, putting in some output as they can be verified by comparing with simulation

    #simulat data with k true predictors, from distribution d and with link l.
	n = 1000
	p = 10000
	k = 10
	d = Poisson
	l = canonicallink(d())

	#set random seed
	Random.seed!(2019)

	#construct snpmatrix, covariate files, and true model b
	x = simulate_random_snparray(n, p, undef)
	xbm = SnpBitMatrix{Float64}(x, model=ADDITIVE_MODEL, center=true, scale=true);
	z = ones(n, 1) # the intercept

	# simulate response, true model b, and the correct non-0 positions of b
	y, true_b, correct_position = simulate_random_response(x, xbm, k, d, l)

	#specify path and folds
	path = collect(1:20)
	num_folds = 3
	folds = rand(1:num_folds, size(x, 1))

	# run threaded IHT
	mses = cv_iht(d(), l, x, z, y, 1, path, num_folds, folds=folds, init=false, use_maf=false, debias=true, parallel=true)

	@test argmin(mses) == 7
	@test all(mses ≈ [ 612.6766519782327; 554.4890824202639; 522.4053768900791;
		485.3207425983446; 445.78015443806163;412.2320772688479; 382.6698340685961;
		385.9384838080234;394.81185628720675; 405.7743050802736; 406.27235883258606;
		413.3816413610929; 423.313949398366; 425.98190340196277;427.5153001291356;
		420.3871855662542; 440.9354230778697;445.96512112259165; 455.5151967326286;
		437.7447856125584])
end

@testset "cv_iht NegativeBinomial" begin
	# Since my code seems to work, putting in some output as they can be verified by comparing with simulation

    #simulat data with k true predictors, from distribution d and with link l.
	n = 1000
	p = 10000
	k = 10
	d = NegativeBinomial
	l = LogLink()

	#set random seed
	Random.seed!(2019)

	#construct snpmatrix, covariate files, and true model b
	x = simulate_random_snparray(n, p, undef)
	xbm = SnpBitMatrix{Float64}(x, model=ADDITIVE_MODEL, center=true, scale=true);
	z = ones(n, 1) # the intercept

	# simulate response, true model b, and the correct non-0 positions of b
	y, true_b, correct_position = simulate_random_response(x, xbm, k, d, l)

	#specify path and folds
	path = collect(1:20)
	num_folds = 3
	folds = rand(1:num_folds, size(x, 1))

	# run threaded IHT
	mses = cv_iht(d(), l, x, z, y, 1, path, num_folds, folds=folds, init=false, use_maf=false, debias=true, parallel=true)

	@test argmin(mses) == 6
	@test all(mses ≈ [ 299.29083418540984; 277.1588834992124; 254.24543899850323;
		239.65111105625778; 234.28391512742712; 226.11206626366229; 231.00815753407852;
		 230.28151617161788; 234.06225621624264; 235.33924496682374; 237.76185304461012;
		 241.4319393553253; 247.51432485217947; 249.26308636719477; 248.54094497173412;
		 253.51544251902232; 259.8047937477371; 253.51719798579745; 263.6992561532014;
		 263.43746844493137])
end

@testset "iht_run_many_models normal" begin
	n = 1000
	p = 10000
	k = 10
	d = Normal
	l = canonicallink(d())

	#set random seed
	Random.seed!(2019)

	#construct snpmatrix, covariate files, and true model b
	x = simulate_random_snparray(n, p, undef)
	xbm = SnpBitMatrix{Float64}(x, model=ADDITIVE_MODEL, center=true, scale=true);
	z = ones(n, 1) # the intercept

	# simulate response, true model b, and the correct non-0 positions of b
	y, true_b, correct_position = simulate_random_response(x, xbm, k, d, l)

	#specify path and run a ton of models
	path = collect(1:20)
	result = iht_run_many_models(d(), l, x, z, y, 1, path, parallel=true);

	@test result[1] ≈ -2272.788840069668
	@test result[2] ≈ -2144.4519706525098
	@test result[3] ≈ -1998.9118082179523
	@test result[4] ≈ -1888.362338848445
	@test result[5] ≈ -1767.649126874247
	@test result[6] ≈ -1617.2385746806704
	@test result[7] ≈ -1484.2709104112216
	@test result[8] ≈ -1447.9560237353137
	@test result[9] ≈ -1421.5794827334514
	@test result[10] ≈ -1406.8807653835697
	@test result[11] ≈ -1396.8628294699147
	@test result[12] ≈ -1390.1763942804469
	@test result[13] ≈ -1383.5471518926422
	@test result[14] ≈ -1376.7041041059629
	@test result[15] ≈ -1371.070140615097
	@test result[16] ≈ -1363.7164318022296
	@test result[17] ≈ -1358.8454592637506
	@test result[18] ≈ -1353.106442197983
	@test result[19] ≈ -1351.709447505716
	@test result[20] ≈ -1343.1979015928894
end


text:
#############################################################################
# Joulia
# A Large-Scale Spatial Power System Model for Julia
# See https://github.com/JuliaEnergy/Joulia.jl
#############################################################################
# This file is the main package file

module Joulia

using DataFrames
using JuMP
using ProgressMeter

# utility functions
include("util.jl")

# power plants
include("powerplants.jl")

# nodes
include("nodes.jl")

# lines
include("lines.jl")

# renewables
include("res.jl")

# storages
include("storages.jl")

# model
include("model.jl")

export
    PowerPlants,
    Lines,
    Nodes,
    RenewableEnergySource,
    Storages,
    week_slices,
    JouliaModel,
    run_model

end # module


text:
using DocStringExtensions
using ProgressMeter
using Gadfly
using Cairo
using RoMEPlotting
Gadfly.set_default_plot_size(35cm,25cm)

# from results: 2020-02-23T03:28:10.097
inittime=[
    13.923553870000001
    # 2.0 # 7.111387647 # precompiling
    # 8.639480731
    # 0.762567549
    # 2.52150559
    14.531640172
    # 0.847736082
    # 0.783871827999999
    # 9.711954108
    # 3.188078154
    49.340871044
    33.347981764
    24.992062225
    22.201434791
    37.045889202
    35.132906422
    33.124455274
    25.782081818
    21.445426782
    38.041981655
    25.667662991
    24.509122983
    18.96485348
    35.120005101
    22.891443964
    19.393177096
    29.566993705
    22.562560487
    20.464169769
    18.329940968
    25.15793886
    19.126782885
    18.384817736
    34.474076872
    25.536059316
    18.814305877
    32.098617146
    23.725097425
    24.661727613
    23.529066873
    20.708167353
    41.850563097
    30.394270757
    27.53120688
    18.184278245
    24.561897334
    28.050592505
    22.277458803
    21.123530868
    25.581882865
    25.839933464
    24.371785428
    28.449811854
    19.297024484
    27.872869139
    17.548716492
    27.280586297
    23.6224749
    19.06120724
    26.650180153
    21.94149314
    25.836919467
    19.094412425
    33.762299562
    38.326195492
    32.017226393
    24.055516607
    30.605293042
    29.898421198
    19.657790709
    25.358863355
    19.702823874
    19.928584292
    25.222084299
    21.744061298
    22.393525834
    46.48569922
    43.689342378
    23.920386841
    40.68907699
    25.069462586
    44.516693615
    21.730817456
    23.7343984
    28.511252744
    28.729642925
    28.958604325
    19.112645443
    35.595642871
    23.72370413
    30.154412078
    27.058608226
    26.685286617
    21.138100832
    21.319813642
    41.143639249
    24.918579089
    40.665687492
    31.076767095
    31.168699379
    31.773255372
    24.11949409
    23.063612059
    24.920672172
    24.633388553
    30.933378561
    30.626149481
    34.718042871
    22.036652509
    19.681987354
    40.114025408
    23.966545466
    21.718494469
    29.748137542
    20.266077557
    27.776297353
    17.636508099
    19.340242204
    20.719850124
    18.43781753
    21.304971816
    19.582438121
    20.097249032
    19.062781597
    19.587003153
    19.031764947
    18.837652548
    19.552649898
    21.096370845
    19.895401358
    18.979098192
    19.719184593
    19.657576037
    19.97702779
    19.884674546
    20.388048656
    24.715450269
    25.178573931
    27.035220654
    21.004084076
    20.69141639
    23.313526345
    21.554464899
    22.589946352
    29.821371524
    26.30997985
    25.350900452
];

solvetime = [
    # 0
    # 0
    # 0
    # 0
    0
    # 0
    # 0
    # 0
    70.0 # 223.298124668 # precompiling is excessive
    127.927759554
    125.413819048
    148.515074234
    114.643111867
    130.067414752
    116.591257116
    125.92697153
    146.91522901
    124.459727457
    112.015100314
    119.873414991
    139.397963179
    125.617490887
    112.820977695
    129.391103947
    123.26573995
    110.883475493
    129.832680594
    136.301901148
    119.74783935
    152.537798002
    126.58402604
    142.098410369
    117.710888204
    146.374041294
    141.393480826
    137.785630885
    154.962325813
    128.055167705
    162.495087517
    120.782532872
    148.97363259
    144.861907273
    165.769649418
    162.056390594
    138.407873697
    140.74692623
    152.207775846
    155.695375416
    152.78733244
    154.058805384
    135.067113018
    142.863716583
    150.918709039
    146.376801785
    156.715065123
    139.197815332
    167.873990906
    146.304335165
    154.659913943
    141.669132261
    157.328337244
    157.663703422
    158.684651059
    161.954268784
    139.224675678
    157.836792744
    139.215477005
    160.034403365
    196.363458843
    177.526854252
    192.500806936
    210.754950187
    186.562554767
    174.188670822
    169.375834258
    171.027302979
    213.349333804
    158.431002099
    183.419584683
    167.103482839
    154.450013071
    166.855410513
    167.144473576
    156.904730049
    177.913141792
    162.445496695
    161.466968961
    159.987589656
    174.016696833
    178.093277683
    153.027023061
    181.590475774
    189.112277624
    186.886697627
    169.787662128
    175.810278451
    143.744676215
    193.088276639
    191.852113533
    184.75878549
    203.878753874
    196.998695861
    183.605287835
    163.128652931
    170.348446623
    190.923487248
    181.001078298
    211.449152831
    201.602216595
    168.809270956
    200.544647688
    210.616895848
    198.639793471
    245.181534504
    179.499840855
    120.849533366
    123.318358076
    118.590452535
    142.919545044
    125.402021766
    125.291528595
    121.632521557
    122.657152356
    128.260106156
    134.983011712
    153.491732159
    122.839969837
    133.198374417
    147.872250012
    143.670244034
    131.515225824
    132.972148121
    130.177216643
    139.176283859
    145.284707366
    141.376052209
    160.664398218
    158.65330199
    135.723196487
    137.019370905
    135.014465904
    149.970985041
    152.706946397
    157.224246642
    139.370237378
    153.155698731
];

pl = Gadfly.plot(
    Gadfly.layer(x=1:10:1390,y=inittime, Geom.line, Theme(default_color=colorant"red",line_width=2pt)),
    Gadfly.layer(x=1:10:1390,y=solvetime, Geom.line, Theme(default_color=colorant"blue",line_width=2pt)),
    Guide.ylabel("time [s]"),
    Guide.xlabel("# poses"),
)

Gadfly.draw(pl, PDF(joinpath("/tmp/caesar/2020-02-23T03:28:10.097", "slamTime.pdf"),15cm,10cm)

pl |> PDF("slamTime.pdf", 12cm, 8cm)

#


text:
# SPD-License-Identifier: MIT

using ArgParse
using HDF5
using JLD2
using Knet
using MarkovModels
using TOML

function main(args)
    indim = h5open(args["trainfea"], "r") do f
        key = iterate(keys(f))[1]
        size(read(f, key), 1)
    end

    pdfid_mapping = load(args["hmms"])["pdfid_mapping"]
    outdim = maximum(values(pdfid_mapping))

    model = buildmodel(TOML.parsefile(args["config"]), indim, outdim)
    save(args["out"], Dict("model" => model))
end

function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table s begin
        "modelfile"
            required = true
            help = "input model file"
        "config"
            required = true
            help = "TOML configuration file"
        "trainfea"
            required = true
            help = "training features"
        "hmms"
            required = true
            help = "hmms saved in JLD2 format"
        "out"
            required = true
            help = "output model"
    end
    s.description = """
    Create a randomly initialized model.
    """
    parse_args(s)
end

args = parse_commandline()
include(args["modelfile"])
main(args)

max_stars_repo_name: MartinKocour/SpeechLab

max_issues_repo_name: MartinKocour/SpeechLab

max_forks_repo_name: MartinKocour/SpeechLab


text:
"""

"""
function average_precipitation_per_pw_bin_dayang(pw,precipitation,max_precip_bin,binspacing)
    precip_bin = Float64[]
    pw_bins = 0:binspacing:max_precip_bin
    n_pw_interval = length(pw_bins)
    for i in 1:n_pw_interval
        pw_i_low = pw_bins[1] + binspacing*(i-1)
        pw_i_up  = pw_bins[1] + binspacing*i
        ind_i    = findall(x -> (x>pw_i_low)&(x<pw_i_up),pw)
        if isempty(ind_i)
            push!(precip_bin,0.0)
        else
            push!(precip_bin,mean(precipitation[ind_i]))
        end
    end
    return precip_bin
end

"""
    average_precipitation_per_pw_bin(pw,precipitation,bins,binspacing)
Computes the mean precipitation as a function of the binned precipitable water following Yang, D., 2018: Boundary Layer Height and Buoyancy Determine the Horizontal Scale of Convective Self-Aggregation. J. Atmos. Sci., 75, 469–478, https://doi.org/10.1175/JAS-D-17-0150.1.
It will be happy if bins is quite large (example from 0 to 300 mm).
"""
function average_precipitation_per_pw_bin(pw,precipitation,pw_bins,binspacing)
    T = eltype(pw)
    tone = one(T)
    tzero = zero(T)
    probability_of_being_in_bin = zeros(T,length(pw_bins))
    average_precipitation_per_bin = zeros(T,length(pw_bins))
    number_of_columns = reduce(*,size(pw))
    factor = number_of_columns * binspacing
    @inbounds @simd for column in CartesianIndices(pw)
        binindex = ceil(Int,pw[column]/binspacing)
        average_precipitation_per_bin[binindex] += precipitation[column]
        probability_of_being_in_bin[binindex] += tone
    end
    @tturbo for ind in eachindex(probability_of_being_in_bin)
        divideby = ifelse(probability_of_being_in_bin[ind] == tzero,tone,probability_of_being_in_bin[ind])
        average_precipitation_per_bin[ind] /= divideby
    end
    @tturbo probability_of_being_in_bin ./= factor
    return probability_of_being_in_bin,average_precipitation_per_bin
end

"""
    velocity_topolar(u,v,index,center)
Take a velocity vector, an origin of said vector and a center and return the tangential and azimuthal velocities with respect to that center.
"""
function velocity_topolar(u,v,index,center)
    pos = index .- center
    theta1 = atan(v,u)
    theta2 = atan(pos[2],pos[1])
    return -hypot(u,v) * cos(theta1 + theta2), hypot(u,v)*sin(theta1 + theta2)
end

"""
    velocity_cartesian_to_polar(u,v,index_of_point = (0,0),index_of_center = (0,0))
    Take matrix of U and V components of velocities as well as a matrix of locations and the location of a center.
    Compute the tangential and radial components of velocity around this center.
"""

function velocity_cartesian_to_polar(u,v,index_of_point = (0,0),index_of_center = (0,0))
    tangential = similar(u)
    radial = similar(u)
    for index in CartesianIndices(u)
        t,r = velocity_cartesian_to_polar(u[index],v[index],Tuple(index)[1:2],index_of_center[1:2])
        radial[index] = r
        tangential[index] = t
    end
    return tangential,radial
end

function velocity_cartesian_to_polar(u :: Number,v :: Number,index_of_point = (0,0),index_of_center = (0,0))
    speed = hypot(u,v)
    x,y = index_of_point .- index_of_center
    angle_of_position_vector = atan(y,x)
    angle_of_velocity_vector = atan(v,u)
    return velocity_polar = speed.*sincos( angle_of_velocity_vector - angle_of_position_vector )
end


text:
# Copyright (c) 2016, Johan Mabille, Sylvain Corlay, Wolf Vollprecht and Martin Renou
# Copyright (c) 2016, QuantStack
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
#
# * Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# * Neither the name of the copyright holder nor the names of its
#   contributors may be used to endorse or promote products derived from
#   this software without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

@inline integer(v::AbstractSIMD{W,Float64}) where {W} =
  _integer(v, has_feature(Val(:x86_64_avx512dq)))
@inline _integer(v::AbstractSIMD{W,Float64}, ::True) where {W} = vconvert(Int64, v)
@inline _integer(v::AbstractSIMD{W,Float64}, ::False) where {W} = vconvert(Int32, v)

@inline function erf_kernel_l9(x::Union{Float64,AbstractSIMD{<:Any,Float64}})
  y = vfmadd(x, 6.49254556481904e-5, 0.00120339380863079)
  z = vfmadd(x, 0.000364915280629351, 0.00849717371168693)
  y = vfmadd(x, y, 0.0403259488531795)
  z = vfmadd(x, z, 0.0869936222615386)
  y = vfmadd(x, y, 0.135894887627278)
  z = vfmadd(x, z, 0.453767041780003)
  y = vfmadd(x, y, 1.12837916709551)
  z = vfmadd(x, z, 1.0)
  res = y / z
  res + vfnmadd(
    8.13035732583580548119176214095147680242299108680658322550212199643608432312984e-16,
    x,
    2.679870718713541577762315771546743935213688171443421284936882991116060314650674e-15,
  )
  # Base.FastMath.div_fast(w * y, z)
end
# erf(x)/x = erf_poly_l9(x*x) # r = x^2; x = sqrt(r)
# erf(sqrt(r))/sqrt(r) = erf_poly_l9(r)
# x in 0...0.65
# r in 0...0.65^2

# using Remez, SpecialFunctions
# erftestsmall(x) = (r = sqrt(x); erf(r)/r)
# ncinit = BigFloat[1.12837916709551, 0.135894887627278, 0.0403259488531795, 0.00120339380863079, 6.49254556481904e-5];
# dcinit = BigFloat[1.0, 0.453767041780003, 0.0869936222615386, 0.00849717371168693, 0.000364915280629351];
# N,D,E,X = ratfn_minimax(erftestsmall, [big"1e-70", big"0.65"^2], 4, 4, (w,x)->BigFloat(1), ncinit, dcinit); @show(E); N

vscalef(b::Bool, x, y, z) = b ? (x * (2^y)) : z
@inline function erf_v56f_kernel(v3f, v4f)
  v29f = v4f * -1.4426950408889634
  v30f = round(v29f)
  v31f = vfmsub(v30f, -0.6931471803691238, v4f)
  v32f = v30f * 1.9082149292705877e-10
  v33f = v31f - v32f
  v34f = v33f * v33f
  v35f = vfmadd(v34f, 4.1381367970572385e-8, -1.6533902205465252e-6)
  v36f = vfmadd(v34f, v35f, 6.613756321437934e-5)
  v37f = vfmadd(v34f, v36f, -0.0027777777777015593)
  v39f = vfnmsub(v34f, v37f, 0.16666666666666602)
  v40f = vfmadd(v34f, v39f, v33f)
  v41f = v40f * v33f
  v42f = 2.0 - v40f
  v43f = Base.FastMath.div_fast(v41f, v42f)
  v44f = 1.0 - v32f
  v45f = v44f + v31f
  v46f = v45f + v43f
  m47 = v4f ≤ 708.3964185322641
  v54f = vscalef(m47, v46f, v30f, zero(v3f))
  # v56f = ifelse(v4f < 709.782712893384, v54f, Inf) # TODO: NaN should return v54f
end
@inline function erf_v97f_kernel(v3f)
  v91f = vfmadd(v3f, 0.0400072964526861, 0.278788439273629)
  v86f = vfmadd(v3f, 0.0225716982919218, 0.157289620742839)
  v92f = vfmadd(v3f, v91f, 1.05074004614827)
  v87f = vfmadd(v3f, v86f, 0.581528574177741)
  v93f = vfmadd(v3f, v92f, 2.38574194785344)
  v88f = vfmadd(v3f, v87f, 1.26739901455873)
  v94f = vfmadd(v3f, v93f, 3.37367334657285)
  v89f = vfmadd(v3f, v88f, 1.62356584489367)
  v95f = vfmadd(v3f, v94f, 2.75143870676376)
  v90f = vfmadd(v3f, v89f, 0.99992114009714)
  v96f = vfmadd(v3f, v95f, 1.0)
  v97f = v90f / v96f
end
@inline function erf_v71_kernel(v3f)
  v65f = vfmadd(v3f, 0.0125304936549413, 0.126579413030178)
  v60f = vfmadd(v3f, 0.00706940843763253, 0.0714193832506776)
  v66f = vfmadd(v3f, v65f, 0.594651311286482)
  v61f = vfmadd(v3f, v60f, 0.331899559578213)
  v67f = vfmadd(v3f, v66f, 1.61876655543871)
  v62f = vfmadd(v3f, v61f, 0.878115804155882)
  v68f = vfmadd(v3f, v67f, 2.65383972869776)
  v63f = vfmadd(v3f, v62f, 1.33154163936765)
  v69f = vfmadd(v3f, v68f, 2.45992070144246)
  v64f = vfmadd(v3f, v63f, 0.999999992049799)
  v70f = vfmadd(v3f, v69f, 1.0)
  v64f, v70f
end
@inline function verf(v0f::Union{Float64,AbstractSIMD{<:Any,Float64}})
  # v1 = reinterpret(UInt64, v)
  # v2 = v1 & 0x7fffffffffffffff
  # v3 = reinterpret(Float64, v2)
  v3f = abs(v0f)
  v4f = v0f * v0f
  # m6 = v3f < 0.65
  # m6 = v3f < 0.68
  m6 = v3f < 0.675
  if vany(collapse_or(m6))
    v19f = v0f * erf_kernel_l9(v4f)
    vall(collapse_and(m6)) && return v19f
  else
    # v19f = zero(v0f)
    v19f = _vundef(v0f)
  end
  m23 = v3f < 2.725
  v56f = erf_v56f_kernel(v3f, v4f)
  if vany(collapse_or(m23 & (~m6))) # any(0.675 < v3f < 2.2) # l58
    v64f, v70f = erf_v71_kernel(v3f)
    v71f = Base.FastMath.div_fast(v56f, v70f)
    v73f = vfnmadd(v71f, v64f, 1.0)
    v78f = copysign(v73f, v0f)
    v84f = ifelse(m6, v19f, v78f)
    vall(collapse_and(m23)) && return v84f
  else
    v84f = v19f
  end
  # l83 # > 2.2
  # v91f = vfmadd(v3f, 0.0400072964526861, 0.278788439273629)
  # v86f = vfmadd(v3f, 0.0225716982919218, 0.157289620742839)
  # v92f = vfmadd(v3f, v91f, 1.05074004614827)
  # v87f = vfmadd(v3f, v86f, 0.581528574177741)
  # v93f = vfmadd(v3f, v92f, 2.38574194785344)
  # v88f = vfmadd(v3f, v87f, 1.26739901455873)
  # v94f = vfmadd(v3f, v93f, 3.37367334657285)
  # v89f = vfmadd(v3f, v88f, 1.62356584489367)
  # v95f = vfmadd(v3f, v94f, 2.75143870676376)
  # v90f = vfmadd(v3f, v89f, 0.99992114009714)
  # v96f = vfmadd(v3f, v95f, 1.0)
  # # v97f = v56f * v90f
  # # v98f = Base.FastMath.div_fast(v97f, v96f)
  # v97f = v90f / v96f
  v97f = erf_v97f_kernel(v3f)
  v99f = vfnmadd(v97f, v56f, 1.0)
  # v99f = sub_ieee(1.0, v98f)
  v104f = copysign(v99f, v0f)
  v105f = ifelse(m23, v84f, v104f)
  return v105f
end
# # erftest(x) = (1 - erf(x)) / VectorizationBase.erf_v56f_kernel(abs(x),x*x)
# # N,D,E,X = ratfn_minimax(erftest, [big"2.2", big"5.9215871960"], 5, 6); @show(E); N
# erf(2.3) = 1.0 - v98f
# (1.0 - erf(2.3) / v56f(2.3)) = v97f
# (1.0 - erf(2.3) / v56f(2.3)) = v90f / v96f
# rational polynomial, degree 5 / degree 6

@inline function verf(v0f::Union{Float32,AbstractSIMD{<:Any,Float32}})
  v3f = abs(v0f)
  m4 = v3f < 0.6666667f0
  v8f = v3f * v3f
  if vany(collapse_or(m4)) # L7
    v9f = vfmadd(v8f, -0.00060952053f0, 0.005013293f0)
    v10f = vfmadd(v8f, v9f, -0.026780106f0)
    v11f = vfmadd(v8f, v10f, 0.11282183f0)
    v12f = vfmadd(v8f, v11f, -0.37612528f0)
    v13f = vfmadd(v8f, v12f, 1.1283791f0)
    v17f = v13f * v0f
    vall(collapse_and(m4)) && return v17f
  else
    v17f = _vundef(v0f)
  end
  v18f = v3f + 1.0f0
  v19f = Base.FastMath.div_fast(v3f, v18f)
  v20f = v19f - 0.4f0
  v23f = -1.442695f0 * v8f
  v24f = round(v23f)
  v25f = vfmsub(v24f, -0.6933594f0, v8f)
  v26f = vfmadd(v24f, 0.00021219444f0, v25f)
  v27f = vfmadd(v26f, 0.0013997796f0, 0.008395563f0)
  v28f = vfmadd(v26f, v27f, 0.0416654f0)
  v31f = v26f * v26f
  v29f = vfmadd(v26f, v28f, 0.16666277f0)
  v30f = vfmadd(v26f, v29f, 0.5f0)
  v32f = vfmadd(v30f, v31f, v26f)
  v33f = v32f + 1.0f0
  m34 = v8f ≤ 88.37626f0
  v42f = vscalef(m34, v33f, v24f, zero(v0f))
  v43f = vfmadd(v20f, -2.6283f0, 6.702909f0)
  v44f = vfmadd(v20f, v43f, -6.4094872f0)
  v45f = vfmadd(v20f, v44f, 3.2607658f0)
  v46f = vfmadd(v20f, v45f, -1.364514f0)
  v47f = vfmadd(v20f, v46f, 0.15627646f0)
  v48f = vfmadd(v20f, v47f, 0.14205085f0)
  v49f = vfmadd(v20f, v48f, 0.38435692f0)
  v50f = vfmadd(v20f, v49f, 0.16037047f0)
  v51f = vfmadd(v20f, v50f, -1.1370356f0)
  v52f = vfmadd(v20f, v51f, 0.5392844f0)
  v53f = v52f * v42f
  v54f = 1.0f0 - v53f
  m55 = v0f < zero(v0f)
  v56f = -v54f
  v57f = ifelse(m55, v56f, v54f)
  v58f = ifelse(m4, v17f, v57f)
  m59 = v3f ≠ Inf32
  m60 = v0f > zero(v0f)
  v61f = ifelse(m60, one(v0f), zero(v0f))
  v62f = ifelse(m55, one(v0f), zero(v0f))
  v63f = v61f - v62f
  v65f = ifelse(v0f == v0f, v63f, NaN32)
  v66f = ifelse(m59, v58f, v65f)
  return v66f
end


text:
export @maybethread

macro maybethread(loop)
    if Threads.nthreads()>1
      quote Threads.@threads $(Expr(loop.head,
                               Expr(loop.args[1].head, esc.(loop.args[1].args)...),
                               esc(loop.args[2]))); end
    else
      # @warn "running single threaded"
      quote $(esc(loop)); end
    end
end


text:
using Luxor

function dump(fig)
    if fig.surfacetype == :png
        Cairo.write_to_png(fig.surface, fig.buffer)
    end

    Cairo.finish(fig.surface)
    Cairo.destroy(fig.surface)

    if fig.filename != ""
        write(fig.filename, fig.bufferdata)
    end
    true
end

function CNOT()
    d = Drawing(400, 400, :png, "cnot")
    origin()

    rect()
    dump(d)
end

@png begin
    fontsize(50)
    circle(Point(0, 0), 150, :stroke)
    text("hello world", halign=:center, valign=:middle)
end

max_stars_repo_name: Roger-luo/QCBM

max_issues_repo_name: Roger-luo/QCBM

max_forks_repo_name: Roger-luo/QCBM


text:
module kNNRecommenders

using CSV
using Random
using SparseArrays

export MovieLens, load_data, split_data, transform_to_rating_matrix
export kNNRecommender, compute_similarity_matrix, predict, evaluate_recall

include("MovieLens.jl")
include("ItemkNN.jl")

end # module


text:
using ITensors

examples_dir = joinpath(ITensors.examples_dir(), "src")
include(joinpath(examples_dir, "trg.jl"))
include(joinpath(examples_dir, "2d_classical_ising.jl"))

function run(; maxdim::Int,
               nsweeps::Int = 20,
               outputlevel::Int = 0,
               cutoff::Float64 = 0.0,
               β::Float64 = 1.001 * βc)
  # Make Ising model MPO
  s = Index(2)
  sₕ = addtags(s, "horiz")
  sᵥ = addtags(s, "vert")
  T = ising_mpo(sₕ, sᵥ, β)
  κ, T = trg(T; χmax = maxdim, cutoff = cutoff, nsteps = nsweeps)
  if outputlevel > 0
    @show nsweeps
    @show cutoff
    @show κ, exp(-β * ising_free_energy(β))
  end
  return ITensors.maxdim(T)
end


text:
################################################################################
#
#     Clgrp.jl : Class group computation of maximal orders in number fields
#
# This file is part of hecke.
#
# Copyright (c) 2015: Claus Fieker, Tommy Hofmann
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# * Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
#
# (C) 2015, 2016 Claus Fieker
#
################################################################################
#
# Todo:
#  - make sure the precision for LLL is high enough (by checking that the
#    resulting elements have a reasonable norm/ length? theory?)
#    done
#  - add reasonable verbose printing
#    done
#  - write hnf from upper_triangular
#  - understand/use profiling information (memory, ...)
#
#  - need different norm function: modular resultant? with a large known
#    factor AND without any proof...
#    otherwise, it takes much too long if the ideal is non-trivial
#  DONE (norm_div)
#
#  - move the various factor, issmooth and similar to a sensible
#    spot. This has nothing to do with class groups
#  - the SingleP version:
#      figure out if a union is better than the current version
#      ie have a version for easy primes
#         and a dumb version as fallback
#      make sure stuff works with fractionals as well!
#      just scaling by the den can affect smoothness (incomplete factor base)
#
#
# Clean up:
#  - write show functions
#  - export
#
# Note: enumerating x,0,0,0 is pointless unless x==1
#
################################################################################

export class_group, FactorBase, issmooth, factor, lll_basis,
       unit_group_fac_elem, unit_group

add_verbose_scope(:ClassGroup)
add_verbose_scope(:ClassGroup_time)
add_verbose_scope(:ClassGroup_gc)

add_assert_scope(:ClassGroup)
set_assert_level(:ClassGroup, 0)
set_assert_level(:LatEnum, 0)

include("Clgp/Ctx.jl")
include("Clgp/FacBase_Euc.jl")
include("Clgp/FacBase_Idl.jl")
include("Clgp/Main_enum.jl")
include("Clgp/Map.jl")
include("Clgp/Proof.jl")
include("Clgp/Rel_add.jl")
include("Clgp/Rel_enum.jl")
include("Clgp/Rel_LLL.jl")
include("Clgp/Main_LLL.jl")
include("Clgp/Rel_Schmettow.jl")

################################################################################
#
#  Main function
#
################################################################################

function class_group_ctx(O::NfOrd; bound::Int = -1, method::Int = 3, large::Int = 1000, redo::Bool = false)

  if !redo
    try
      c = _get_ClassGrpCtx_of_order(O)::ClassGrpCtx{SMat{fmpz}}
      return c
    end
  end

  if bound == -1
    bound = Int(ceil(log(abs(discriminant(O)))^2*0.3))
    (bound == 0) && (bound = 1)
  end

  c = class_group_init(O, bound, complete = false)::ClassGrpCtx{SMat{fmpz}}

  _set_ClassGrpCtx_of_order(O, c)

  c.B2 = bound * large

  if method == 2
    class_group_find_relations2(c)
  else
    d = root(abs(discriminant(O)), 2)
    c.expect = class_group_expected(c, 100)
    class_group_via_lll(c)
  end

  return c
end

################################################################################
#
#  Verification
#
################################################################################

# think of a sensible function name

function _validate_class_unit_group(c::ClassGrpCtx, U::UnitGrpCtx)
  @vprint :UnitGroup 1 "Validating unit group and class group ... \n"
  O = U.order

  # The residue of the zeta function cannot be computed for degree 1 (K = Q),
  # so we shortcircuit it.

  if degree(O) == 1
    if c.h == 1 && U.tentative_regulator == 1
      return fmpz(1)
    else
      error("Something odd for K = Q")
    end
  end

  if !isdefined(U, :torsion_units)
    @vprint :UnitGroup 1 "Computing torsion structure ... \n"
    U.torsion_units = torsion_units(O)
    U.torsion_units_order = length(U.torsion_units)
    U.torsion_units_gen = torsion_units_gen(O)
  end

  w = U.torsion_units_order

  r1, r2 = signature(O)

  if !isdefined(U, :residue)
    @vprint :UnitGroup 1 "Computing residue of Dedekind zeta function ... \n"
    U.residue = zeta_log_residue(O, 0.6931/2)  #log(2)/2
  end
  residue = U.residue

  pre = prec(parent(residue))

  Ar = ArbField(pre, false)

  loghRtrue = Ar(residue) + log(Ar(w)*sqrt(abs(Ar(discriminant(O))))//(Ar(2)^(r1+r2) * const_pi(Ar)^r2))

  # I should assert that radius(loghRtrue) < log(2)

  @assert isfinite(loghRtrue)

  while true
    loghRapprox = log(c.h* abs(U.tentative_regulator))

    @assert isfinite(loghRapprox)

    if contains(loghRtrue, loghRapprox)
      return fmpz(1)
    elseif !overlaps(loghRtrue, loghRapprox)
      e = exp(loghRapprox - loghRtrue)
      return abs_upper_bound(e, fmpz)
    end

    error("Not yet implemented")
  end
end

function _class_unit_group(O::NfOrd; bound::Int = -1, method::Int = 3, large::Int = 1000, redo::Bool = false, unit_method::Int = 1)

  @vprint :UnitGroup 1 "Computing tentative class and unit group ... \n"

  @v_do :UnitGroup 1 pushindent()
  c = class_group_ctx(O, bound = bound, method = method, large = large, redo = redo)
  @v_do :UnitGroup 1 popindent()

  if c.finished
    U = _get_UnitGrpCtx_of_order(O)
    @assert U.finished
    @vprint :UnitGroup 1 "... done (retrieved).\n"
    return c, U, 1
  end

  @vprint :UnitGroup 1 "Tentative class number is now $(c.h)\n"

  U = UnitGrpCtx{FacElem{nf_elem, AnticNumberField}}(O)

  need_more = true

  hnftime = 0.0

  do_units = true
  while true
    @v_do :UnitGroup 1 pushindent()
    if do_units
      if unit_method == 1
        @vtime_add_elapsed :UnitGroup 1 c :unit_time r = _unit_group_find_units(U, c)
      else
        @vtime_add_elapsed :UnitGroup 1 c :unit_hnf_time module_trafo_assure(c.M)
        @vtime_add_elapsed :UnitGroup 1 c :unit_time r = _unit_group_find_units_with_trafo(U, c)
      end

      @v_do :UnitGroup 1 popindent()
    end
    if r == 1  # use saturation!!!!
      if _validate_class_unit_group(c, U) == 1
        break
      end
    end
    #TODO: use LLL?
    if need_more
      d = root(abs(discriminant(O)), 2)
      c.expect = class_group_expected(d, degree(O), Int(norm(c.FB.ideals[1])), 100)
      need_more = false
    end
    class_group_new_relations_via_lll(c, extra = unit_rank(O) - length(U.units) +1)
    h_old = c.h
    class_group_get_pivot_info(c)
    if h_old == c.h
      do_units = true
    else
      do_units = false
    end
  end
  @assert U.full_rank
  _set_UnitGrpCtx_of_order(O, U)

  c.finished = true
  U.finished = true

  @vprint :ClassGroup 1 "hnftime $(c.time[:hnf_time])\n"

  return c, U, _validate_class_unit_group(c, U)
end

function unit_group_ctx(c::ClassGrpCtx; redo::Bool = false)
  O = order(c.FB.ideals[1])
  if !redo
    try
      U = _get_UnitGrpCtx_of_order(O)::UnitGrpCtx
      return U
    end
  end

  U = UnitGrpCtx{FacElem{nf_elem, AnticNumberField}}(O)
  need_more = true
  while true
    r = _unit_group_find_units(U, c)
    if r == 0
      if need_more
        d = root(abs(discriminant(O)), 2)
        c.expect = class_group_expected(d, degree(O), Int(norm(c.FB.ideals[1])), 100)
        need_more = false
      end
      class_group_new_relations_via_lll(c, E)
    else
      _set_UnitGrpCtx_of_order(O, U)
      return U
    end
  end
end

function unit_group(c::ClassGrpCtx)
  U = unit_group_ctx(c)
  return unit_group(c, U)
end

function unit_group(c::ClassGrpCtx, U::UnitGrpCtx)
  O = order(c.FB.ideals[1])
  K = nf(O)
  U, mU = unit_group_fac_elem(c, U)

  r = MapUnitGrp{typeof(U), typeof(O)}()
  r.header = Hecke.MapHeader(U, O,
    x->O(evaluate(image(mU, x))),
    x->preimage(mU, FacElem([K(x)], fmpz[1])))
  return U, r
end

doc"""
***
    class_group(O::NfOrd; bound = -1, method = 3, redo = false, large = 1000) -> GrpAbFinGen, Map

> Returns a group $A$ and a map $f$ from $A$ to the set of ideals of $O$.
> The inverse of the map is the projection onto the group of ideals modulo the
> group of principal ideals.
> \texttt{redo} allows to trigger a re-computation, thus avoiding the cache.
> \texttt{bound}, when given, is the bound for the factor base.
"""
function class_group(O::NfOrd; bound::Int = -1, method::Int = 3, redo::Bool = false, unit_method::Int = 1, large::Int = 1000)
  c, U, b = _class_unit_group(O, bound = bound, method = method, redo = redo, unit_method = unit_method, large = large)
  @assert b==1
  return class_group(c)
end

doc"""
***
    unit_group(O::NfOrd) -> GrpAbFinGen, Map

> Returns a group $U$ and an isomorphism map $f \colon U \to \mathcal O^\times$.
> A set of fundamental units of $\mathcal O$ can be
> obtained via `[ f(U[1+i]) for i in 1:unit_rank(O) ]`.
> `f(U[1])` will give a generator for the torsion subgroup.
"""
function unit_group(O::NfOrd; method::Int = 3, unit_method::Int = 1)
  c, U, b = _class_unit_group(O, method = method, unit_method = unit_method)
  @assert b==1
  return unit_group(c, U)
end

doc"""
***
    unit_group_fac_elem(O::NfOrd) -> GrpAbFinGen, Map

> Returns a group $U$ and an isomorphism map $f \colon U \to \mathcal O^\times$.
> A set of fundamental units of $\mathcal O$ can be
> obtained via `[ f(U[1+i]) for i in 1:unit_rank(O) ]`.
> `f(U[1])` will give a generator for the torsion subgroup.
> All elements will be returned in factored form.
"""
function unit_group_fac_elem(O::NfOrd; method::Int = 3, unit_method::Int = 1)
  c, U, b = _class_unit_group(O, method = method, unit_method = unit_method)
  @assert b==1
  return unit_group_fac_elem(c, U)
end


text:
#!/usr/bin/env julia
# Tim Sterne-Weiler 2015

using Pkg

const dir = abspath( splitdir(@__FILE__)[1] )
const ver = readline(open(dir * "/VERSION"))

start = time_ns()
println( stderr, "Whippet $ver loading... " )

Pkg.activate(dir * "/..")

using Whippet
using Serialization
using ArgParse

function parse_cmd()
  s = ArgParseSettings()
  # TODO finish options...
  @add_arg_table s begin
    "filename.fastq[.gz]"
      arg_type = String
      required = true
    "paired_mate.fastq[.gz]"
      arg_type = String
      required = false
    "--index", "-x"
      help     = "Output prefix for saving index 'dir/prefix' (default Whippet/index/graph[.jls])"
      arg_type = String
      default  = fixpath( "$dir/../index/graph" )
    "--out", "-o"
      help     = "Where should the gzipped output go 'dir/prefix'?"
      arg_type = String
      default  = fixpath( "./output" )
    "--seed-len", "-L"
      help     = "Seed length"
      arg_type = Int
      default  = 18
    "--seed-try", "-M"
      help     = "Number of failed seeds to try before giving up"
      arg_type = Int
      default  = 3
    "--seed-tol", "-T"
      help     = "Number of seed hits to tolerate"
      arg_type = Int
      default  = 4
    "--seed-buf", "-B"
      help     = "Ignore this many bases from beginning and end of read for seed"
      arg_type = Int
      default  = 5
    "--seed-inc", "-I"
      help     = "Number of bases to increment seed each iteration"
      arg_type = Int
      default  = 18
    "--pair-range", "-P"
      help     = "Seeds for paired end reads must match within _ bases of one another"
      arg_type = Int
      default  = 5000
    "--mismatches", "-X"
      help     = "Allowable number of mismatches in alignment (counted as 1-10^(-phred/10))"
      arg_type = Int
      default  = 3
    "--score-min", "-S"
      help     = "Minimum percent matching (matches - mismatches) / read_length"
      arg_type = Float64
      default  = 0.7
    "--sam"
      help     = "Should SAM format be sent to stdout?"
      action   = :store_true
    "--biascorrect"
      help     = "Apply fragment GC-content and 5' Sequence bias correction methods for more stable PSI values"
      action   = :store_true
    "--stranded"
      help     = "Is the data strand specific in fwd orientation? If so, increase speed with this flag"
      action   = :store_true
    "--pair-same-strand"
      help     = "Whippet by default tries to align fwd/rev pairs, if your data is fwd/fwd or rev/rev set this flag"
      action   = :store_true
    "--phred-33"
      help     = "Qual string is encoded in Phred+33 integers (default)"
      action   = :store_true
    "--phred-64"
      help     = "Qual string is encoded in Phred+64 integers"
      action   = :store_true
    "--circ"
      help     = "Allow back/circular splicing, this will allow output of `BS`-type lines"
      action   = :store_true
    "--force-gz"
      help     = "Regardless of suffix, consider read input as gzipped"
      action   = :store_true
  end
  return parse_args(s)
end

function main()

   args = parse_cmd()

   println(stderr, " $( round( (time_ns()-start)/1e9, digits=6 ) ) seconds." )

   indexpath = fixpath( args["index"] )
   indexname = hasextension(indexpath, ".jls") ? indexpath : indexpath * ".jls"
   println(stderr, "Loading splice graph index... $indexname")
   @timer lib = open(deserialize, indexname)

   ispaired = args["paired_mate.fastq[.gz]"] != nothing

    args["filename.fastq[.gz]"] = fixpath(args["filename.fastq[.gz]"])
    if ispaired
       args["paired_mate.fastq[.gz]"] = fixpath(args["paired_mate.fastq[.gz]"])
    end

   ContainerType = ispaired ? SGAlignPaired : SGAlignSingle

   if args["biascorrect"]
      CounterType   = JointBiasCounter
      BiasModelType = JointBiasMod
   else
      CounterType   = DefaultCounter
      BiasModelType = DefaultBiasMod
   end

   param = AlignParam( args, ispaired, kmer=lib.kmer )
   quant = GraphLibQuant{ContainerType,CounterType}( lib )
   multi = MultiMapping{ContainerType,CounterType}()
   mod   = BiasModelType()

   enc_offset = args["phred-64"] ? 64 : 33

   parser = make_fqparser( args["filename.fastq[.gz]"], forcegzip=args["force-gz"] )

   if ispaired
      mate_parser = make_fqparser( args["paired_mate.fastq[.gz]"], forcegzip=args["force-gz"] )
   end

   #TODO: first implementation was too slow, ie too much communication overhead
   # try speeding up by using Distributed.nprocs > 1
   println(stderr, "Processing reads from file...")
    if ispaired
       println(stderr, "FASTQ_1: " * args["filename.fastq[.gz]"])
       println(stderr, "FASTQ_2: " * args["paired_mate.fastq[.gz]"])
       @timer mapped,totreads,readlen = process_paired_reads!( parser, mate_parser, param, lib, quant, multi, mod,
                                                           sam=args["sam"], qualoffset=enc_offset )
       readlen = Int(floor(readlen))
       println(stderr, "Finished mapping $mapped paired-end reads of length $readlen each out of a total of $totreads mate-pairs...")
    else
       println(stderr, "FASTQ: " * args["filename.fastq[.gz]"])
       @timer mapped,totreads,readlen = process_reads!( parser, param, lib, quant, multi, mod,
                                                     sam=args["sam"], qualoffset=enc_offset )
       readlen = Int(floor(readlen))
       println(stderr, "Finished mapping $mapped single-end reads of length $readlen out of a total of $totreads reads...")
    end

   # TPM_EM
   println(stderr, "Calculating expression values and MLE of equivalence classes with EM:")
   primer_normalize!( mod )
   primer_adjust!( quant, mod )
   primer_adjust!( multi, mod )
   println(stderr, "- $( length(multi.map) ) multi-gene mapping read equivalence classes...")
   build_equivalence_classes!( quant, lib, assign_long=true )
   println(stderr, "- $( length(quant.classes) ) multi-isoform equivalence classes...")
   calculate_tpm!( quant, readlen=readlen )
   @timer iter = gene_em!( quant, multi, sig=1, readlen=readlen, maxit=10000 )
   println(stderr, "Finished calculating transcripts per million (TpM) after $iter iterations of EM...")
   set_gene_tpm!( quant, lib )
   gc_normalize!( mod, lib, quant )
   gc_adjust!( quant, mod )
   gc_adjust!( multi, mod )

   output_tpm( args["out"], lib, quant )
   output_stats( args["out"] * ".map.gz", lib, quant, param, indexpath, totreads, mapped, Int(round(total_multi(multi))), readlen, ver )
   output_junctions( args["out"] * ".jnc.gz", lib, quant )

   println(stderr, "Assigning multi-mapping reads based on maximum likelihood estimate..")
   # Now assign multi to edges.
   @timer assign_ambig!( quant, lib, multi )

   effective_lengths!( lib, quant, 1, 0)
   println(stderr, "Calculating maximum likelihood estimate of events..." )
   @timer process_events( args["out"] * ".psi.gz" , lib, quant, isnodeok=false, iscircok=args["circ"], readlen=readlen )
   println(stderr, "Whippet $ver done." )
end

@timer main()


text:
append!(empty!(LOAD_PATH), Base.DEFAULT_LOAD_PATH)
using Pkg
exampledir = joinpath(@__DIR__, "..", "examples")
Pkg.activate(exampledir)
Pkg.instantiate()
include(joinpath(exampledir, "generate_notebooks.jl"))


text:
@stage_testset "stan" "stan-interface.jl" begin
    using CmdStan

    chn = sample(gdemo_default, 2000, 1000, false, 1,
        CmdStan.Adapt(), CmdStan.Hmc(stepsize_jitter=0))
    check_gdemo(chn)
end


text:
export CG_PR

function CG_PR(nlp :: AbstractNLPModel,
               stp :: NLPStopping;
               verbose :: Bool=false,
               verboseLS :: Bool = false,
               linesearch :: Function = armijo_ls,
               scaling :: Bool = true,
               kwargs...)

    return CG_generic(nlp,
                      stp;
                      verbose = verbose,
                      verboseLS = verboseLS,
                      linesearch  = linesearch,
                      CG_formula  = formula_PR,
                      scaling  = scaling,
                      kwargs...)
end


text:
module TestDataset

using Base.Test
using EconDatasets
using Dates

dates = Date(1960, 11, 26):Date(2014, 7, 23)
ticker = "^GDAXI"

## read in single ticker
kk = readYahooFinance(dates, ticker, :w)

## read multiple ticker, complete data
ticker = ["^GDAXI", "^GSPC", "BMW.DE"]
kk = map((x) -> readYahooFinance(dates, x), ticker)

## read adjusted closing prices for multiple ticker symbols
dates = Date(1960, 11, 26):Date(2014, 7, 23)
ticker = ["^GDAXI", "^GSPC", "BMW.DE"]
vals = readYahooAdjClose(dates, ticker)

end