text:
# # This  function fits Kernel Density and evaluates density at v (Needs KernelDensity.jl)
# function useKDE2fit(x::AbstractArray{Float64,3},v::AbstractMatrix)
#     q,w,e=size(x)
#     y = similar(x,q,w)

#     for a=1:q, s=1:w
#         z = kde( x[a,s,:] )
#         y[a,s] = pdf(z, v[a,s])
#     end

#     return y
# end

# This  function fits Normal pdf and evaluates density at v (Needs Distributions.jl)
function useN2fit(x::AbstractArray{Float64,3}, v::AbstractMatrix)
    q, w, e = size(x)
    y = similar(x, q, w)

    for a in 1:q, s in 1:w
        z = fit(Normal, x[a,s,:])
        y[a,s] = pdf.(z, v[a,s])
    end

    return y
end

# fitting the multivariate version
function useMvN2fit(x::AbstractArray{Float64,3}, v::AbstractMatrix)
    # w is no of variables
    # q is no of forecast horizon
    # e are number of draws
    q, w, e = size(x)
    y = similar(x, q)

    for a = 1:q
        z = fit(MvNormal, x[a,:,:])
        y[a] = pdf(z, v[a,:])
    end

    return y
end

#############################################
## fitting PIT
"
PSPIT Calculate Probability Integral Transform time series.
   [Z, ZHIST] = PSPIT(FORECAST, OBS) obtains the PIT
   z-series for a given vector of observations in OBS, with associated
   empirical probability forecasts given in the matrix FORECASTS. Each
   row of FORECASTS contains an empirical distribution for each element
   of OBS. [...] = PSPIT(..., ZBINS) specifies the number of histogram bins
   on the interval [0,1] (default is 10 bins). [...] =
   PSPIT(..., ZBINS, ZSIGNIF) specifies a significance probability for the
   null hypothesis test of the z-series bin counts (default 0.95 - 95%
   significance level).
   Output parameter Z is the z-series, and ZHIST is a structure
   with the following elements:
    'h.counts'   - Count of the number of z-series elements falling within
                 each histogram bin on the interval [0,1].
    'h.centres'  - Z-series histogram bin centres.
    'confmean' - Expected z-series histogram count, under the null
                 hypothesis of z-series being iid uniform on [0,1].
    'confpos'  - Upper confidence interval for the z-series histogram
                 counts, under the null hypothesis of z-series being iid
                 uniform on [0,1].
    'confneg'  - Lower confidence interval for the z-series histogram
                 counts, under the null hypothesis of z-series being iid
                 uniform on [0,1].
(cc) Max Little, 2008. This software is licensed under the
Attribution-Share Alike 2.5 Generic Creative Commons license:
http://creativecommons.org/licenses/by-sa/2.5/
If you use this work, please cite:
Little MA et al. (2008), Parsimonious Modeling of UK Daily Rainfall for
Density Forecasting, in Geophysical Research Abstracts, EGU General
Assembly, Vienna 2008, Volume 10.
"
function pspit(forecasts::Matrix, obs::Vector; zsignif::Float64 = 0.95, nbins::Int = 10, returnZ::Bool=false)
    N = length(obs)
    M, S = size(forecasts)

    @assert M == N "Length of OBS must match number of rows in FORECASTS"

    z = similar(obs)
    freqs = (1. : S) ./ S
    scale = similar(forecasts, S)
    for i = 1:N
        scale .= sort(forecasts[i,:])

        @views obsI = obs[i]
        if obsI < scale[1]
            z[i] = 0.
        elseif obsI > scale[end]
            z[i] = 1.
        else
            idx = findlast(scale .<= obsI)
            if idx > 1
                z[i] = freqs[idx - 1] + rand(Uniform()) * (freqs[idx] - freqs[idx - 1])
            else
                z[i] = 0.
            end
        end
    end

    if returnZ
        return z
    end

    # Obtaining a histogram
    bins = LinRange(0., 1., nbins+1)
    h = fit(Histogram, z, bins)

    # Find the approximate confintv confidence interval for
    # binomially-distributed PIT histogram bin counts, under
    # the assumption of iid U(0,1) PIT z-series

    bprob = 1. / nbins
    meanheight = convert(Int, floor(N/nbins))::Int
    binocumul(x) = cdf(Binomial(N,bprob), x)
    for confwidth = 1:meanheight
        cip = meanheight + confwidth
        cin = meanheight - confwidth
        ci = binocumul(cip) - binocumul(cin - 1.)
        if ci >= zsignif
            return (z=z, zhist=h, confmean= meanheight, confpos = cip, confneg = cin)
        end
    end

    @warn "the confidence interval could not be compute due to few observation"
    return (z=z, zhist=h, confmean= meanheight)

end

"PSCRPS Calculate Continuous-Ranked Probability Score.
   CRPS = PSCRPS(FORECAST, OBS) obtains the CRPS for a given vector of
   observations in OBS, with associated empirical probability forecasts
   given in the matrix FORECASTS. Each row of FORECASTS contains an
   empirical distribution for each element of OBS. Output parameter is the
   CRPS value."
function pscrps(forecasts::Matrix, obs::Vector)
    N = length(obs)
    M, S = size(forecasts)

    @assert M == N "Length of OBS must match number of rows in FORECASTS"

    crps_individual = similar(obs)
    for i = 1:N
        crps1 = mean(abs.(forecasts[i, :] .- obs[i]))
        crps2 = mean(abs.(diff(forecasts[i, randperm(S)])))
        crps_individual[i] = crps1 - 0.5 * crps2
    end
    crps = mean(crps_individual)
    return (crps = crps, crps_individual = crps_individual)
end

max_stars_repo_name: zymbuzz/myJuliaUtils

max_issues_repo_name: zymbuzz/myJuliaUtils

max_forks_repo_name: zymbuzz/myJuliaUtils


text:
module FrostHelperBadelineChaserBlock

using ..Ahorn, Maple

@mapdef Entity "FrostHelper/BadelineChaserBlock" BadelineChaserBlock(x::Integer, y::Integer, width::Integer=Maple.defaultBlockWidth, height::Integer=Maple.defaultBlockWidth, reversed::Bool=false)

const placements = Ahorn.PlacementDict(
    "Badeline Chaser Block (Frost Helper)" => Ahorn.EntityPlacement(
        BadelineChaserBlock,
        "rectangle",
    )
)

Ahorn.minimumSize(entity::BadelineChaserBlock) = 16, 16
Ahorn.resizable(entity::BadelineChaserBlock) = true, true

function Ahorn.selection(entity::BadelineChaserBlock)
    x, y = Ahorn.position(entity)

    width = Int(get(entity.data, "width", 8))
    height = Int(get(entity.data, "height", 8))

    return Ahorn.Rectangle(x, y, width, height)
end

function getTextures(entity::BadelineChaserBlock)
    return "objects/FrostHelper/badelineChaserBlock/solid", "objects/FrostHelper/badelineChaserBlock/emblemsolid"
end

function renderSwapBlock(ctx::Ahorn.Cairo.CairoContext, x::Number, y::Number, width::Number, height::Number, midResource::String, frame::String)
    midSprite = Ahorn.getSprite(midResource, "Gameplay")

    tilesWidth = div(width, 8)
    tilesHeight = div(height, 8)

    for i in 2:tilesWidth - 1
        Ahorn.drawImage(ctx, frame, x + (i - 1) * 8, y, 8, 0, 8, 8)
        Ahorn.drawImage(ctx, frame, x + (i - 1) * 8, y + height - 8, 8, 16, 8, 8)
    end

    for i in 2:tilesHeight - 1
        Ahorn.drawImage(ctx, frame, x, y + (i - 1) * 8, 0, 8, 8, 8)
        Ahorn.drawImage(ctx, frame, x + width - 8, y + (i - 1) * 8, 16, 8, 8, 8)
    end

    for i in 2:tilesWidth - 1, j in 2:tilesHeight - 1
        Ahorn.drawImage(ctx, frame, x + (i - 1) * 8, y + (j - 1) * 8, 8, 8, 8, 8)
    end

    Ahorn.drawImage(ctx, frame, x, y, 0, 0, 8, 8)
    Ahorn.drawImage(ctx, frame, x + width - 8, y, 16, 0, 8, 8)
    Ahorn.drawImage(ctx, frame, x, y + height - 8, 0, 16, 8, 8)
    Ahorn.drawImage(ctx, frame, x + width - 8, y + height - 8, 16, 16, 8, 8)

    Ahorn.drawImage(ctx, midSprite, x + div(width - midSprite.width, 2), y + div(height - midSprite.height, 2))
end

function Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::BadelineChaserBlock, room::Maple.Room)
    startX, startY = Int(entity.data["x"]), Int(entity.data["y"])

    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    frame, mid = getTextures(entity)

    renderSwapBlock(ctx, startX, startY, width, height, mid, frame)
end

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::BadelineChaserBlock, room::Maple.Room)
    startX, startY = Int(entity.data["x"]), Int(entity.data["y"])

    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    frame, mid = getTextures(entity)

    renderSwapBlock(ctx, startX, startY, width, height, mid, frame)
end

end

ext: jl

max_stars_repo_name: JaThePlayer/FrostHelper

max_issues_repo_name: JaThePlayer/FrostHelper

max_forks_repo_name: JaThePlayer/FrostHelper


text:
################################################################################
# Visuals
################################################################################
function visualize(env::Environment, traj::Vector{Vector{T}}) where T
	@assert size(traj[1]) == size(env.state)
    storage = generate_storage(env.mechanism, [env.representation == :minimal ? minimal_to_maximal(env.mechanism, x) : x for x in traj])
    visualize(env.mechanism, storage,
        vis=env.vis)
end

################################################################################
# Miscellaneous
################################################################################
type2symbol(H) = Symbol(lowercase(String(H.name.name)))


text:
module TestWater

using Test
using LinearAlgebra
using Canopy.Constants: T_0
using Canopy.Water

rho_w_validat = hcat(
    [0., 4., 5., 10., 15., 20., 25., 30.] .+ T_0,
    [999.8395, 999.9720, 999.96, 999.7026, 999.1026,
     998.2071, 997.0479, 995.6502])
pK_w_validat = hcat(
    collect(0:25:100) .+ T_0,
    [14.95, 13.99, 13.26, 12.70, 12.25])

@time @testset "Canopy.Water" begin
    @test norm(water_density.(rho_w_validat[:, 1]) -
               rho_w_validat[:, 2]) / norm(rho_w_validat[:, 2]) < 1e-4
    @test norm(water_dissoc.(pK_w_validat[:, 1]) - pK_w_validat[:, 2]) /
        norm(pK_w_validat[:, 2]) < 1e-3
end

end  # module


text:
using LinearAlgebra: eigen, cholesky, diag, diagm, UpperTriangular, LowerTriangular
using Statistics: mean

"""
  Functions for mvrnorm, covriance to correlation conversion (cov2cor)
  and vice versa cor2cov.
"""
function generateRandomMatrix(p::Int64, n::Int64)
  x = rand(n, p)
  return x' * x
end

"""
  Function to convert covariance to correlation matrix
"""
function cov2cor(E::Array{T, 2}) where {T <: AbstractFloat}
  p = size(E)[1]
  d::Array{T, 1} = diag(E)
  d .^= -0.5
  C = (d * d') .* E
  C = 0.5.*(UpperTriangular(C) + LowerTriangular(C)')
  C = C + C'
  for i in 1:p
    C[i, i] = T(1)
  end
  return C
end

"""
  Function to convert correlation to covariance matrix
  d is the variance vector of the covariance matrix
"""
function cor2cov(C::Array{T, 2}, d::Array{T, 1}) where {T <: AbstractFloat}
  p = size(C)[1]
  d1 = d .^0.5
  E = (d1 * d1') .* C
  E = 0.5.*(UpperTriangular(E) + LowerTriangular(E)')
  E = E + E'
  for i in 1:p
    E[i, i] /= 2
  end
  return E
end

"""
  Function that carries samples from a multivariate normal
  distribution when supplied with the number of samples (n),
  the mean vector (mu), and the covariance matrix (sigma).

  Example:
  mvrnorm(100, zeros(10), cov2cor(generateRandomMatrix(10, 10000)))
"""
function mvrnorm(n::Int64, mu::Array{T, 1}, sigma::Array{T, 2}) where {T <: AbstractFloat}
  p = size(sigma)[1]
  A = cholesky(sigma).L
  output = zeros(T, (n, p))
  for i = 1:n
    output[i, :] = mu + A * randn(p)
  end
  return output
end

#=======================================================================================================#
"""
  Random Number Generator for Beta (and Uniform) Distribution
"""
abstract type AbstractSampleDistribution end
struct BetaSampleDistribution{T} <: AbstractSampleDistribution
  alpha::T
  beta::T
  function BetaSampleDistribution(alpha::T, beta::T) where {T <: AbstractFloat}
    return new{T}(alpha, beta)
  end
  function BetaSampleDistribution(alpha::T, beta::T) where {T <: Integer}
    return BetaSampleDistribution(Float64(alpha), Float64(beta))
  end
end

function calcSample(ualpha::T, vbeta::T) where {T <: AbstractFloat}
  return ualpha/(ualpha + vbeta)
end
function condition(ualpha::T, vbeta::T) where {T <: AbstractFloat}
  return (ualpha + vbeta) > T(1)
end

"""
  Version 1 of the random sample from beta distribution
  Reference: C. P. Robert, G. Casella,
            Monte Carlo Statistical Methods, Example 2.11 p44
"""
function sample(distrib::BetaSampleDistribution{T}, shape::Tuple{Vararg{Int64}}) where {T <: AbstractFloat}
  n = prod(shape)
  U = rand(T, n); V = rand(T, n)
  Y::Array{T} = zeros(T, shape)
  ialpha = 1/distrib.alpha; ibeta = 1/distrib.beta
  for i in 1:n
    u::T = U[i]; v::T = V[i]
    ualpha = u .^ialpha; vbeta = v .^ibeta;
    while condition(ualpha, vbeta)
      u = rand(T, 1)[1]; v = rand(T, 1)[1]
      ualpha = u .^ialpha; vbeta = v .^ibeta;
    end
    Y[i] = calcSample(ualpha, vbeta)
  end
  return Y
end

function sample(distrib::BetaSampleDistribution{T}, shape::Int64) where {T <: AbstractFloat}
  return sample(distrib, (shape,))
end

#=======================================================================================================#

struct UniformSampleDistribution{T} <: AbstractSampleDistribution
  min::T
  max::T
  function UniformSampleDistribution(min::T, max::T) where {T <: AbstractFloat}
    @assert(min < max, "Minimum value is not less than maximum value")
    return new{T}(min, max)
  end
  function UniformSampleDistribution(min::T, max::T) where {T <: Integer}
    return UniformSampleDistribution(Float64(min), Float64(max))
  end
end

function sample(distrib::UniformSampleDistribution{T}, shape::Tuple{Vararg{Int64}}) where {T <: AbstractFloat}
  rsample::Array{T} = rand(T, shape)
  return (rsample .* (distrib.max - distrib.min)) .+ distrib.min
end
function sample(distrib::UniformSampleDistribution{T}, shape::Int64) where {T <: AbstractFloat}
  return sample(distrib, (shape,))
end

function Base.min(x::Array{T}) where {T <: AbstractFloat}
  n = length(x)
  @assert(n > 0, "Length of array is zero.")
  if n == 1
    return x[1]
  end
  ret::T = x[1]
  for i in 2:n
    ret = ret > x[i] ? x[i] : ret
  end
  return ret
end

function Base.max(x::Array{T}) where {T <: AbstractFloat}
  n = length(x)
  @assert(n > 0, "Length of array is zero.")
  if n == 1
    return x[1]
  end
  ret::T = x[1]
  for i in 2:n
    ret = ret < x[i] ? x[i] : ret
  end
  return ret
end

function Base.range(x::Array{T}) where {T <: AbstractFloat}
  return [min(x), max(x)]
end

#=======================================================================================================#

"""
  Identity Matrix
"""
function I(::Type{T}, p::Int64) where {T <: Number}
  x = zeros(T, (p, p))
  for i in 1:p
    x[i, i] = T(1)
  end
  return x
end
"""
  Convenience overload
"""
function I(p::Int64)
  return I(Float64, p)
end

"""
  Types for generating random correlation matrices
"""
abstract type AbstractRandomCorrelationMatrix end

struct BetaGenerator <: AbstractRandomCorrelationMatrix end
struct OnionGenerator <: AbstractRandomCorrelationMatrix end
struct UniformGenerator <: AbstractRandomCorrelationMatrix end
struct VineGenerator <: AbstractRandomCorrelationMatrix end

"""
  Vine method for generating a random correlation matrix
  Source:
    https://stats.stackexchange.com/questions/2746/how-to-efficiently-generate-random-positive-semidefinite-correlation-matrices/125017#125017

    Check it with to make sure the solutions are correct:
    Journal of Multivariate Analysis 100 (2009) 1989–2001
    Generating random correlation matrices based on vines and extended
    onion method. Daniel Lewandowski(a)*, Dorota Kurowicka (a),
    Harry Joe (b).

  randomCorrelationMatrix(VineGenerator(), 10, 2.0)
"""
function randomCorrelationMatrix(::VineGenerator, d::Int64, eta::T) where {T <: AbstractFloat}

  beta = eta + (d - 1)/2
  P = zeros(T, (d, d))
  S = I(T, d)

  for k in 1:(d - 1)
    beta = beta - (1/2)
    distrib = BetaSampleDistribution(beta, beta)

    for i in (k + 1):d
      P[k, i] = sample(distrib, 1)[1]
      P[k, i] = (P[k, i] - 0.5)*2
      p = P[k, i]

      for l in (k - 1):-1:1
        p = p * sqrt((1 - P[l, i]^2)*(1 - P[l, k]^2)) + P[l, i]*P[l, k]
      end

      S[k, i] = p
      S[i, k] = p
    end
  end

  return S
end

"""
  Random correlation matrix using the Onion Generator
"""
function randomCorrelationMatrix(::OnionGenerator, d::Int64, eta::T) where {T <: AbstractFloat}
  beta::T = eta + (d - 2)/2
  distrib = BetaSampleDistribution(beta, beta)
  u = sample(distrib, 1)[1]
  r = I(T, 2)
  r[1, 2] = r[2, 1] = 2*u - 1
  for k in 2:(d - 1)
    beta -= T(1/2)
    distrib = BetaSampleDistribution(T(k/2), beta)
    y = sample(distrib, 1)[1]
    U = rand(T, k)
    w = sqrt(y) .* U

    ev = eigen(r)
    A = ev.vectors * diagm(abs.(ev.values).^(0.5)) * ev.vectors'
    z = A * w
    r = [r z; z' T(1)]
  end
  return r
end

"""
  Random correlation matrix by sampling from Beta Distribution
  # Example:
  randomCorrelationMatrix(BetaGenerator(), 10, (1.0, 1.0))
"""
function randomCorrelationMatrix(::BetaGenerator, d::Int64, (alpha, beta)::Tuple{T, T}) where {T <: AbstractFloat}
  distrib = BetaSampleDistribution(alpha, beta)
  r = sample(distrib, (d, d))
  # Change range to (-1, 1)
  r .= (r .* 2) .+ T(-1)
  # Symmetry
  r .= T(0.5) .* (r + r')
  for i in 1:d
    r[i, i] = T(1)
  end
  ev = eigen(r)
  r = ev.vectors * diagm(sort(abs.(ev.values))) * ev.vectors'
  maxR = max(r)
  r ./= maxR
  r .= T(0.5) .* (r + r')
  for i in 1:d
    r[i, i] = T(1)
  end
  return r
end

"""
  Random correlation matrix by sampling from Uniform Distribution
  # Example:
  randomCorrelationMatrix(Float64, UniformGenerator(), 10)
"""
function randomCorrelationMatrix(::Type{T}, ::UniformGenerator, d::Int64) where {T <: AbstractFloat}
  distrib = UniformSampleDistribution(T(-1), T(1))
  r = sample(distrib, (d, d))
  r .= T(0.5) .* (r + r')
  for i in 1:d
    r[i, i] = T(1)
  end
  ev = eigen(r)
  r = ev.vectors * diagm(sort(abs.(ev.values))) * ev.vectors'
  maxR = max(abs.(r)...)
  r ./= maxR
  r .= T(0.5) .* (r + r')
  for i in 1:d
    r[i, i] = T(1)
  end
  return r
end

#=======================================================================================================#

"""
  Function to simulate X and eta, p = number of parameters
  (including intercept), n = number of samples.

  # Example:
    using Random: seed!
    seed!(0)
    simulateData(Float64, 10, 1000)
"""
function simulateData(::Type{T}, p::Int64, n::Int64, delta::T = T(0)) where {T <: AbstractFloat}

  corr = randomCorrelationMatrix(T, UniformGenerator(), p)
  mu = zeros(T, p)
  X = mvrnorm(n, mu, corr)

  #= The intercept =#
  X[:, 1] .= T(1)

  b = zeros(T, p)
  idist = UniformSampleDistribution(T(0), T(0.3))
  b[1] = sample(idist, 1)[1]
  if length(b) > 1
    distrib = UniformSampleDistribution(T(-0.1), T(0.1))
    b[2:p] = sample(distrib, p - 1)
  end

  eta = X*b
  sd = 0.5*abs(mean(eta))
  eta = delta .+ eta .+ sd .* randn(n)

  return (X = X, eta = eta)
end

max_stars_repo_name: dataPulverizer/glmSolver

max_issues_repo_name: dataPulverizer/glmSolver

max_forks_repo_name: dataPulverizer/glmSolver


text:
using Submodular
using FactCheck

Tol = 1e-3

facts("Algorithms") do

  context("card-inc-fix") do
    n = 20
    g = rand(n)                      # the vector used to generate the cardinality-based function
    g = sort(g, rev=true)
    g = (g + 0.1)/1.01
    for i = 2: n
      g[i] = g[i] + g[i-1]
    end
    S = SetVariable(n)
    gg = card(g, S)                  # F(S) = g(|S|)

    y = rand(n)                      # the point to be projected on the base polytope of F(S) = g(|S|)
    euclidean_proj = card_inc_fix(gg, y, "euclidean")

    # sanity checks:
    # Is the sorted order of indices in y and the projection the same?
    # this is a known property of projections under uniform divergences over cardinality-based polytopes
    @fact sortperm(y) == sortperm(euclidean_proj[1:n]) --> true

    # Is the constraint x(E) = F(E) = g(n) satisfied upto an error of Tol^3?
    @fact sum(euclidean_proj) - g[n] --> roughly(0, Tol^3)
  end

  context("frank-wolfe with away steps") do
    n = 4
    x = Variable(n)
    g = norm(x - [3, 2, 5, 1])
    S = SetVariable(n)
    F = card(S)
    p(z) = -0.5*z^2 + n*z + 0.5 * z
    perm_func = compose(p, F)
    P = BasePoly(perm_func)
    prob = SCOPEminimize(g, x in P)
    @fact frank_wolfe_away(prob, verbose = false) --> roughly([3, 2, 4, 1], Tol)
  end

end


text:
# JuMP Objectives/constraints #
# ========================== #
function JuMP.normalized_coefficient(
    con_ref::ConstraintRef{Model, JuMP._MOICON{F, S}}, decision
    ) where {S, T, F <: AffineDecisionFunction{T}}
    con = JuMP.constraint_object(con_ref)
    return JuMP._affine_coefficient(con.func, decision)
end

function JuMP.set_objective_coefficient(model::Model, decision_or_known::Union{DecisionRef, KnownRef}, coeff::Real)
    if model.nlp_data !== nothing && _nlp_objective_function(model) !== nothing
        error("A nonlinear objective is already set in the model")
    end
    obj_fct_type = objective_function_type(model)
    if obj_fct_type == VariableRef || obj_fct_type == AffExpr || obj_fct_type == QuadExpr
        current_obj = objective_function(model)
        set_objective_function(model, add_to_expression!(coeff * decision_or_known, current_obj))
    elseif obj_fct_type == typeof(decision_or_known)
        current_obj = objective_function(model)
        if index(current_obj) == index(decision_or_known)
            set_objective_function(model, coeff * decision_or_known)
        else
            set_objective_function(model, add_to_expression!(coeff * decision_or_known, current_obj))
        end
    elseif obj_fct_type == DecisionAffExpr{Float64} && decision_or_known isa DecisionRef
        MOI.modify(
            backend(model),
            MOI.ObjectiveFunction{moi_function_type(obj_fct_type)}(),
            DecisionCoefficientChange(index(decision_or_known), coeff))
    elseif obj_fct_type == DecisionAffExpr{Float64} && decision_or_known isa KnownRef
        MOI.modify(
            backend(model),
            MOI.ObjectiveFunction{moi_function_type(obj_fct_type)}(),
            KnownCoefficientChange(index(decision_or_known), coeff, value(decision_or_known)))
    else
        error("Objective function type not supported: $(obj_fct_type)")
    end
    return nothing
end

function JuMP.normalized_rhs(con_ref::ConstraintRef{Model, JuMP._MOICON{F, S}}) where {
        T, S <: Union{MOI.LessThan{T}, MOI.GreaterThan{T}, MOI.EqualTo{T}},
        F <: AffineDecisionFunction{T}}
    con = constraint_object(con_ref)
    return MOI.constant(con.set)
end

function JuMP.set_normalized_coefficient(
    con_ref::ConstraintRef{Model, JuMP._MOICON{F, S}}, variable::VariableRef, coeff
    ) where {S, T, F <: AffineDecisionFunction{T}}
    MOI.modify(backend(owner_model(con_ref)), index(con_ref),
               MOI.ScalarCoefficientChange(index(variable), convert(T, coeff)))
    return nothing
end

function JuMP.set_normalized_coefficient(
    con_ref::ConstraintRef{Model, JuMP._MOICON{F, S}}, decision::DecisionRef, coeff
    ) where {S, T, F <: AffineDecisionFunction{T}}
    MOI.modify(backend(owner_model(con_ref)), index(con_ref),
               DecisionCoefficientChange(index(decision), convert(T, coeff)))
    return nothing
end

function JuMP.set_normalized_coefficient(
    con_ref::ConstraintRef{Model, JuMP._MOICON{F, S}}, known::KnownRef, coeff
    ) where {S, T, F <: AffineDecisionFunction{T}}
    MOI.modify(backend(owner_model(con_ref)), index(con_ref),
               KnownCoefficientChange(index(known), convert(T, coeff), convert(T, value(known))))
    return nothing
end

# Internal update functions #
# ========================== #
function update_decisions!(model::JuMP.Model, change::Union{DecisionModification, KnownModification})
    update_decision_objective!(model, objective_function_type(model), change)
    update_decision_variable_constraints!(model, change)
    update_decision_constraints!(model, change)
end

function update_decision_objective!(::JuMP.Model, ::DataType, ::Union{DecisionModification, KnownModification})
    # Nothing to do if objective does not have decisions
    return nothing
end

function update_decision_objective!(model::JuMP.Model, func_type::Type{<:DecisionQuadExpr}, change::DecisionModification)
    MOI.modify(backend(model),
               MOI.ObjectiveFunction{JuMP.moi_function_type(func_type)}(),
               change)
    return nothing
end

function update_decision_objective!(model::JuMP.Model, func_type::Type{F}, change::KnownModification) where F <: Union{DecisionAffExpr, DecisionQuadExpr}
    MOI.modify(backend(model),
               MOI.ObjectiveFunction{JuMP.moi_function_type(func_type)}(),
               change)
    return nothing
end

function update_decision_variable_constraints!(::JuMP.Model, ::Union{DecisionModification, KnownModification})
    # Nothing to do in most cases
    return nothing
end

function update_decision_variable_constraints!(model::JuMP.Model, change::DecisionStateChange)
    for F in [DecisionRef, ]
        for S in [MOI.EqualTo{Float64}, MOI.LessThan{Float64}, MOI.GreaterThan{Float64}, FreeDecision]
            for cref in all_constraints(model, F, S)
                update_decision_constraint!(cref, change)
            end
        end
    end
    return nothing
end

function update_decision_variable_constraints!(model::JuMP.Model, ::DecisionsStateChange)
    for F in [DecisionRef, ]
        for S in [MOI.EqualTo{Float64}, MOI.LessThan{Float64}, MOI.GreaterThan{Float64}, FreeDecision]
            for cref in all_constraints(model, F, S)
                # Fetch the decision
                dref = JuMP.jump_function(model, MOI.get(model, MOI.ConstraintFunction(), cref))
                # Perform specific decision state change
                change = DecisionStateChange(index(dref), state(dref), 0.0)
                update_decision_constraint!(cref, change)
            end
        end
    end
    return nothing
end

function update_decision_constraints!(model::JuMP.Model, change::Union{DecisionModification, KnownModification})
    for F in [DecisionAffExpr{Float64}, DecisionQuadExpr{Float64}]
        for S in [MOI.EqualTo{Float64}, MOI.LessThan{Float64}, MOI.GreaterThan{Float64}]
            for cref in all_constraints(model, F, S)
                update_decision_constraint!(cref, change)
            end
        end
    end
    for F in [Vector{DecisionAffExpr{Float64}}]
        for S in [MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives]
            for cref in all_constraints(model, F, S)
                update_decision_constraint!(cref, change)
            end
        end
    end
    return nothing
end

function update_decision_constraint!(cref::ConstraintRef, change::Union{DecisionModification, KnownModification})
    update_decision_constraint!(backend(owner_model(cref)), cref.index, change)
    return nothing
end

function update_decision_constraint!(model::MOI.ModelLike, ci::MOI.ConstraintIndex{SingleDecision, S}, change::DecisionModification) where {T,S}
    MOI.modify(model, ci, change)
    return nothing
end

function update_decision_constraint!(model::MOI.ModelLike, ci::MOI.ConstraintIndex{VectorOfDecisions, S}, change::VectorDecisionModification) where {T,S}
    MOI.modify(model, ci, change)
    return nothing
end

function update_decision_constraint!(model::MOI.ModelLike, ci::MOI.ConstraintIndex{AffineDecisionFunction{T}, S}, change::Union{DecisionCoefficientChange, KnownModification}) where {T,S}
    MOI.modify(model, ci, change)
    return nothing
end

function update_decision_constraint!(model::MOI.ModelLike, ci::MOI.ConstraintIndex{AffineDecisionFunction{T}, S}, change::Union{DecisionStateChange, DecisionsStateChange}) where {T,S}
    return nothing
end

function update_decision_constraint!(model::MOI.ModelLike, ci::MOI.ConstraintIndex{VectorAffineDecisionFunction{T}, S}, change::Union{DecisionMultirowChange, KnownModification}) where {T,S}
    MOI.modify(model, ci, change)
    return nothing
end

function update_decision_constraint!(model::MOI.ModelLike, ci::MOI.ConstraintIndex{VectorAffineDecisionFunction{T}, S}, change::Union{DecisionStateChange, DecisionsStateChange}) where {T,S}
    return nothing
end

function update_decision_constraint!(model::MOI.ModelLike, ci::MOI.ConstraintIndex{<:QuadraticDecisionFunction{T}, S}, change::Union{DecisionModification, KnownModification}) where {T,S}
    MOI.modify(model, ci, change)
    return nothing
end


text:
module Distance

include("dist.jl")

export
    #   Funtion
    GetDistance,

    #   Types of distances
    Euclidean,
    CityBlock,
    TotalVariation,
    Chebyshev,
    Jaccard,
    BrayCurtis,
    CosineDist,
    SpanNormDist
end


text:
using FromFile
using Random: randperm
using LossFunctions
@from "Core.jl" import Options, Dataset, Node
@from "EquationUtils.jl" import countNodes
@from "EvaluateEquation.jl" import evalTreeArray, differentiableEvalTreeArray

function Loss(x::AbstractArray{T}, y::AbstractArray{T}, options::Options{A,B,C})::T where {T<:Real,A,B,C<:SupervisedLoss}
    value(options.loss, y, x, AggMode.Mean())
end
function Loss(x::AbstractArray{T}, y::AbstractArray{T}, options::Options{A,B,C})::T where {T<:Real,A,B,C<:Function}
    sum(options.loss.(x, y))/length(y)
end

function Loss(x::AbstractArray{T}, y::AbstractArray{T}, w::AbstractArray{T}, options::Options{A,B,C})::T where {T<:Real,A,B,C<:SupervisedLoss}
    value(options.loss, y, x, AggMode.WeightedMean(w))
end
function Loss(x::AbstractArray{T}, y::AbstractArray{T}, w::AbstractArray{T}, options::Options{A,B,C})::T where {T<:Real,A,B,C<:Function}
    sum(options.loss.(x, y, w))/sum(w)
end

# Loss function. Only MSE implemented right now. TODO
# Also need to put actual loss function in scoreFuncBatch!
function EvalLoss(tree::Node, dataset::Dataset{T}, options::Options;
                  allow_diff=false)::T where {T<:Real}
    if !allow_diff
        (prediction, completion) = evalTreeArray(tree, dataset.X, options)
    else
        (prediction, completion) = differentiableEvalTreeArray(tree, dataset.X, options)
    end
    if !completion
        return T(1000000000)
    end

    if dataset.weighted
        return Loss(prediction, dataset.y, dataset.weights, options)
    else
        return Loss(prediction, dataset.y, options)
    end
end

# Score an equation
function scoreFunc(dataset::Dataset{T},
                   baseline::T, tree::Node,
                   options::Options; allow_diff=false)::T where {T<:Real}
    mse = EvalLoss(tree, dataset, options; allow_diff=allow_diff)
    return mse / baseline + countNodes(tree)*options.parsimony
end

# Score an equation with a small batch
function scoreFuncBatch(dataset::Dataset{T}, baseline::T,
                        tree::Node, options::Options)::T where {T<:Real}
    batchSize = options.batchSize
    batch_idx = randperm(dataset.n)[1:options.batchSize]
    batch_X = dataset.X[:, batch_idx]
    batch_y = dataset.y[batch_idx]
    (prediction, completion) = evalTreeArray(tree, batch_X, options)
    if !completion
        return T(1000000000)
    end

    if !dataset.weighted
        mse = Loss(prediction, batch_y, options)
    else
        batch_w = dataset.weights[batch_idx]
        mse = Loss(prediction, batch_y, batch_w, options)
    end
    return mse / baseline + countNodes(tree) * options.parsimony
end


text:
ENV["GKSwstype"]="100"

using Test
using WaterWaves1D

include("./testmodels.jl")

# @testset "LoadSave" begin
#
#     dump  = convert( ProblemSave, problem1 )
#     pload = convert( Problem, dump )
#
#     save(problem1, "testsave")
#     pload = loadpb("testsave")
#
#     @test pload.model.kwargs == problem1.model.kwargs
#     @test pload.solver.Uhat == problem1.solver.Uhat
#
# end

param = ( ϵ  = 1/2, μ = 1)
paramX= ( N  = 2^8, L  = 10)
paramT= ( T  = 5, dt = 0.1)

@testset "Parameters" begin

    @test param.ϵ  == 0.5
    @test param.μ  == 1
    @test paramX.N  == 256
    @test paramX.L  == 10
    @test paramT.T  == 5
    @test paramT.dt == 0.1

end


text:
function calcJacobian!(rr::RadauIntegrator{T_object,NX,NC,NR,NSM}, x0::Vector{Float64}, t::Float64) where {T_object,NX,NC,NR,NSM}
    N_Loop, n_rem = divrem(NX, NC)
    (n_rem != 0) && (N_Loop += 1)
    i_now = 1:NC
    for k = 1:N_Loop
        i_clamp = i_now[1]:min(i_now[end], NX)
        seed_indices!(rr.cv.x_dual, x0, i_clamp, rr.cv.seed)
        rr.de_object.de(rr.cv.xx_dual, rr.cv.x_dual, rr.de_object, t)
        write_indices!(rr, i_clamp)
        i_now = i_now .+ NC
    end
    return nothing
end

function seed_indices!(duals::Vector{ForwardDiff.Dual{T,V,N}}, x::Vector{Float64}, index::UnitRange{Int64},
    seed::NTuple{N,ForwardDiff.Partials{N,V}}) where {T,V,N}

    duals .= x
    i = 1
    for k = index
        duals[k] = ForwardDiff.Dual{T,V,N}(x[k], seed[i])
        i += 1
    end
    return nothing
end

function write_indices!(rr::RadauIntegrator{T_object,NX,NC,NR,NSM}, index::UnitRange{Int64}) where {T_object,NX,NC,NR,NSM}
    for i = 1:NX
        xx_dual_i = rr.cv.xx_dual[i]
        rr.cv.xx_0[i] = ForwardDiff.value(xx_dual_i)
        the_partials = ForwardDiff.partials(xx_dual_i)
        k = 1
        for j = index  # for each x
            rr.cv.neg_J[i, j] = -the_partials[k]  ### the first dual is the partial of the first element wrt all partials ###
            k += 1
        end
    end
    return nothing
end

function zeroFill!(tup_vec_in::NTuple{N,Vector{T}}, table::RadauTable{NS}) where {N,T,NS}
    for i = 1:NS
        fill!(tup_vec_in[i], zero(T))
    end
    return nothing
end

function initialize_X!(rr::RadauIntegrator{T_object,NX,NC,NR,NSM}, table::RadauTable{NS}, x0::Vector{Float64}) where {T_object,NX,NC,NR,NSM,NS}
    if rr.dense.is_has_X
        initialize_X_with_interp!(rr, table)
    else
        initialize_X_with_X0!(rr, table, x0)
    end
end

function initialize_X_with_X0!(rr::RadauIntegrator{T_object,NX,NC,NR,NSM}, table::RadauTable{NS}, x0::Vector{Float64}) where {T_object,NX,NC,NR,NSM,NS}
    for i = 1:NS
        LinearAlgebra.BLAS.blascopy!(NX, x0, 1, rr.ct.X_stage[i], 1)
    end
    return nothing
end

function updateFX!(rr::RadauIntegrator{T_object,NX,NC,NR,NSM}, table::RadauTable{NS}, x0::Vector{Float64}, t::Float64) where {T_object,NX,NC,NR,NSM,NS}
    for i = 1:NS
        time_stage = table.c[i] * rr.step.h + t
        rr.de_object.de(rr.ct.F_X_stage[i], rr.ct.X_stage[i], rr.de_object, time_stage)
    end
    return nothing
end

function calcEw!(rr::RadauIntegrator{T_object,NX,NC,NR,NSM}, table::RadauTable{NS}, x0::Vector{Float64}) where {T_object,NX,NC,NR,NSM,NS}
    residual = 0.0
    for i = 1:NS
        # rr.cv.store_float .= rr.ct.X_stage[i] - x0 - rr.step.h * sum( rr.A[i, j] * rr.ct.F_X_stage[j])
        LinearAlgebra.BLAS.blascopy!(NX, rr.ct.X_stage[i], 1, rr.cv.store_float, 1)
        LinearAlgebra.BLAS.axpy!(-1.0, x0, rr.cv.store_float)
        for j = 1:NS
            coefficient = -rr.step.h * table.A[i, j]
            # rr.cv.store_float .-= (rr.h * rr.A[i, j]) * rr.ct.F_X_stage[j]
            LinearAlgebra.BLAS.axpy!(coefficient, rr.ct.F_X_stage[j], rr.cv.store_float)
        end
        residual += dot(rr.cv.store_float, rr.cv.store_float)
        for j = 1:NS
            # rr.ct.Ew_stage[j] .+= (rr.step.h⁻¹ * rr.λ[j] * rr.T⁻¹[j, i]) * rr.cv.store_float
            coefficient = rr.step.h⁻¹[1] * table.λ[j] * table.T⁻¹[j, i]
            LinearAlgebra.BLAS.axpy!(coefficient, rr.cv.store_float, rr.ct.Ew_stage[j])
        end
    end
    return residual
end

function updateInvC!(rr::RadauIntegrator{T_object,NX,NC,NR,NSM}, table::RadauTable{NS}) where {T_object,NX,NC,NR,NSM,NS}
    for i = 1:NS
        rr.ct.inv_C_stage[i] .= rr.cv.neg_J
        for k = 1:NX
            rr.ct.inv_C_stage[i][k, k] += rr.step.h⁻¹[1] * table.λ[i]
        end
        ### NOTE: Negative sign taken care of when update X_stage ###
        _, ipiv, info = LinearAlgebra.LAPACK.getrf!(rr.ct.inv_C_stage[i])  # rr.cv.store_complex .= - inv(rr.cv.C) * rr.ct.Ew_stage[i]
        LinearAlgebra.LAPACK.getri!(rr.ct.inv_C_stage[i], ipiv)
    end
    return nothing
end

function updateStageX!(rr::RadauIntegrator{T_object,NX,NC,NR,NSM}, table::RadauTable{NS}) where {T_object,NX,NC,NR,NSM,NS}
    for i = 1:NS
        ### NOTE: Negative sign taken care of when update X_stage ###
        # rr.cv.store_complex .= rr.ct.inv_C_stage[i] * rr.ct.Ew_stage[i]
        LinearAlgebra.BLAS.gemv!('N', one(ComplexF64), rr.ct.inv_C_stage[i], rr.ct.Ew_stage[i], zero(ComplexF64), rr.cv.store_complex)
        for j = 1:NS
            # rr.ct.delta_Z_stage[j] .+= rr.T[j, i] * rr.cv.store_complex
            LinearAlgebra.BLAS.axpy!(table.T[j, i], rr.cv.store_complex, rr.ct.delta_Z_stage[j])
        end
    end
    for i = 1:NS  # Update X_stage
        # LinearAlgebra.BLAS.axpy!(1.0, real.(rr.ct.delta_Z_stage[i]), rr.ct.X_stage[i])
        rr.ct.X_stage[i] .-= real.(rr.ct.delta_Z_stage[i])
    end
    return nothing
end


text:
using LinearAlgebra
using SparseArrays

function degree_selector(t, M, U, p)
    C = ceil.(abs.(t)*M)'*U
    C = zero_to_inf.(C)

    # idx is a CarthesianIndex if C' is a Matrix, or a scalar if C' is a row
    # vector. idx[1] extract the first index, i.e. row, of the CarthesianIndex
    cost, idx = findmin(C')
    m = idx[1]

    if cost == Inf
        cost = 0
    end
    s = max(cost/m,1);
    return (m, s)
end

function zero_to_inf(x::Number)
    if x == 0.
        Inf
    else
        x
    end
end


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ 0abe9cd0-7a26-11eb-0d6d-f986ed77c036
using SpecialFunctions, HypergeometricFunctions, StaticArrays

# ╔═╡ c009caa0-7a27-11eb-28ef-bb4011cdcd50
using ClenshawCurtisBessel

# ╔═╡ 2176cb22-7a63-11eb-0e79-377be1de3061
using DoubleFloats

# ╔═╡ d311f31a-7a5f-11eb-3f9b-c3fd0c2655da
using LinearAlgebra

# ╔═╡ 386fbd2c-7a2a-11eb-051b-13b127ecfc30
struct BrandersPiessensProblem{T,N,M} <: ClenshawCurtisBessel.OliverProblem{T,N,M}
	a::T
	ν::T
end

# ╔═╡ bc00d5ee-7a5d-11eb-2eb5-5b938e7df938
md"""

Computes the modifed moments from Branders & Piessens 1987.

```math
\begin{aligned}
\frac{a^2}{16} &M_{k}(a, \nu) + \left[ (k+1)^2 - \nu^2 - \frac{a^2}{4} \right] M_{k + 2}(a, \nu) + \left[ 4 \nu^2 - 2(k+4) + 4 \right] M_{k+3}(a,\nu) \\
&- \left[ 2(k+4)^2 - 6 + 6\nu^2 - \frac{3a^2}{8} \right] M_{k+4}(a,\nu) + (4\nu^2 + 2(k+4) + 4)M_{k+5}(a,\nu) \\
&+ \left[ (k+7)^2 - \nu^2 - \frac{a^2}{4} \right]M_{k+6} + \frac{a^2}{16} M_{k+8}(a,\nu) \quad = 0
\end{aligned}
```
"""

# ╔═╡ d9adcfb6-7a4e-11eb-2ada-fbdffca3d2bc

function ClenshawCurtisBessel.OliverP(
		OP::BrandersPiessensProblem{T}, s, k) where T
	a, ν = OP.a, OP.ν
	if s == 8
		return a^2/16
	elseif s == 7
		return zero(T)
	elseif s == 6
		return (k+1)^2 - ν^2 - a^2/4
	elseif s == 5
		return 4ν^2 - 2(k+4) + 4
	elseif s == 4
		return -(2 * (k+4)^2 - 6 + 6ν^2 - 3 * a^2 / 8)
	elseif s == 3
		return 4ν^2 + 2(k+4) + 4
	elseif s == 2
		return (k+7)^2 - ν^2 - a^2/4
	elseif s == 1
		return zero(T)
	elseif s == 0
		return a^2/16
	end
end

# ╔═╡ 29239542-7a56-11eb-152d-b386e0588444
"""
	OliverR(OP::BrandersPiessensProblem{T}, i)

Computes ``R(i)``, the right side of the recurrence relation.
"""
ClenshawCurtisBessel.OliverR(
	OP::BrandersPiessensProblem{T}, i) where T = zero(T)

# ╔═╡ dbf33dbe-7a54-11eb-0410-9b91fb2d22c8
import ClenshawCurtisBessel: M₀, M₁, M₂, M₃, Mk_asymptotic

# ╔═╡ 97550070-7a54-11eb-3df8-6ff173e946fc
function generate_BC(OP::BrandersPiessensProblem{T,8,2}, maxorder) where T
	a, ν = OP.a, OP.ν
	MBC = Dict{Int,T}()
	MBC[-3] = M₃(a, ν)
	MBC[-2] = M₂(a, ν)
	MBC[-1] = M₁(a, ν)
	MBC[ 0] = M₀(a, ν)
	MBC[ 1] = M₁(a, ν)
	MBC[ 2] = M₂(a, ν)
	MBC[ 3] = M₃(a, ν)

	# Branders & Piessens call the max moment N, we call it maxorder
	ℓ = max(maxorder, Int(ceil(a + 10)))
	print(typeof(a), " ", a, " ", ν, " ", ℓ, "\n")
	MBC[ℓ] = Mk_asymptotic(a, ν, ℓ)
	MBC[ℓ+1] = Mk_asymptotic(a, ν, ℓ+1)
	return ℓ+1, MBC
end

# ╔═╡ 9740e914-7a54-11eb-1a68-0f136447f214
import ClenshawCurtisBessel: assembleP, assembleρ
# using BandedMatrices

# ╔═╡ 3be2f3f0-7a63-11eb-248f-7727012b9b60
HypergeometricFunctions.logabsgamma(x::Double64) =
	DoubleFloats.loggamma(abs(x)), sign(gamma(x))

# ╔═╡ 6bbd6c44-7a70-11eb-3859-bd57ba21c45d
# SpecialFunctions.besselj(x::Double64, y::AbstractFloat) =
# 	SpecialFunctions.besselj(Float64(x), Float64(x))

# ╔═╡ 972d5cdc-7a54-11eb-0619-85688bb73fcf
@time begin
	bpp = BrandersPiessensProblem{DoubleFloat,8,2}(10.0, 3.5)
	indexBC, MBC = generate_BC(bpp, 200)
	ρ = assembleρ(bpp, -3, indexBC, MBC)
	𝐏 = assembleP(bpp, -3, indexBC)
	sol = qr(𝐏) \ ρ
end

# ╔═╡ ae66d5ca-7a58-11eb-3001-bf4a310c41d4
begin
	ref = big"0.0002069511037367724863632484164263304887654"
	(sol[11] .- ref) ./ ref
end

# ╔═╡ Cell order:
# ╠═0abe9cd0-7a26-11eb-0d6d-f986ed77c036
# ╠═c009caa0-7a27-11eb-28ef-bb4011cdcd50
# ╠═386fbd2c-7a2a-11eb-051b-13b127ecfc30
# ╟─bc00d5ee-7a5d-11eb-2eb5-5b938e7df938
# ╠═d9adcfb6-7a4e-11eb-2ada-fbdffca3d2bc
# ╠═29239542-7a56-11eb-152d-b386e0588444
# ╠═dbf33dbe-7a54-11eb-0410-9b91fb2d22c8
# ╠═97550070-7a54-11eb-3df8-6ff173e946fc
# ╠═9740e914-7a54-11eb-1a68-0f136447f214
# ╠═2176cb22-7a63-11eb-0e79-377be1de3061
# ╠═d311f31a-7a5f-11eb-3f9b-c3fd0c2655da
# ╠═3be2f3f0-7a63-11eb-248f-7727012b9b60
# ╠═6bbd6c44-7a70-11eb-3859-bd57ba21c45d
# ╠═972d5cdc-7a54-11eb-0619-85688bb73fcf
# ╠═ae66d5ca-7a58-11eb-3001-bf4a310c41d4


text:
module Syntax
export @traits, @traits_show_implementation

using ExprParsers
import WhereTraits
using WhereTraits: CONFIG
using WhereTraits.Utils
using WhereTraits.InternalState
using Suppressor
include("Lowering.jl")
using .Lowering
include("Parsing.jl")
using .Parsing
include("Rendering.jl")
using .Rendering

"""
@traits f(a, b) where {!isempty(a), !isempty(b)} = (a[1], b[1])
"""
macro traits(expr_original)
  expr_expanded = macroexpand(__module__, expr_original)
  expr_traits = _traits(@MacroEnv, expr_expanded, expr_original)
  expr_traits = esc(expr_traits)
  if CONFIG.suppress_on_traits_definitions
    expr_traits = :(@suppress $expr_traits)
  end
  expr_traits
end

function _traits(env, expr_expanded::Expr, expr_original::Expr)
  parser = EP.AnyOf(EP.Function(), EP.anything)
  _traits_parsed(env, parse_expr(parser, expr_expanded), expr_original)
end

function _traits_parsed(env, func_parsed::EP.Function_Parsed, expr_original::Expr)
  basefunc, lowerings = lower_args_default(func_parsed)

  basefunc_outer, basefunc_inner = parse_traitsfunction(env, basefunc, expr_original)
  exprs = merge_and_render_update(env, basefunc_outer, basefunc_inner, doc = true)

  for lowering in lowerings
    # As lowering dropped variables, also traits may need to be dropped. Do this silently.
    lowered_outer, lowered_inner = parse_traitsfunction(env, lowering, expr_original, on_traits_dropped = msg -> nothing)
    # we don't document lowerings
    lowered_exprs = merge_and_render_update(env, lowered_outer, lowered_inner, doc = false)
    append!(exprs, lowered_exprs)
  end
  # return nothing in order to not return implementation detail
  flatten_blocks(Expr(:block, exprs..., nothing))
end

function _traits_parsed(env, parsed, expr_original)
  throw(ArgumentError("@traits macro expects function expression, got `$expr_original`"))
end

end # module


text:
using AbnormalReturns
using Documenter

Documenter.makedocs(
    modules = [AbnormalReturns],
    sitename = "AbnormalReturns.jl",
    pages = [
        "Introduction" => "index.md",
        "Example" => "example.md",
        "API" => "api.md"
    ]
)

deploydocs(
    repo = "github.com/junder873/AbnormalReturns.jl.git",
    target = "build",
)


text:
# phonopy_project.jl
# Demonstration / scratchpad for using JuliaPhonons module

push!(LOAD_PATH,"../src") # Temporary versions of modules in PWD
using JuliaPhonons

# The REPO comes with a basic MAPI cubic unit cell Phonopy calculation
P=read_POSCAR(open("POSCAR"),expansion=[2,2,2])
eigenvectors,eigenmodes=read_meshyaml(open("mesh.yaml"),P)

myf=open("mode_decomposition.dat","w")
gnuplot_header(P,f=myf)

for (count,(eigenvector,eigenmode)) in enumerate(zip(eigenvectors,eigenmodes))
    output_animated_xyz(P,count,eigenvector,eigenmode) # generates files anim_{count}.xyz

    decompose_eigenmode_atomtype(P,count,eigenvector,eigenmode)
    # Output again, to a file for later printing
    decompose_eigenmode_atomtype(P,count,eigenvector,eigenmode,f=myf)

    decompose_eigenmode_atom_contributions(P,count,eigenvector)
end

close(myf)

max_stars_repo_name: jarvist/Julia-Phonons

max_issues_repo_name: lucydot/Julia-Phonons

max_forks_repo_name: lucydot/Julia-Phonons


text:
using SafeTestsets

@safetestset "Probability Objects" begin
    using InformationGeometry, Test, LinearAlgebra, Distributions

    DS = DataSet([0,0.5,1,1.5],[1.,3.,7.,8.1],[1.2,2.,0.6,1.])
    model(x,θ) = θ[1] * x + θ[2]
    DM = DataModel(DS,model)
    p = rand(2)

    @test IsLinear(DM)
    Dist = DataDist(ydata(DM),ysigma(DM))
    @test abs(loglikelihood(DM,p) - logpdf(Dist,EmbeddingMap(DM,p))) < 1e-13
    @test Score(DM,p) ≈ transpose(EmbeddingMatrix(DM,p)) * gradlogpdf(Dist,EmbeddingMap(DM,p))
    @test FisherMetric(DM,p) ≈ transpose(EmbeddingMatrix(DM,p)) * inv(cov(Dist)) * EmbeddingMatrix(DM,p)

    # Test AD vs manual derivative
    @test norm(AutoScore(DM,p) - Score(DM,p)) < 2e-13
    @test norm(AutoMetric(DM,p) .- FisherMetric(DM,p), 1) < 2e-9

    # Do these tests in higher dimensions, check that OrthVF(PL) IsOnPlane....
    # @test OrthVF(DM,XYPlane,p) == OrthVF(DM,p)

    @test dot(OrthVF(DM,p),Score(DM,p)) < 1e-14
    @test norm(FindMLE(DM) - [5.01511545953636, 1.4629658803705]) < 5e-10
end

@safetestset "Confidence Regions" begin
    using InformationGeometry, Test, Plots

    DS = DataSet([0,0.5,1,1.5],[1.,3.,7.,8.1],[1.2,2.,0.6,1.])
    model(x,θ) = θ[1] * x + θ[2];    DM = DataModel(DS,model)
    DME = DataModel(DataSetExact([0,0.5,1,1.5],0.1ones(4),[1.,3.,7.,8.1],[1.2,2.,0.6,1.]), model)

    sols = ConfidenceRegions(DM,1:2; tol=1e-6)
    @test IsStructurallyIdentifiable(DM,sols[1]) == true
    @test size(SaveConfidence(sols,50)) == (50,4)
    @test size(SaveGeodesics(sols,50)) == (50,2)
    @test size(SaveDataSet(DM)) == (4,3)
    @test ConfidenceRegionVolume(DM,sols[1];N=5000) < ConfidenceRegionVolume(DM,sols[2];N=5000,WE=true)

    @test size(ConfidenceBands(DM,sols[1]; N=50, plot=false)) == (50,3)
    @test size(PlotMatrix(inv(FisherMetric(DM,MLE(DM))),MLE(DM); N=50,plot=false)) == (50,2)
    @test typeof(FittedPlot(DM)) <: Plots.Plot
    @test typeof(FittedPlot(DME)) <: Plots.Plot
    @test typeof(ResidualPlot(DM)) <: Plots.Plot

    @test typeof(VisualizeGeos([MBAM(DM)])) <: Plots.Plot
    simplermodel(x,p) = p[1]*x;    DMSimp = DataModel(DS,simplermodel)
    @test length(ConfidenceRegion(DMSimp,1.)) == 2
    @test ModelComparison(DM,DMSimp)[2] > 0.

    @test FindFBoundary(DM,1)[1] - FindConfBoundary(DM,1)[1] > 0
end

@safetestset "More Boundary tests" begin
    using InformationGeometry, Test, Random, Distributions, OrdinaryDiffEq, LinearAlgebra

    Random.seed!(31415);    normerr(sig::Number) = rand(Normal(0,sig));     quarticlin(x,p) = p[1]*x.^4 .+ p[2]
    X = collect(0:0.2:3);   err = 2. .+ 2sqrt.(X);      Y = quarticlin(X,[1,8.]) + normerr.(err)
    ToyDME = DataModel(DataSetExact(X,0.1ones(length(X)),Y,err), (x,p) -> 15p[1]^3 * x.^4 .+ p[2]^5)

    @test InterruptedConfidenceRegion(BigFloat(ToyDME), 8.5; tol=1e-5) isa ODESolution

    NewX, NewP = TotalLeastSquares(ToyDME)
    @test LogLike(Data(ToyDME), NewX, EmbeddingMap(Data(ToyDME),Predictor(ToyDME),NewP,NewX)) > loglikelihood(ToyDME, MLE(ToyDME))

    @test ModelMap(Predictor(ToyDME), PositiveDomain(2)) isa ModelMap

    sol = ConfidenceRegion(ToyDME,1; tol=1e-6)
    @test ApproxInRegion(sol, MLE(ToyDME)) && !ApproxInRegion(sol, sol.u[1] + 1e-5BasisVector(1,2))

    #Check that bounding box from ProfileLikelihood coincides roughly with exact box.
    Mats = ProfileLikelihood(ToyDME,2; plot=false)
    ProfBox = ProfileBox(ToyDME, InterpolatedProfiles(Mats),1)
    ExactBox = ConstructCube(ConfidenceRegion(ToyDME,1; tol=1e-6))
    @test norm(Center(ProfBox) - Center(ExactBox)) < 3e-5
    @test norm(CubeWidths(ProfBox) - CubeWidths(ExactBox)) < 3e-4
    @test 0 < PracticallyIdentifiable(Mats) < 2

    # Method for general cost functions on 2D domains
    sol = GenerateBoundary(x->-norm(x,1.5), [1., 0])
    @test 0.23 ≤ length(GenerateBoundary(x->-norm(x,1.5), [1., 0]; Boundaries=(u,t,int)->u[1]<0.).u) / length(sol.u) ≤ 0.27
end

@safetestset "ODE-based models" begin
    using InformationGeometry, Test, OrdinaryDiffEq, LinearAlgebra

    function SIR!(du,u,p,t)
        S, I, R = u
        β, γ = p
        du[1] = -β * I * S
        du[2] = +β * I * S - γ * I
        du[3] = +γ * I
        nothing
    end
    SIRsys = ODEFunction(SIR!)
    infected = [3, 8, 28, 75, 221, 291, 255, 235, 190, 126, 70, 28, 12, 5]
    SIRDS = InformNames(DataSet(collect(1:14), infected, 5ones(14)), ["Days"], ["Infected"])

    SIRinitial = X->([763.0-X[1], X[1], 0.0], X[2:3])

    # Use SplitterFunction SIRinitial to infer initial condition I₀ as first parameter
    SIRDM = DataModel(SIRDS, SIRsys, SIRinitial, x->x[2], [0.6,0.0023,0.46]; tol=1e-6)
    @test SIRDM isa DataModel
    @test DataModel(SIRDS, SIRsys, [762, 1, 0.], [2], [0.0022,0.45], true; meth=Tsit5(), tol=1e-6) isa DataModel
    @test norm(2*EmbeddingMap(Data(SIRDM), Predictor(SIRDM), MLE(SIRDM)) - EmbeddingMap(Data(SIRDM), ModifyODEmodel(SIRDM, x->2*x[2]), MLE(SIRDM))) < 2e-4
end

@safetestset "Model Transformations" begin
    using InformationGeometry, Test

    PiDM = DataModel(DataSet([0,1], [0.5π,1.5π], [0.5,0.5]), ModelMap((x,p)->p[1], θ->θ[1]-1, HyperCube([[0,5]])))
    @test !IsInDomain(Predictor(PiDM), [0.9]) && IsInDomain(Predictor(PiDM), [1.1])

    # Translation
    PiDM2 = DataModel(Data(PiDM), TranslationTransform(Predictor(PiDM),[1.]))
    @test !IsInDomain(Predictor(PiDM2), [-0.1]) && IsInDomain(Predictor(PiDM2), [0.1])

    # LogTransform
    PiDM3 = DataModel(Data(PiDM), LogTransform(Predictor(PiDM),trues(1)))
    @test !IsInDomain(Predictor(PiDM3), exp.([1])-[0.1]) && IsInDomain(Predictor(PiDM3), exp.([1])+[0.1])

    DS = DataSet([0,0.5,1,1.5],[1.,3.,7.,8.1],[1.2,2.,0.6,1.])
    @test FisherMetric(LinearDecorrelation(DataModel(DS, (x,θ)->θ[1] * x + θ[2])), zeros(2)) ≈ [1 0; 0 1]

    # TranstrumModel = ModelMap((x::Real,p::AbstractVector)->exp(-p[1]*x) + exp(-p[2]*x), θ::AbstractVector -> θ[1]>θ[2], PositiveDomain(2, 1e2), (1,1,2))
    # TranstrumDM = DataModel(DataSet([0.33, 1, 3], [0.88,0.5,0.35], [0.1,0.3,0.2]), TranstrumModel)
    # linTranstrum = LogTransform(TranstrumDM)
    # RicciScalar(TranstrumDM, MLE(TranstrumDM)), RicciScalar(linTranstrum, MLE(linTranstrum))
    # loglikelihood(TranstrumDM, MLE(TranstrumDM)), loglikelihood(linTranstrum, MLE(linTranstrum))

    # Try with normal functions too, not only ModelMaps.
    # Try Ricci in particular, maybe as BigFloat.

    # Does Score / FisherMetric and AutoDiff still work?
end

@safetestset "In-place ModelMaps" begin
    using InformationGeometry, Test, LinearAlgebra

    DM = DataModel(DataSet([1,2,3],[4,1,5,2,6.5,3.5],[0.5,0.5,0.45,0.45,0.6,0.6], (3,1,2)), (x,p)-> [p[1]^3*x, p[2]^2*x])
    dm = InplaceDM(DM)

    @test EmbeddingMap(DM, MLE(DM)) ≈ EmbeddingMap(dm, MLE(dm))
    @test EmbeddingMatrix(DM,MLE(DM)) ≈ EmbeddingMatrix(dm,MLE(dm))
    @test Score(DM, MLE(DM)) ≈ Score(dm, MLE(dm))
    @test FisherMetric(DM, MLE(DM)) ≈ FisherMetric(dm, MLE(dm))
end

@safetestset "Inputting Datasets of various shapes" begin
    using InformationGeometry, Test, LinearAlgebra, Random, Distributions, StaticArrays, Plots

    ycovtrue = Diagonal([1,2,3]) |> x->convert(Matrix,x)
    ptrue = [1.,π,-5.];        ErrorDistTrue = MvNormal(zeros(3),ycovtrue)

    model(x::AbstractVector{<:Number},p::AbstractVector{<:Number}) = SA[p[1] * x[1]^2 + p[3]^3 * x[2],
                                                        sinh(p[2]) * (x[1] + x[2]), exp(p[1]*x[1] + p[1]*x[2])]
    Gen(t) = float.([t,0.5t^2]);    Xdata = Gen.(0.5:0.1:3)
    Ydata = [model(x,ptrue) + rand(ErrorDistTrue) for x in Xdata]
    Sig = BlockMatrix(ycovtrue,length(Ydata));    DS = DataSet(Xdata,Ydata,Sig)
    DM = DataModel(DS,model)
    @test norm(MLE(DM) - ptrue) < 5e-2
    DME = DataModel(DataSetExact(DS), model)
    P = MLE(DM) + 0.5rand(length(MLE(DM)))
    @test loglikelihood(DM,P) ≈ loglikelihood(DME,P)
    @test Score(DM,P) ≈ Score(DME,P)

    Planes, sols = ConfidenceRegion(DM,1)
    @test typeof(VisualizeSols(Planes,sols)) <: Plots.Plot

    ODM = OptimizedDM(DME)
    @test norm(EmbeddingMatrix(DME,MLE(DME)) .- EmbeddingMatrix(ODM,MLE(DME)), 1) < 1e-9

    CDM = DataModel(CompositeDataSet(Data(ODM)), Predictor(ODM), dPredictor(ODM), MLE(ODM))
    @test abs(loglikelihood(ODM, P) - loglikelihood(CDM, P)) < 5e-6
    @test norm(Score(ODM, P) - Score(CDM, P)) < 2e-6
    @test norm(FisherMetric(ODM, P) - FisherMetric(CDM, P)) < 2e-6
    @test norm(InformationGeometry.ResidualStandardError(ODM) - InformationGeometry.ResidualStandardError(CDM)) < 1e-10

    lastDS = Data(Data(CDM))[3]
    newCDS = vcat(Data(Data(CDM))[1:end-1], [SubDataSet(lastDS, 1:2:Npoints(lastDS))], [SubDataSet(lastDS, 2:2:Npoints(lastDS))]) |> CompositeDataSet
    # repeat last component
    newmodel(x::AbstractVector{<:Number},p::AbstractVector{<:Number}) = SA[p[1] * x[1]^2 + p[3]^3 * x[2], sinh(p[2]) * (x[1] + x[2]),
                                                        exp(p[1]*x[1] + p[1]*x[2]), exp(p[1]*x[1] + p[1]*x[2])]
    splitCDM = DataModel(newCDS, newmodel, MLE(CDM))
    @test abs(loglikelihood(splitCDM, P) - loglikelihood(CDM, P)) < 1e-5
    @test norm(Score(splitCDM, P) - Score(CDM, P)) < 2e-4
    @test norm(FisherMetric(splitCDM, P) - FisherMetric(CDM, P)) < 2e-3
end

@safetestset "Priors" begin
    using InformationGeometry, Test, LinearAlgebra, Distributions

    DS1 = DataSet([0,0.5],[1.,3.],[1.2,2.]);    DS2 = DataSet([1,1.5],[7.,8.1],[0.6,1.])
    DS = join(DS1, DS2);    model(x,θ) = θ[1] * x + θ[2];
    DM1 = DataModel(DS1,model);     DM = DataModel(DS,model);

    logprior(X) = loglikelihood(DM1, X)

    DM12 = DataModel(DS2, model, MLE(DM1), logprior)
    @test norm(MLE(DM) - MLE(DM12)) < 1e-6
    @test loglikelihood(DM, MLE(DM)) ≈ loglikelihood(DM12, MLE(DM))
    @test norm(Score(DM, MLE(DM)) - Score(DM12, MLE(DM))) < 1e-12
    @test FisherMetric(DM, MLE(DM)) ≈ FisherMetric(DM12, MLE(DM))

    # Intentionally not giving information about LogPrior input here
    dm = DataModel(DS, model, MLE(DM), x->0.0)
    @test loglikelihood(DM, MLE(DM)) ≈ loglikelihood(dm, MLE(DM))
    @test Score(DM, MLE(DM)) ≈ Score(dm, MLE(DM))
    @test FisherMetric(DM, MLE(DM)) ≈ FisherMetric(dm, MLE(DM))
end

@safetestset "Kullback-Leibler Divergences" begin
    using InformationGeometry, Test, LinearAlgebra, Distributions

    # Analytical divergences via types defined in Distributions.jl
    @test KullbackLeibler(MvNormal([1,2,3.],diagm([2,4,5.])),MvNormal([1,-2,-3.],diagm([1,5,2.]))) ≈ 11.056852819440055
    @test KullbackLeibler(Normal(1,2),Normal(5,1)) ≈ 8.806852819440055
    @test KullbackLeibler(Cauchy(1,3),Cauchy(5,2)) ≈ 0.5355182363563621
    @test KullbackLeibler(Exponential(5),Exponential(1)) ≈ 2.3905620875658995
    @test KullbackLeibler(Weibull(12,2),Weibull(3,4)) ≈ 2.146124463755512
    @test KullbackLeibler(Gamma(5,15),Gamma(20,1)) ≈ 29.409061308330323

    # Numerically calculated for via arbitrary types defined in Distributions.jl
    # ALSO ADD TESTS FOR DISCRETE DISTRIBUTIONS, DISTRIBUTIONS WITH LIMITED DOMAIN
    @test abs(KullbackLeibler(Cauchy(1,2),Normal(5,1),HyperCube([-20,20]); tol=1e-8) - 16.77645704773449) < 1e-5
    @test abs(KullbackLeibler(Cauchy(1,2),Normal(5,1),HyperCube([-20,20]); Carlo=true, N=Int(3e6)) - 16.7764) < 5e-2
    @test abs(KullbackLeibler(MvTDist(1,[3,2,1.],diagm([1.,2.,3.])),MvNormal([1,2,3.],diagm([2,4,5.])),HyperCube([[-10,10.] for i in 1:3]); Carlo=true, N=Int(3e6)) - 1.6559288) < 3e-1

    # Product distributions, particularly Normal and Cauchy
    P = [Normal(0,1), Cauchy(1,2)] |> product_distribution
    Q = [Cauchy(1,1), Cauchy(2,4)] |> product_distribution
    R = [Normal(2,4), Normal(-1,0.5)] |> product_distribution
    @test abs(KullbackLeibler(P, Q, HyperCube([[-20,20] for i in 1:2]); tol=1e-7) - 0.719771180) < 1e-8
    @test abs(KullbackLeibler(R, P, HyperCube([[-20,20] for i in 1:2]); tol=1e-7) - 9.920379769) < 1e-8
    @test abs(KullbackLeibler(P, R, HyperCube([[-20,20] for i in 1:2]); tol=1e-7) - 48.99179438) < 1e-8

    # Via any positive (hopefully normalized) functions
    @test abs(KullbackLeibler(x->pdf(Normal(1,3),x),y->pdf(Normal(5,2),y),HyperCube([-20,20]); Carlo=true, N=Int(3e6)) - KullbackLeibler(Normal(1,3),Normal(5,2))) < 2e-2
    @test abs(KullbackLeibler(x->pdf(Normal(1,3),x),y->pdf(Normal(5,2),y),HyperCube([-20,20]); tol=1e-8) - KullbackLeibler(Normal(1,3),Normal(5,2))) < 1e-5
    P = MvNormal([1,2,3.],diagm([1,2,1.5]));    Q = MvNormal([1,-2,-3.],diagm([2,1.5,1.]));     Cube = HyperCube([[-15,15] for i in 1:3])
    @test abs(KullbackLeibler(x->pdf(P,x),y->pdf(Q,y),Cube; Carlo=true, N=Int(3e6)) - KullbackLeibler(x->pdf(P,x),y->pdf(Q,y),Cube; tol=1e-8)) < 0.8
end

@safetestset "Differential Geometry" begin
    using InformationGeometry, Test, LinearAlgebra, StaticArrays

    S2metric((θ,ϕ)) = [1.0 0; 0 sin(θ)^2]
    function S2Christoffel((θ,ϕ))
        Symbol = zeros(suff(ϕ),2,2,2);    Symbol[1,2,2] = -sin(θ)*cos(θ)
        Symbol[2,1,2] = Symbol[2,2,1] = cos(θ)/sin(θ);  Symbol
    end
    # Calculation by hand works out such that in this special case:
    S2Ricci(x) = S2metric(x)
    ConstMetric(x) = Diagonal(ones(2))

    # Test Numeric Christoffel Symbols, Riemann and Ricci tensors, Ricci Scalar
    # Test WITH AND WITHOUT BIGFLOAT
    x = rand(2)
    @test norm(ChristoffelSymbol(S2metric,x) .- S2Christoffel(x), 1) < 5e-9
    @test norm(ChristoffelSymbol(S2metric,BigFloat.(x)) .- S2Christoffel(BigFloat.(x)), 1) < 1e-40

    @test abs(RicciScalar(S2metric,x) - 2) < 5e-4
    @test abs(RicciScalar(S2metric,BigFloat.(x)) - 2) < 2e-22

    # Use wilder metric and test AutoDiff vs Finite
    import InformationGeometry: MetricPartials, ChristoffelPartials
    Y = rand(3)
    Metric3(x) = [sinh(x[3]) exp(x[1])*sin(x[2]) 0; 0 cosh(x[2]) cos(x[2])*x[3]*x[2]; exp(x[2]) cos(x[3])*x[1]*x[2] 0.]
    @test MetricPartials(Metric3, Y; ADmode=Val(true)) ≈ MetricPartials(Metric3, Y; ADmode=Val(false))
    @test ChristoffelSymbol(Metric3, Y; ADmode=Val(true)) ≈ ChristoffelSymbol(Metric3, Y; ADmode=Val(false))
    @test maximum(abs.(ChristoffelPartials(Metric3, Y; ADmode=Val(true)) - ChristoffelPartials(Metric3, Y; ADmode=Val(false), BigCalc=true))) < 1e-11
    @test maximum(abs.(Riemann(Metric3, Y; ADmode=Val(true)) - Riemann(Metric3, Y; ADmode=Val(false), BigCalc=true))) < 1e-11
    # Test with static arrays
    Metric3SA(x) = SA[sinh(x[3]) exp(x[1])*sin(x[2]) 0; 0 cosh(x[2]) cos(x[2])*x[3]*x[2]; exp(x[2]) cos(x[3])*x[1]*x[2] 0.]
    @test MetricPartials(Metric3SA, Y; ADmode=Val(true)) ≈ MetricPartials(Metric3SA, Y; ADmode=Val(false))
    @test ChristoffelSymbol(Metric3SA, Y; ADmode=Val(true)) ≈ ChristoffelSymbol(Metric3SA, Y; ADmode=Val(false))
    @test maximum(abs.(ChristoffelPartials(Metric3SA, Y; ADmode=Val(true)) - ChristoffelPartials(Metric3SA, Y; ADmode=Val(false), BigCalc=true))) < 1e-11
    @test maximum(abs.(Riemann(Metric3SA, Y; ADmode=Val(true)) - Riemann(Metric3SA, Y; ADmode=Val(false), BigCalc=true))) < 1e-11
    # Test with BigFloat
    @test -45 > MetricPartials(Metric3SA, BigFloat.(Y); ADmode=Val(true)) - MetricPartials(Metric3SA, BigFloat.(Y); ADmode=Val(false)) |> maximum |> log10 |> Float64
    @test -45 > ChristoffelSymbol(Metric3SA, BigFloat.(Y); ADmode=Val(true)) - ChristoffelSymbol(Metric3SA, BigFloat.(Y); ADmode=Val(false)) |> maximum |> log10 |> Float64
    @test -20 > ChristoffelPartials(Metric3SA, BigFloat.(Y); ADmode=Val(true)) - ChristoffelPartials(Metric3SA, BigFloat.(Y); ADmode=Val(false)) |> maximum |> log10 |> Float64
    @test -20 > Riemann(Metric3SA, BigFloat.(Y); ADmode=Val(true)) - Riemann(Metric3SA, BigFloat.(Y); ADmode=Val(false)) |> maximum |> log10 |> Float64

    @test abs(GeodesicDistance(ConstMetric,[0,0],[1,1]) - sqrt(2)) < 2e-8
    @test abs(GeodesicDistance(S2metric,[π/4,1],[3π/4,1]) - π/2) < 1e-9
    @test abs(GeodesicDistance(S2metric,[π/2,0],[π/2,π/2]) - π/2) < 1e-8

    DS = DataSet([0,0.5,1],[1.,3.,7.],[1.2,2.,0.6]);    DM = DataModel(DS, (x,p) -> p[1]^3 *x + p[2]^3)
    y = MLE(DM) + 0.2(rand(2) .- 0.5)
    geo = GeodesicBetween(DM, MLE(DM), y; tol=1e-11)
    @test norm(MLE(DM) - [1.829289173660125,0.942865200406147]) < 1e-7

    Len = GeodesicLength(DM,geo)
    @test abs(InformationGeometry.ParamVol(geo) * InformationGeometry.GeodesicEnergy(DM,geo) - Len^2) < 1e-8
    Confnum = InvConfVol(ChisqCDF(pdim(DM), 2*(LogLikeMLE(DM) - loglikelihood(DM, y))))
    @test InformationGeometry.GeodesicRadius(DM, Confnum) - Len < 1e-5

    # Apply logarithmic map first since it is typically multi-valued for positively curved manifolds.
    @test norm(ExponentialMap(FisherMetric(DM), MLE(DM), LogarithmicMap(FisherMetric(DM), MLE(DM), y)) - y) < 1

end

@safetestset "Numerical Helper Functions" begin
    using InformationGeometry, Test, BenchmarkTools, ForwardDiff

    # Compare Integrate1D and IntegrateND

    # Test integration, differentiation, Monte Carlo, GeodesicLength
    # TEST WITH AND WITHOUT BIGFLOAT
    @test abs(InformationGeometry.MonteCarloArea(x->((x[1]^2 + x[2]^2) < 1), HyperCube([[-1,1],[-1,1]])) - π) < 1.5e-3
    @test abs(Integrate1D(cos, (0,π/2); tol=1e-12) - IntegrateND(cos, (0,π/2); tol=1e-12)) < 1e-10
    z = 3rand()
    @test abs(Integrate1D(x->2/sqrt(π) * exp(-x^2), [0,z/sqrt(2)]) - ConfVol(z)) < 1e-12
    @test abs(LineSearch(x->(x < BigFloat(π))) - π) < 1e-14
    @test abs(LineSearch(x->(x < BigFloat(π)), BigFloat(1e-14); tol=1e-30) - BigFloat(π)) < 1e-25
    @test abs(CubeVol(TranslateCube(HyperCube([[0,1],[0,π],[-sqrt(2),0]]),rand(3))) - sqrt(2)*π) < 3e-15

    k = rand(1:20);     r = 10rand()
    @test InvChisqCDF(k,Float64(ChisqCDF(k,r))) ≈ r
    @test abs(InvChisqCDF(k,ChisqCDF(k,BigFloat(r)); tol=1e-20) - r) < 1e-18
end


text:
#!/usr/bin/env julia
using Printf
using Distributed
using TimerOutputs

function ζ(N,a)
  s = 0.0
  e = 1.0*a

  VLEN      = 512   #vector length

  __P_SUM   = zeros(Float64,VLEN)

  __DEC     = zeros(Float64,VLEN)
  __D_POW	  = zeros(Float64,VLEN)
  __L 		  = zeros(Float64,VLEN)

  start_index = N - (N % VLEN);
  end_index = VLEN;
  #@printf("start = %i, end = %i\n",start_index,end_index)
  #@printf("N  = %i\n",N)

  for k in 1:VLEN
    __DEC[k]	= convert(Float64,VLEN)
    __L[k]		= convert(Float64,start_index - (k - 1));
  end

  #@printf("(pre warmup) s=%f\n",s)
  # warm up loop to handle non-vector size clean up before main loop
  if N == start_index
    s = 0.0
  else
   for i ∈ N:-1:start_index+1
    s += (1.0/i)^e
   end
  end
  #@printf("(post warmup) s=%f\n",s)

  # main loop
  for k in start_index:-VLEN:end_index
    __D_POW	  = (1.0 ./ __L).^e
    __P_SUM	 += __D_POW
    __L      -= __DEC
  end

  # clean up loop to handle post vector loop calc

  #for i in end_index:-1:1
  # s += (1.0/i)^e
  #end

  # sum reduction
  for i in 1:VLEN
       s+=__P_SUM[i]
  end

  s
end

to = TimerOutput()

N = 16000000000

st = 1
s =  zeros(Float64,32)
pd = zeros(Float64,32)
ndx = 1
#ENV["JULIA_DEBUG"] = "CUDAnative"

s[1] = ζ(N,2)

for i=1:1
  @printf("s[%i] = %f\n",i,s[i])
end

@timeit to "ζ"       pd[ndx] = ζ(N,2)
ndx+=1

for i=1:1
  @printf("pd[%i] = %f\n",i,pd[i])
end

show(to)
@printf("\n");

max_stars_repo_name: joelandman/rzf

max_issues_repo_name: joelandman/rzf

max_forks_repo_name: joelandman/rzf


text:
using BinDeps
using CMakeWrapper

@BinDeps.setup

function cflags_validator(pkg_names...)
    return (name, handle) -> begin
        for pkg_name in pkg_names
            try
                run(`pkg-config --cflags $(pkg_name)`)
                return true
            catch ErrorException
            end
        end
        false
    end
end

"""
Director's libraries are unversioned, so there's no possible way
to know if a system-installed version of Director is compatible
with this interface. Instead, we have to restrict ourselves to
only versions of Director which are built locally.
"""
function is_local_build(name, handle)
    startswith(relpath(name, @__DIR__), "usr/")
end

basedir = dirname(@__FILE__)
director_version = "0.1.0-234-g74cea84"
director_sha = "02c2ef65f8d1d9f3de1d56d129351cd43846d70b"

@static if Sys.islinux()
    python = library_dependency("python", aliases=["libpython2.7.so",], validate=cflags_validator("python", "python2"))
    qt4 = library_dependency("QtCore", aliases=["libQtCore.so", "libQtCore.so.4.8"], depends=[python])
    qt4_opengl = library_dependency("QtOpenGL", aliases=["libQtOpenGL.so", "libQtOpenGL.so.4.8"], depends=[qt4])
    director = library_dependency("ddApp", aliases=["libddApp"], depends=[python, qt4, qt4_opengl], validate=is_local_build)

    linux_distributor = strip(read(`lsb_release -i -s`, String))
    linux_version = try
        VersionNumber(strip(read(`lsb_release -r -s`, String)))
    catch e
        if isa(e, ArgumentError)
            v"0.0.0"
        else
            rethrow(e)
        end
    end

    provides(AptGet, Dict("libqt4-dev"=>qt4, "libqt4-opengl-dev"=>qt4_opengl, "python-dev"=>python))

    force_source_build = lowercase(get(ENV, "DIRECTOR_BUILD_FROM_SOURCE", "")) in ["true", "1"]

    director_binary = nothing
    if !force_source_build
        if linux_distributor == "Ubuntu"
            if linux_version >= v"16.04"
                director_binary = "ubuntu-16.04"
            elseif linux_version >= v"14.04"
                director_binary = "ubuntu-14.04"
            end
        elseif linux_distributor == "Debian"
            if linux_version >= v"8.7"
                director_binary = "ubuntu-14.04"
            end
        end
    end

    if director_binary !== nothing
        provides(BuildProcess, (@build_steps begin
            FileDownloader("https://dl.bintray.com/patmarion/director/director-$(director_version)-$(director_binary).tar.gz",
                           joinpath(basedir, "downloads", "director.tar.gz"))
            CreateDirectory(joinpath(basedir, "usr"))
            (`tar xzf $(joinpath(basedir, "downloads", "director.tar.gz")) --directory=usr --strip-components=1`)
        end), director)
    else
        provides(Sources,
                 URI("https://github.com/RobotLocomotion/director/archive/$(director_sha).zip"),
                 director,
                 unpacked_dir="director-$(director_sha)")
        provides(CMakeProcess(srcdir=joinpath(basedir, "src",
                                              "director-$(director_sha)", "distro", "superbuild"),
                              cmake_args=["-DUSE_LCM=ON",
                                          "-DUSE_EXTERNAL_INSTALL:BOOL=ON"],
                              targetname=""),
                 director)
    end

elseif Sys.isapple()
    # Use the libvtkDRCFilters library instead of libddApp
    # to work around weird segfault when dlclose()-ing libddApp
    director = library_dependency("vtkDRCFilters", aliases=["libvtkDRCFilters.dylib"], validate=is_local_build)
    provides(BuildProcess, (@build_steps begin
        FileDownloader("https://dl.bintray.com/patmarion/director/director-$(director_version)-mac.tar.gz",
                       joinpath(basedir, "downloads", "director.tar.gz"))
        CreateDirectory(joinpath(basedir, "usr"))
        (`tar xzf $(joinpath(basedir, "downloads", "director.tar.gz")) --directory=usr --strip-components=1`)
    end), director)

end

@BinDeps.install Dict(:ddApp => :libddApp)


text:
@doc raw"""
    cronbach(Σ::AbstractMatrix)

Calculate [Cronbach's alpha](https://en.wikipedia.org/wiki/Cronbach%27s_alpha).
This is also known as tau-equivalent reliability (``\rho_T``).
`Σ` is a covariance matrix (which means it must be square).

# Examples
```jldoctest
julia> C = [10  6  6  6;  # fictitious data
             6 11  6  6;
             6  6 12  6;
             6  6  6 13];

julia> cronbach(C)
0.8135593220338984
```
"""
function cronbach(Σ::AbstractMatrix)
	k, m = size(Σ)
    k == m || throw(DimensionMismatch("dimensions must match"))

    σᵢⱼ = zero(eltype(Σ))
    @inbounds for i in 2:k, j in 1:i-1
        σᵢⱼ += Σ[i, j]
    end
    σₓ² = 2σᵢⱼ + sum(diag(Σ))
    σᵢⱼ /= (k - 1) * k / 2

    return k^2 * σᵢⱼ / σₓ²
end


text:
nodes(g::Graph) = keys(fadj(g))

out_neighbor(g::Graph, u) = keys(fadj(g)[u])
in_neighbor(g::Graph, u) = keys(badj(g)[u])

has_node(g::Graph, u) = u in nodes(g)


text:
#!/usr/bin/env julia

using Test

# write your own tests here
@test 1 == 1

using CORBITS

function test_corbits()
  #@compat is_windows() ? error("# CORBITS won't work on windows") : nothing
  a =  Cdouble[0.05, 0.15]
  r_star = convert(Cdouble,0.005)
  r = Cdouble[0.0001,0.0001]
  ecc = Cdouble[0.02, 0.1]
  Omega = Cdouble[0.0, 0.0]
  omega = Cdouble[ 0.0, 0.5]
  inc = Cdouble[pi/2, pi/2]
  use_pl = Cint[1,1]
  prob_of_transits_approx(a, r_star, r, ecc, Omega, omega, inc, use_pl)
end

isapprox(test_corbits(), 0.03367003367003367, atol= 0.0001)


text:
using Test
import ParallelStencil
using ParallelStencil.ParallelKernel
import ParallelStencil.ParallelKernel: @reset_parallel_kernel, @is_initialized, @get_package, @get_numbertype, SUPPORTED_PACKAGES, PKG_CUDA, PKG_NONE, NUMBERTYPE_NONE
import ParallelStencil.ParallelKernel: @require, @symbols, longnameof
TEST_PACKAGES = SUPPORTED_PACKAGES
@static if PKG_CUDA in TEST_PACKAGES
    import CUDA
    if !CUDA.functional() TEST_PACKAGES = filter!(x->x≠PKG_CUDA, TEST_PACKAGES) end
end

@static for package in TEST_PACKAGES  eval(:(
    @testset "$(basename(@__FILE__)) (package: $(nameof($package)))" begin
        @testset "1. Reset of ParallelKernel" begin
            @testset "Reset if not initialized" begin
                @require !@is_initialized()
                @reset_parallel_kernel()
                @test !@is_initialized()
                @test @get_package() == $PKG_NONE
                @test @get_numbertype() == $NUMBERTYPE_NONE
            end;
            @testset "Reset if initialized" begin
                @require !@is_initialized()
                @init_parallel_kernel($package, Float64)
                @require @is_initialized() && @get_package() == $package
                @reset_parallel_kernel()
                @test length(@symbols($(@__MODULE__), Data)) == 1
                @test !@is_initialized()
                @test @get_package() == $PKG_NONE
                @test @get_numbertype() == $NUMBERTYPE_NONE
            end;
        end;
    end;
)) end == nothing || true;


text:
# Starting point for conditional example for Learning Julia

# Create an if-else conditional
x = 17
if x < 10
    println("x is small")
else
    println("x is big")
end

# multiple conditions can be specified with if-elseif-else
if x < 10
    println("x is small")
elseif x >= 10 && x < 25
    println("x is medium")
else
    println("x is big")
end

# The ternary operator can condense a comparison
println(x < 10 ? "x is less than 10" : "x is 10 or greater")
# ternary operator condense if-else conditional statement
# if the condition is true, execute the portion between ? and :
# if the condition is false, execute the portion after :

max_stars_repo_name: scaotravis/Julia-Course

max_issues_repo_name: scaotravis/Julia-Course

max_forks_repo_name: scaotravis/Julia-Course


text:
import ClausenFunctions

n = 10_000_000
x_min = 0.0
x_max = pi
data = (x_max - x_min)*rand(Float64, n) + x_min*ones(n)

println("Benchmarking cl1::Float64")

time_cl1(data) = @time map(ClausenFunctions.cl1, data)

map(ClausenFunctions.cl1, data) # trigger compilation
time_cl1(data)                  # trigger compilation
time_cl1(data)
time_cl1(data)

println("Benchmarking cl2::Float64")

time_cl2(data) = @time map(ClausenFunctions.cl2, data)

map(ClausenFunctions.cl2, data) # trigger compilation
time_cl2(data)                  # trigger compilation
time_cl2(data)
time_cl2(data)

println("Benchmarking cl3::Float64")

time_cl3(data) = @time map(ClausenFunctions.cl3, data)

map(ClausenFunctions.cl3, data) # trigger compilation
time_cl3(data)                  # trigger compilation
time_cl3(data)
time_cl3(data)

println("Benchmarking cl4::Float64")

time_cl4(data) = @time map(ClausenFunctions.cl4, data)

map(ClausenFunctions.cl4, data) # trigger compilation
time_cl4(data)                  # trigger compilation
time_cl4(data)
time_cl4(data)

println("Benchmarking cl5::Float64")

time_cl5(data) = @time map(ClausenFunctions.cl5, data)

map(ClausenFunctions.cl5, data) # trigger compilation
time_cl5(data)                  # trigger compilation
time_cl5(data)
time_cl5(data)

println("Benchmarking cl6::Float64")

time_cl6(data) = @time map(ClausenFunctions.cl6, data)

map(ClausenFunctions.cl6, data) # trigger compilation
time_cl6(data)                  # trigger compilation
time_cl6(data)
time_cl6(data)

println("Benchmarking cl::Float64")

time_cl(k, data) = @time map(x -> ClausenFunctions.cl(k, x), data)

n = 1_000_000
data = (x_max - x_min)*rand(Float64, n) + x_min*ones(n)

for k in vcat(collect(1:16), [1000, 1001, 1_000_000])
    println("Benchmarking cl($(k),x)::Float64")
    map(x -> ClausenFunctions.cl(k, x), data) # trigger compilation
    time_cl(k, data)            # trigger compilation
    time_cl(k, data)
    time_cl(k, data)
end

println("Benchmarking sl::Float64")

time_sl(k, data) = @time map(x -> ClausenFunctions.sl(k, x), data)

n = 1_000_000
data = (x_max - x_min)*rand(Float64, n) + x_min*ones(n)

for k in vcat(collect(1:31), [1000, 1001, 1_000_000])
    println("Benchmarking sl($(k),x)::Float64")
    map(x -> ClausenFunctions.sl(k, x), data) # trigger compilation
    time_sl(k, data)            # trigger compilation
    time_sl(k, data)
    time_sl(k, data)
end


text:
module ASTInterface

__precompile__(false)

using Distributions
using CommonRLInterface
import Base: rand

export ASTMDP,
	   Simulation,
	   Environment,
	   EnvironmentValue,
	   reset!,
	   environment,
	   observe,
	   step!,
	   isterminal,
	   isevent,
	   distance,
	   flatten,
	   unflatten

include("AST.jl")
include("interface.jl")
include("RL.jl")

end


text:
"""
$SIGNATURES

Resolve arguments for `LxObj` by interpolating `#k` appropriately.
"""
function resolve_args(base::AS, braces::Vector{OCBlock})
    res = base
    for (i, brace) in enumerate(braces)
        brace_content = stent(brace)
        # space-sensitive 'unsafe' one
        res = replace(res, "!#$i" => brace_content)
        # space-insensitive 'safe' one (e.g. for things like `\mathbb#1`)
        res = replace(res, "#$i" => " " * brace_content)
    end
    return res
end

"""
$(SIGNATURES)

Take a `<: LxObj` and try to resolve it by looking up the appropriate definition, applying it then
reparsing the result.
"""
function resolve_lxobj(lxo::LxObj, lxdefs::Vector{LxDef};
                       inmath::Bool=false)::String
    # retrieve the definition the environment points to
    lxd = getdef(lxo)
    env = lxo isa LxEnv

    # in case it's defined in Utils or in Franklin
    name = getname(lxo)
    fun  = Symbol(ifelse(env, "env_", "lx_") * name)
    # it will be `nothing` in math mode or when defined in utils
    if isnothing(lxd)
        # check if it's defined in Utils and act accordingly
        if isdefined(Main, utils_symb()) && isdefined(utils_module(), fun)
            raw = Core.eval(utils_module(), :($fun($lxo, $lxdefs)))
            return reprocess(raw, lxdefs)
        else
            # let the math backend deal with the string
            return lxo.ss
        end
    end

    # the definition can be empty (which can be on purpose, for internal defs)
    if (!env && isempty(lxd)) || (env && isempty(lxd.first) && isempty(lxd.second))
        name = getname(lxo)
        isdefined(Franklin, fun) && return eval(:($fun($lxo, $lxdefs)))
        return ""
    end

    # non-empty cases
    if !env
        partial = resolve_args(lxd, lxo.braces)
    else
        partial  = resolve_args(lxd.first, lxo.braces)
        partial *= content(lxo)
        partial *= resolve_args(lxd.second, lxo.braces)
    end

    # if 'inmath' surround accordingly so that this information is preserved
    inmath && (partial = mathenv(partial))

    return reprocess(partial, lxdefs)
end

"""
$SIGNATURES

Convenience function to take a markdown string (e.g. produced by a latex command) and re-parse it.
"""
function reprocess(s::AS, lxdefs::Vector{<:LxDef}; nostripp=false) where T
    r = convert_md(s, lxdefs;
                   isrecursive=true, isconfig=false, has_mddefs=false,
                   nostripp=nostripp)
    return r
end


text:
export PermList
export randpermlist, leastmoved, greatestmoved, supportsize, support, fixed
export pow2
export pivtopermlist

import DataStructures: list

immutable PermList{T<:Real} <: AbstractPerm{T}
    data::Vector{T}
end

setindex!(p::PermList, i::Int, k::Integer) = p.data[k] = i
getindex{T}(p::PermList{T}, k::Real) = k > length(p.data) ? convert(T,k) : (p.data)[k]
length(p::PermList) = length(p.data)

randcyclelist(n::Integer) = PermList(randcycle(n))

pivtopermlist{T<:Real}(piv::AbstractArray{T}) = PermList(PermPlain.ipiv2perm(piv))
pivtopermlist{T<:Real}(piv::AbstractArray{T},n) = PermPlain.ipiv2perm(piv,n)

## Apply permutation, and permutation operations ##

getindex{T}(v::Array{T,1},p::PermList{Bool}) = error("Silence compiler. You don't want Bool, anyway")
getindex(v::Array, p::PermList) = v[p.data] # How to define this for everything?
getindex(v::String, p::PermList) = v[p.data] # How to define this for everything?
*(p::PermList, k::Real) = k > length(p) ? k : p[k]
*{T}(p::PermList, a::AbstractVector{T}) = PermPlain.permapply(p.data,a)
# updating ops are "syntactic-operators". Don't know how to define method for them
# *=(p::PermList, q::PermList) = (permcompose!(p.data,q.data), p)
/(p::PermList, q::PermList) = PermList(PermPlain.permcompose(p.data,invperm(q.data)))
# preimage of k under p
/(k::Int, p::PermList) = PermPlain.preimage(p.data,k)

## Output ##

show(io::IO, p::PermList) = print(io,p)
show(p::PermList) = print(p)
# This is needed to avoid trying to print PermList with showarray and failing in 1000 ways
writemime(io::IO, ::MIME"text/plain", p::PermList) = print(io,p)


text:
"""
    OrbitInfo

Struct containing orbit, start time, and observer location info.

Fields:

- `start_time`: start time in UTC in the format
                (year, month, day, hour, minute, second)
- `start_time_julian_date`: start time in Julian days
- `site_loc_lla`: site latitude, longitude, and height in (deg, deg, meters)
- `site_loc_ecef`: site ECEF position in (meters, meters, meters)
- `tle_file_name`: file name of satellite TLE file
- `tle`: `TLE` struct (data loaded from TLE file)
- `orb`: single or vector of `OrbitPropagator` structs
- `eop`: Earth Orientation Parameters object
"""
mutable struct OrbitInfo{T1,T2,T3,T4,T5,T6,T7,T8}
    start_time::T1
    start_time_julian_date::T2
    site_loc_lla::T3
    site_loc_ecef::T4
    tle_file_name::T5
    tle::T6
    orb::T7
    eop::T8
end

"""
    initorbitinfo(source_tle::String, target_tle::String, start_time,
                  site_loc_lla)

Initialize the struct OrbitInfo for multiple satellites. Provide the file names
of the individual TLE files.

Required Arguments:

- `source_tle::String`: TLE string for satellite that is trasmitting signal
- `target_tle::String`: TLE string for satellite that reflectes transmitted
                        signal from source satellite
- `start_time`: start time in UTC in the format
                (year, month, day, hour, minute, second)
- `site_loc_lla`: site latitude, longitude, and height in (deg, deg, meters)

Returns:

- `OrbitInfo` struct
"""
function initorbitinfo(source_tle::String, target_tle::String, start_time,
                       site_loc_lla)
    start_time_julian_date = DatetoJD(start_time...)
    lat = deg2rad(site_loc_lla[1])
    lon = deg2rad(site_loc_lla[2])
    height = site_loc_lla[3]
    site_loc_ecef = GeodetictoECEF(lat, lon, height)
    tle = [read_tle(source_tle)[1], read_tle(target_tle)[1]]
    orb = [init_orbit_propagator(Val(:sgp4), tle[1], sgp4_gc_wgs84),
           init_orbit_propagator(Val(:sgp4), tle[2], sgp4_gc_wgs84)]
    eop = get_iers_eop(:IAU1980)
    tle_names = [source_tle, target_tle]
    return OrbitInfo(start_time, start_time_julian_date,
                     site_loc_lla, site_loc_ecef,
                     tle_names, tle, orb, eop)
end

"""
    initorbitinfo(source_tle::TLE, target_tle::TLE, start_time,
                  site_loc_lla)

Initialize the struct OrbitInfo for multiple satellites. Provide the already
loaded TLE files as `SatelliteToolbox` `TLE` structs.

Required Arguments:

- `source_tle::TLE`: TLE struct for satellite that is trasmitting signal
- `target_tle::TLE`: TLE struct for satellite that reflectes transmitted
                        signal from source satellite
- `start_time`: start time in UTC in the format
                (year, month, day, hour, minute, second)
- `site_loc_lla`: site latitude, longitude, and height in (deg, deg, meters)

Returns:

- `OrbitInfo` struct
"""
function initorbitinfo(source_tle::TLE, target_tle::TLE, start_time,
                       site_loc_lla)
    start_time_julian_date = DatetoJD(start_time...)
    lat = deg2rad(site_loc_lla[1])
    lon = deg2rad(site_loc_lla[2])
    height = site_loc_lla[3]
    site_loc_ecef = GeodetictoECEF(lat, lon, height)
    tle = [source_tle, target_tle]
    orb = [init_orbit_propagator(Val(:sgp4), tle[1], sgp4_gc_wgs84),
           init_orbit_propagator(Val(:sgp4), tle[2], sgp4_gc_wgs84)]
    eop = get_iers_eop(:IAU1980)
    tle_names = [source_tle.name, target_tle.name]
    return OrbitInfo(start_time, start_time_julian_date,
                     site_loc_lla, site_loc_ecef,
                     tle_names, tle, orb, eop)
end

"""
    initorbitinfo(source_tle::String, start_time, site_loc_lla)

Initialize the struct OrbitInfo for single satellite. Provide the file name of
the individual TLE file.

Required Arguments:

- `source_tle::String`: TLE string for satellite that is trasmitting signal
- `start_time`: start time in UTC in the format
                (year, month, day, hour, minute, second)
- `site_loc_lla`: site latitude, longitude, and height in (deg, deg, meters)

Returns:

- `OrbitInfo` struct
"""
function initorbitinfo(source_tle::String, start_time, site_loc_lla)
    start_time_julian_date = DatetoJD(start_time...)
    lat = deg2rad(site_loc_lla[1])
    lon = deg2rad(site_loc_lla[2])
    height = site_loc_lla[3]
    site_loc_ecef = GeodetictoECEF(lat, lon, height)
    tle = read_tle(source_tle)[1]
    orb = init_orbit_propagator(Val(:sgp4), tle, sgp4_gc_wgs84)
    eop = get_iers_eop(:IAU1980)
    return OrbitInfo(start_time, start_time_julian_date,
                     site_loc_lla, site_loc_ecef,
                     source_tle, tle, orb, eop)
end

"""
    initorbitinfo(source_tle::TLE, start_time, site_loc_lla)

Initialize the struct OrbitInfo for single satellite. Provide the already loaded
TLE file as a `SatelliteToolbox` `TLE` struct.

Required Arguments:

- `source_tle::TLE`: TLE struct for satellite that is trasmitting signal
- `start_time`: start time in UTC in the format
                (year, month, day, hour, minute, second)
- `site_loc_lla`: site latitude, longitude, and height in (deg, deg, meters)

Returns:

- `OrbitInfo` struct
"""
function initorbitinfo(source_tle::TLE, start_time, site_loc_lla)
    start_time_julian_date = DatetoJD(start_time...)
    lat = deg2rad(site_loc_lla[1])
    lon = deg2rad(site_loc_lla[2])
    height = site_loc_lla[3]
    site_loc_ecef = GeodetictoECEF(lat, lon, height)
    tle = source_tle
    orb = init_orbit_propagator(Val(:sgp4), tle, sgp4_gc_wgs84)
    eop = get_iers_eop(:IAU1980)
    return OrbitInfo(start_time, start_time_julian_date,
                     site_loc_lla, site_loc_ecef,
                     source_tle.name, tle, orb, eop)
end

"""
    SatelliteRAE

Holds information on a given satellite range, azimuth, and elevation from a
observer position.

Fields:

- `name`: name of satellite
- `sat_tle`: `Satellite` or `TLE` struct
- `julian_date_range`: two element vector containing start and end times
                       in Julian days
- `obs_lla`: observer latitude, longitude, and height
- `obs_ecef`: observer ECEF coordinates
- `Δt`: step size of data
- `ts`: vector of time corresponding to data
- `sat_range`: distance between satellite and observer
- `sat_azimuth`: vector of satellite azimuths
- `sat_elevation`: vector of satellite elevations
- `sat_ecef`: vector of satellite ECEF coordinates
"""
struct SatelliteRAE{A1,A2,A3,A4,A5,A6,A7,A8,A9,A10,A11}
    name::A1
    sat_tle::A2
    julian_date_range::A3
    obs_lla::A4
    obs_ecef::A5
    Δt::A6
    ts::A7
    sat_range::A8
    sat_azimuth::A9
    sat_elevation::A10
    sat_ecef::A11
end

"""
    calcenumatrix(obs_lla)

Calculate the ECEF to ENU transformation matrix using the observer's position in
LLA.

**NOTE:** Latitudes and longitudes are in radians.

Required Arguments:

- `obs_lla`: observer latitude, longitude, and height in (deg, deg, meters)

Returns:

- ECEF to ENU transformation matrix
"""
function calcenumatrix(obs_lla)
    lat = deg2rad(obs_lla[1])  # rad
    lon = deg2rad(obs_lla[2])  # rad
    h = obs_lla[3]  # meters
    return [         -sin(lon)           cos(lon)        0;
            -sin(lat)*cos(lon) -sin(lat)*sin(lon) cos(lat);
             cos(lat)*cos(lon)  cos(lat)*sin(lon) sin(lat)]
end

"""
    calcelevation(sat_tle::TLE, julian_date_range, eop, obs_lla;
                  name="Satellite", Δt=1/60/60/24)

Calculates the elevation of a given satellite relative to the observer for every
second between the range specified in `julian_date_range`.

Required Arguments:

- `sat_tle::TLE`: satellite TLE struct
- `julian_date_range`: two element vector containing start and end times
                       in Julian days
- `eop`: Earth Orientation Parameters object
- `obs_lla`: observer latitude, longitude, and height in (deg, deg, meters)

Optional Arguments:

- `name`: name of satellite `(default = "Satellite"`
- `Δt`: step size in seconds `(default = 1/60/60/24)`

Returns:

- `SatelliteRAE` struct
"""
function calcelevation(sat_tle::TLE, julian_date_range, eop, obs_lla;
                       name="Satellite", Δt=1/60/60/24)
    obs_ecef = GeodetictoECEF(deg2rad(obs_lla[1]), deg2rad(obs_lla[2]),
                              obs_lla[3])
    sat_orb = init_orbit_propagator(Val{:sgp4}, sat_tle)
    ts = Array(julian_date_range[1]:Δt:julian_date_range[2])
    # Propagate orbit to ts
    sat_orb, rs, vs = propagate_to_epoch!(sat_orb, ts)
    # Allocate space for storage
    sat_ranges = Array{Float64}(undef, length(ts))
    azs = Array{Float64}(undef, length(ts))
    els = Array{Float64}(undef, length(ts))
    sat_ecefs = Array{Float64}(undef, length(ts), 3)
    # Calculate ENU transformation matrix
    R_ENU = calcenumatrix(obs_lla)
    for i in 1:length(ts)
        # Convert orbits to state vectors
        sat_teme = kepler_to_sv(sat_orb[i])
        # Transform TEME to ECEF frame and extract position
        sat_ecef = svECItoECEF(sat_teme, TEME(), ITRF(), sat_teme.t, eop).r
        # Calculate user-to-sat vector
        user_to_sat = sat_ecef - obs_ecef
        # Caluclate satellite range
        sat_range = norm(user_to_sat)
        # Normalize `user_to_sat`
        user_to_sat_norm = user_to_sat./sat_range
        # Transform normalized user-to-sat vector to ENU
        enu = R_ENU*user_to_sat_norm  # [East, North, South]
        # Calculate satellite azimuth
        az = atan(enu[1], enu[2])
        # Calculate satellite elevation
        el = asin(enu[3]/norm(enu))
        # Save values
        sat_ranges[i] = sat_range
        azs[i] = rad2deg(az)
        els[i] = rad2deg(el)
        sat_ecefs[i,:] = sat_ecef
    end
    return SatelliteRAE(name, sat_tle, julian_date_range,
                        obs_lla, obs_ecef, Δt, ts, sat_ranges,
                        azs, els, sat_ecefs)
end

"""
    calcelevation(satellite::Satellite, obs_lla; name=string(satellite.id))

Calculates the elevation of a given satellite relative to the observer for all
time specified in `satellite.t`.

Required Arguments:

- `satellite::Satellite`: struct containing satellite orbit information
- `obs_lla`: observer latitude, longitude, and height in (deg, deg, meters)

Optional Arguments:

- `name`: name of satellite `(default = string(satellite.id))`

Returns:

- `SatelliteRAE` struct
"""
function calcelevation(satellite::Satellite, obs_lla; name=string(satellite.id))
    obs_ecef = GeodetictoECEF(deg2rad(obs_lla[1]), deg2rad(obs_lla[2]),
                              obs_lla[3])
    ts = satellite.t
    Δt = ts[2] - ts[1]
    julian_date_range = (ts[1], ts[end])
    # Allocate space for storage
    sat_ranges = Array{Float64}(undef, length(ts))
    azs = Array{Float64}(undef, length(ts))
    els = Array{Float64}(undef, length(ts))
    # Calculate ENU transformation matrix
    R_ENU = calcenumatrix(obs_lla)
    for i in 1:length(ts)
        sat_ecef = satellite.r_ecef[i,:]
        # Calculate user-to-sat vector
        user_to_sat = sat_ecef - obs_ecef
        # Caluclate satellite range
        sat_range = norm(user_to_sat)
        # Normalize `user_to_sat`
        user_to_sat_norm = user_to_sat./sat_range
        # Transform normalized user-to-sat vector to ENU
        enu = R_ENU*user_to_sat_norm  # [East, North, South]
        # Calculate satellite azimuth
        az = atan(enu[1], enu[2])
        # Calculate satellite elevation
        el = asin(enu[3]/norm(enu))
        # Save values
        sat_ranges[i] = sat_range
        azs[i] = rad2deg(az)
        els[i] = rad2deg(el)
    end
    return SatelliteRAE(name, satellite, julian_date_range,
                        obs_lla, obs_ecef, Δt, ts, sat_ranges,
                        azs, els, satellite.r_ecef)
end

"""
    calcelevation(sat_ecef, obs_lla)

Calculates the elevation of an object at a given ECEF coordinate specified by
`sat_ecef`. Returns only `(sat_range, az, el)` instead of `SatelliteRAE` struct.
Format is `(meters, deg, deg)`.

Required Arguments:

- `sat_ecef`: satellite ECEF coordinates with all components in meters
- `obs_lla`: observer latitude, longitude, and height in (deg, deg, meters)

Returns:

- `sat_range`: range between satellite and observer in meters
- `az`: satellite azimuth in degrees
- `el`: satellite elevation in degrees
"""
function calcelevation(sat_ecef, obs_lla)
    # Calculate ENU transformation matrix
    R_ENU = calcenumatrix(obs_lla)
    obs_ecef = GeodetictoECEF(deg2rad(obs_lla[1]), deg2rad(obs_lla[2]),
                              obs_lla[3])
    # Calculate user-to-sat vector
    user_to_sat = sat_ecef - obs_ecef
    # Caluclate satellite range
    sat_range = norm(user_to_sat)
    # Normalize `user_to_sat`
    user_to_sat_norm = user_to_sat./sat_range
    # Transform normalized user-to-sat vector to ENU
    enu = R_ENU*user_to_sat_norm  # [East, North, South]
    # Calculate satellite azimuth
    az = atan(enu[1], enu[2])
    # Calculate satellite elevation
    el = asin(enu[3]/norm(enu))
    # Convert from degrees to radians
    az = rad2deg(az)
    el = rad2deg(el)
    return (sat_range, az, el)
end

"""
    getGPSSatnums(obs_time_JD, prns)

Get the NORAD ID for each PRN for a given observationn time in Julian Days. Read
from the `GPSData` dictionary. Returns dictionary with keys being the prns.

Required Arguments:

- `obs_time_JD`: observation time in Julian Days
- `prns`: PRNs to get NORAD ID for

Returns:

- `gps_satnums::Dict`: dictionary containing NORAD IDs for each PRN
    * dictionary keys are the PRN numbers
"""
function getGPSSatnums(obs_time_JD, prns)
    gps_satnums = Dict{Int,Int}()
    for prn in prns
        for sv in collect(keys(GPSData[prn]))
            if GPSData[prn][sv]["active"]
                if obs_time_JD >= GPSData[prn][sv]["start_julian_date"]
                    gps_satnums[prn] = GPSData[prn][sv]["satnum"]
                end
            else
                if ((obs_time_JD > GPSData[prn][sv]["start_julian_date"]) &&
                   (obs_time_JD < GPSData[prn][sv]["end_julian_date"]))
                    gps_satnums[prn] = GPSData[prn][sv]["satnum"]
                end
            end
        end
    end
    return gps_satnums
end

"""
    getTLEs(obs_time_JD, satnums; Δdays=5)

Query Space-Track.org for TLEs matching time and satnum criteria. Parse and
determine TLEs for eac satnum that is closest but before the observation time.

Required Arguments:

- `obs_time_JD`: observation time in Julian Days
- `satnums`: vector of satelleite NORAD IDs

Optional Arguments:

- `Δdays`: Number of days to look before and up to `obs_time_JD`

Returns:

- `filtered_tles::Dict{Int,TLE}`: dictionary containing the TLEs for each
                                  satellite NORAD ID
    * dictionary keys are the NORAD IDs
"""
function getTLEs(obs_time_JD, satnums; Δdays=5)
    obs_time_JD_begin = obs_time_JD - Δdays
    obs_time_begin = JDtoDate(obs_time_JD_begin)
    obs_time = JDtoDate(obs_time_JD+1)
    satnum_list = string(satnums[1])
    for i in 2:length(satnums)
        satnum_list = string(satnum_list, ",", string(satnums[i]))
    end
    tle_file = string(homedir(), "/.GNSSTools/tles.tle")
    date_range = string(string(obs_time_begin[1], pad=4), "-",
                        string(obs_time_begin[2], pad=2), "-",
                        string(obs_time_begin[3], pad=2), "--",
                        string(obs_time[1], pad=4), "-",
                        string(obs_time[2], pad=2), "-",
                        string(obs_time[3], pad=2))
    login_url = "https://www.space-track.org/ajaxauth/login"
    base_query = "query=https://www.space-track.org/basicspacedata/query/class/tle/EPOCH/"
    query_tail = "/orderby/NORAD_CAT_ID/format/3le"
    norad_cat_id = "/NORAD_CAT_ID/"
    # Get username from user
    print("Provide Space-Track.org username: ")
    username = readline(stdin)
    secret_pass = Base.getpass("Provide Space-Track.org user password")
    password = read(secret_pass, String)
    run(`curl -o $tle_file $login_url -d identity=$username"""&"""password=$password"""&"""$base_query$date_range$norad_cat_id$satnum_list$query_tail`;
        wait=false);
    Base.shred!(secret_pass)
    tles = read_tle(tle_file)
    filtered_tles = Dict{Int,TLE}()
    for satnum in satnums
        Δts = []
        idxs = []
        for i in 1:length(tles)
            tle = tles[i]
            if tle.sat_num == satnum
                append!(Δts, abs(tle.epoch-obs_time_JD))
                append!(idxs, i)
            end
        end
        if ~isempty(Δts)
            min_Δt_idx = argmin(Δts)
            filtered_tles[satnum] = tles[idxs[min_Δt_idx]]
        else
            @warn "$satnum had no TLE results. Increase Δdays (currently set to $Δdays)."
        end
    end
    return filtered_tles
end


text:
#=
Contains code derived from the Python module "Newick" Copyright 2003-2008 by Thomas Mailund, released under the GPL v.2
=#

struct Token
    name::Symbol
    value::Any

    # Constructor that preprocesses `:string` tokens
    # TODO: Do the preprocessing elsewhere, this is too specific!
    function Token(name::Symbol, str::String)
        name == :string ? value = strip(str, ['"', '\'', ' ', '\n']) : value = str
        value = String(value)
        new(name, value)
    end
end

mutable struct Lexer
    input::String
    next_token::Union{Nothing, Token}
    tokendefs::Vector{Tuple{Symbol,Regex}}
end

"""
Return the current token in the lexer without taking it out from the lexer.
"""
function peektoken(lxr::Lexer)::Union{Token,Nothing}
    isnothing(lxr.next_token) || return lxr.next_token
    length(lxr.input) == 0 && return nothing
    m = nothing

    for (defname, defregex) in lxr.tokendefs
        m = match(defregex, lxr.input)
        if ! isnothing(m)
            lxr.next_token = Token(defname, convert(String, m.match))
            return lxr.next_token
        end
    end

    error("Invalid token received")
end

"""
Take out the current token in the lexer and return it
"""
function next!(lxr::Lexer)::Token
    token = lxr.next_token

    if isnothing(token)
        token = peektoken(lxr)
    end
    if ! isnothing(token)
        # lxr.input = lxr.input[length(token.value) + 1 : end]
        lxr.input = lxr.input[nextind(lxr.input, lastindex(token.value)):end]
        lxr.next_token = nothing
    end

    return token
end

function readtoken!(lxr::Lexer, name::Symbol)::Token
    token = peektoken(lxr)
    if token.name ≠ name
        error("Expected a `$(name)` token, but received a `$(token.name)` token. String: $(lxr.input)")
    end
    next!(lxr)

    return token
end

function skiptoken!(lxr::Lexer, name::Symbol)::Lexer
    while peektoken(lxr).name == name
        readtoken!(lxr, name)
    end

    return lxr
end


text:
using FTCTests
using Test

@testset "FTCTests.jl" begin
    include("privileged_name.jl")
    include("run_sim.jl")
    include("evaluate.jl")
end


text:
"""
LowRankTransform
```
    P = rand(10,5)
    tr = LowRankTransform(P)
```
Apply the low-rank projection realised by the matrix `P`
The second dimension of `P` must match the number of features of the target.
"""
struct LowRankTransform{T<:AbstractMatrix{<:Real}} <: Transform
    proj::T
end

function set!(t::LowRankTransform{<:AbstractMatrix{T}},M::AbstractMatrix{T}) where {T<:Real}
    @assert size(t) == size(M) "Size of the given matrix $(size(M)) and the projection matrix $(size(t)) are not the same"
    t.proj .= M
end

params(t::LowRankTransform) = t.proj

Base.size(tr::LowRankTransform,i::Int) = size(tr.proj,i)
Base.size(tr::LowRankTransform) = size(tr.proj) #  TODO Add test

function transform(t::LowRankTransform,X::AbstractMatrix{<:Real},obsdim::Int=defaultobs)
    @boundscheck size(t,2) != size(X,feature_dim(obsdim)) ?
        throw(DimensionMismatch("The projection matrix has size $(size(t)) and cannot be used on X with dimensions $(size(X))")) : nothing
    @inbounds _transform(t,X,obsdim)
end

function transform(t::LowRankTransform,x::AbstractVector{<:Real},obsdim::Int=defaultobs) #TODO Add test
    @assert size(t,2) == length(x) "Vector has wrong dimensions $(length(x)) compared to projection matrix"
    t.proj*x
end

_transform(t::LowRankTransform,X::AbstractVecOrMat{<:Real},obsdim::Int=defaultobs) = obsdim == 2 ? t.proj * X : X * t.proj'


text:
@testset "paretoset" begin
    p = hcat(([x,y] for x in 0:0.2:1 for y in 0:0.2:1 if (x-1)^2+(y-1)^2<=1)...)
    @test p[:,ParetoSet(p)] == [0.0 0.2 0.4 1.0; 1.0 0.4 0.2 0.0]
    @test p[:,ParetoSet(p,[:Min, :Min])] == [0.0 0.2 0.4 1.0; 1.0 0.4 0.2 0.0]
    @test p[:,ParetoSet(p,[EfficientGlobalOptimization.Min, EfficientGlobalOptimization.Min])] == [0.0 0.2 0.4 1.0; 1.0 0.4 0.2 0.0]
    p = hcat(([x,y] for x in 0:0.2:1 for y in 0:0.2:1 if x^2+y^2<=1)...)
    @test p[:,ParetoSet(p, [:Max, :Max])] == [0.0 0.6 0.8 1.0; 1.0 0.8 0.6 0.0]
    @test p[:,ParetoSet(p, [EfficientGlobalOptimization.Max, EfficientGlobalOptimization.Max])] == [0.0 0.6 0.8 1.0; 1.0 0.8 0.6 0.0]
end