text:
abstract type AbstractNormalizationLayer{affine, track_stats} <: AbstractExplicitLayer end

"""
    BatchNorm(chs::Integer, activation=identity; init_bias=zeros32, init_scale=ones32, affine=true, track_stats=true, epsilon=1f-5, momentum=0.1f0)

[Batch Normalization](https://arxiv.org/abs/1502.03167) layer.

`BatchNorm` computes the mean and variance for each `D_1×...×D_{N-2}×1×D_N` input slice and normalises the input accordingly.

## Arguments

* `chs`: Size of the channel dimension in your data. Given an array with `N` dimensions, call the `N-1`th the channel dimension. For a batch of feature vectors this is just the data dimension, for `WHCN` images it's the usual channel dimension.
* `activation`: After normalisation, elementwise activation `activation` is applied.

## Keyword Arguments

* If `affine=true`, it also applies  a shift and a rescale to the input through to learnable per-channel bias bias and scale scale parameters.
    - `init_bias`: Controls how the `bias` is initiliazed
    - `init_scale`: Controls how the `scale` is initiliazed
* If `track_stats=true`, accumulates mean and variance statistics in training phase that will be used to renormalize the input in test phase.
* `epsilon`: a value added to the denominator for numerical stability
* `momentum`:  the value used for the `running_mean` and `running_var` computation

## Inputs

* `x`: Array where `size(x, N - 1) = chs` and `ndims(x) > 2`

## Returns

* `y`: Normalized Array
* Update model state

## Parameters

* `affine=true`
    * `bias`: Bias of shape `(chs,)`
    * `scale`: Scale of shape `(chs,)`
* `affine=false` - Empty `NamedTuple()`

## States

* Statistics if `track_stats=true`
    - `running_mean`: Running mean of shape `(chs,)`
    - `running_var`: Running variance of shape `(chs,)`
* Statistics if `track_stats=false`
    - `running_mean`: nothing
    - `running_var`: nothing
* `training`: Used to check if training/inference mode

Use [`Lux.testmode`](@ref) during inference.

## Example

```julia
m = Chain(
    Dense(784 => 64),
    BatchNorm(64, relu),
    Dense(64 => 10),
    BatchNorm(10)
)
```

See also [`GroupNorm`](@ref)
"""
struct BatchNorm{affine, track_stats, F1, F2, F3, N} <:
       AbstractNormalizationLayer{affine, track_stats}
    activation::F1
    epsilon::N
    momentum::N
    chs::Int
    init_bias::F2
    init_scale::F3
end

function BatchNorm(chs::Int,
                   activation=identity;
                   init_bias=zeros32,
                   init_scale=ones32,
                   affine::Bool=true,
                   track_stats::Bool=true,
                   epsilon=1.0f-5,
                   momentum=0.1f0)
    activation = NNlib.fast_act(activation)
    return BatchNorm{affine, track_stats, typeof(activation), typeof(init_bias),
                     typeof(init_scale), typeof(epsilon)}(activation, epsilon, momentum,
                                                          chs, init_bias, init_scale)
end

function initialparameters(rng::AbstractRNG, l::BatchNorm{affine}) where {affine}
    return affine ? (scale=l.init_scale(rng, l.chs), bias=l.init_bias(rng, l.chs)) :
           NamedTuple()
end
function initialstates(rng::AbstractRNG,
                       l::BatchNorm{affine, track_stats}) where {affine, track_stats}
    return if track_stats
        (running_mean=zeros32(rng, l.chs), running_var=ones32(rng, l.chs),
         training=Val(true))
    else
        (running_mean=nothing, running_var=nothing, training=Val(true))
    end
end

parameterlength(l::BatchNorm{affine}) where {affine} = affine ? (l.chs * 2) : 0
function statelength(l::BatchNorm{affine, track_stats}) where {affine, track_stats}
    (track_stats ? 2 * l.chs : 0) + 1
end

function (BN::BatchNorm)(x::AbstractArray{T, N}, ps, st::NamedTuple) where {T, N}
    @assert size(x, N - 1) == BN.chs
    @assert !istraining(st)||size(x, N) > 1 "During `training`, `BatchNorm` can't handle Batch Size == 1"

    x_normalized, xmean, xvar = normalization(x,
                                              st.running_mean,
                                              st.running_var,
                                              ps.scale,
                                              ps.bias,
                                              BN.activation,
                                              collect([1:(N - 2); N]),
                                              st.training,
                                              BN.momentum,
                                              BN.epsilon)

    st = merge(st, (running_mean=xmean, running_var=xvar))

    return x_normalized, st
end

function (BN::BatchNorm{affine, track_stats})(x::Union{CuArray{T, 2}, CuArray{T, 4},
                                                       CuArray{T, 5}}, ps,
                                              st::NamedTuple) where {
                                                                     T <:
                                                                     Union{Float32, Float64
                                                                           }, affine,
                                                                     track_stats}
    # NNlibCUDA silently updates running_mean and running_var so copying them
    if istraining(st)
        running_mean2 = track_stats ? copy(st.running_mean) : nothing
        running_var2 = track_stats ? copy(st.running_var) : nothing
    else
        if track_stats
            running_mean2 = copy(st.running_mean)
            running_var2 = copy(st.running_var)
        else
            N = ndims(x)
            reduce_dims = collect([1:(N - 2); N])
            running_mean2 = mean(x; dims=reduce_dims)
            running_var2 = var(x; mean=running_mean2, dims=reduce_dims,
                               corrected=false)
        end
    end
    res = applyactivation(BN.activation,
                          batchnorm(affine ? ps.scale : nothing,
                                    affine ? ps.bias : nothing,
                                    x,
                                    running_mean2,
                                    running_var2,
                                    BN.momentum;
                                    eps=BN.epsilon,
                                    training=istraining(st)))
    if track_stats
        st = merge(st, (running_mean=running_mean2, running_var=running_var2))
    end
    return res, st
end

function Base.show(io::IO, l::BatchNorm{affine, track_stats}) where {affine, track_stats}
    print(io, "BatchNorm($(l.chs)")
    (l.activation == identity) || print(io, ", $(l.activation)")
    affine || print(io, ", affine=false")
    track_stats || print(io, ", track_stats=false")
    return print(io, ")")
end

"""
    GroupNorm(chs::Integer, groups::Integer, activation=identity; init_bias=zeros32, init_scale=ones32, affine=true, track_stats=false, epsilon=1f-5, momentum=0.1f0)

[Group Normalization](https://arxiv.org/abs/1803.08494) layer.

## Arguments

* `chs`: Size of the channel dimension in your data. Given an array with `N` dimensions, call the `N-1`th the channel dimension. For a batch of feature vectors this is just the data dimension, for `WHCN` images it's the usual channel dimension.
* `groups` is the number of groups along which the statistics are computed. The number of channels must be an integer multiple of the number of groups.
* `activation`: After normalisation, elementwise activation `activation` is applied.

## Keyword Arguments

* If `affine=true`, it also applies  a shift and a rescale to the input through to learnable per-channel bias bias and scale scale parameters.
    - `init_bias`: Controls how the `bias` is initiliazed
    - `init_scale`: Controls how the `scale` is initiliazed
* If `track_stats=true`, accumulates mean and variance statistics in training phase that will be used to renormalize the input in test phase.
* `epsilon`: a value added to the denominator for numerical stability
* `momentum`:  the value used for the `running_mean` and `running_var` computation

## Inputs

* `x`: Array where `size(x, N - 1) = chs` and `ndims(x) > 2`

## Returns

* `y`: Normalized Array
* Update model state

## Parameters

* `affine=true`
    * `bias`: Bias of shape `(chs,)`
    * `scale`: Scale of shape `(chs,)`
* `affine=false` - Empty `NamedTuple()`

## States

* Statistics if `track_stats=true`
    - `running_mean`: Running mean of shape `(groups,)`
    - `running_var`: Running variance of shape `(groups,)`
* Statistics if `track_stats=false`
    - `running_mean`: nothing
    - `running_var`: nothing
* `training`: Used to check if training/inference mode

Use [`Lux.testmode`](@ref) during inference.

## Example

```julia
m = Chain(
    Dense(784 => 64),
    GroupNorm(64, 4, relu),
    Dense(64 => 10),
    GroupNorm(10, 5)
)
```

!!! warning
    GroupNorm doesn't have CUDNN support. The GPU fallback is not very efficient.

See also [`BatchNorm`](@ref)
"""
struct GroupNorm{affine, track_stats, F1, F2, F3, N} <:
       AbstractNormalizationLayer{affine, track_stats}
    activation::F1
    epsilon::N
    momentum::N
    chs::Int
    init_bias::F2
    init_scale::F3
    groups::Int
end

function GroupNorm(chs::Int,
                   groups::Int,
                   activation=identity;
                   init_bias=zeros32,
                   init_scale=ones32,
                   affine::Bool=true,
                   track_stats::Bool=true,
                   epsilon=1.0f-5,
                   momentum=0.1f0)
    @assert chs % groups==0 "The number of groups ($(groups)) must divide the number of channels ($chs)"
    activation = NNlib.fast_act(activation)
    return GroupNorm{affine, track_stats, typeof(activation), typeof(init_bias),
                     typeof(init_scale), typeof(epsilon)}(activation, epsilon, momentum,
                                                          chs, init_bias, init_scale,
                                                          groups)
end

function initialparameters(rng::AbstractRNG, l::GroupNorm{affine}) where {affine}
    return affine ? (scale=l.init_scale(rng, l.chs), bias=l.init_bias(rng, l.chs)) :
           NamedTuple()
end
function initialstates(rng::AbstractRNG,
                       l::GroupNorm{affine, track_stats}) where {affine, track_stats}
    return if track_stats
        (running_mean=zeros32(rng, l.groups), running_var=ones32(rng, l.groups),
         training=Val(true))
    else
        (running_mean=nothing, running_var=nothing, training=Val(true))
    end
end

parameterlength(l::GroupNorm{affine}) where {affine} = affine ? (l.chs * 2) : 0
function statelength(l::GroupNorm{affine, track_stats}) where {affine, track_stats}
    (track_stats ? 2 * l.groups : 0) + 1
end

function (GN::GroupNorm)(x::AbstractArray{T, N}, ps, st::NamedTuple) where {T, N}
    sz = size(x)
    @assert N > 2
    @assert sz[N - 1] == GN.chs

    x_ = reshape(x, sz[1:(N - 2)]..., sz[N - 1] ÷ GN.groups, GN.groups, sz[N])

    x_normalized, xmean, xvar = normalization(x_,
                                              st.running_mean,
                                              st.running_var,
                                              ps.scale,
                                              ps.bias,
                                              GN.activation,
                                              collect(1:(N - 1)),
                                              st.training,
                                              GN.momentum,
                                              GN.epsilon)

    st = merge(st, (running_mean=xmean, running_var=xvar))

    return reshape(x_normalized, sz), st
end

function Base.show(io::IO, l::GroupNorm{affine, track_stats}) where {affine, track_stats}
    print(io, "GroupNorm($(l.chs), $(l.groups)")
    (l.activation == identity) || print(io, ", $(l.activation)")
    affine || print(io, ", affine=false")
    track_stats || print(io, ", track_stats=false")
    return print(io, ")")
end

"""
    WeightNorm(layer::AbstractExplicitLayer, which_params::NTuple{N,Symbol}, dims::Union{Tuple,Nothing}=nothing)

Applies [weight normalization](https://arxiv.org/abs/1602.07868) to a parameter in the given layer.

``w = g\\frac{v}{\\|v\\|}``

Weight normalization is a reparameterization that decouples the magnitude of a weight tensor from its direction. This updates the parameters in `which_params` (e.g. `weight`) using two parameters: one specifying the magnitude (e.g. `weight_g`) and one specifying the direction (e.g. `weight_v`).

## Arguments

* `layer` whose parameters are being reparameterized
* `which_params`: parameter names for the parameters being reparameterized
* By default, a norm over the entire array is computed. Pass `dims` to modify the dimension.

## Inputs

* `x`: Should be of valid type for input to `layer`

## Returns

* Output from `layer`
* Updated model state of `layer`

## Parameters

* `normalized`: Parameters of `layer` that are being normalized
* `unnormalized`: Parameters of `layer` that are not being normalized

## States

* Same as that of `layer`
"""
struct WeightNorm{which_params, L <: AbstractExplicitLayer, D} <: AbstractExplicitLayer
    layer::L
    dims::D
end

function WeightNorm(layer::AbstractExplicitLayer, which_params::NTuple{N, Symbol},
                    dims::Union{Tuple, Nothing}=nothing) where {N}
    return WeightNorm{Val{which_params}, typeof(layer), typeof(dims)}(layer, dims)
end

function initialparameters(rng::AbstractRNG,
                           wn::WeightNorm{Val{which_params}}) where {which_params}
    ps_layer = initialparameters(rng, wn.layer)
    ps_normalized = []
    ps_unnormalized = []
    i = 1
    for k in propertynames(ps_layer)
        v = ps_layer[k]
        if k ∈ which_params
            dim = wn.dims === nothing ? ndims(v) : wn.dims[i]
            push!(ps_normalized, Symbol(string(k) * "_g") => _norm_except(v, dim))
            push!(ps_normalized, Symbol(string(k) * "_v") => v)
            i += 1
        else
            push!(ps_unnormalized, k => v)
        end
    end
    ps_unnormalized = length(ps_unnormalized) == 0 ? NamedTuple() : (; ps_unnormalized...)
    return (normalized=(; ps_normalized...), unnormalized=ps_unnormalized)
end

initialstates(rng::AbstractRNG, wn::WeightNorm) = initialstates(rng, wn.layer)

function (wn::WeightNorm)(x, ps::Union{ComponentArray, NamedTuple}, s::NamedTuple)
    _ps = get_normalized_parameters(wn, wn.dims, ps.normalized)
    return wn.layer(x, merge(_ps, ps.unnormalized), s)
end

@inbounds @generated function get_normalized_parameters(::WeightNorm{Val{which_params}},
                                                        dims::T,
                                                        ps::Union{ComponentArray, NamedTuple
                                                                  }) where {T, which_params}
    parameter_names = string.(which_params)
    v_parameter_names = Symbol.(parameter_names .* "_v")
    g_parameter_names = Symbol.(parameter_names .* "_g")
    normalized_params_symbol = [gensym(p) for p in parameter_names]

    function get_norm_except_invoke(i)
        return if T <: Tuple
            :(_norm_except(ps.$(v_parameter_names[i]), dims[$i]))
        else
            :(_norm_except(ps.$(v_parameter_names[i])))
        end
    end

    calls = []
    for i in 1:length(parameter_names)
        push!(calls,
              :($(normalized_params_symbol[i]) = ps.$(v_parameter_names[i]) .*
                                                 (ps.$(g_parameter_names[i]) ./
                                                  $(get_norm_except_invoke(i)))))
    end
    push!(calls,
          :(return NamedTuple{$(which_params)}(tuple($(Tuple(normalized_params_symbol)...)))))

    return Expr(:block, calls...)
end

function Base.show(io::IO, w::WeightNorm{Val{which_params}}) where {which_params}
    return print(io, "WeightNorm{", which_params, "}(", w.layer, ")")
end


text:
using SmallParser
using Test

@testset "SmallParser.jl" begin
    # Write your tests here.
end

ext: jl


text:
# these are all custom diagnostic functions. May help with debugging
using Statistics
using LinearAlgebra

"""
	est_d(f, input; dif=1e-7)

Estimate the derivative/gradient of `f` at `inputs` with forward differences

# Keyword Arguments
- `dif::Real=1e-7`: The forward differencing step size
"""
est_d(f::Function, input; dif::Real=1e-7) = (f(input + dif) - f(input)) ./ dif

"""
	est_∇nlogL_kep

Estimate the gradient of `nlogL_kep(data, times, covariance, ks)` with forward
differences

# Keyword Arguments
- `data_unit::Unitful.Velocity=1u"m/s"`: The units of the input data vector
- `include_priors::Bool=true`: Whether to add a prior evaluation to the output
- `use_hk::Bool=true`: Whether to use h and k or e and ω parameterization
- `kwargs...`: keyword arguments that are passed to GLOM.est_∇()
"""
function est_∇nlogL_kep(
    data::Vector{T},
    times::Vector{T2} where T2<:Unitful.Time,
    covariance::Union{Cholesky,Diagonal},
    ks::kep_signal;
    data_unit::Unitful.Velocity=1u"m/s",
    include_priors::Bool=true,
    use_hk::Bool=true,
    kwargs...
    ) where T<:Real

    K_u = unit(ks.K)
    P_u = unit(ks.P)
    γ_u = unit(ks.γ)
    og = ustrip.([ks.K, ks.P, ks.M0, ks.h, ks.k, ks.γ])
    function f2(vector)
        ks_internal = ks_from_vec(vector, K_u, P_u, γ_u; use_hk=true)
        return GLOM.nlogL(covariance, remove_kepler(data, times, ks_internal; data_unit=data_unit)) - (include_priors * logprior_kepler(ks_internal; use_hk=use_hk))
    end
    return GLOM.est_∇(f2, og; kwargs...)
end
function est_∇nlogL_kep(
    data::Vector{T},
    times::Vector{T2} where T2<:Unitful.Time,
    covariance::Union{Cholesky,Diagonal},
    ks::kep_signal_wright;
    data_unit::Unitful.Velocity=1u"m/s",
    include_priors::Bool=true,
    use_hk::Bool=false,
    kwargs...
    ) where T<:Real

    P_u = unit(ks.P)
    og = ustrip.([ks.P, ks.M0, ks.e])
    function f2(vector)
        ks_internal = fit_kepler_wright_linear_step(data, times, covariance, vector[1] * P_u, vector[2], vector[3]; data_unit=data_unit)
        return GLOM.nlogL(covariance, remove_kepler(data, times, ks_internal; data_unit=data_unit)) - (include_priors * logprior_kepler(ks_internal; use_hk=use_hk))
    end
    return GLOM.est_∇(f2, og; kwargs...)
end
est_∇nlogL_kep(
    glo_rv::GLO_RV,
    covariance::Union{Cholesky,Diagonal},
    ks::Union{kep_signal, kep_signal_wright};
    include_priors::Bool=true,
    use_hk::Bool=false,
    kwargs...
    ) where T<:Real = est_∇nlogL_kep(glo_rv.GLO.y_obs, glo_rv.time, covariance, ks; data_unit=glo_rv.rv_factor, include_priors=include_priors, use_hk=use_hk, kwargs...)

function test_∇nlogL_kep(
    data::Vector{T},
    times::Vector{T2} where T2<:Unitful.Time,
    covariance::Union{Cholesky,Diagonal},
    ks::Union{kep_signal, kep_signal_wright};
    data_unit::Unitful.Velocity=1u"m/s",
    include_priors::Bool=true
    ) where T<:Real

    grad_est = est_∇nlogL_kep(data, times, covariance, ks; data_unit=data_unit, include_priors=include_priors)
    grad = ∇nlogL_kep(data, times, covariance, ks; data_unit=data_unit, include_priors=include_priors)
    return GLOM.test_∇(grad_est, grad; function_name="∇nlogL_kep")
end
test_∇nlogL_kep(
    glo_rv::GLO_RV,
    covariance::Union{Cholesky,Diagonal},
    ks::Union{kep_signal, kep_signal_wright};
    include_priors::Bool=true
    ) where T<:Real = test_∇nlogL_kep(glo_rv.GLO.y_obs, glo_rv.time, covariance, ks; data_unit=glo_rv.rv_factor, include_priors=include_priors)

function est_∇∇nlogL_kep(
    data::Vector{T},
    times::Vector{T2} where T2<:Unitful.Time,
    covariance::Union{Cholesky,Diagonal},
    ks::kep_signal;
    kwargs...
    ) where T<:Real

    K_u = unit(ks.K)
    P_u = unit(ks.P)
    γ_u = unit(ks.γ)
    og = ustrip.([ks.K, ks.P, ks.M0, ks.h, ks.k, ks.γ])
    g2(vector) = ∇nlogL_kep(data, times, covariance, ks_from_vec(vector, K_u, P_u, γ_u; use_hk=true); kwargs...)
    return GLOM.est_∇∇(g2, og)
end
est_∇∇nlogL_kep(
    glo_rv::GLO_RV,
    covariance::Union{Cholesky,Diagonal},
    ks::kep_signal;
    kwargs...
    ) where T<:Real = est_∇∇nlogL_kep(glo_rv.GLO.y_obs, glo_rv.time, covariance, ks; data_unit=glo_rv.rv_factor, kwargs...)

function test_∇∇nlogL_kep(
    data::Vector{T},
    times::Vector{T2} where T2<:Unitful.Time,
    covariance::Union{Cholesky,Diagonal},
    ks::kep_signal;
    kwargs...
    ) where T<:Real

    hess_est = est_∇∇nlogL_kep(data, times, covariance, ks; kwargs...)
    hess = ∇∇nlogL_kep(data, times, covariance, ks; kwargs...)
    return GLOM.test_∇∇(hess_est, hess; function_name="∇∇nlogL_kep")
end
test_∇∇nlogL_kep(
    glo_rv::GLO_RV,
    covariance::Union{Cholesky,Diagonal},
    ks::kep_signal;
    kwargs...
    ) where T<:Real = test_∇∇nlogL_kep(glo_rv.GLO.y_obs, glo_rv.time, covariance, ks; data_unit=glo_rv.rv_factor, kwargs...)


text:
using MPI
using StaticArrays
using CLIMA
using CLIMA.VariableTemplates
using CLIMA.Mesh.Topologies
using CLIMA.Mesh.Grids
using CLIMA.MPIStateArrays
using CLIMA.DGmethods
using CLIMA.DGmethods.NumericalFluxes
using Printf
using LinearAlgebra
using Logging
using GPUifyLoops

import CLIMA.DGmethods:
    BalanceLaw,
    vars_aux,
    vars_state,
    vars_gradient,
    vars_diffusive,
    vars_integrals,
    integral_load_aux!,
    flux_nondiffusive!,
    flux_diffusive!,
    source!,
    wavespeed,
    update_aux!,
    indefinite_stack_integral!,
    reverse_indefinite_stack_integral!,
    boundary_state!,
    gradvariables!,
    init_aux!,
    init_state!,
    init_ode_state,
    LocalGeometry,
    integral_set_aux!,
    vars_reverse_integrals,
    reverse_integral_load_aux!,
    reverse_integral_set_aux!

struct IntegralTestSphereModel{T} <: BalanceLaw
    Rinner::T
    Router::T
end

function update_aux!(
    dg::DGModel,
    m::IntegralTestSphereModel,
    Q::MPIStateArray,
    t::Real,
)
    indefinite_stack_integral!(dg, m, Q, dg.auxstate, t)
    reverse_indefinite_stack_integral!(dg, m, Q, dg.auxstate, t)

    return true
end

vars_integrals(::IntegralTestSphereModel, T) = @vars(v::T)
vars_reverse_integrals(::IntegralTestSphereModel, T) = @vars(v::T)
vars_aux(m::IntegralTestSphereModel, T) =
    @vars(int::vars_integrals(m, T), rev_int::vars_integrals(m, T), r::T, a::T)

vars_state(::IntegralTestSphereModel, T) = @vars()
vars_diffusive(::IntegralTestSphereModel, T) = @vars()

flux_nondiffusive!(::IntegralTestSphereModel, _...) = nothing
flux_diffusive!(::IntegralTestSphereModel, _...) = nothing
source!(::IntegralTestSphereModel, _...) = nothing
boundary_state!(_, ::IntegralTestSphereModel, _...) = nothing
init_state!(::IntegralTestSphereModel, _...) = nothing
wavespeed(::IntegralTestSphereModel, _...) = 1

function init_aux!(m::IntegralTestSphereModel, aux::Vars, g::LocalGeometry)

    x, y, z = g.coord
    aux.r = hypot(x, y, z)
    θ = atan(y, x)
    ϕ = asin(z / aux.r)
    # Exact integral
    aux.a = 1 + cos(ϕ)^2 * sin(θ)^2 + sin(ϕ)^2
    aux.int.v = exp(-aux.a * aux.r^2) - exp(-aux.a * m.Rinner^2)
    aux.rev_int.v = exp(-aux.a * m.Router^2) - exp(-aux.a * aux.r^2)
end

@inline function integral_load_aux!(
    m::IntegralTestSphereModel,
    integrand::Vars,
    state::Vars,
    aux::Vars,
)
    integrand.v = -2 * aux.r * aux.a * exp(-aux.a * aux.r^2)
end

@inline function integral_set_aux!(
    m::IntegralTestSphereModel,
    aux::Vars,
    integral::Vars,
)
    aux.int.v = integral.v
end

@inline function reverse_integral_load_aux!(
    m::IntegralTestSphereModel,
    integral::Vars,
    state::Vars,
    aux::Vars,
)
    integral.v = aux.int.v
end

@inline function reverse_integral_set_aux!(
    m::IntegralTestSphereModel,
    aux::Vars,
    integral::Vars,
)
    aux.rev_int.v = integral.v
end

if !@isdefined integration_testing
    const integration_testing = parse(
        Bool,
        lowercase(get(ENV, "JULIA_CLIMA_INTEGRATION_TESTING", "false")),
    )
end

using Test
function run(mpicomm, topl, ArrayType, N, FT, Rinner, Router)
    grid = DiscontinuousSpectralElementGrid(
        topl,
        FloatType = FT,
        DeviceArray = ArrayType,
        polynomialorder = N,
        meshwarp = Topologies.cubedshellwarp,
    )
    dg = DGModel(
        IntegralTestSphereModel(Rinner, Router),
        grid,
        Rusanov(),
        CentralNumericalFluxDiffusive(),
        CentralNumericalFluxGradient(),
    )

    Q = init_ode_state(dg, FT(0))
    dQdt = similar(Q)

    exact_aux = copy(dg.auxstate)
    dg(dQdt, Q, nothing, 0.0)
    euclidean_distance(exact_aux, dg.auxstate)
end

let
    CLIMA.init()
    ArrayType = CLIMA.array_type()

    mpicomm = MPI.COMM_WORLD

    polynomialorder = 4

    base_Nhorz = 4
    base_Nvert = 2
    Rinner = 1 // 2
    Router = 1

    polynomialorder = 4

    expected_result = [
        4.662884229467401e-7,
        7.218989778540723e-9,
        1.1258613174916711e-10,
        1.7587739986848968e-12,
    ]
    lvls = integration_testing ? length(expected_result) : 1

    for FT in (Float64,) # Float32)
        err = zeros(FT, lvls)
        for l in 1:lvls
            @info (ArrayType, FT, "sphere", l)
            Nhorz = 2^(l - 1) * base_Nhorz
            Nvert = 2^(l - 1) * base_Nvert
            Rrange = grid1d(FT(Rinner), FT(Router); nelem = Nvert)
            topl = StackedCubedSphereTopology(mpicomm, Nhorz, Rrange)
            err[l] = run(
                mpicomm,
                topl,
                ArrayType,
                polynomialorder,
                FT,
                FT(Rinner),
                FT(Router),
            )
            @test expected_result[l] ≈ err[l] rtol = 1e-3 atol = eps(FT)
        end
        if integration_testing
            @info begin
                msg = ""
                for l in 1:(lvls - 1)
                    rate = log2(err[l]) - log2(err[l + 1])
                    msg *= @sprintf("\n  rate for level %d = %e\n", l, rate)
                end
                msg
            end
        end
    end
end

nothing

max_stars_repo_name: trontrytel/CLIMA

max_issues_repo_name: trontrytel/CLIMA

max_forks_repo_name: trontrytel/CLIMA


text:
using DelimitedFiles

"""
    DataSplit

A basic struct for encapsulating the four components of supervised training.
"""
struct DataSplit
    train_x::AdaptiveResonance.RealMatrix
    test_x::AdaptiveResonance.RealMatrix
    train_y::AdaptiveResonance.IntegerVector
    test_y::AdaptiveResonance.IntegerVector
    DataSplit(train_x, test_x, train_y, test_y) = new(train_x, test_x, train_y, test_y)
end # DataSplit

"""
    DataSplit(data_x::Array, data_y::Array, ratio::Float)

Return a DataSplit struct that is split by the ratio (e.g. 0.8).
"""
function DataSplit(data_x::Array, data_y::Array, ratio::Real)
    dim, n_data = size(data_x)
    split_ind = Integer(floor(n_data*ratio))

    train_x = data_x[:, 1:split_ind]
    test_x = data_x[:, split_ind+1:end]
    train_y = data_y[1:split_ind]
    test_y = data_y[split_ind+1:end]

    return DataSplit(train_x, test_x, train_y, test_y)
end # DataSplit(data_x::Array, data_y::Array, ratio::Real)

"""
    train_test_art(art::ARTModule, data::DataSplit; supervised::Bool=false, art_opts...)

Train and test an ART module.
"""
function train_test_art(art::ARTModule, data::DataSplit; supervised::Bool=false, train_opts::NamedTuple=NamedTuple(), test_opts::NamedTuple=NamedTuple())
    # Default performance to undefined
    perf = NaN
    # If the module is unsupervised by default
    if art isa ART
        # Check if training with lazy supervision
        if supervised
            # Train with the supervised ART syntax
            train!(art, data.train_x, y=data.train_y; train_opts...)
            y_hat = classify(art, data.test_x; test_opts...)

            # Calculate performance
            perf = performance(y_hat, data.test_y)
        # Otherwise, train in an unsupervised fashion
        else
            train!(art, data.train_x; train_opts...)
        end
    # Otherwise, necessarily train on a supervised model
    elseif art isa ARTMAP
        # Train and classify
        train!(art, data.train_x, data.train_y; train_opts...)
        y_hat = classify(art, data.test_x; test_opts...)

        # Calculate performance
        perf = performance(y_hat, data.test_y)
    else
        error("Incompatible ART module passed for testing")
    end

    # If the performance is not a NaN (potentially unsupervsied), then display perf
    if !isnan(perf)
        @info "$(typeof(art)): performance is $perf"
    end

    return perf
end

"""
    showtypetree(T, level=0)

Show the tree of subtypes for a type.
```julia
showtypetree(Number)
```
"""
function showtypetree(T, level=0)
    println("\t" ^ level, T)
    for t in subtypes(T)
        showtypetree(t, level+1)
    end
end # showtypetree(T, level=0)

"""
    load_iris(data_path::String ; split_ratio::Real = 0.8)

Loads the iris dataset for testing and examples.
"""
function load_iris(data_path::String ; split_ratio::Real = 0.8)
    raw_data = readdlm(data_path,',')
    labels = ["Iris-setosa", "Iris-versicolor", "Iris-virginica"]
    raw_x = Array{AdaptiveResonance.RealFP}(raw_data[2:end, 2:5])
    raw_y_labels = raw_data[2:end, 6]
    raw_y = Array{Integer}(undef, 0)
    for ix = 1:length(raw_y_labels)
        for jx = 1:length(labels)
            if raw_y_labels[ix] == labels[jx]
                push!(raw_y, jx)
            end
        end
    end
    n_samples, n_features = size(raw_x)

    # Julia is column-major, so use columns for features
    raw_x = permutedims(raw_x)

    # Create the datasplit object
    data = DataSplit(raw_x, raw_y, split_ratio)

    return data
end # load_iris(data_path::String ; split_ratio::Real = 0.8)


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

# declaring the Geometry module. in the future we might add or move existing code to it. currently it contains the basic types: Vec3, Vec4 and Transform
module Geometry

include("Transform.jl")

end # module Geometry
export Geometry

using .Geometry

include("Ray.jl")

include("Surface.jl")

# not ideal having this like this, but we need the types for Intersection
include("../Optical/OpticalInterface.jl")

include("CSG/Intersection.jl")
include("CSG/Interval.jl")

include("Primitives/NonCSG/Triangle.jl")

include("BoundingBox.jl")

include("Primitives/Plane.jl")
include("Primitives/Cylinder.jl")
include("Primitives/Sphere.jl")
include("Primitives/SphericalCap.jl")

include("Primitives/NonCSG/Rectangle.jl")
include("Primitives/NonCSG/Hexagon.jl")
include("Primitives/NonCSG/Ellipse.jl")
include("Primitives/NonCSG/Stop.jl")

# include("BoundingVolumeHierarchy.jl")

include("AccelSurface.jl")

include("Primitives/Curves/Knots.jl")
include("Primitives/Curves/Spline.jl")
include("Primitives/Curves/BSpline.jl")
include("Primitives/Curves/Bezier.jl")
include("Primitives/Curves/PowerBasis.jl")

include("CSG/CSG.jl")

include("Primitives/Zernike.jl")
include("Primitives/Qtype.jl")
include("Primitives/Chebyshev.jl")
include("Primitives/GridSag.jl")

include("AnalyticIntersection.jl")

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

"""
    BoundedCylinder(radius::T, height::T; interface::NullOrFresnel{T} = nullinterface(T)) -> CSGGenerator{T}

Create a cylinder with planar caps on both ends centred at `(0, 0, 0)` with axis `(0, 0, 1)`.
"""
function BoundedCylinder(radius::T, height::T; interface::NullOrFresnel{T} = NullInterface(T)) where {T<:Real}
    barrel = Cylinder(radius, height, interface = interface)
    top = Plane(SVector{3,T}(0.0, 0.0, 1.0), SVector{3,T}(0.0, 0.0, height / 2), vishalfsizeu = radius, vishalfsizev = radius)
    bottom = Plane(SVector{3,T}(0.0, 0.0, -1.0), SVector{3,T}(0.0, 0.0, -height / 2), vishalfsizeu = radius, vishalfsizev = radius)
    return csgintersection(leaf(barrel), csgintersection(top, bottom))
end
export BoundedCylinder

"""
    Cuboid(halfsizex::T, halfsizey::T, halfsizez::T; interface::NullOrFresnel{T} = nullinterface(T)) -> CSGGenerator{T}

Create a cuboid centred at `(0, 0, 0)`.
"""
function Cuboid(halfsizex::T, halfsizey::T, halfsizez::T; interface::NullOrFresnel{T} = NullInterface(T)) where {T<:Real}
    xmin = Plane(SVector{3,T}(-1.0, 0.0, 0.0), SVector{3,T}(-halfsizex, 0.0, 0.0), vishalfsizeu = halfsizez, vishalfsizev = halfsizey, interface = interface)
    xmax = Plane(SVector{3,T}(1.0, 0.0, 0.0), SVector{3,T}(halfsizex, 0.0, 0.0), vishalfsizeu = halfsizez, vishalfsizev = halfsizey, interface = interface)
    ymin = Plane(SVector{3,T}(0.0, -1.0, 0.0), SVector{3,T}(0.0, -halfsizey, 0.0), vishalfsizeu = halfsizez, vishalfsizev = halfsizex, interface = interface)
    ymax = Plane(SVector{3,T}(0.0, 1.0, 0.0), SVector{3,T}(0.0, halfsizey, 0.0), vishalfsizeu = halfsizez, vishalfsizev = halfsizex, interface = interface)
    zmin = Plane(SVector{3,T}(0.0, 0.0, -1.0), SVector{3,T}(0.0, 0.0, -halfsizez), vishalfsizeu = halfsizex, vishalfsizev = halfsizey, interface = interface)
    zmax = Plane(SVector{3,T}(0.0, 0.0, 1.0), SVector{3,T}(0.0, 0.0, halfsizez), vishalfsizeu = halfsizex, vishalfsizev = halfsizey, interface = interface)
    xint = csgintersection(xmin, xmax)
    yint = csgintersection(ymin, ymax)
    zint = csgintersection(zmin, zmax)
    return csgintersection(xint, csgintersection(yint, zint))
end
export Cuboid

"""
    HexagonalPrism(side_length::T, visheight::T = 2.0; interface::NullOrFresnel{T} = nullinterface(T)) -> CSGGenerator{T}

Create an infinitely tall hexagonal prism with axis `(0, 0, 1)`, the longer hexagon diameter is along the x axis.
For visualization `visheight` is used, **note that this does not fully represent the surface**.
"""
function HexagonalPrism(side_length::T, visheight::T = 2.0; interface::NullOrFresnel{T} = NullInterface(T)) where {T<:Real}
    h = side_length * sqrt(3.0) / 2.0
    hcos = h * sqrt(3.0) / 2.0
    hsin = h / 2.0
    side_length *= 1.0001
    s1 = Plane(SVector{3,T}(0.0, 1.0, 0.0), SVector{3,T}(0.0, h, 0.0), vishalfsizeu = visheight / 2.0, vishalfsizev = side_length / 2.0, interface = interface)
    s2 = Plane(SVector{3,T}(0.0, -1.0, 0.0), SVector{3,T}(0.0, -h, 0.0), vishalfsizeu = visheight / 2.0, vishalfsizev = side_length / 2.0, interface = interface)
    s3 = Plane(SVector{3,T}(hcos, hsin, 0.0), SVector{3,T}(hcos, hsin, 0.0), vishalfsizeu = visheight / 2.0, vishalfsizev = side_length / 2.0, interface = interface)
    s4 = Plane(SVector{3,T}(-hcos, -hsin, 0.0), SVector{3,T}(-hcos, -hsin, 0.0), vishalfsizeu = visheight / 2.0, vishalfsizev = side_length / 2.0, interface = interface)
    s5 = Plane(SVector{3,T}(hcos, -hsin, 0.0), SVector{3,T}(hcos, -hsin, 0.0), vishalfsizeu = visheight / 2.0, vishalfsizev = side_length / 2.0, interface = interface)
    s6 = Plane(SVector{3,T}(-hcos, hsin, 0.0), SVector{3,T}(-hcos, hsin, 0.0), vishalfsizeu = visheight / 2.0, vishalfsizev = side_length / 2.0, interface = interface)
    i1 = csgintersection(s1, s2)
    i2 = csgintersection(s3, s4)
    i3 = csgintersection(s5, s6)
    return csgintersection(i1, csgintersection(i2, i3))
end
export HexagonalPrism

"""
    RectangularPrism(halfsizex::T, halfsizey::T, visheight::T=2.0; interface::NullOrFresnel{T} = nullinterface(T)) -> CSGGenerator{T}

Create an infinitely tall rectangular prism with axis `(0, 0, 1)`.
For visualization `visheight` is used, **note that this does not fully represent the surface**.
"""
function RectangularPrism(halfsizex::T, halfsizey::T, visheight::T = 2.0; interface::NullOrFresnel{T} = NullInterface(T)) where {T<:Real}
    s1 = Plane(SVector{3,T}(0.0, 1.0, 0.0), SVector{3,T}(0.0, halfsizey, 0.0), vishalfsizeu = visheight / 2.0, vishalfsizev = halfsizex, interface = interface)
    s2 = Plane(SVector{3,T}(0.0, -1.0, 0.0), SVector{3,T}(0.0, -halfsizey, 0.0), vishalfsizeu = visheight / 2.0, vishalfsizev = halfsizex, interface = interface)
    s3 = Plane(SVector{3,T}(1.0, 0.0, 0.0), SVector{3,T}(halfsizex, 0.0, 0.0), vishalfsizeu = visheight / 2.0, vishalfsizev = halfsizey, interface = interface)
    s4 = Plane(SVector{3,T}(-1.0, 0.0, 0.0), SVector{3,T}(-halfsizex, 0.0, 0.0), vishalfsizeu = visheight / 2.0, vishalfsizev = halfsizey, interface = interface)
    i1 = csgintersection(s1, s2)
    i2 = csgintersection(s3, s4)
    return csgintersection(i1, i2)
end
export RectangularPrism

"""
    TriangularPrism(side_length::T, visheight::T = 2.0; interface::NullOrFresnel{T} = nullinterface(T)) -> CSGGenerator{T}

Create an infinitely tall triangular prism with axis `(0, 0, 1)`.
For visualization `visheight` is used, **note that this does not fully represent the surface**.
"""
function TriangularPrism(side_length::T, visheight::T = 2.0; interface::NullOrFresnel{T} = NullInterface(T)) where {T<:Real}
    side_length = side_length / 2
    p1 = Plane(SVector{3,T}(-1.0, 0.0, 0.0), SVector{3,T}(-side_length / 2, 0.0, 0.0), interface = interface, vishalfsizev = side_length, vishalfsizeu = visheight / 2)
    p2 = Plane(SVector{3,T}(sind(30), cosd(30), 0.0), SVector{3,T}(side_length / 2 * sind(30), side_length / 2 * cosd(30), 0.0), interface = interface, vishalfsizev = side_length, vishalfsizeu = visheight / 2)
    p3 = Plane(SVector{3,T}(sind(30), -cosd(30), 0.0), SVector{3,T}(side_length / 2 * sind(30), -side_length / 2 * cosd(30), 0.0), interface = interface, vishalfsizev = side_length, vishalfsizeu = visheight / 2)
    return csgintersection(leaf(p1), csgintersection(p2, p3))
end
export TriangularPrism

"""
    Spider(narms::Int, armwidth::T, radius::T, origin::SVector{3,T} = SVector{3,T}(0.0, 0.0, 0.0), normal::SVector{3,T} = SVector{3,T}(0.0, 0.0, 1.0)) -> Vector{Rectangle{T}}

Creates a 'spider' obscuration with `narms` rectangular arms evenly spaced around a circle defined by `origin` and `normal`.
Each arm is a rectangle `armwidth`×`radius`.

e.g. for 3 and 4 arms we get:
```
   |         _|_
  / \\         |
```
"""
function Spider(narms::Int, armwidth::T, radius::T, origin::SVector{3,T} = SVector{3,T}(0.0, 0.0, 0.0), normal::SVector{3,T} = SVector{3,T}(0.0, 0.0, 1.0)) where {T<:Real}
    dθ = 2π / narms
    rects = Vector{Rectangle{T}}(undef, 0)
    for i in 0:(narms - 1)
        θ = i * dθ
        r = Rectangle(armwidth / 2, radius / 2, normal, origin + SVector(radius / 2 * cos(θ), radius / 2 * sin(θ), 0.0), rotationvec = SVector(cos(θ), sin(θ), 0.0), interface = opaqueinterface(T))
        push!(rects, r)
    end
    return rects
end
export Spider


text:
function relorientation(face, simplex)

    v = setdiff(simplex, face)
    length(v) == 1 || return 0

    # find the position of the missing vertex
    v = v[1]
    #i = Base.findfirst(simplex, v)
    i = something(findfirst(isequal(v), simplex),0)
    s = (-1)^(i-1)

    # remove that vertex from the simplex
    face2 = Array{Int}(undef,length(simplex)-1)
    for j in 1 : i-1
        face2[j] = simplex[j]
    end
    for j in i : length(simplex)-1
        face2[j] = simplex[j+1]
    end

    # get the permutation that maps face to face2
    #p = indexin(face, face2)
    #p = [ findfirst(face2,v) for v in face ]
    p = [something(findfirst(isequal(v),face2),0) for v in face]

    return s * levicivita(p) * i
end


text:
"""
`get_data(b, series [; startyear, endyear, catalog])`
Request one or multiple series using the BLS API.

Arguments
---------
* `b`: A `Bls` connection
* `series`: A string, or array of strings, identifying the time series
* `startyear`: A four-digit year identifying the start of the data request. Defaults to
    9 or 19 years before `endyear`, depending on the API version used.
* `endyear`: A four-digit year identifying the end of the data request. Defaults to
    9 or 19 years after `endyear`, depending on the API version used; or, this year, if
    neither `startyear` nor `endyear` is provided.
* `catalog`: Whether to return any available metadata about the series. Defaults to `false`.

Returns
-------
A `BlsSeries`, or an array of `BlsSeries`.
"""
function get_data(b::Bls, series::AbstractString; kwargs...)
   return get_data(b, [series]; kwargs...)
end
function get_data(b::Bls, series::Array{T, 1};
                  startyear::Int = typemin(Int),
                  endyear::Int   = typemin(Int),
                  catalog::Bool  = false) where {T<:AbstractString}
    # Resolve default and user-specified date ranges
    if startyear == endyear == typemin(Int)
        # If neither startyear nor endyear is specified
        endyear = Int(Dates.year(Dates.now()))
        startyear = endyear - (LIMIT_YEARS_PER_QUERY[get_api_version(b)]-1)
    elseif startyear == typemin(Int) && endyear ≠ typemin(Int)
        # If only endyear is specified
        startyear = endyear - (LIMIT_YEARS_PER_QUERY[get_api_version(b)]-1)
    elseif startyear ≠ typemin(Int) && endyear == typemin(Int)
        # If only startyear is specified
        endyear = startyear + (LIMIT_YEARS_PER_QUERY[get_api_version(b)]-1)
    end
    @assert endyear > startyear
    @assert startyear ≤ Int(Dates.year(Dates.now()))

    # Make multiple requests for year range greater than limit
    limit = LIMIT_YEARS_PER_QUERY[get_api_version(b)]
    nrequests = div(endyear-startyear, limit) + 1
    if nrequests > requests_remaining(b)
        println("Insufficient number of requests remaining ", "(", nrequests, " needed ",
            requests_remaining(b), " remaining).")
        return if length(series) > 1
            [EMPTY_RESPONSE() for i in 1:length(series)]
        else
            EMPTY_RESPONSE()
        end
    end

    # Actually make requests to server
    data = []
    for i=1:nrequests
        t0 = startyear + limit * i - limit
        t1 = startyear + limit * i - 1
        if t1 > endyear
            t1 = endyear
        end
        result = _get_data(b, series, t0, t1, catalog)

        # Check that we got a non-empty result
        # TODO this check fails if series metadata has been added to the result
        if isempty(result)
            continue
        end

        # Append to existing results
        if isempty(data)
            data = result
        else
            append_result!(data, result)
        end

    end

    if length(data) == 1
        return data[1]
    else
        return data
    end
end

# Append data frame and catalog information for each series.
function append_result!(data, result)
    for i=1:length(result)
        # Simply don't attempt to append empty results
        if isempty(result[i])
            println("Empty response from server ignored.")
            continue
        end

        @assert result[i].id==data[i].id

        # It's possible that the first df, from 'data', has no rows but differently typed
        # columns from the second df. In this case, we ditch the first df entirely.
        if nrow(data[i].data) > 0
            append!(data[i].data, result[i].data)
        else
            data[i].data = result[i].data
        end

        # Add non-empty catalog strings
        if !isempty(result[i].catalog)
            data[i].catalog = join(data[i].catalog, ". ")
        end
    end
end

function is_catalog_okay(catalog, message)
    # Catalog okay?
    catalog &&
        !isempty(message) &&
        !any(s->occursin(BLS_RESPONSE_CATALOG_FAIL1, s), message) &&
        !any(s->occursin(BLS_RESPONSE_CATALOG_FAIL2, s), message)
end

# Worker method for a single request
function _get_data(b::Bls, series::Array{T,1},
                   startyear::Int, endyear::Int, catalog::Bool) where {T<:AbstractString}
    n_series = length(series)

    # Setup payload
    url     = get_api_url(b);
    headers = Dict("Content-Type" => "application/json")
    payload = Dict("seriesid"     => series,
                   "startyear"    => startyear,
                   "endyear"      => endyear,
                   "catalog"      => catalog)
    key     = get_api_key(b);
    if !isempty(key)
        payload["registrationKey"] = key
    end

    # Submit POST request to BLS
    body = JSON.json(payload)
    response = HTTP.request("POST", url, headers, body)

    # Check if request succeeded
    status = response.status
    if response.status == 200
        response_json = JSON.parse(String(copy(response.body)))
        increment_requests!(b)
    elseif response.status == 202
        # A 202 status code can be returned to note that "Your request is processing.". It's
        # unclear how BLS treats this, so we just try to process anyway. Exceptions are
        # expected.
        response_json = JSON.parse(String(copy(response.body)))
        increment_requests!(b)
    elseif haskey(BLS_STATUS_CODE_REASONS, response.status)
        reason = BLS_STATUS_CODE_REASONS[status]
        reason_http = HttpCommon.STATUS_CODES[status]
        error("API request failed with status $(status) ($(reason_http)): $(reason)")
    else
        if DEBUG log_error(url, payload, headers, response) end
        error("API request failed unexpectedly with status $(status)")
    end

    # Response okay?
    if response_json["status"] ≠ BLS_RESPONSE_SUCCESS
        if DEBUG log_error(url, payload, headers, response) end
        status = response_json["status"]
        message = if haskey(response_json, "message")
            join(hcat(response_json["message"]), ";")
        else
            "<no message returned>"
        end
        println("API request failed with status $(status): $(message)")

        # Return empty response for each series
        return [EMPTY_RESPONSE() for i in 1:n_series]
    end

    catalog_okay = is_catalog_okay(catalog, response_json["message"])

    # Parse response into DataFrames, one for each series
    @assert n_series == length(response_json["Results"]["series"])
    out = Array{BlsSeries,1}(undef, n_series)
    for (i, series) in enumerate(response_json["Results"]["series"])
        seriesID = series["seriesID"]
        catalog_out = if catalog_okay
            series["catalog"]
        else
            ""
        end
        catalog_out = vcat(catalog_out)
        catalog_out = join(catalog_out, ". ")

        data   = map(parse_period_dict, series["data"])
        dates  = reverse([x[1] for x in data],1)
        values = reverse([x[2] for x in data],1)
        df = DataFrame(date=dates, value=values)

        # Data may not be returned in order, for some reason.
        sort!(df)

        out[i] = BlsSeries(seriesID, df, catalog_out)
    end

    return out
end

function parse_period_dict(dict::Dict{T,Any}) where {T<:AbstractString}
    value = parse(Float64, dict["value"])
    year  = parse(Int, dict["year"])

    period = dict["period"]
    # Monthly data
    if occursin(r"M\d\d", period) && period ≠ "M13"
        month = parse(Int, period[2:3])
        date = Dates.Date(year, month, 1)

    # Quarterly data
    elseif occursin(r"Q\d\d", period)
        quarter = parse(Int, period[3])
        date = Dates.Date(year, 3*quarter-2, 1)

    # Annual data
    elseif occursin(r"A\d\d", period)
        date = Dates.Date(year, 1, 1)

    # Not implemented
    else
        error("Data of frequency ", period, " not implemented")
    end

    return (date, value)
end


text:
# ---
# title: 1114. Print in Order
# id: problem1114
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories:
# link: <https://leetcode.com/problems/print-in-order/description/>
# hidden: true
# ---
#
# Suppose we have a class:
#
#
#
#     public class Foo {
#       public void first() { print("first"); }
#       public void second() { print("second"); }
#       public void third() { print("third"); }
#     }
#
#
# The same instance of `Foo` will be passed to three different threads. Thread A
# will call `first()`, thread B will call `second()`, and thread C will call
# `third()`. Design a mechanism and modify the program to ensure that `second()`
# is executed after `first()`, and `third()` is executed after `second()`.
#
#
#
# **Example 1:**
#
#
#
#     Input: [1,2,3]
#     Output: "firstsecondthird"
#     Explanation: There are three threads being fired asynchronously. The input [1,2,3] means thread A calls first(), thread B calls second(), and thread C calls third(). "firstsecondthird" is the correct output.
#
#
# **Example 2:**
#
#
#
#     Input: [1,3,2]
#     Output: "firstsecondthird"
#     Explanation: The input [1,3,2] means thread A calls first(), thread B calls third(), and thread C calls second(). "firstsecondthird" is the correct output.
#
#
#
# **Note:**
#
# We do not know how the threads will be scheduled in the operating system, even
# though the numbers in the input seems to imply the ordering. The input format
# you see is mainly to ensure our tests' comprehensiveness.
#
#
## @lc code=start
using LeetCode

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


text:
export fill_sac_evh!

"""
    writesac(W::SeisEvent[, v=0])

Write all data in SeisEvent structure `W` to auto-generated SAC files. Event
header information is written from W.hdr; W.source is not used as there is no
standard header position for event source information.
"""
function writesac(S::SeisEvent;
  nvhdr::Integer=6,
  v::Integer=KW.v)

  tdata = Array{Float32, 1}(undef, 0)
  reset_sacbuf()

  ev_id  = codeunits(S.hdr.id == "" ? "-12345  " : S.hdr.id)
  ev_id  = ev_id[1:min(length(ev_id),16)]
  BUF.sac_cv[9:length(ev_id)+8] .= ev_id

  # Values from event header
  if S.hdr.loc.lat != 0.0
    setindex!(BUF.sac_fv, Float32(S.hdr.loc.lat), 36)
    setindex!(BUF.sac_dv, S.hdr.loc.lat, 18)
  end
  if S.hdr.loc.lon != 0.0
    setindex!(BUF.sac_fv, Float32(S.hdr.loc.lon), 37)
    setindex!(BUF.sac_dv, S.hdr.loc.lon, 17)
  end
  S.hdr.loc.dep == 0.0 || setindex!(BUF.sac_fv, Float32(S.hdr.loc.dep), 39)
  S.hdr.mag.val == -5.0f0 || setindex!(BUF.sac_fv, Float32(S.hdr.mag.val), 40)
  BUF.sac_cv[9:length(ev_id)+8] .= ev_id
  t_evt = d2u(S.hdr.ot)

  # Ints
  BUF.sac_iv[7] = Int32(7)
  try
    BUF.sac_iv[9] = parse(Int32, S.hdr.id)
  catch err
    @warn(string("Can't write non-integer event ID ", S.hdr.id, " to SAC."))
  end

  for i = 1:S.data.n
    BUF.sac_fv[8] = Float32(t_evt - S.data.t[i][1,2]*μs)
    BUF.sac_dv[2] = t_evt - S.data.t[i][1,2]*μs
    write_sac_channel(S.data, i, nvhdr, "", v)
  end
  return nothing
end

mk_ot!(Ev::SeisEvent, i::Int, o::T) where T<:AbstractFloat = (Ev.hdr.ot =
  u2d(o + Ev.data.t[i][1,2]*μs))

"""
    fill_sac_evh!(Ev::SeisEvent, fname::String; k=i)

Fill (overwrite) values in `Ev.hdr` with data from SAC file `fname`. Keyword
`k=i` specifies the reference channel `i` from which the absolute origin time
`Ev.hdr.ot` is set. Potentially affects header fields `:id`, `:loc` (subfields
.lat, .lon, .dep), and `:ot`.
"""
function fill_sac_evh!(Ev::SeisEvent, fname::String; k::Int=1)
  reset_sacbuf()
  io = open(fname, "r")
  swap = should_bswap(io)
  fv = BUF.sac_fv
  iv = BUF.sac_iv
  cv = BUF.sac_cv

  # read
  seekstart(io)
  fastread!(io, fv)
  fastread!(io, iv)
  fastread!(io, cv)
  if swap == true
    fv .= bswap.(fv)
    iv .= bswap.(iv)
  end
  sac_v = getindex(iv, 7)
  (iv[9] == sac_nul_i)  || (Ev.hdr.id = string(iv[9]))                # id
  (fv[8] == sac_nul_f)  || (mk_ot!(Ev, k, fv[8]))                     # ot
  (fv[36] == sac_nul_f) || (Ev.hdr.loc.lat = Float64(fv[36]))         # lat
  (fv[37] == sac_nul_f) || (Ev.hdr.loc.lon = Float64(fv[37]))         # lon
  (fv[39] == sac_nul_f) || (Ev.hdr.loc.dep = Float64(fv[39]))         # dep
  (fv[40] == sac_nul_f) || (Ev.hdr.mag.val = fv[40])                  # mag

  if sac_v > 6
    fastskip(io, 4*getindex(iv, 10))
    dv = BUF.sac_dv
    fastread!(io, dv)
    swap && (dv .= bswap.(dv))

    # parse doubles 4 (o), 17 (evlo), 18 (evla)
    (dv[4] == sac_nul_d)  || mk_ot!(Ev, k, dv[4])                     # ot
    (dv[17] == sac_nul_d) || (Ev.hdr.loc.lon = Float64(dv[17]))       # lon
    (dv[18] == sac_nul_d) || (Ev.hdr.loc.lat = Float64(dv[18]))       # lat
  end

  reset_sacbuf()
  close(io)
  return nothing
end


text:
## load in the data
ENV["OMP_NUM_THREADS"] = 6
using PowerSpectra
using Healpix
using PyCall, PyPlot
using CSV, DataFrames, LinearAlgebra
using BenchmarkTools
hp = pyimport("healpy")
nmt = pyimport("pymaster")

data_dir = "/home/zequnl/.julia/dev/PowerSpectra/notebooks/data/"
# mask = readMapFromFITS(data_dir * "mask.fits", 1, Float64)
nside = 256
lmax = 3 * nside - 1

flat_beam = SpectralVector(ones(3*nside))
flat_mask = HealpixMap{Float64, RingOrder}(ones(nside2npix(nside)) )

mask = readMapFromFITS("test/example_mask_1.fits", 1, Float64)
clf()
hp.mollview(mask.pixels, title="Simple Mask")
gcf()

##
m1 = Field("143_hm1", mask, flat_mask, flat_beam)
m2 = Field("143_hm2", mask, flat_mask, flat_beam)
workspace = SpectralWorkspace(m1, m2, m1, m2)
@time mcm12 = compute_mcm_TT(workspace, "143_hm1", "143_hm2")

##
m1_P = CovField("143_hm1", mask, mask, flat_mask, flat_mask, flat_mask, flat_beam, flat_beam)
m2_P = CovField("143_hm2", mask, mask, flat_mask, flat_mask, flat_mask, flat_beam, flat_beam)
workspace_P = PolarizedSpectralWorkspace(m1, m2, m1, m2)
@time mcm12 = compute_mcm_EE(workspace, "143_hm1", "143_hm2")
# @time factorized_mcm12 = lu(mcm12.parent);

##
m1_P = CovField("143_hm1", mask, mask, flat_mask, flat_mask, flat_mask, flat_beam, flat_beam)
m2_P = CovField("143_hm2", mask, mask, flat_mask, flat_mask, flat_mask, flat_beam, flat_beam)
workspace_P = PolarizedSpectralWorkspace(m1_P, m2_P, m1_P, m2_P)
@time mcm12 = compute_mcm_TE(workspace_P, "143_hm1", "143_hm2")
# @time factorized_mcm12 = lu(mcm12.parent);

##

using DelimitedFiles
f_0 = nmt.NmtField(mask.pixels, [flat_mask.pixels])
f_2 = nmt.NmtField(mask.pixels, [flat_mask.pixels, flat_mask.pixels])
b = nmt.NmtBin.from_nside_linear(nside, 1)

##
w = nmt.NmtWorkspace()
@time w.compute_coupling_matrix(f_2, f_2, b)
writedlm("test/mcm_EE_diag.txt", diag(w.get_coupling_matrix()[1:4:4*lmax, 1:4:4*lmax])[3:767])

##
w = nmt.NmtWorkspace()
@time w.compute_coupling_matrix(f_0, f_0, b)
writedlm("test/mcm_TT_diag.txt", diag(w.get_coupling_matrix()[1:lmax, 1:lmax])[3:767])

##
w = nmt.NmtWorkspace()
@time w.compute_coupling_matrix(f_0, f_2, b)
writedlm("test/mcm_TE_diag.txt", diag(w.get_coupling_matrix()[1:2:2*lmax, 1:2:2*lmax])[3:767])

# @time w.compute_coupling_matrix(f_0, f_0, b)

##
clf()
plt.plot(diag(w.get_coupling_matrix()[1:2:2*lmax, 1:2:2*lmax]), "-")
plt.plot(diag(mcm12.parent), "--")
# plt.ylim(0.3,0.4)
plt.xlim(0,20)
# plt.yscale("log")
gcf()

##

clf()
plt.plot(diag(w.get_coupling_matrix()[1:4:4*lmax, 1:4:4*lmax])[3:767]  .- diag(mcm12.parent)[3:767], "-")
# plt.xlim(0,10)
gcf()

##

map1 = readMapFromFITS("test/example_map.fits", 1, Float64)

using DelimitedFiles
f_0 = nmt.NmtField(mask.pixels, [map1.pixels])
f_2 = nmt.NmtField(mask.pixels, [flat_mask.pixels, flat_mask.pixels])
b = nmt.NmtBin.from_nside_linear(nside, 1)
w = nmt.NmtWorkspace()
@time w.compute_coupling_matrix(f_0, f_0, b)
writedlm("test/mcm_TT_diag.txt", diag(w.get_coupling_matrix()[1:lmax, 1:lmax])[3:767])

##
b = nmt.NmtBin.from_nside_linear(nside, 1)
cl_00 = nmt.compute_full_master(f_0, f_0, b)

writedlm("test/example_TT_spectrum.txt", cl_00[1,:])

## OLD STUFF

##
#######
flat_mask = HealpixMap{Float64, RingOrder}(ones(nside2npix(nside)) )
theory = CSV.read(data_dir * "theory.csv"; limit=2000)

clf()
plot(theory.cltt .* theory.ell.^2, "-")
plot(1e-2 .* theory.ell.^2, "-")
gcf()

##
m0 = hp.synfast(theory.cltt, nside=nside, verbose=false, pixwin=true, new=true)
m = HealpixMap{Float64, RingOrder}(nside)
m.pixels .= m0;
# map = readMapFromFITS(data_dir * "map.fits", 1, Float64)

wl = SpectralVector(hp.pixwin(nside))
m1 = Field("143_hm1", mask, flat_mask, wl)
m2 = Field("143_hm2", mask, flat_mask, wl)
workspace = SpectralWorkspace(m1, m2, m1, m2)

@time mcm12 = mcm(workspace, "143_hm1", "143_hm1")
@time factorized_mcm12 = cholesky(Hermitian(mcm12.parent));

##
import PowerSpectra: TT

cltt = SpectralVector(convert(Vector, theory.cltt))
spectra = Dict{PowerSpectra.SpectrumName, SpectralVector{Float64, Vector{Float64}}}(
    (TT, "143_hm1", "143_hm1") => cltt,
    (TT, "143_hm2", "143_hm2") => cltt,
    (TT, "143_hm1", "143_hm2") => cltt,
    (TT, "143_hm2", "143_hm1") => cltt)
@time C = compute_covmat(workspace, spectra, factorized_mcm, factorized_mcm,
                         m1, m2, m1, m2);

##
@time Clhat = spectra_from_masked_maps(m * mask, m * mask, factorized_mcm, wl, wl)
ells = collect(0:(length(Clhat)-1));

##
using JLD2
function generate_sim_array(spec, nsims, nside, mask_, fact_mcm)

    result = Array{Float64, 2}(undef, (3 * nside, nsims))
    for i in 1:nsims
        m0 = hp.synfast(spec, nside=nside, verbose=false, pixwin=false, new=true)
        norm = mean(m0)
        m0 .-= norm
        m = HealpixMap{Float64, RingOrder}(nside)
        m.pixels .= m0
        result[:,i] .= spectra_from_masked_maps(m * mask_, m * mask_, factorized_mcm, wl, wl)
    end

    return result
end
# sims = generate_sim_array(theory.cltt, 100, nside, mask, factorized_mcm)
# @save "sims.jld2" sims
@load "sims.jld2" sims

##

clf()
plot(theory.cltt[1:(3*nside)] ./ Clhat, "-")
ylim(0,2)
xlim(0, 2nside)
ylabel(raw"$\hat{C}_{\ell} / C_{\ell}^{\mathrm{theory}}$")
xlabel(raw"Multipole moment, $\ell$")
# plot(1e-2 .* theory.ell.^2, "-")
gcf()

##

##

using Statistics
clf()
plot( (Statistics.var(sims, dims=2) ./ diag(parent(C)))[2:2nside] )
# plot( 1 ./ wl.parent .^ 4)
ylabel(raw"$\mathrm{Var}^{\mathrm{sim}}(C_{\ell}) / \mathrm{Var}^{\mathrm{analytic}}(C_{\ell})$")
xlabel(raw"Multipole moment, $\ell$")
# yscale("log")
gcf()

##
nmt = pyimport("pymaster")
f_0 = nmt.NmtField(mask.pixels, [mask.pixels])
b = nmt.NmtBin.from_nside_linear(nside, 1)
w = nmt.NmtWorkspace()
##

@time w.compute_coupling_matrix(f_0, f_0, b)
##

##

cw = nmt.NmtCovarianceWorkspace()
cw.compute_coupling_coefficients(f_0, f_0, f_0, f_0)
n_ell = length(b.get_effective_ells())
covar_00_00 = nmt.gaussian_covariance(cw,
                                      0, 0, 0, 0,  # Spins of the 4 fields
                                      [cltt],  # TT
                                      [cltt],  # TT
                                      [cltt],  # TT
                                      [cltt],  # TT
                                      w, wb=w)

##

using Statistics
clf()

# plot( ( diag(covar_00_00) ./ diag(parent(C))[3:(3nside)])[2:2nside],
#     label=raw"$\mathrm{Var}^{\mathrm{nmt}}(C_{\ell}) / \mathrm{Var}^{\mathrm{planck}}(C_{\ell})$")

plot( (Statistics.var(sims, dims=2) ./ diag(parent(C)))[1:2nside],
    label=raw"$\mathrm{Var}^{\mathrm{sim}}(C_{\ell}) / \mathrm{Var}^{\mathrm{planck}}(C_{\ell})$")
plot( 2:(2nside+1), (Statistics.var(sims, dims=2)[3:(3nside)] ./ diag(covar_00_00))[1:2nside],
    label=raw"$\mathrm{Var}^{\mathrm{sim}}(C_{\ell}) / \mathrm{Var}^{\mathrm{nmt}}(C_{\ell})$")
legend()
xlabel(raw"Multipole moment, $\ell$")
ylim(0.5,1.5)
# xlim(0,50)
# yscale("log")
gcf()

##
# M .*= sqrt( sum(mask.pixels) / nside2npix(nside) )

# using LinearAlgebra

# m.pixels .*= mask.pixels
# Cl_hat = (alm2cl(map2alm(m; niter=3)))
# Cl_hat = inv(Symmetric(parent(M))) * Cl_hat

# println(sum(Cl_hat[100:300]) / 200)
##
using PyPlot
clf()
axhline(1)
plot(Cl_hat .* (1:length(Cl_hat)).^2 , alpha=0.5, color="C2")
ylim(0,2)
gcf()

##
clf()
hp.mollview(m.pixels)
gcf()

##
using LinearAlgebra
using PyPlot
plt.clf()
plt.plot(diag(parent(M),4))
plt.gcf()
##

PowerSpectra.effective_weights!(workspace, m1, m2, m1, m2)
PowerSpectra.W_spectra!(workspace)

using BenchmarkTools
# @btime cov($workspace, $m1, $m2, $m1, $m2)
##

max_stars_repo_name: xzackli/PlanckCov

max_issues_repo_name: xzackli/PlanckCov


text:
module Graphsearch

using Graphtypes

export
    depth_first_search

function depth_first_search(g::Graph, starting_point::Int64, max_depth::Int64,
    f::Function, blacklist::Array{Int64,1}=Array((Int64),0))

    path = [starting_point]
    visited = [starting_point, blacklist]

    while length(path) != 0
        current = path[end]
        neighbours = list_outgoing(g, current)

        if f(g, current)

            #ignore match if it's a blacklisted starting point
            if !((current == starting_point) && (starting_point in blacklist))
                #expected condition is met
                return path
        end end

        #find next place to be
        next_index = findfirst(
                    n->(!(n in visited)),list_outgoing(g, current))

        if next_index != 0
            #move forward
            push!(path, neighbours[next_index])
            push!(visited, neighbours[next_index])
        else
            #move backward
            pop!(path)
        end
    end

    #search failed
    return []
end

end

max_stars_repo_name: AaltoScienceInst/pebblegame

max_issues_repo_name: AaltoScienceInst/pebblegame

max_forks_repo_name: AaltoScienceInst/pebblegame


text:
include("./Cholesky.jl")
using LinearAlgebra
using IntervalSets
using SpectralGaussianSimulation
using DirectGaussianSimulation

function exactMatrix(A::NoiseCov)
  return  Matrix(A.L * A.L') \ Matrix((A.L * A.L' + spdiagm(0 => A.invσ)) * spdiagm(0 => (1 ./ A.invσ)))
end

function appMatrix(A::NoiseCov)
  return  Matrix(A.L * A.L') \ Matrix(A.U' * A.U * spdiagm(0 => 1 ./ A.invσ))
end

function KL(M,N)
  N = Hermitian(N)
  M = Hermitian(M)
  try
    return 1/2 * (tr(N \ M) + logdet(N) - logdet(M) - size(M, 1))
  catch err
    if isa(err, DomainError)
      return typemax(eltype(M))
    end
  end
end

function KL(M,N)
  N = Hermitian(N)
  M = Hermitian(M)
  try
    return 1/2 * (tr(N \ M) + logdet(N) - logdet(M) - size(M, 1))
  catch err
    if isa(err, DomainError)
      return typemax(eltype(M))
    end
  end
end

# Samples from a given variogram using the FFT, on the unit square. Leaves a padding of length n on either side to decrease artifacts due to fft sampling
function specSampleGrid2D(n, NSamples, padFact::Int, var)
  x = zeros(2, n^2)
  Y = zeros(n^2, NSamples)
  for k in Iterators.partition(1:NSamples, 10)
    simProblem = SimulationProblem(RegularGrid((-padFact * 1.0, -padFact * 1.0), (1.0 + padFact * 1.0, 1.0 + padFact * 1.0), dims=((1 + 2 * padFact) * n, (1 + 2 * padFact) * n)), :z => Float64, length(k))

    simSolver = SpecGaussSim(:z => (variogram=var,))
    simSolution = solve(simProblem, simSolver)

    xRaw = coordinates(domain(simSolution))
    YRaw = reduce(hcat, simSolution.realizations[:z])

    xInds = xRaw .∈ [0.0 .. 1.0]
    xInds = findall(vec(xInds[1, :] .* xInds[2, :]))

    @assert  length(xInds) == n^2

    # only for first iteration
    if k[1] == 1
      x .= xRaw[:, xInds]
      Y[:, k] .= YRaw[xInds, :]
    else
      @assert x == xRaw[:, xInds]
      Y[:, k] .= YRaw[xInds, :]
      GC.gc()
    end
  end
  return x, Y
end


text:
function  plot_models(df::DataFrame, type = :waic;
    fig = "", title)

    models = size(df, 1)
    colors = [:blue, :red, :green, :darkred, :black]

    pars = ["dSE", "WAIC"]

    xmin = 0.95min(minimum(df.WAIC), minimum(df.lppd))
    xmax = 1.1max(maximum(df.WAIC), maximum(df.lppd))
    for r in eachrow(df)
        oos = r.WAIC
        xmin = min(xmin, oos - r.SE)
        xmax = max(xmax, oos + r.SE)
    end

    ylabs = String[]
    for j in 1:models
        l = length(df.models[j])
        str = df.models[j] * repeat(" ", 9-l)
        append!(ylabs, [str])
        for i in 1:2
            l = length(String(pars[i]))
            str = repeat(" ", 9-l) * String(pars[i])
            append!(ylabs, [str])
        end
    end
    reverse!(ylabs)

    ys = [string(ylabs[i]) for i = 1:length(ylabs)]
    p = plot(xlims=(xmin, xmax), leg=false, framestyle=:grid)
    title!(title)
    yran = range(1, stop=length(ylabs), length=length(ys))
    yticks!(yran, ys)

    line = 0
    for r in reverse(eachrow(df))
        line += 1
        hline!([line] .+ length(pars), color=:darkgrey, line=(2, :dash))
        for i in 1:2
            line += 1
            if i == 1
                syms = Symbol.(pars[i])
                ypos = line - 1
                oos = r.WAIC
                lower = oos - r.SE
                upper = oos + r.SE
                plot!([lower, upper], [ypos, ypos], leg=false, color=:black)
                scatter!([oos], [ypos], color=:darkred, markersize=8)
                scatter!([r.lppd], [ypos], color=:darkblue)
            else
                syms = Symbol.(pars[i])
                ypos = line - 1
                is = r.WAIC
                lower = is - r.dSE
                upper = is + r.dSE
                if r.dSE > 0
                    plot!([lower, upper], [ypos, ypos], leg=false, color=:grey)
                    scatter!([is], [ypos], color=:lightgrey,
                        marker = ([:utriangle], 8, 0.8, Plots.stroke(3, :gray)))
                end
            end
        end
     end
    p
end

function plot_models(s::Vector{NamedTuple}, mnames::Vector{String}, pars::Vector{Symbol};
   fig="", title="")

  levels = length(mnames) * (length(pars) + 1)
  colors = [:blue, :red, :green, :darkred, :black]

  xmin = 0; xmax = 0.0
  for i in 1:length(s)
    for par in pars
      syms = Symbol.(keys(s[i]))
      if Symbol(par) in syms
        mp = s[i][par].mean
        xmin = min(xmin, s[i][par].lower)
        xmax = max(xmax, s[i][par].upper)
      end
    end
  end

  ylabs = String[]
  for j in 1:length(mnames)
    for i in 1:length(pars)
      l = length(String(pars[i]))
      str = repeat(" ", 9-l) * String(pars[i])
      append!(ylabs, [str])
    end
    l = length(mnames[j])
    str = mnames[j] * repeat(" ", 9-l)
    append!(ylabs, [str])
  end

  ys = [string(ylabs[i]) for i = 1:length(ylabs)]
  p = plot(xlims=(xmin, xmax), leg=false, framestyle=:grid)
  title!(title)
  yran = range(1, stop=length(ylabs), length=length(ys))
  yticks!(yran, ys)

  line = 0
  for (mindx, model) in enumerate(mnames)
    line += 1
    hline!([line] .+ length(pars), color=:darkgrey, line=(2, :dash))
    for (pindx, par) in enumerate(pars)
      line += 1
      syms = Symbol.(keys(s[mindx]))
      if Symbol(par) in syms
        ypos = (line - 1)
        mp = s[mindx][Symbol(par)].mean
        lower = s[mindx][Symbol(par)].lower
        upper = s[mindx][Symbol(par)].upper
        plot!([lower, upper], [ypos, ypos], leg=false, color=colors[pindx])
        scatter!([mp], [ypos], color=colors[pindx])
      end
    end
  end
  if length(fig) > 0
    savefig(p, fig)
    return s
  end
  p
end

export
    plot_models


text:
################################################################################
#
#  NfAbs/NonSimple.jl : non-simple absolute number fields
#
# This file is part of Hecke.
#
# Copyright (c) 2015, 2016, 2017, 2018: Claus Fieker, Tommy Hofmann
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# * Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
#
#  Copyright (C) 2018 Tommy Hofmann, Claus Fieker
#
################################################################################

export NfAbsNS, NfAbsNSElem

export issimple, simple_extension

@inline base_ring(K::NfAbsNS) = FlintQQ

@inline base_field(K::NfAbsNS) = FlintQQ

@inline degree(K::NfAbsNS) = K.degree

@inline degrees(K::NfAbsNS) = K.degrees

@inline ngens(K::NfAbsNS) = length(K.pol)

function ismaximal_order_known(K::NfAbsNS)
  return has_attribute(K, :maximal_order)
end

################################################################################
#
#  Copy
#
################################################################################

function Base.deepcopy_internal(a::NfAbsNSElem, dict::IdDict)
  # TODO: Fix this once deepcopy is fixed for fmpq_mpoly
  # z = NfAbsNSElem(Base.deepcopy_internal(data(a), dict))
  z = NfAbsNSElem(parent(a), Base.deepcopy(data(a)))
  return z
end

#julia's a^i needs copy
function Base.copy(a::NfAbsNSElem)
  return parent(a)(a.data)
end

################################################################################
#
#  Comply with Nemo ring interface
#
################################################################################

order_type(::NfAbsNS) = NfAbsOrd{NfAbsNS, NfAbsNSElem}

order_type(::Type{NfAbsNS}) = NfAbsOrd{NfAbsNS, NfAbsNSElem}

function iszero(a::NfAbsNSElem)
  reduce!(a)
  return iszero(data(a))
end

function isone(a::NfAbsNSElem)
  reduce!(a)
  return isone(data(a))
end

Nemo.zero(K::NfAbsNS) = K(Nemo.zero(parent(K.pol[1])))

Nemo.one(K::NfAbsNS) = K(Nemo.one(parent(K.pol[1])))

Nemo.one(a::NfAbsNSElem) = one(a.parent)

function Nemo.zero!(a::NfAbsNSElem)
  a.data = zero(a.data)
  return a
end

function Nemo.one!(a::NfAbsNSElem)
  a.data = one(a.data)
  return a
end

function Nemo.one!(a::fmpq_mpoly)
  ccall((:fmpq_mpoly_one, libflint), Nothing,
      (Ref{fmpq_mpoly}, Ref{FmpqMPolyRing}), a, parent(a))
  return a
end

################################################################################
#
#  Random
#
################################################################################

RandomExtensions.maketype(K::NfAbsNS, r) = elem_type(K)

function rand(rng::AbstractRNG, sp::SamplerTrivial{<:Make2{NfAbsNSElem,NfAbsNS,<:UnitRange}})
  K, r = sp[][1:end]
  # TODO: This is super slow
  b = basis(K, copy = false)
  z::Random.gentype(sp) = K() # type-assert to help inference on Julia 1.0 and 1.1
  for i in 1:degree(K)
    z += rand(rng, r) * b[i]
  end
  return z
end

rand(K::NfAbsNS, r::UnitRange) = rand(GLOBAL_RNG, K, r)
rand(rng::AbstractRNG, K::NfAbsNS, r::UnitRange) = rand(rng, make(K, r))

################################################################################
#
#  Basis matrix
#
################################################################################

function basis_matrix(A::Array{NfAbsNSElem})
  @assert length(A) > 0
  n = length(A)
  d = degree(parent(A[1]))

  MM = zero_matrix(FlintQQ, n, d)
  for i in 1:n
    elem_to_mat_row!(MM, i, A[i])
  end
  return MM
end

function basis_matrix(A::Vector{NfAbsNSElem}, ::Type{FakeFmpqMat})
  return FakeFmpqMat(basis_matrix(A))
end

################################################################################
#
#  Field access
#
################################################################################

@inline Nemo.data(a::NfAbsNSElem) = a.data

@inline Nemo.parent(a::NfAbsNSElem) = a.parent::NfAbsNS

issimple(a::NfAbsNS) = false

issimple(::Type{NfAbsNS}) = false

function basis(K::NfAbsNS; copy::Bool = true)
  if isdefined(K, :basis)
    if copy
      return deepcopy(K.basis)::Vector{NfAbsNSElem}
    else
      return K.basis::Vector{NfAbsNSElem}
    end
  end
  Rx = parent(K.pol[1])
  b = Vector{NfAbsNSElem}(undef, degree(K))
  ind = 1
  d = degrees(K)
  it = cartesian_product_iterator([0:d[i]-1 for i = 1:length(d)], inplace = true)
  for i in it
    el = Rx()
    setcoeff!(el, i, fmpq(1))
    b[ind] = K(el, false)
    ind += 1
  end
  K.basis = b
  if copy
    return deepcopy(b)::Vector{NfAbsNSElem}
  else
    return b::Vector{NfAbsNSElem}
  end
end

# Given an exponent vector b, the following function returns the index of
# the basis element corresponding to b.
function monomial_to_index(K::NfAbsNS, b::Vector{T}) where {T}
  n = ngens(K)
  idx = b[n]
  d = degrees(K)
  for j in n-1:-1:1
    idx *= d[j]
    idx += b[j]
  end
  return Int(idx + 1)
end

################################################################################
#
#  Reduction
#
################################################################################

function reduce!(a::NfAbsNSElem)
  q, a.data = divrem(a.data, parent(a).pol)
  return a
end

################################################################################
#
#  Denominator
#
################################################################################

denominator(a::NfAbsNSElem) = denominator(a.data)

################################################################################
#
#  String I/O
#
################################################################################

function Base.show(io::IO, a::NfAbsNS)
  @show_name(io, a)
  @show_special(io, a)
  print(io, "Non-simple number field with defining polynomials ", a.pol)
end

function Base.show(io::IO, a::NfAbsNSElem)
  print(io, AbstractAlgebra.obj_to_string(a, context = io))
end

function AbstractAlgebra.expressify(x::NfAbsNSElem; context = nothing)
  return AbstractAlgebra.expressify(data(x), symbols(parent(x)), context = context)
end

################################################################################
#
#  Unary operators
#
################################################################################

function Base.:(-)(a::NfAbsNSElem)
  return NfAbsNSElem(parent(a), -data(a))
end

################################################################################
#
#  Binary operators
#
################################################################################

function Base.:(+)(a::NfAbsNSElem, b::NfAbsNSElem)
  parent(a) == parent(b) || force_op(+, a, b)::NfAbsNSElem
  return NfAbsNSElem(parent(a), data(a) + data(b))
end

function Base.:(-)(a::NfAbsNSElem, b::NfAbsNSElem)
  parent(a) == parent(b) || force_op(-, a, b)::NfAbsNSElem
  return NfAbsNSElem(parent(a), data(a) - data(b))
end

function Base.:(*)(a::NfAbsNSElem, b::NfAbsNSElem)
  parent(a) == parent(b) || force_op(*, a, b)::NfAbsNSElem
  return parent(a)(data(a) * data(b))
end

function Base.:(//)(a::NfAbsNSElem, b::NfAbsNSElem)
  return div(a, b)
end

function Base.div(a::NfAbsNSElem, b::NfAbsNSElem)
  parent(a) == parent(b) || force_op(div, a, b)::NfAbsNSElem
  return a * inv(b)
end

Nemo.divexact(a::NfAbsNSElem, b::NfAbsNSElem; check::Bool = false) = div(a, b)

################################################################################
#
#  Powering
#
################################################################################

function Base.:(^)(a::NfAbsNSElem, b::Integer)
  if b < 0
    return inv(a)^(-b)
  elseif b == 0
    return one(parent(a))
  elseif b == 1
    return deepcopy(a)
  elseif mod(b, 2) == 0
    c = a^(div(b, 2))
    return c*c
  else#if mod(b, 2) == 1
    return a^(b - 1)*a
  end
end

function Base.:(^)(a::NfAbsNSElem, b::fmpz)
  if b < 0
    return inv(a)^(-b)
  elseif b == 0
    return one(parent(a))
  elseif b == 1
    return deepcopy(a)
  elseif mod(b, 2) == 0
    c = a^(div(b, 2))
    return c*c
  else# mod(b, 2) == 1
    return a^(b - 1)*a
  end
end

################################################################################
#
#  Comparison
#
################################################################################

function Base.:(==)(a::NfAbsNSElem, b::NfAbsNSElem)
  parent(a) == parent(b) || force_op(==, a, b)::Bool
  return data(a) == data(b)
end

################################################################################
#
#  Unsafe operations
#
################################################################################

function Nemo.mul!(c::NfAbsNSElem, a::NfAbsNSElem, b::NfAbsNSElem)
  mul!(c.data, a.data, b.data)
  c = reduce!(c)
  return c
end

function Nemo.add!(c::NfAbsNSElem, a::NfAbsNSElem, b::NfAbsNSElem)
  add!(c.data, a.data, b.data)
  return c
end

function Nemo.add!(c::NfAbsNSElem, a::NfAbsNSElem, b::fmpz)
  add!(c.data, a.data, parent(c.data)(b))
  return c
end

function Nemo.add!(c::NfAbsNSElem, a::NfAbsNSElem, b::Integer)
  add!(c.data, a.data, parent(c.data)(b))
  return c
end

function Nemo.addeq!(b::NfAbsNSElem, a::NfAbsNSElem)
  addeq!(b.data, a.data)
  b = reduce!(b)
  return b
end

function Nemo.mul!(c::NfAbsNSElem, a::NfAbsNSElem, b::fmpz)
  mul!(c.data, a.data, b)
  return c
end

function Nemo.mul!(c::NfAbsNSElem, a::NfAbsNSElem, b::Integer)
  mul!(c.data, a.data, parent(c.data)(b))
  return c
end

################################################################################
#
#  Conversion to matrix
#
################################################################################

function elem_to_mat_row!(M::fmpz_mat, i::Int, d::fmpz, a::NfAbsNSElem)
  K = parent(a)
  # TODO: This is super bad
  # Proper implementation needs access to the content of the underlying
  # fmpq_mpoly

  for j in 1:ncols(M)
    M[i, j] = zero(FlintZZ)
  end

  one!(d)

  if length(data(a)) == 0
    return nothing
  end

  z = zero_matrix(FlintQQ, 1, ncols(M))
  elem_to_mat_row!(z, 1, a)
  z_q = FakeFmpqMat(z)

  for j in 1:ncols(M)
    M[i, j] = z_q.num[1, j]
  end

  ccall((:fmpz_set, libflint), Nothing, (Ref{fmpz}, Ref{fmpz}), d, z_q.den)

  return nothing
end

function elem_to_mat_row!(M::fmpq_mat, i::Int, a::NfAbsNSElem)
  K = parent(a)
  for j in 1:ncols(M)
    M[i, j] = zero(FlintQQ)
  end
  adata = data(a)
  for j in 1:length(adata)
    exps = exponent_vector(adata, j)
    k = monomial_to_index(K, exps)
    M[i, k] = coeff(adata, j)
  end
  return M
end

function elem_from_mat_row(K::NfAbsNS, M::fmpq_mat, i::Int)
  a = K()
  b = basis(K, copy = false)
  for c = 1:ncols(M)
    a += M[i, c]*b[c]
  end
  return a
end

function elem_from_mat_row(K::NfAbsNS, M::fmpz_mat, i::Int, d::fmpz)
  b = basis(K, copy = false)
  Qxy = parent(b[1].data)
  a = Qxy()
  tmp = Qxy()
  for c = 1:ncols(M)
    mul!(tmp, b[c].data, M[i, c])
    add!(a, a, tmp)
    #a += M[i, c]*b[c]
  end
  return divexact(K(a), d)
end

function SRow(a::NfAbsNSElem)
  sr = SRow(FlintQQ)
  adata = data(a)
  for i=1:length(adata)
    # TODO: Do this inplace with preallocated exps array
    exps = exponent_vector(adata, i)
    push!(sr.pos, monomial_to_index(parent(a), exps))
    push!(sr.values, coeff(adata, i))
  end
  p = sortperm(sr.pos)
  sr.pos = sr.pos[p]
  sr.values = sr.values[p]
  return sr
end

################################################################################
#
#  Discriminant
#
################################################################################

function discriminant(K::NfAbsNS)
  Qx = FlintQQ["x"][1]
  d = fmpq(1)
  for i = 1:length(K.pol)
    d *= discriminant(Qx(K.pol[i]))^(div(degree(K), total_degree(K.pol[i])))
  end
  return d
end

################################################################################
#
#  Minimal polynomial
#
################################################################################

function minpoly_dense(a::NfAbsNSElem)
  K = parent(a)
  n = degree(K)
  M = zero_matrix(FlintQQ, degree(K)+1, degree(K))
  z = a^0
  elem_to_mat_row!(M, 1, z)
  z *= a
  elem_to_mat_row!(M, 2, z)
  i = 2
  Qt, _ = PolynomialRing(FlintQQ,"t", cached=false)
  while true
    if n % (i-1) == 0 && rank(M) < i
      N = nullspace(transpose(sub(M, 1:i, 1:ncols(M))))
      @assert N[1] == 1
      v = Vector{fmpq}(undef, i)
      for j in 1:i
        v[j] = N[2][j, 1]
      end
      #f = Qt([N[2][j, 1] for j=1:i])
      f = Qt(v)
      return f*inv(leading_coefficient(f))
    end
    z *= a
    elem_to_mat_row!(M, i+1, z)
    i += 1
  end
end

function minpoly_sparse(a::NfAbsNSElem)
  K = parent(a)
  n = degree(K)
  M = sparse_matrix(FlintQQ)
  z = a^0
  sz = SRow(z)
  i = 0
  push!(sz.values, FlintQQ(1))
  push!(sz.pos, n+i+1)
  push!(M, sz)
  z *= a
  sz = SRow(z)
  i = 1
  Qt, t = PolynomialRing(FlintQQ, "x", cached = false)
  while true
    if n % i == 0
      echelon!(M)
      fl, so = can_solve_ut(sub(M, 1:i, 1:n), sz)
      if fl
        so = mul(so, sub(M, 1:i, n+1:ncols(M)))
        # TH: If so is the zero vector, we cannot use the iteration,
        # so we do it by hand.
        if length(so.pos) == 0
          f = t^i
        else
          f = t^i - sum(c*t^(k-1) for (k,c) = so)
        end
        return f
      end
    end
    push!(sz.values, FlintQQ(1))
    push!(sz.pos, n+i+1)
    push!(M, sz)
    z *= a
    sz = SRow(z)
    i += 1
    if i > degree(parent(a))
      error("too large")
    end
  end
end

function minpoly(a::NfAbsNSElem)
  return minpoly_via_trace(a)::fmpq_poly
end

function minpoly(Qx::FmpqPolyRing, a::NfAbsNSElem)
  return Qx(minpoly(a))
end

function minpoly(Rx::FmpzPolyRing, a::NfAbsNSElem)
  f = minpoly(a)
  if !isone(denominator(f))
    error("element is not integral")
  end
  return Rx(denominator(f)*f)
end

function minpoly(a::NfAbsNSElem, R::FlintIntegerRing)
  return minpoly(PolynomialRing(R, cached = false)[1], a)
end

function minpoly(a::NfAbsNSElem, ::FlintRationalField)
  return minpoly(a)
end

################################################################################
#
#  Characteristic polynomial
#
################################################################################

function charpoly(a::NfAbsNSElem)
  f = minpoly(a)
  return f^div(degree(parent(a)), degree(f))
end

function charpoly(Rx::FmpqPolyRing, a::NfAbsNSElem)
  return Qx(charpoly(a))
end

function charpoly(Rx::FmpzPolyRing, a::NfAbsNSElem)
  f = charpoly(a)
  if !isone(denominator(f))
    error("element is not integral")
  end
  return Rx(denominator(f)*f)
end

function charpoly(a::NfAbsNSElem, R::FlintIntegerRing)
  return charpoly(PolynomialRing(R, cached = false)[1], a)
end

function charpoly(a::NfAbsNSElem, ::FlintRationalField)
  return charpoly(a)
end

################################################################################
#
#  Inverse
#
################################################################################

function inv(a::NfAbsNSElem)
  if iszero(a)
    error("division by zero")
  end
  f = minpoly(a)
  z = parent(a)(coeff(f, degree(f)))
  for i=degree(f)-1:-1:1
    z = z*a + coeff(f, i)
  end
  return -z*inv(coeff(f, 0))
end

################################################################################
#
#  Norm
#
################################################################################

function norm(a::NfAbsNSElem)
  f = minpoly(a)
  return (-1)^degree(parent(a)) * coeff(f, 0)^div(degree(parent(a)), degree(f))
end

################################################################################
#
#  Representation matrix
#
################################################################################

function representation_matrix(a::NfAbsNSElem)
  K = parent(a)
  b = basis(K, copy = false)
  M = zero_matrix(FlintQQ, degree(K), degree(K))
  for i=1:degree(K)
    elem_to_mat_row!(M, i, a*b[i])
  end
  return M
end

function representation_matrix_q(a::NfAbsNSElem)
  M = representation_matrix(a)
  return _fmpq_mat_to_fmpz_mat_den(M)
end

################################################################################
#
#  mod function
#
################################################################################

function mod(a::NfAbsNSElem, p::fmpz)
  b = copy(a)
  mod!(b, p)
  return b
end

# TODO: Dan says that it is better to use a BuilderCtx if the result has
# denominator 1
function mod!(b::NfAbsNSElem, p::fmpz)
  for i=1:length(b.data)
    el = coeff(b.data, i)
    dnew, cp = ppio(denominator(el), p)
    el *= cp
    n = mod(numerator(el), dnew * p)
    setcoeff!(b.data, i, fmpq(n, dnew))
  end
  combine_like_terms!(b.data)
  return b
end

################################################################################
#
#  Substitution
#
################################################################################

#function isunivariate(f::fmpq_mpoly)
#  deg = 0
#  var = 0
#  for i = 1:length(f)
#    exps = exponent_vector(f, i)
#    for j = 1:length(exps)
#      if !iszero(exps[j])
#        if iszero(var)
#          var = j
#          deg = exps[j]
#        elseif var != j
#          return false, fmpq_poly()
#        elseif deg < exps[j]
#          deg = exps[j]
#        end
#      end
#    end
#  end
#
#  Qx = PolynomialRing(FlintQQ, "x")[1]
#  coeffs = Vector{fmpq}(undef, deg+1)
#  if iszero(deg)
#    if iszero(f)
#      coeffs[1] = 0
#      return true, Qx(coeffs)
#    end
#    #f is a constant
#    coeffs[1] = coeff(f, 1)
#    return true, Qx(coeffs)
#  end
#  for i = 1:length(f)
#    exps = exponent_vector(f, i)
#    coeffs[exps[var]+1] = coeff(f, i)
#  end
#  for i = 1:length(coeffs)
#    if !isassigned(coeffs, i)
#      coeffs[i] = fmpq(0)
#    end
#  end
#  return true, Qx(coeffs)
#
#end

# TODO: - Preallocate the exps array
#       - Do we still need this?
function msubst(f::fmpq_mpoly, v::Vector{T}) where {T}
  n = length(v)
  @assert n == nvars(parent(f))
  variables = vars(f)
  if length(f) == 0
    return zero(fmpq) * one(parent(v[1]))
  end
  if length(variables) == 1
    fl = isunivariate(f)
    p = to_univariate(Globals.Qx, f)
    @assert fl
    #I need the variable. Awful
    vect_exp = exponent_vector(variables[1], 1)
    i = 1
    while iszero(vect_exp[i])
      i += 1
    end
    return evaluate(p, v[i])
  end
  powers = Dict{Int, Dict{Int, T}}()
  for i = 1:n
    powers[i] = Dict{Int, T}()
  end
  exps = exponent_vector(f, length(f))
  powers[1][exps[1]] = v[1]^exps[1]
  r = coeff(f, length(f)) * powers[1][exps[1]]
  for j = 2:n
    if iszero(exps[j])
      continue
    end
    el = v[j]^exps[j]
    powers[j][exps[j]] = el
    r *= el
  end
  R = parent(r)
  for i = length(f)-1:-1:1
    exps = exponent_vector(f, i)
    #r += coeff(f, i) * prod(v[j]^exps[j] for j=1:n)
    r1 = coeff(f, i)*one(R)
    for j = 1:n
      if iszero(exps[j])
        continue
      end
      if haskey(powers[j], exps[j])
        r1 *= powers[j][exps[j]]
      else
        el = v[j]^exps[j]
        powers[j][exps[j]] = el
        r1 *= el
      end
    end
    r += r1
  end
  return r
end

################################################################################
#
#  Simple extensions
#
################################################################################

function simple_extension(K::NfAbsNS; cached::Bool = true, check = true, simplified::Bool = false)
  if simplified
    return simplified_simple_extension(K, cached = cached)
  end
  n = ngens(K)
  g = gens(K)
  if n == 1
    #The extension is already simple
    f = to_unvariate(Globals.Qx, K.pol[1])
    Ka, a = NumberField(f, "a", cached = cached, check = check)
    mp = NfAbsToNfAbsNS(Ka, K, g[1], [a])
    return Ka, mp
  end
  pe = g[1]
  i = 1
  ind = Int[1]
  f = minpoly(pe)
  #TODO: use resultants rather than minpoly??
  while i < n
    i += 1
    j = 1
    f = minpoly(pe + j * g[i])
    while degree(f) < prod(total_degree(K.pol[k]) for k in 1:i)
      j += 1
      f = minpoly(pe + j * g[i])
    end
    push!(ind, j)
    pe += j * g[i]
  end
  Ka, a = number_field(f, check = check, cached = cached)
  k = base_ring(K)
  M = zero_matrix(k, degree(K), degree(K))
  z = one(K)
  elem_to_mat_row!(M, 1, z)
  if degree(K) > 1
    elem_to_mat_row!(M, 2, pe)
    z = mul!(z, z, pe)
    for i=3:degree(K)
      z = mul!(z, z, pe)
      elem_to_mat_row!(M, i, z)
    end
  end
  N = zero_matrix(k, n, degree(K))
  for i = 1:n
    elem_to_mat_row!(N, i, g[i])
  end
  s = solve(transpose(M), transpose(N))
  b = basis(Ka)
  emb = Vector{nf_elem}(undef, n)
  for i = 1:n
    emb[i] = zero(Ka)
    for j = 1:degree(Ka)
      emb[i] += b[j] * s[j, i]
    end
  end
  h = hom(Ka, K, pe, inverse = emb)
  embed(h)
  embed(MapFromFunc(x->preimage(h, x), K, Ka))
  return Ka, h
end

function NumberField(K1::AnticNumberField, K2::AnticNumberField; cached::Bool = false, check::Bool = false)
  K , l = number_field([K1.pol, K2.pol], "_\$", check = check, cached = cached)
  mp1 = hom(K1, K, l[1], check = false)
  mp2 = hom(K2, K, l[2], check = false)
  embed(mp1)
  embed(mp2)
  return K, mp1, mp2
end

function NumberField(fields::Vector{AnticNumberField}; cached::Bool = true, check::Bool = true)
  pols = Vector{fmpq_poly}(undef, length(fields))
  for i = 1:length(fields)
    pols[i] = fields[i].pol
  end
  K, gK = number_field(pols, "\$", check = check, cached = cached)
  mps = Vector{NfAbsToNfAbsNS}(undef, length(fields))
  for i = 1:length(fields)
    mps[i] = hom(fields[i], K, gK[i])
    if cached
      embed(mps[i])
    end
  end
  return K, mps
end

################################################################################
#
#  Constructors and parent object overloading
#
################################################################################

@doc Markdown.doc"""
    number_field(f::Vector{fmpq_poly}, s::String="_\$") -> NfAbsNS

Let $f = (f_1, \ldots, f_n)$ be univariate rational polynomials, then
we construct
 $$K = Q[t_1, \ldots, t_n]/\langle f_1(t_1), \ldots, f_n(t_n)\rangle .$$
The ideal must be maximal, however, this is not tested.
"""
function NumberField(f::Vector{fmpq_poly}, s::String="_\$"; cached::Bool = false, check::Bool = true)
  n = length(f)
  if occursin('#', s)
    lS = Symbol[Symbol(replace(s, "#"=>"$i")) for i=1:n]
  else
    lS = Symbol[Symbol("$s$i") for i=1:n]
  end
  return NumberField(f, lS, cached = cached, check = check)
end

function NumberField(f::Vector{fmpq_poly}, s::Vector{String}; cached::Bool = false, check::Bool = true)
  lS = Symbol[Symbol(x) for x=s]
  return NumberField(f, lS, cached = cached, check = check)
end

function NumberField(f::Vector{fmpq_poly}, S::Vector{Symbol}; cached::Bool = false, check::Bool = true)
  length(S) == length(f) || error("number of names must match the number of polynomials")
  n = length(S)
  s = var(parent(f[1]))
  Qx, x = PolynomialRing(FlintQQ, ["$s$i" for i=1:n], cached = false)
  K = NfAbsNS(f, fmpq_mpoly[f[i](x[i]) for i=1:n], S, cached)
  K.degrees = [degree(f[i]) for i in 1:n]
  K.degree = prod(K.degrees)
  if check
    if !_check_consistency(K)
      error("The fields are not linearly disjoint!")
    end
  end
  return K, gens(K)
end

function NumberField(f::Vector{fmpz_poly}, s::String="_\$"; cached::Bool = false, check::Bool = true)
  Qx, _ = PolynomialRing(FlintQQ, var(parent(f[1])), cached = false)
  return NumberField(fmpq_poly[Qx(x) for x = f], s, cached = cached, check = check)
end

function NumberField(f::Vector{fmpz_poly}, s::Vector{String}; cached::Bool = false, check::Bool = true)
  Qx, _ = PolynomialRing(FlintQQ, var(parent(f[1])), cached = false)
  return NumberField(fmpq_poly[Qx(x) for x = f], s, cached = cached, check = check)
end

function NumberField(f::Vector{fmpz_poly}, S::Vector{Symbol}; cached::Bool = false, check::Bool = true)
  Qx, _ = PolynomialRing(FlintQQ, var(parent(f[1])), cached = false)
  return NumberField(fmpq_poly[Qx(x) for x = f], S, cached = cached, check = check)
end

function gens(K::NfAbsNS)
  l = Vector{NfAbsNSElem}(undef, ngens(K))
  degs = degrees(K)
  gQxy = gens(parent(K.pol[1]))
  for i = 1:length(gQxy)
    if isone(degs[i])
      l[i] = K(gQxy[i])
    else
      l[i] = NfAbsNSElem(K, gQxy[i])
    end
  end
  return l
end

function vars(E::NfAbsNS)
  return E.S
end
function symbols(E::NfAbsNS)
  return vars(E)
end

function Base.names(E::NfAbsNS)
  v = vars(E)
  res = Vector{String}(undef, length(v))
  for i = 1:length(res)
    res[i] = string(v[i])
  end
  return res
end

function (K::NfAbsNS)(a::fmpq_mpoly, red::Bool = true)
  if red
    q, a = divrem(a, K.pol)
  end
  z = NfAbsNSElem(K, a)
  return z
end

function (K::NfAbsNS)(a::Vector{fmpq})
  return dot(a, basis(K))
end

(K::NfAbsNS)(a::Integer) = K(parent(K.pol[1])(a))

(K::NfAbsNS)(a::Rational{T}) where {T <: Integer} = K(parent(K.pol[1])(a))

(K::NfAbsNS)(a::fmpz) = K(parent(K.pol[1])(a))

(K::NfAbsNS)(a::fmpq) = K(parent(K.pol[1])(a))

(K::NfAbsNS)() = zero(K)

(K::NfAbsNS)(a::NumFieldElem) = force_coerce(K, a)

function (K::NfAbsNS)(a::NfAbsNSElem)
  if parent(a) === K
    return deepcopy(a)
  end
  error("not compatible")
end

function show_sparse_cyclo(io::IO, a::NfAbsNS)
  print(io, "Sparse cyclotomic field of order $(get_attribute(a, :cyclo))")
end

function cyclotomic_field(::Type{NonSimpleNumField}, n::Int; cached::Bool = false)
  lf = factor(n)
  x = gen(Hecke.Globals.Zx)
  lp = [cyclotomic(Int(p^k), x) for (p,k) = lf.fac]
  ls = ["z($n)_$(p^k)" for (p,k) = lf.fac]
  C, g = number_field(lp, ls, cached = cached, check = false)
  set_attribute!(C, :show => show_sparse_cyclo, :cyclo => n)
  return C, g
end

function trace_assure(K::NfAbsNS)
  if isdefined(K, :traces)
    return
  end
  Qx, x = PolynomialRing(FlintQQ, cached = false)
  K.traces = [polynomial_to_power_sums(Qx(f), total_degree(f)-1) for f = K.pol]
end

#= Idea
  if k = Q[x,y]/<f, g>
    then
      tr(x^i) = power_sums(f)
      tr(y^i) = power_sums(g)
      tr(x^i y^j) = tr(x^i) tr(y^j):
        in the tower of fields:
          tr_<f,g>(xy) = tr_<f>(x (tr_<g> y)) = tr_f x * tr_g y
  so trace_assure computes trace(x^i)
  and tr assembles....
=#

function tr(a::NfAbsNSElem)
  k = parent(a)
  if iszero(a)
    return fmpq()
  end
  trace_assure(k)
  t = fmpq()
  for trm = terms(a.data)
    c = coeff(trm, 1)::fmpq
    e = exponent_vector(trm, 1)
    tt = fmpq(1)
    for i=1:length(e)
      if e[i] != 0
        tt *= k.traces[i][e[i]]
      else
        tt *= total_degree(k.pol[i])
      end
    end
    t += c*tt
  end
  return t
end

#TODO:
#  test f mod p first
#  if all polys are monic, the test if traces have non-trivial gcd
function minpoly_via_trace(a::NfAbsNSElem)
  k = parent(a)
  d = degree(k)
  b = a
  l = [tr(b)]
  i = 1
  while i <= d
    while d % i != 0
      b *= a
      push!(l, tr(b))
      i += 1
    end
    q = fmpq(1, div(d, i))
    f = power_sums_to_polynomial([x*q for x = l])
    if iszero(subst(f, a))  #TODO: to checks first...
      return f::fmpq_poly
    end
    b *= a
    push!(l, tr(b))
    i += 1
  end
  error("cannot happen")
end

function isnorm_divisible(a::NfAbsNSElem, n::fmpz)
  return iszero(mod(norm(a), n))
end

function valuation(a::NfAbsOrdElem, p::NfAbsOrdIdl)
  i = 1
  q = p
  while true
    if !(a in q)
      return i - 1
    end
    q = q * p
    i = i + 1
  end
end

#TODO: find a better algo.
function degree(a::NfAbsNSElem)
  return degree(minpoly(a))
end

function degree(a::nf_elem)
  return degree(minpoly(a))
end

#TODO: Improve the algorithm
function primitive_element(K::NfAbsNS)
  g = gens(K)
  pe = g[1]
  d = total_degree(K.pol[1])
  i = 1
  while i < length(g)
    i += 1
    d *= total_degree(K.pol[i])
    while true
      pe += g[i]
      f = minpoly(pe)
      degree(f) == d && break
    end
  end
  return pe
end

@doc Markdown.doc"""
    factor(f::PolyElem{NfAbsNSElem}) -> Fac{Generic.Poly{NfAbsNSElem}}

The factorisation of $f$ (using Trager's method).
"""
function factor(f::PolyElem{NfAbsNSElem})
  Kx = parent(f)
  K = base_ring(f)

  iszero(f) && error("poly is zero")

  if degree(f) == 0
    r = Fac{typeof(f)}()
    r.fac = Dict{typeof(f), Int}()
    r.unit = Kx(leading_coefficient(f))
    return r
  end

  f_orig = deepcopy(f)
  @vprint :PolyFactor 1 "Factoring $f\n"
  @vtime :PolyFactor 2 g = gcd(f, derivative(f))
  if degree(g) > 0
    f = div(f, g)
  end

  if degree(f) == 1
    multip = div(degree(f_orig), degree(f))
    r = Fac{typeof(f)}()
    r.fac = Dict{typeof(f), Int}(f*(1//leading_coefficient(f)) => multip)
    r.unit = one(Kx) * leading_coefficient(f_orig)
    return r
  end

  f = f*(1//leading_coefficient(f))

  k = 0
  g = f
  @vtime :PolyFactor 2 N = norm(g)

  pe = K()
  while isconstant(N) || !issquarefree(N)
    k = k + 1
    if k == 1
      pe = primitive_element(K)
    end
    g = compose(f, gen(Kx) - k*pe)
    @vtime :PolyFactor 2 N = norm(g)
  end
  @vtime :PolyFactor 2 fac = factor(N)

  res = Dict{PolyElem{NfAbsNSElem}, Int64}()

  for i in keys(fac.fac)
    t = change_base_ring(K, i, parent = Kx)
    t = compose(t, gen(Kx) + k*pe)
    @vtime :PolyFactor 2 t = gcd(f, t)
    res[t] = 1
  end

  r = Fac{typeof(f)}()
  r.fac = res
  r.unit = Kx(1)

  if f != f_orig
    error("factoring with mult not implemented")
  end
  r.unit = one(Kx)* leading_coefficient(f_orig)//prod((leading_coefficient(p) for (p, e) in r))
  return r
end

################################################################################
#
#  Hashing
#
################################################################################

function Base.hash(a::NfAbsNSElem, h::UInt)
  return Base.hash(a.data, h)
end


text:
## Loading FCB test data

# must be imported first to avoid unrelated libz error
using ImageMagick

using HDF5, MAT

datadir = joinpath(ENV["HOME"], "data", "fcb")

file = matopen(joinpath(datadir, "test_fcb_col_notes.mat"))
fcbCols = read(file, "fcb_data_col_notes")
close(file)

file = matopen(joinpath(datadir, "test_fcb_class.mat"))
fcbClass = read(file, "fcb_class")
close(file)

file = matopen(joinpath(datadir, "test_fcb_time.mat"))
fcbTime = read(file, "time")
close(file)

file = matopen(joinpath(datadir, "test_fcb_data.mat"))
fcbData = read(file, "fcb_data")
close(file)

## Basic plotting of data
using Makie

# don't open a separate Makie plot window -- use atom plot pane instead
AbstractPlotting.inline!(true)

## Makie basic plotting example
# generate some data
# x = rand(10)
# y = rand(10)
# colors = rand(10)
# scene = scatter(x, y, color = colors, markersize = 0.05*ones(length(x)))
# scene = scatter(x, y, markersize = 0.05*ones(length(x)))

frame = 1

# plot(fcb_data{1}(:,4),fcb_data{1}(:,1),'.') %PE vs SSC
x = log10.(fcbData[frame][:,4])
y = log10.(fcbData[frame][:,1])

# prep color information
colo = fcbClass[frame][:]
colo[isnan.(colo)] .= 7.0
@assert sum(isnan.(colo)) == 0
colo ./= 8.0

scene = scatter(x, y, color=colo, markersize = 0.03*ones(length(x)))

## subselect two classes for development

frame = 1

msk1 = fcbClass[frame][:] .== 1
msk2 = fcbClass[frame][:] .== 6

msk = msk1 .| msk2

x = log10.(fcbData[frame][msk,4])
y1 = log10.(fcbData[frame][msk,1])
y3 = log10.(fcbData[frame][msk,3])

colo = fcbClass[frame][msk]
colo ./= 7.0

scene = scatter(x, y1, color=colo, markersize = 0.03*ones(length(x)))

scene = scatter(x, y3, color=colo, markersize = 0.03*ones(length(x)))

#


text:
test("Ranks and permutations") do

  test("ranks") do

    @t FeldtLib.ranks([1]) == [1]

    @t FeldtLib.ranks([1, 2]) == [1, 2]
    @t FeldtLib.ranks([2, 1]) == [2, 1]

    @t FeldtLib.ranks([3, 2, 1]) == [3, 2, 1]
    @t FeldtLib.ranks([3, 1, 2]) == [3, 1, 2]
    @t FeldtLib.ranks([2, 1, 3]) == [2, 1, 3]
    @t FeldtLib.ranks([1, 2, 3]) == [1, 2, 3]
    @t FeldtLib.ranks([1, 3, 2]) == [1, 3, 2]

    @t FeldtLib.ranks([4.0, 2.5, 1.2, 5.6]) == [3, 2, 1, 4]

    @t FeldtLib.ranks([7.0, 3.5, 2.2, 10.6, 456.0]) == [3, 2, 1, 4, 5]

  end

  test("factorial") do
    @t factorial(1) == 1
    @t factorial(2) == 2
    @t factorial(3) == 6
    @t factorial(4) == 24
    @t factorial(5) == 120
  end

  test("combinations") do
    @t combinations(7, 14) == 3432
  end
end


text:
import Base.@deprecate

import Base: real, abs, size

@deprecate real(fp::AbstractPath) canonicalize(fp)
@deprecate abs(fp::AbstractPath) absolute(fp)
@deprecate isabs(fp::AbstractPath) isabsolute(fp)
@deprecate move(src::AbstractPath, dest::AbstractPath; kwargs...) mv(src, dest; kwargs...)
@deprecate remove(fp::AbstractPath; kwargs...) rm(fp; kwargs...)
@deprecate size(fp::AbstractPath) filesize(fp)