text:
####
#### Output Structure
####
"""
    struct RNNOutput
        y
        hidden
        memory
        indices
    end

Outputs of the RNN models are always `RNNOutput`
`hidden`,`memory` and `indices` may be nothing depending on the kwargs you used in forward.

`y` is last hidden states of each layer. `size(y)=(H/2H,[B,T])`.
If you use unequal length instances in a batch input, `y` becomes 2D array `size(y)=(H/2H,sum_of_sequence_lengths)`.
See `indices` and `PadRNNOutput` to get correct time outputs for a specific instance or to pad whole output.

`h` is the hidden states in each timesstep. `size(h) = h,B,L/2L`

`c` is the hidden states in each timesstep. `size(h) = h,B,L/2L`

`indices` is corresponding instace indices for your `RNNOutput.y`. You may call `yi = y[:,indices[i]]`.
"""
struct RNNOutput{T,V,Z}
    y::T
    hidden::V
    memory::Z
    indices::Union{Vector{Vector{Int}},Nothing}
end

####
#### RNN Types
####
"""
    SRNN(;input=inputSize, hidden=hiddenSize, activation=:relu, options...)
    LSTM(;input=inputSize, hidden=hiddenSize, options...)
    GRU(;input=inputSize, hidden=hiddenSize, options...)

    (1) (l::T)(x; kwargs...) where T<:AbstractRNN
    (2) (l::T)(x::Array{Int}; batchSizes=nothing, kwargs...) where T<:AbstractRNN
    (3) (l::T)(x::Vector{Vector{Int}}; sorted=false, kwargs...) where T<:AbstractRNN

All RNN layers has above forward run(1,2,3) functionalities.

(1) `x` is an input array with size equals d,[B,T]

(2) For this You need to have an RNN with embedding layer.
`x` is an integer array and inputs coressponds one hot vector indices.
You can give 2D array for minibatching as rows corresponds to one instance.
You can give 1D array with minibatching by specifying batch batchSizes argument.
Checkout `Knet.rnnforw` for this.

(3) For this You need to have an RNN with embedding layer.
`x` is an vector of integer vectors. Every integer vector corresponds to an
instance. It automatically batches inputs. It is better to give inputs as sorted.
If your inputs sorted you can make `sorted` argument true to increase performance.

see RNNOutput

# options

* `embed=nothing`: embedding size or and embedding layer
* `numLayers=1`: Number of RNN layers.
* `bidirectional=false`: Create a bidirectional RNN if `true`.
* `dropout=0`: Dropout probability. Ignored if `numLayers==1`.
* `skipInput=false`: Do not multiply the input with a matrix if `true`.
* `dataType=eltype(KnetLayers.arrtype)`: Data type to use for weights. Default is Float32.
* `algo=0`: Algorithm to use, see CUDNN docs for details.
* `seed=0`: Random number seed for dropout. Uses `time()` if 0.
* `winit=xavier`: Weight initialization method for matrices.
* `binit=zeros`: Weight initialization method for bias vectors.
* `usegpu=(KnetLayers.arrtype <: KnetArray)`: GPU used by default if one exists.

# Keywords

* hx=nothing : initial hidden states
* cx=nothing : initial memory cells
* hy=false   : if true returns h
* cy=false   : if true returns c

"""
AbstractRNN

for layer in (:SRNN, :LSTM, :GRU)
    layername=string(layer)
    @eval begin
        mutable struct $layer <: AbstractRNN
            embedding::Union{Nothing,Embed}
            params
            specs::Knet.RNN
            gatesview::Dict
        end

        (m::$layer)(x,h...;o...) = RNNOutput(_forw(m,x,h...;o...)...)

        function $layer(;input::Integer, hidden::Integer, embed=nothing, activation=:relu,
                         usegpu=(arrtype <: KnetArray), dataType=eltype(arrtype), o...)
            embedding,inputSize = _getEmbed(input,embed)
            rnnType = $layer==SRNN ? activation : Symbol(lowercase($layername))
            r = Knet.RNN(inputSize, hidden; rnnType=rnnType, usegpu=usegpu, dataType=dataType, o...)
            gatesview  = Dict{Symbol,Any}()
            for l in 1:r.numLayers, d in 0:r.direction
                for (g,id) in gate_mappings($layer)
                    for (ih,ihid) in input_mappings, (ty,param) in param_mappings
                         gatesview[Symbol(ty,ih,g,l,d)] = rnnparam(r, (r.direction+1)*(l-1)+d+1, id[ihid], param; useview=true)
                    end
                end
            end
            $layer(embedding,r.w,r,gatesview)
        end
    end
end
gate_mappings(::Type{SRNN}) = Dict(:h=>(1,2))
gate_mappings(::Type{GRU})  = Dict(:r=>(1,4),:u=>(2,5),:n=>(3,6))
gate_mappings(::Type{LSTM}) = Dict(:i=>(1,5),:f=>(2,6),:n=>(3,7),:o=>(4,8))
const input_mappings = Dict(:i=>1,:h=>2)
const param_mappings = Dict(:w=>1,:b=>2)

####
#### Utils
####
"""

    PadSequenceArray(batch::Vector{Vector{T}}) where T<:Integer

Pads a batch of integer arrays with zeros

```
julia> PadSequenceArray([[1,2,3],[1,2],[1]])
3×3 Array{Int64,2}:
 1  2  3
 1  2  0
 1  0  0
 ```

"""
function PadSequenceArray(batch::Vector{Vector{T}}; pad=0) where T<:Integer
    B      = length(batch)
    lngths = length.(batch)
    Tmax   = maximum(lngths)
    padded = Array{T}(undef,B,Tmax)
    @inbounds for n = 1:B
        padded[n,1:lngths[n]] = batch[n]
        padded[n,lngths[n]+1:end] .= pad
    end
    return padded
end

"""

    PadRNNOutput(s::RNNOutput)
Pads a rnn output if it is produces by unequal length batches
`size(s.y)=(H/2H,sum_of_sequence_lengths)` becomes `(H/2H,B,Tmax)`

"""
function PadRNNOutput(s::RNNOutput)
    s.indices == nothing && return s,nothing
    d = size(s.y,1)
    B = length(s.indices)
    lngths = length.(s.indices)
    Tmax = maximum(lngths)
    mask = trues(d,B,Tmax)
    cw = []
    @inbounds for i=1:B
        y1 = s.y[:,s.indices[i]]
        df = Tmax-lngths[i]
        if df > 0
            mask[:,:,end-df+1:end] .= false
            cpad = zeros(Float32,d*df) # zeros(Float32,2d,df)
            kpad = atype(cpad)
            ypad = reshape(cat1d(y1,kpad),d,Tmax) # hcat(y1,kpad)
            push!(cw,ypad)
        else
            push!(cw,y1)
        end
    end
    RNNOutput(reshape(vcat(cw...),d,B,Tmax),s.hidden,s.memory,nothing),mask
end

function _pack_sequence(batch::Vector{Vector{T}}) where T<:Integer
    tokens = Int[]
    bsizes = Int[]
    B      = length(batch)
    Lmax   = length(first(batch))
    @inbounds for t = 1:Lmax
        bs = 0
        @inbounds for k = 1:B
            if t<=length(batch[k])
                push!(tokens,batch[k][t])
                bs += 1
            end
        end
        push!(bsizes,bs)
    end
    return tokens,bsizes
end

function _batchSizes2indices(batchSizes)
    B = batchSizes[1]
    inds = Vector{Int}[]
    for i=1:B
        ind = i.+cumsum(filter(x->(x>=i),batchSizes)[1:end-1])
        push!(inds,append!(Int[i],ind))
    end
    return inds
end

_getEmbed(input::Int,embed::Nothing) = (nothing,input)
_getEmbed(input::Int,embed::Embed)   = size(embed.w,2) == input ? (embed,input) : error("dimension mismatch in your embedding")
_getEmbed(input::Int,embed::Integer) = (Embed(input=input,output=embed),embed)

function _forw(rnn::AbstractRNN,seq::Array{T},h...;batchSizes=nothing,o...) where T<:Integer
    rnn.embedding === nothing && error("rnn has no embedding!")
    ndims(seq) == 1 && batchSizes === nothing && (seq = reshape(seq,1,length(seq)))
    y,hidden,memory,_ = rnnforw(rnn.specs,rnn.params,rnn.embedding(seq),h...;batchSizes=batchSizes,o...)
    return y,hidden,memory,nothing
end

function _forw(rnn::AbstractRNN,batch::Vector{Vector{T}},h...;sorted=true,o...) where T<:Integer
    rnn.embedding === nothing && error("rnn has no embedding!")
    if all(length.(batch).==length(batch[1]))
        return _forw(rnn,cat(batch...;dims=2)',h...;o...)
    end
    if !sorted
        v   = sortperm(batch;by=length,rev=true,alg=MergeSort)
        rev = sortperm(v;alg=MergeSort)
        batch= batch[v]
    end
    tokens,bsizes = _pack_sequence(batch)
    y,hidden,memory,_  = _forw(rnn,tokens,h...;batchSizes=bsizes,o...)
    inds    = _batchSizes2indices(bsizes);
    if !sorted
        inds=inds[rev];
        hidden!=nothing && (hidden=_sort3D(hidden,rev))
        memory!=nothing && (memory=_sort3D(memory,rev))
    end
    y,hidden,memory,inds
end

function _forw(rnn::AbstractRNN,x,h...;o...)
    if rnn.embedding === nothing
        y,hidden,memory,_ = rnnforw(rnn.specs,rnn.params,x,h...;o...)
    else
        y,hidden,memory,_ = rnnforw(rnn.specs,rnn.params,rnn.embedding(x),h...;o...)
    end
    y,hidden,memory,nothing
end

_sort3D(hidden::Array,inds) = hidden[:,inds,:]
function _sort3D(h::KnetArray,inds)
    container = [];
    for i=1:size(h,3)
        push!(container,h[:,:,i][:,inds])
    end
    reshape(cat1d(container...),size(h))
end


text:
using MuseInference, Soss, Turing, Distributions, Random,
    Zygote, MeasureTheory, Test, StableRNGs, LinearAlgebra

##

rng = StableRNG(0)

@testset "MuseInference" begin

    @testset "Turing w/ Distributions" begin

        Turing.@model function turing_funnel()
            θ ~ Distributions.Normal(0, 3)
            z ~ Distributions.MvNormal(zeros(512), exp(θ) * I)
            x ~ Distributions.MvNormal(z, I)
        end

        @testset "$name" for (name,autodiff) in [
            ("ForwardDiff", MuseInference.ForwardDiffBackend()),
            ("Zygote", MuseInference.ZygoteBackend())
        ]

            (;x) = rand(copy(rng), turing_funnel() | (θ=0,))
            prob = TuringMuseProblem(turing_funnel() | (;x); autodiff)
            MuseInference.check_self_consistency(prob, (θ=1,), has_volume_factor=true, rng=copy(rng))
            result = muse(prob, (θ=1,); rng=copy(rng), get_covariance=true)
            @test result.dist.μ / result.dist.σ < 2

        end

    end

    ##

    @testset "Soss w/ Distributions" begin

        soss_funnel = Soss.@model (σ) begin
            θ ~ Distributions.Normal(0, σ)
            z ~ Distributions.MvNormal(zeros(512), exp(θ) * I)
            x ~ Distributions.MvNormal(z, I)
        end

        @testset "$name" for (name,autodiff) in [
            ("ForwardDiff", MuseInference.ForwardDiffBackend()),
            ("Zygote", MuseInference.ZygoteBackend())
        ]

            (;x) = Soss.predict(copy(rng), soss_funnel(3), θ=0)
            prob = SossMuseProblem(soss_funnel(3) | (;x); autodiff)
            MuseInference.check_self_consistency(prob, (θ=1,), has_volume_factor=false, rng=copy(rng))
            result = muse(prob, (θ=1,); rng=copy(rng), get_covariance=true)
            @test result.dist.μ / result.dist.σ < 2

        end

    end

    @testset "Soss w/ MeasureTheory" begin

        soss_funnel = Soss.@model (σ) begin
            θ ~ MeasureTheory.Normal(0, σ)
            z ~ MeasureTheory.Normal(0, exp(θ/2)) ^ 512
            x ~ For(z) do zᵢ
                MeasureTheory.Normal(zᵢ, 1)
            end
        end

        @testset "$name" for (name,autodiff) in [
            ("ForwardDiff", MuseInference.ForwardDiffBackend()),
            # ("Zygote", MuseInference.ZygoteBackend()) # broken
        ]

            (;x) = Soss.predict(copy(rng), soss_funnel(3), θ=0)
            prob = SossMuseProblem(soss_funnel(3) | (;x); autodiff)
            MuseInference.check_self_consistency(prob, (θ=1,), has_volume_factor=false, rng=copy(rng))
            result = muse(prob, (θ=1,); rng=copy(rng), get_covariance=true)
            @test result.dist.μ / result.dist.σ < 2

        end

    end

end

##

if get(ENV, "BUILD_DOCS", false) != false
    include(joinpath(@__DIR__, "../docs/make.jl"))
end

ext: jl


text:
@testset "Paths" begin
  grid = RegularGrid{Float64}(100,100)

  # make sure iteration works
  for path in [SimplePath(grid), RandomPath(grid)]
    count = 0
    for location in path
      count += 1
    end
    @test count == length(path)
  end
end


text:
for func in [:clblasSrotg, :clblasDrotg, :clblasCrotg, :clblasZrotg]

    @eval @blasfun $func(CA::cl.CL_mem, offCA::Csize_t,
                              CB::cl.CL_mem, offCB::Csize_t,
                              C::cl.CL_mem, offC::Csize_t,
                              S::cl.CL_mem, offS::Csize_t,
                              n_queues::cl.CL_uint,
                              queues::Ptr{cl.CL_command_queue},
                              n_events_in_wait_list::cl.CL_uint,
                              event_wait_list::Ptr{cl.CL_event},
                              events::Ptr{cl.CL_event})

    @eval @blasfun2 $func(CA::cl.CL_mem, offCA::Csize_t,
                               CB::cl.CL_mem, offCB::Csize_t,
                               C::cl.CL_mem, offC::Csize_t,
                               S::cl.CL_mem, offS::Csize_t)

end


text:
export conductor, isconductor

########################################################################################
#
#  Tools for conductor
#
########################################################################################

function _norm_group_gens_small(C::ClassField)

  mp=C.mq

  mR = C.rayclassgroupmap
  mS = C.quotientmap

  R=domain(mR)
  cond=mR.modulus_fin
  inf_plc1=mR.modulus_inf
  O=parent(cond).order
  E=order(domain(mp))
  expo=Int(exponent(domain(mp)))
  K=O.nf

  mS=inv(mS)
  dom=domain(mS)
  M=zero_matrix(FlintZZ,ngens(dom), ngens(codomain(mS)))
  for i=1:ngens(dom)
    elem=mS(dom[i]).coeff
    for j=1:ngens(codomain(mS))
      M[i,j]=elem[1,j]
    end
  end
  S1=Hecke.GrpAbFinGenMap(domain(mS),codomain(mS),M)
  T,mT=Hecke.kernel(S1)

  Sgens=find_gens_sub(mR,mT)

  return Sgens

end

#
#  Find small primes generating a subgroup of the ray class group
#

function find_gens_sub(mR::MapRayClassGrp, mT::GrpAbFinGenMap)

  O = order(codomain(mR))
  R = domain(mR)
  T = domain(mT)
  m = Hecke._modulus(mR)
  l = minimum(m)
  lp = NfOrdIdl[]
  sR = GrpAbFinGenElem[]

  if isdefined(mR, :prime_ideal_cache)
    S = mR.prime_ideal_cache
  else
    S = prime_ideals_up_to(O, 1000)
    mR.prime_ideal_cache = S
  end
  q, mq = quo(T, sR, false)
  for (i,P) in enumerate(S)
    if divisible(l,P.minimum)
      continue
    end
    if haskey(mR.prime_ideal_preimage_cache, P)
      f = mR.prime_ideal_preimage_cache[P]
    else
      f = mR\P
      mR.prime_ideal_preimage_cache[P] = f
    end
    bool, pre = haspreimage(mT, f)
    if !bool
      continue
    end
    if iszero(mq(pre))
      continue
    end
    push!(sR, pre)
    push!(lp, P)
    q, mq = quo(T, sR, false)
    if order(q) == 1
      break
    end
  end
  if order(q) == 1
    return lp
  else
    error("Not enough primes")
  end
end

#
#  This functions constructs generators for 1+p^u/1+p^u+1
#

function _1pluspk_1pluspk1(K::AnticNumberField, p::NfOrdIdl, pk::NfOrdIdl, pv::NfOrdIdl, lp::Dict{NfOrdIdl, Int}, prime_power::Dict{NfOrdIdl, NfOrdIdl}, a::Int, n::Int)

  O=maximal_order(K)
  b=basis(pk)
  N = basis_mat(pv, Val{false})*basis_mat_inv(pk, Val{false})
  G=AbelianGroup(N.num)
  S,mS=snf(G)
  #Generators
  gens=Array{NfOrdElem,1}(ngens(S))
  for i=1:ngens(S)
    gens[i]=O(0)
    for j=1:ngens(G)
      gens[i]+=mod(mS.map[i,j], S.snf[end])*b[j]
    end
  end
  for i=1:ngens(S)
    gens[i]+=1
  end
  if length(lp) > 1
    i_without_p = ideal(O,1)
    for (p2,vp2) in lp
      (p != p2) && (i_without_p *= prime_power[p2])
    end

    alpha, beta = idempotents(prime_power[p],i_without_p)
    for i in 1:length(gens)
      gens[i] = beta*gens[i] + alpha
    end
  end
  if mod(n,2)==0
    for i=1:length(gens)
      gens[i]=make_positive(gens[i],a)
    end
  end
  return gens
end

#######################################################################################
#
#  Conductor functions
#
#######################################################################################

doc"""
***
  conductor(C::Hecke.ClassField) -> NfOrdIdl, Array{InfPlc,1}

> Return the conductor of the abelian extension corresponding to C
***
"""
function conductor(C::Hecke.ClassField)

  if isdefined(C,:conductor)
    return C.conductor
  end
  mR = C.rayclassgroupmap
  mS = C.quotientmap
  mp = inv(mS)* mR
  G = domain(mp)
  #
  #  First, we need to find the subgroup
  #

  R = domain(mR)
  cond = mR.modulus_fin
  inf_plc = mR.modulus_inf
  O = parent(cond).order
  if isone(cond) && isempty(inf_plc)
    return ideal(O,1), InfPlc[]
  end
  E=order(G)
  expo=Int(exponent(G))
  K=O.nf

  #
  #  Some of the factors of the modulus are unnecessary for order reasons:
  #
  L=deepcopy(mR.fact_mod)
  for (p,vp) in L
    if !divisible(E,p.minimum)
      if gcd(E, norm(p)-1)==1
        Base.delete!(L,p)
      else
        L[p]=1
      end
    else
      if L[p]==1
        Base.delete!(L,p)
      end
    end
  end

  prime_power=Dict{NfOrdIdl, NfOrdIdl}()
  for (p,v) in mR.fact_mod
    prime_power[p]=p^v
  end

  if !isempty(inf_plc)
    totally_positive_generators(mR,Int(cond.gen_one))
  end

  #Finite part of the modulus
  if !isempty(L)
    tmg=mR.tame_mult_grp
    wild=mR.wild_mult_grp
  end
  for (p,v) in L
    if v==1
      Q,mQ=quo(G,[preimage(mp, ideal(O,tmg[p][1]))],false)
      if order(Q)==E
        Base.delete!(L,p)
      end
    else
      k1=v-1
      k2=v
      gens=GrpAbFinGenElem[]
      Q=DiagonalGroup(Int[])
      while k1>=1
        multg=_1pluspk_1pluspk1(K, p, p^k1, p^k2, mR.fact_mod, prime_power, Int(cond.gen_one),Int(E))
        for i=1:length(multg)
          push!(gens, preimage(mp, ideal(O,multg[i])))
        end
        Q,mQ=quo(G,gens,false)
        if order(Q)!=E
          L[p]=k2
          break
        end
        k1-=1
        k2-=1
      end
      if order(Q)==E
        if haskey(tmg, p)
          push!(gens, mp\ideal(O,tmg[p][1]))
          Q,mQ=quo(G,gens,false)
          if order(Q)==E
            delete!(L,p)
          end
        else
          delete!(L,p)
        end
      end
    end
  end
  cond=ideal(O,1)
  for (p,vp) in L
    cond*=p^vp
  end
  #Infinite part of the modulus
  cond_inf=InfPlc[]
  if !isempty(inf_plc)
    S, ex, lo=carlos_units(O)
    for i=1:length(inf_plc)
      pl=inf_plc[i]
      j=1
      while true
        if !ispositive(ex(S[j]),pl)
          break
        end
        i+=1
      end
      delta=minimum(cond)*ex(S[j])
      el=1+delta
      con=abs_upper_bound(1/real(conjugates_arb(delta))[j], fmpz)
      el+=con*delta
      Q,mQ=quo(G,mp\ideal(O,el),false)
      if order(Q)!=e
        push!(cond_inf, pl)
      end
    end
  end

  return cond, cond_inf

end

###############################################################################
#
#  isconductor function
#
###############################################################################

doc"""
***
  isconductor(C::Hecke.ClassField, m::NfOrdIdl, inf_plc::Array{InfPlc,1}=InfPlc[]; check) -> NfOrdIdl, Array{InfPlc,1}

> Checks if m, inf_plc is the conductor of the abelian extension corresponding to C. If check is false, it assumes that the
> given modulus is a multiple of the conductor.

***
"""
function isconductor(C::Hecke.ClassField, m::NfOrdIdl, inf_plc::Array{InfPlc,1}=InfPlc[]; check::Bool=true)
  #
  #  First, we need to find the subgroup
  #

  mR = C.rayclassgroupmap
  mS = C.quotientmap
  G = codomain(mS)
  mp = inv(mS) * mR

  R = domain(mR)
  cond = mR.modulus_fin
  inf_plc2 = mR.modulus_inf
  O=parent(cond).order
  E=order(G)
  expo=Int(exponent(G))
  K=O.nf
  tmg=mR.tame_mult_grp
  wild=mR.wild_mult_grp

  if check
    mS=inv(mS)
    dom=domain(mS)
    M=zero_matrix(FlintZZ,ngens(dom), ngens(codomain(mS)))
    for i=1:ngens(dom)
      elem=mS(dom[i]).coeff
      for j=1:ngens(codomain(mS))
        M[i,j]=elem[1,j]
      end
    end
    S1=Hecke.GrpAbFinGenMap(domain(mS),codomain(mS),M)
    T,mT=Hecke.kernel(S1)

    Sgens=find_gens_sub(mR,mT)

    r,mr=ray_class_group(m,inf_plc, n_quo= expo)
    quot=GrpAbFinGenElem[mr\s for s in Sgens]
    s,ms=quo(r, quot, false)
    if order(s)!=E
      return false
    end
  end

  #
  #  Some of the factors of the modulus are unnecessary for order reasons:
  #

  L=factor(m)
  for (p,vp) in L
    if !haskey(mR.fact_mod, p) || vp>mR.fact_mod[p]
      return false
    end
    if !divisible(E,minimum(p))
      if gcd(E, norm(p)-1)==1
        return false
      elseif vp>1
        return false
      end
    elseif vp==1
      return false
    end
  end

  if isodd(E) && !isempty(inf_plc)
    return false
  end
  for pl in inf_plc
    if !(pl in inf_plc2)
      return false
    end
  end

  prime_power=Dict{NfOrdIdl, NfOrdIdl}()
  for (p,v) in L
    prime_power[p]=p^v
  end

  if !isempty(inf_plc2)
    totally_positive_generators(mR,Int(cond.gen_one))
  end

  #Finite part of the modulus
  for (p,v) in L
    if v==1
      Q,mQ=quo(G,[preimage(mp, ideal(O,tmg[p][1]))])
      if order(Q)==E
        return false
      end
    else
      multg=_1pluspk_1pluspk1(K, p, p^(v-1), p^v, mR.fact_mod, prime_power, Int(cond.gen_one),Int(E))
      gens=Array{GrpAbFinGenElem,1}(length(multg))
      for i=1:length(multg)
        gens[i]= preimage(mp, ideal(O,multg[i]))
      end
      Q,mQ=quo(G,gens, false)
      if order(Q)==E
        return false
      end
    end
  end

  #Infinite part of the modulus
  if !isempty(inf_plc2)
    S, ex, lo=carlos_units(O)
    for i=1:length(inf_plc)
      i=1
      while true
        if !ispositive(ex(S[i]),pl)
          break
        end
        i+=1
      end
      el=1+minimum(cond)*ex(S[i])
      while !ispositive(el, pl)
        el+=minimum(cond)*ex(S[i])
      end
      Q,mQ=quo(G,preimage(mp, ideal(O,el)), false)
      if order(Q)==e
        return false
      end
    end
  end

  return true

end

####################################################################################
#
#  Discriminant function
#
####################################################################################

function discriminant(C::ClassField)

  if isdefined(C,:conductor)
    m=C.conductor[1]
    inf_plc=C.conductor[2]
  else
    C.conductor=conductor(C)
    m=C.conductor[1]
    inf_plc=C.conductor[2]
  end
  @assert typeof(m)==NfOrdIdl

  mp = inv(C.quotientmap) * C.rayclassgroupmap
  R = domain(mp)
  n = order(R)

  mR = C.rayclassgroupmap

  cond=mR.modulus_fin
  inf_plc=mR.modulus_inf
  O=parent(cond).order
  E=order(R)
  expo=Int(exponent(R))
  K=O.nf
  a=Int(minimum(m))
  relative_disc=Dict{NfOrdIdl,Int}()
  abs_disc=Dict{fmpz,fmpz}()
  lp=mR.fact_mod
  if isempty(lp)
    C.relative_discriminant=relative_disc
  end
  tmg=mR.tame_mult_grp
  wldg=mR.wild_mult_grp
  prime_power=Dict{NfOrdIdl, NfOrdIdl}()
  for (p,v) in lp
    prime_power[p]=p^v
  end
  fact=factor(m)

  for (p,v) in fact
    if v==1
      ap=n
      if isprime(n)
        ap-=1
      else
        el=mp\ideal(O,tmg[p][1]) #The generator is totally positive, we modified it before
        q,mq=quo(R,[el])
        ap-= order(q)
      end
      qw=divexact(degree(O),prime_decomposition_type(O,Int(p.minimum))[1][2])*ap
      abs_disc[p.minimum] = qw
      relative_disc[p] = ap
    else
      ap=n*v
      if isprime(n)
        ap-=v
      else
        if length(lp) > 1
          i_without_p = ideal(O,1)
          for (p2,vp2) in lp
            (p != p2) && (i_without_p *= prime_power[p2])
          end
          alpha, beta = idempotents(prime_power[p],i_without_p)
        end
        s=v
        @hassert :QuadraticExt 1 s>=2
        els=GrpAbFinGenElem[]
        for k=2:v
          s=s-1
          pk=p^s
          pv=pk*p
          gens=_1pluspk_1pluspk1(K, p, pk, pv, lp, prime_power, a,Int(n))
          for i=1:length(gens)
            push!(els,mp\ideal(O,gens[i]))
          end
          ap-=order(quo(R,els)[1])
          @hassert :QuadraticExt 1 ap>0
        end
        if haskey(tmg, p)
          push!(els, mp\ideal(O,tmg[p][1]))
        end
        ap-=order(quo(R,els)[1])
        @hassert :QuadraticExt 1 ap>0
      end
      td=prime_decomposition_type(O,Int(p.minimum))
      np=fmpz(p.minimum)^(td[1][1]*length(td)*ap)
      abs_disc[p.minimum]=td[1][1]*length(td)*ap
      relative_disc[p]=ap
    end
  end
  C.relative_discriminant=relative_disc
  C.absolute_discriminant=abs_disc
  return C.relative_discriminant

end

###############################################################################
#
#  conductor function for abelian extension function
#
###############################################################################
#
#  For this function, we assume the base field to be normal over Q and the conductor of the extension we are considering to be invariant
#  The input must be a multiple of the minimum of the conductor, we don't check for consistency.
#

function _is_conductor_min_normal(C::Hecke.ClassField, a::Int)
# a is a positive integer in the modulus

  mr = C.rayclassgroupmap

  mp = inv(C.quotientmap) * mr
  R = domain(mp)

  lp=mr.fact_mod
  if isempty(lp)
    return true
  end
  O=order(first(keys(lp)))
  K=nf(O)
  tmg=mr.tame_mult_grp
  #first, tame part
  primes_done=fmpz[]
  for p in keys(tmg)
    if p.minimum in primes_done
      continue
    end
    push!(primes_done, p.minimum)
    el=mp\ideal(O,tmg[p][1])
    if iszero(el)
      return false
    end
  end

  #wild part
  if !isempty(mr.wild_mult_grp)
    prime_power=Dict{NfOrdIdl, NfOrdIdl}()
    for (p,v) in lp
      prime_power[p]=p^v
    end
    wldg=mr.wild_mult_grp
    primes_done=fmpz[]
    for p in keys(wldg)
      if p.minimum in primes_done
        continue
      end
      push!(primes_done, p.minimum)
      @assert lp[p]>=2
      k=lp[p]-1
      pk=p^k
      pv=prime_power[p]
      gens=_1pluspk_1pluspk1(K, p, pk, pv, lp, prime_power, a, Int(order(R)))
      iscond=false
      for i in 1:length(gens)
        if !iszero(mp\ideal(O,gens[i]))
          iscond=true
          break
        end
      end
      if !iscond
        return false
      end
    end
  end
  return true

end

#
#  Same function as above, but in the assumption that the map comes from a ray class group over QQ
#

function _is_conductor_minQQ(C::Hecke.ClassField, n::Int)

  mr = C.rayclassgroupmap
  mp = inv(C.quotientmap) * mr
  m = mr.modulus_fin
  mm = Int(minimum(m))
  lp = factor(m.minimum)

  O=order(m)
  K=nf(O)

  R=ResidueRing(FlintZZ, mm, cached=false)
  for (p,v) in lp.fac
    if isodd(p)
      if v==1
        x=_unit_grp_residue_field_mod_n(Int(p), n)[1]
        s=divexact(mm,Int(p))
        d,a,b=gcdx(s,Int(p))
        l=a*s*R(x)+p*b
        if iszero(mp\(ideal(O,Int(l.data))))
          return false
        end
      else
        s=divexact(mm,p^v)
        d,a,b=gcdx(s,p^v)
        x=a*s*R(1+p)^((p-1)*p^(v-2))+p^v*b
        if iszero(mp\(ideal(O,Int(x.data))))
          return false
        end
      end
    else
      if v==1
        return false
      end
      if v==2
        s=divexact(mm,4)
        d,a,b=gcdx(s,4)
        l=a*s*R(-1)+4*b
        if iszero(mp\(ideal(O,Int(l.data))))
          return false
        end
      else
        s=divexact(mm,2^v)
        d,a,b=gcdx(s, 2^v)
        l=a*s*R(5)^(2^(v-3))+2^v*b
        if iszero(mp\(ideal(O,Int(l.data))))
          return false
        end
      end
    end
  end
  return true

end

#######################################################################################
#
#  relative discriminant for abelian extension function
#
#######################################################################################

function discriminant_conductor(O::NfOrd, C::ClassField, a::Int, mr::MapRayClassGrp, bound::fmpz, n::Int)

  lp=mr.fact_mod
  abs_disc=factor(discriminant(O)^n).fac
  if isempty(lp)
    C.absolute_discriminant=abs_disc
    return true
  end
  K=nf(O)

  discr=fmpz(1)
  mp=inv(C.quotientmap) * C.rayclassgroupmap
  R=domain(mp)

  cyc_prime= isprime(n)==true

  #first, tamely ramified part
  tmg=mr.tame_mult_grp
  primes_done=fmpz[]
  for p in keys(tmg)
    if p.minimum in primes_done || haskey(mr.wild_mult_grp, p)
      continue
    end
    ap=n
    push!(primes_done, p.minimum)
    if cyc_prime
      ap-=1
    else
      el=mp\ideal(O,tmg[p][1]) #The generator is totally positive, we modified it before
      q,mq=quo(R,[el])
      ap-= order(q)
    end
    qw=divexact(degree(O),prime_decomposition_type(O,Int(p.minimum))[1][2])*ap
    discr*=fmpz(p.minimum)^qw
    if discr>bound
      @vprint :QuadraticExt 2 "too large\n"
      return false
    else
      if haskey(abs_disc, p.minimum)
        abs_disc[p.minimum]+=qw
      else
        abs_disc[p.minimum]=qw
      end
      #for q in keys(tmg)
      #  if minimum(q)==minimum(p)
      #    relative_disc[q]=ap
      #  end
      #end
    end
  end

  #now, wild ramification
  if !isempty(mr.wild_mult_grp)
    prime_power=Dict{NfOrdIdl, NfOrdIdl}()
    for (p,v) in lp
      prime_power[p]=p^v
    end
    wldg=mr.wild_mult_grp
    primes_done=fmpz[]
    for p in keys(wldg)
      if p.minimum in primes_done
        continue
      end
      push!(primes_done, p.minimum)
      ap=n*lp[p]
      if cyc_prime
        ap-=lp[p]
      else
        if length(lp) > 1
          i_without_p = ideal(O,1)
          for (p2,vp2) in lp
            (p != p2) && (i_without_p *= prime_power[p2])
          end

          alpha, beta = idempotents(prime_power[p],i_without_p)
        end
        s=lp[p]
        @hassert :QuadraticExt 1 s>=2
        els=GrpAbFinGenElem[]
        for k=2:lp[p]
          s=s-1
          pk=p^s
          pv=pk*p
          gens=_1pluspk_1pluspk1(K, p, pk, pv, lp, prime_power, a,n)
          for i=1:length(gens)
            push!(els,mp\ideal(O,gens[i]))
          end
          ap-=order(quo(R,els)[1])
          @hassert :QuadraticExt 1 ap>0
        end
        if haskey(tmg, p)
          push!(els, mp\ideal(O,tmg[p][1]))
        end
        ap-=order(quo(R,els)[1])
        @hassert :QuadraticExt 1 ap>0
      end
      td=divexact(degree(O),prime_decomposition_type(O,Int(p.minimum))[1][2])*ap
      np=p.minimum^td
      discr*=np
      if discr>bound
        @vprint :QuadraticExt 2 "too large\n"
        return false
      else
        if haskey(abs_disc, p.minimum)
          abs_disc[p.minimum]+=td
        else
          abs_disc[p.minimum]=td
        end
      #  for q in keys(tmg)
      #    if minimum(q)==minimum(p)
      #      relative_disc[q]=ap
      #    end
      #  end
      end
    end
  end
  #C.relative_discriminant=relative_disc
  C.absolute_discriminant=abs_disc
  return true

end

#same function but for ray class groups over QQ

function discriminant_conductorQQ(O::NfOrd, C::ClassField, m::Int, bound::fmpz, n::Int)

  discr=fmpz(1)
  mp = inv(C.quotientmap) * C.rayclassgroupmap
  G=domain(mp)

  cyc_prime= isprime(n)==true

  lp=factor(m).fac
  abs_disc=Dict{fmpz,Int}()

  R=ResidueRing(FlintZZ, m, cached=false)

  for (p,v) in lp
    if v==1
      ap=n
      if cyc_prime
        ap-=1
      else
        x=_unit_grp_residue_field_mod_n(Int(p),n)[1]
        s=divexact(m,Int(p))
        d,a,b=gcdx(s, Int(p))
        l=Int((R(x)*a*s+b*Int(p)).data)
        el=mp\ideal(O,l)
        q,mq=quo(G,[el])
        ap-= order(q)
      end
      discr*=p^ap
      if discr>bound
        @vprint :QuadraticExt 2 "too large\n"
        return false
      else
        abs_disc[p]=ap
      end
    else
      ap=n*v
      pow=Int(p)^Int(v)
      el=R(1)
      if cyc_prime
        ap-=v
      else
        if isodd(p)
          s=divexact(m,pow)
          d,a,b=gcdx(pow,s)
          s1=R(1+p)^(p-1)
          el=G[1]
          if v==2
            el=mp\ideal(O,Int((b*s*R(s1)+a*pow).data))
            ap-=order(quo(G,[el])[1])
          else
            for k=0:v-2
              el=mp\ideal(O,Int((b*s*R(s1)^(p^k)+a*pow).data))
              ap-=order(quo(G,[el])[1])
              @hassert :QuadraticExt 1 ap>0
            end
          end
          if gcd(n,p-1)==1
            ap-=order(quo(G,[mp\(ideal(O,fmpz((b*s*R(s1)+a*pow).data)))])[1])
          else
            x=_unit_grp_residue_field_mod_n(Int(p),n)[1]
            el1=mp\ideal(O,Int((R(x)*b*s+a*pow).data))
            ap-=order(quo(G,[mp\(ideal(O,Int((b*s*R(s1)+a*pow).data))), el1])[1])
          end
        else
          s=divexact(m,2^v)
          d,a,b=gcdx(2^v,s)
          el=0*G[1]
          for k=v-3:-1:0
            el=mp\ideal(O,Int((R(5)^(2^k)*b*s+a*2^v).data))
            ap-=order(quo(G,[el])[1])
          end
          el1=mp\ideal(O,Int((R(-1)*b*s+a*p^v).data))
          ap-=2*order(quo(G,[el, el1])[1])
        end
      end
      discr*=p^ap
      if discr>bound
        @vprint :QuadraticExt 2 "too large\n"
        return false
      else
        abs_disc[p]=ap
      end
    end
  end
  C.absolute_discriminant=abs_disc
  return true

end

##############################################################################
#
#  Is Abelian function
#
##############################################################################

doc"""
***
  function isabelian(f::Nemo.Generic.Poly, K::Nemo.AnticNumberField) -> Bool

 > Check if the extension generated by a root of the irreducible polynomial $f$ over a number field $K$ is abelian
 > The function is probabilistic.

***
"""

function isabelian(f::Nemo.PolyElem, K::Nemo.AnticNumberField)

  O=maximal_order(K)
  d=discriminant(f)
  N=num(norm(K(d)))
  n=degree(f)

  inf_plc=real_places(K)
  m=ideal(O,O(d))
  lp=collect(keys(factor(n)))
  M=zero_matrix(FlintZZ,0,0)
  Grps=Any[]
  R=AbelianGroup(fmpz[])
  for i=1:length(lp)
    T,mT=ray_class_group_p_part(Int(lp[i]),m,inf_plc)
    if valuation(order(T),lp[i])<valuation(n,lp[i])
      return false
    end
    push!(Grps, [T,mT])
  end
  for i=1:length(lp)
    R=direct_product(R,Grps[i][1])
  end
  function mR(J::NfOrdIdl)
    x=(Grps[1][2]\J).coeff
    for i=2:length(lp)
      hcat!(x,((Grps[i][2])\J).coeff)
    end
    return R(x)
  end

  S,mS=snf(R)
  M=rels(S)

  p=1
  Ox,x=PolynomialRing(O,"y", cached=false)
  f1=Ox([O(coeff(f,i)) for i=0:n])

  determinant=order(S)
  new_mat=M

  B=log(abs(discriminant(O)))*degree(f)+log(N)
  B=4*B+2.5*degree(f)*degree(O)+5
  B=B^2

  #
  # Adding small primes until they generate the norm group
  #

  while determinant > n
    p=next_prime(p)
    if p>B
      return false #Bach bound says that the norm group must be generated by primes $\leq B$
    end
    if !divisible(N,p)
      L=prime_decomposition(O,p)
      for i=1:length(L)
        F,mF=ResidueField(O,L[i][1])
        Fz,z= PolynomialRing(F,"z", cached=false)
        g=mF(f1)
        D=factor_shape(g)
        if length(D)>1
          return false
        end
        candidate=mR(((L[i][1])^first(keys(D))))
        new_mat=vcat(new_mat,(mS(candidate)).coeff)
        new_mat=hnf(new_mat)
        new_mat=sub(new_mat,1:ngens(S), 1:ngens(S))
        determinant=abs(det(new_mat))
      end
    end
  end
  if determinant==n
    return true
  else
    return false
  end

end

################################################################################
#
#  Norm group function
#
################################################################################

doc"""
***
  function norm_group(f::Nemo.PolyElem, mR::Hecke.MapRayClassGrp) -> Hecke.FinGenGrpAb, Hecke.FinGenGrpAbMap

 > Computes the subgroup of the Ray Class Group $R$ given by the norm of the extension generated by the roots of $f$

***
"""

function norm_group(f::Nemo.PolyElem, mR::Hecke.MapRayClassGrp)

  R=mR.header.domain
  O=mR.modulus_fin.parent.order
  K=O.nf
  d=discriminant(f)
  N=num(norm(K(d)))
  N1=fmpz(norm(mR.modulus_fin))
  n=degree(f)

  Q,mQ=quo(R,n, false)

  S,mS=snf(Q)
  M=rels(S)

  p=1
  Ox,x=PolynomialRing(O,"y", cached=false)
  f=Ox([O(coeff(f,i)) for i=0:n])

  determinant=abs(det(M))
  listprimes=typeof(R[1])[]
  new_mat=M

  #
  # Adding small primes until they generate the norm group
  #

  while determinant!= n
    p=next_prime(p)
    if !divisible(N,p) && !divisible(N1,p)
      L=prime_decomposition(O,p)
      for i=1:length(L)
        F,mF=ResidueField(O,L[i][1])
        Fz,z= PolynomialRing(F, "z", cached=false)
        g=mF(f)
        D=factor_shape(g)
        E=collect(keys(D))[1]
        candidate=mR\(((L[i][1]))^E)
        new_mat=vcat(new_mat,((mS*mQ)(candidate)).coeff)
        new_mat=hnf(new_mat)
        new_mat=sub(new_mat,1:ngens(S), 1:ngens(S))
        new_det=abs(det(new_mat))
        if determinant!=new_det
          push!(listprimes, candidate)
          determinant=new_det
        end
      end
    end
  end

  #
  # Computing the Hermite normal form of the subgroup
  #
  subgrp=[el for el in listprimes]
  for i=1:ngens(R)
    push!(subgrp, n*R[i])
  end
  return sub(R, subgrp, false)

end


text:
module RAVE

using Reexport

@reexport using Knet, MLToolkit, Distributions

import AutoGrad

include("data.jl")
export get_data
include("plotting.jl")
export write_vae, make_Z_plot

function flatstr(args::Dict)
    TR_LIST = [:n_epochs, :batch_size, :lr, :beta1, :beta2]
    IGNORE_LIST = vcat([:continue, :store, :gpu_id, :check_every, :eval_every,
                        :bnp, :inf, :is_deep, :obs], TR_LIST)
    args_str_flat = flatten_dict(args; exclude=IGNORE_LIST, equal_sym="-")
    args_str_flat_tr = flatten_dict(args; include=TR_LIST, equal_sym="-")
    # Method
    args_str_flat_method = "$(args[:bnp])-$(args[:inf])-$(args[:is_deep] ? "deep" : "linear")-$(args[:obs])"
    return args_str_flat, args_str_flat_method, args_str_flat_tr
end

export flatstr

# IBP
include("ibpvae.jl")
export IBPDeepEncoder, IBPLinearEncoder, BernoulliIBPLinearDecoder, GaussianIBPLinearDecoder
include("mfibpvae.jl")
include("sibpvae.jl")
include("rrsibpvae.jl")

function IBPVAE(x_dim, args; sig_X=FT(0.1))
    N = args[:tr_sz]
    prior = IBP(FT(args[:alpha]))
    h_dim = args[:h_dim]
    k_max = args[:k_max]
    obs_type = Symbol(args[:obs])
    inf = Symbol(args[:inf])

    ibp_encoder = if inf == :mf
        MeanFieldSBC(h_dim, k_max)
    elseif inf == :s || inf == :rrs
        StructuredSBC(prior, h_dim, k_max)
    else
        throw("[IBPVAE] Unkown inference method: $(args[:inference])")
    end

    if args[:is_deep]
        encoder = IBPDeepEncoder(k_max, x_dim, h_dim, ibp_encoder)
        obs_layer = obs_type == :gauss ? GaussianLogVarDense(h_dim, x_dim) : BernoulliLogitDense(h_dim, x_dim)
        decoder = Chain(LazyDense(k_max, h_dim; f=relu), obs_layer)
    else
        encoder = IBPLinearEncoder(k_max, x_dim, h_dim, ibp_encoder)
        decoder = obs_type == :gauss ? GaussianIBPLinearDecoder(sig_X, k_max, x_dim) : BernoulliIBPLinearDecoder(k_max, x_dim)
    end

    return if inf == :mf
        MFIBPVAE(prior, N, k_max, x_dim, h_dim, encoder, decoder)
    elseif inf == :s
        SIBPVAE(prior, N, k_max, x_dim, h_dim, encoder, decoder)
    elseif inf == :rrs
        RRSIBPVAE(prior, N, k_max, x_dim, h_dim, encoder, decoder, Param(Rho(1)))
    end
end

export IBPVAE

end # module


text:
const su = 1 #Sun's index
const ea = 4 #Earth's index
const mo = 5 #Moon's index

# vector of G*m values in au^3/day^2
GM1 = 0.491248045036476000E-10 # Mercury
GM2 = 0.724345233264412000E-09 # Venus
GM3 = 0.888769244512563400E-09 # Earth
GM4 = 0.954954869555077000E-10 # Mars
GM5 = 0.282534584083387000E-06 # Jupiter
GM6 = 0.845970607324503000E-07 # Saturn
GM7 = 0.129202482578296000E-07 # Uranus
GM8 = 0.152435734788511000E-07 # Neptune
GM9 = 0.217844105197418000E-11 # Pluto
GMS = 0.295912208285591100E-03 # Sun
GMM = 0.109318945074237400E-10 # Moon
GMB = 0.899701139019987100E-09 # Earth-Moon barycenter

const μ = [GMS, GM1, GM2, GM3, GMM, GM4, GM5, GM6, GM7, GM8, GM9,

1.4004765561723440E-13, # 1 Ceres
3.8547501878088100E-14, # 4 Vesta
3.1044481989387130E-14, # 2 Pallas
1.2358007872941250E-14, # 10 Hygiea
6.3432804736486020E-15, # 31 Euphrosyne
5.2561686784936620E-15, # 704 Interamnia
5.1981269794574980E-15, # 511 Davida
4.6783074183509050E-15, # 15 Eunomia
3.6175383171479370E-15, # 3 Juno
3.4115868261938120E-15, # 16 Psyche
3.1806592826525410E-15, # 65 Cybele
2.5771141273110470E-15, # 88 Thisbe
2.5310917260150680E-15, # 48 Doris
2.4767881012558670E-15, # 52 Europa
2.2955593906374620E-15, # 451 Patientia
2.1992951735740730E-15, # 87 Sylvia
2.1364344425714070E-15, # 7 Iris
2.1124383605999520E-15, # 423 Diotima
1.9758423651245200E-15, # 29 Amphitrite
1.8939016675253820E-15, # 24 Themis
1.7970048945074460E-15, # 13 Egeria
1.7558991833247000E-15, # 790 Pretoria
1.7459557262705000E-15, # 372 Palma
1.6717209917006440E-15, # 107 Camilla
1.5850986571596890E-15, # 354 Eleonora
1.5465676956243250E-15, # 96 Aegle
1.5079333711965190E-15, # 386 Siegena
1.4975196825567010E-15, # 39 Laetitia
1.3886265898561990E-15, # 324 Bamberga
1.3315362554599750E-15, # 11 Parthenope
1.2792300000000000E-15, # 94 Aurora
1.2026244434834600E-15, # 22 Kalliope
1.1889849499520080E-15, # 120 Lachesis
1.1355858944839220E-15, # 185 Eunike
1.1006456795750680E-15, # 14 Irene
1.0975631032822510E-15, # 536 Merapi
1.0778410042407300E-15, # 9 Metis
1.0684599027046930E-15, # 165 Loreley
1.0356448401311940E-15, # 19 Fortuna
9.9366295459092490E-16, # 130 Elektra
9.6501295105487520E-16, # 128 Nemesis
9.4193644635953100E-16, # 187 Lamberta
9.3242237621988690E-16, # 41 Daphne
9.3159485940656200E-16, # 532 Herculina
9.2540854530185370E-16, # 85 Io
9.0708048441145050E-16, # 444 Gyptis
8.8950672849270450E-16, # 702 Alauda
8.5612605995538920E-16, # 137 Meliboea
8.4594307289596830E-16, # 471 Papagena
8.4256780185679340E-16, # 45 Eugenia
8.3124192126733720E-16, # 6 Hebe
8.3122000000000010E-16, # 76 Freia
8.2983887671636940E-16, # 154 Bertha
8.2144999999999990E-16, # 409 Aspasia
7.5582329262288920E-16, # 145 Adeona
7.2405478852581390E-16, # 308 Polyxo
7.0079069220413430E-16, # 121 Hermione
7.0078739271302920E-16, # 349 Dembowska
6.9951633549830870E-16, # 144 Vibilia
6.9079712474674250E-16, # 216 Kleopatra
6.3394427275876510E-16, # 59 Elpis
6.2808585539363830E-16, # 259 Aletheia
6.2674000000000000E-16, # 566 Stereoskopia
6.2157460662366880E-16, # 747 Winchester
5.8942565297069080E-16, # 8 Flora
5.6477307179764760E-16, # 93 Minerva
5.6241736501924590E-16, # 54 Alexandra
5.5764804768085390E-16, # 405 Thia
5.5435235615988890E-16, # 47 Aglaja
5.4846724291131710E-16, # 489 Comacina
5.3973999999999980E-16, # 106 Dione
5.3689097042583350E-16, # 344 Desiderata
5.3470999999999980E-16, # 168 Sibylla
5.2966693580739110E-16, # 238 Hypatia
5.2669011092543480E-16, # 268 Adorea
5.1446100208767350E-16, # 69 Hesperia
5.0204242212925650E-16, # 712 Boliviana
4.8352000000000000E-16, # 420 Bertholda
4.8122396678018730E-16, # 104 Klymene
4.7637535316542470E-16, # 690 Wratislavia
4.6542452473979750E-16, # 129 Antigone
4.6386880160793900E-16, # 618 Elfriede
4.5585992488069250E-16, # 375 Ursula
4.5559595663772370E-16, # 814 Tauris
4.5336468487968120E-16, # 150 Nuwa
4.5014136561179280E-16, # 196 Philomela
4.4713680178417890E-16, # 117 Lomia
4.2242882143774450E-16, # 139 Juewa
4.1832602031724900E-16, # 164 Eva
4.1456571725352780E-16, # 225 Henrietta
4.0681729177164310E-16, # 147 Protogeneia
4.0369435176860730E-16, # 92 Undina
3.9416000000000000E-16, # 146 Lucina
3.9124285470831010E-16, # 173 Ino
3.8880038985457880E-16, # 27 Euterpe
3.8674996577516790E-16, # 212 Medea
3.8350891708780020E-16, # 596 Scheila
3.7864750330894840E-16, # 895 Helio
3.7661406016708330E-16, # 141 Lumen
3.7487362845520320E-16, # 5 Astraea
3.7154667975341450E-16, # 105 Artemis
3.6992883127021260E-16, # 56 Melete
3.6806019206396510E-16, # 57 Mnemosyne
3.6782000000000000E-16, # 419 Aurelia
3.6611682430617230E-16, # 127 Johanna
3.6414973977832910E-16, # 490 Veritas
3.6186546970297340E-16, # 410 Chloris
3.5953000000000000E-16, # 654 Zelinda
3.5934934807617610E-16, # 381 Myrrha
3.5073744512956140E-16, # 74 Galatea
3.4333444128591690E-16, # 388 Charybdis
3.4310265912379690E-16, # 68 Leto
3.4176985017472390E-16, # 505 Cava
3.4106769797534870E-16, # 508 Princetonia
3.4020311574394290E-16, # 89 Julia
3.3677605743530070E-16, # 360 Carlova
3.3629754927610390E-16, # 909 Ulla
3.3620465420887160E-16, # 134 Sophrosyne
3.3407005297045190E-16, # 481 Emita
3.2728000000000000E-16, # 46 Hestia
3.2624933794345620E-16, # 334 Chicago
3.2284000000000000E-16, # 469 Argentina
3.1316732532224060E-16, # 140 Siwa
3.0567119446636530E-16, # 776 Berbericia
3.0465115564904100E-16, # 211 Isolda
3.0054836259460050E-16, # 241 Germania
2.9524140803084220E-16, # 40 Harmonia
2.9445412915212860E-16, # 34 Circe
2.9376520795314820E-16, # 514 Armida
2.9262727442945280E-16, # 28 Bellona
2.9175003506482500E-16, # 705 Erminia
2.9127000000000000E-16, # 328 Gudrun
2.8990971962115820E-16, # 175 Andromache
2.8848684039510110E-16, # 303 Josephina
2.8824066183805000E-16, # 159 Aemilia
2.8716404826701960E-16, # 1093 Freda
2.7688888401578460E-16, # 70 Panopaea
2.7653166434743810E-16, # 42 Isis
2.7485668803401500E-16, # 554 Peraga
2.7230587289059820E-16, # 194 Prokne
2.7166197083932590E-16, # 95 Arethusa
2.7097270724162130E-16, # 247 Eukrate
2.7071416736527810E-16, # 466 Tisiphone
2.6816131961351820E-16, # 356 Liguria
2.6336000000000000E-16, # 156 Xanthippe
2.5931073098872930E-16, # 209 Dido
2.5705491133531450E-16, # 51 Nemausa
2.5294428720409990E-16, # 21 Lutetia
2.5206144925338950E-16, # 772 Tanete
2.5119412656578330E-16, # 192 Nausikaa
2.4428317417320690E-16, # 98 Ianthe
2.4404616777010060E-16, # 91 Aegina
2.4193160564640660E-16, # 476 Hedwig
2.3729075379342010E-16, # 162 Laurentia
2.3522561732418410E-16, # 35 Leukothea
2.2979758127971450E-16, # 426 Hippo
2.2861593339578100E-16, # 455 Bruchsalia
2.2671033414599050E-16, # 804 Hispania
2.2481550488436830E-16, # 276 Adelheid
2.2211133615828450E-16, # 595 Polyxena
2.2099166077176560E-16, # 346 Hermentaria
2.1856205771130560E-16, # 37 Fides
2.1581556817636010E-16, # 788 Hohensteina
2.1523995570228910E-16, # 86 Semele
2.1501467433612250E-16, # 602 Marianna
2.0818012830905660E-16, # 762 Pulcova
2.0815063964697380E-16, # 17 Thetis
2.0528101779869580E-16, # 283 Emma
2.0089277366511320E-16, # 18 Melpomene
1.9735313911048410E-16, # 506 Marion
1.9725153985693760E-16, # 769 Tatjana
1.9715919666254550E-16, # 233 Asterope
1.9688550185992420E-16, # 250 Bettina
1.9576160906423330E-16, # 773 Irmintraud
1.9484986115571250E-16, # 545 Messalina
1.9317757851829200E-16, # 12 Victoria
1.9151562798850780E-16, # 488 Kreusa
1.8953317604197830E-16, # 23 Thalia
1.8915247466562090E-16, # 326 Tamara
1.8654000000000000E-16, # 275 Sapientia
1.8492472080300540E-16, # 203 Pompeja
1.8413301875237820E-16, # 266 Aline
1.8168406201733030E-16, # 221 Eos
1.8124493964458610E-16, # 521 Brixia
1.7989374811476500E-16, # 751 Faina
1.7619560261002570E-16, # 357 Ninina
1.7608987155151350E-16, # 230 Athamantis
1.7554616782526140E-16, # 200 Dynamene
1.7252555236382890E-16, # 191 Kolga
1.7050000000000000E-16, # 114 Kassandra
1.6995669526331390E-16, # 635 Vundtia
1.6970600184097090E-16, # 36 Atalante
1.6796766833548250E-16, # 181 Eucharis
1.6713475520093610E-16, # 176 Iduna
1.6643778725882150E-16, # 626 Notburga
1.6595259516345300E-16, # 148 Gallia
1.6373439522610840E-16, # 26 Proserpina
1.6333263911175180E-16, # 50 Virginia
1.6308230966598090E-16, # 784 Pickeringia
1.5923470884492270E-16, # 675 Ludmilla
1.5830198804699920E-16, # 786 Bredichina
1.5678504140474830E-16, # 407 Arachne
1.5586000000000000E-16, # 393 Lampetia
1.5456008538659910E-16, # 709 Fringilla
1.5352978885566200E-16, # 980 Anacostia
1.5109884909293550E-16, # 171 Ophelia
1.4938668030133560E-16, # 694 Ekard
1.4869856293449610E-16, # 416 Vaticana
1.4858970838252890E-16, # 491 Carina
1.4857405347656310E-16, # 780 Armenia
1.4820190164375290E-16, # 30 Urania
1.4756769258889110E-16, # 335 Roberta
1.4655734975631910E-16, # 412 Elisabetha
1.4521678954842190E-16, # 404 Arsinoe
1.4383748004467890E-16, # 674 Rachele
1.4268157848333010E-16, # 713 Luscinia
1.4244927463509560E-16, # 71 Niobe
1.4144228920898080E-16, # 377 Campania
1.4100078405763090E-16, # 350 Ornamenta
1.4076985722105040E-16, # 110 Lydia
1.3906039685198460E-16, # 223 Rosa
1.3720617596137790E-16, # 373 Melusina
1.3659771964686970E-16, # 100 Hekate
1.3558171306734880E-16, # 449 Hamburga
1.3232859647467680E-16, # 38 Leda
1.3145378148802100E-16, # 210 Isabella
1.2987929239170220E-16, # 498 Tokio
1.2812662256605980E-16, # 102 Miriam
1.2583765926450050E-16, # 1015 Christa
1.2573126556318860E-16, # 84 Klio
1.2351963628284910E-16, # 90 Antiope
1.2317521169884130E-16, # 345 Tercidina
1.2196487517410750E-16, # 791 Ani
1.2178097995831790E-16, # 358 Apollonia
1.2049549488578180E-16, # 663 Gerlinde
1.2048459545465020E-16, # 227 Philosophia
1.1995850162334400E-16, # 32 Pomona
1.1892814439317690E-16, # 740 Cantabia
1.1789669849357310E-16, # 206 Hersilia
1.1614439541131080E-16, # 80 Sappho
1.1584591337407960E-16, # 313 Chaldaea
1.1563954713645830E-16, # 366 Vincentina
1.1534210849315850E-16, # 445 Edna
1.1363293901133810E-16, # 236 Honoria
1.1260861147488830E-16, # 143 Adria
1.1187183314500760E-16, # 503 Evelyn
1.1177734423974850E-16, # 696 Leonora
1.1152801330348170E-16, # 1467 Mashona
1.1048622528732650E-16, # 240 Vanadis
1.1031851822968280E-16, # 385 Ilmatar
1.0968348906260350E-16, # 83 Beatrix
1.0890481919600570E-16, # 62 Erato
1.0886266508800200E-16, # 517 Edith
1.0826185861581930E-16, # 109 Felicitas
1.0436103612506560E-16, # 683 Lanzia
1.0354972501527590E-16, # 849 Ara
1.0314956358376310E-16, # 97 Klotho
1.0248487058744050E-16, # 160 Una
1.0223675545561340E-16, # 81 Terpsichore
1.0147296895716560E-16, # 201 Penelope
1.0044659839630940E-16, # 387 Aquitania
1.0011180168586460E-16, # 103 Hera
9.9581620100255190E-17, # 680 Genoveva
9.5491282158876470E-17, # 1107 Lictoria
9.5223891872172390E-17, # 213 Lilaea
9.5156050418462080E-17, # 135 Hertha
9.4537052258201020E-17, # 691 Lehigh
9.1998074776309120E-17, # 20 Massalia
9.1514999225958940E-17, # 205 Martha
9.0942944157907600E-17, # 667 Denise
8.8772708265623380E-17, # 124 Alkeste
8.6580000000000000E-17, # 163 Erigone
8.5978021881831310E-17, # 568 Cheruskia
8.5976202637173240E-17, # 735 Marghanna
8.4811739114660020E-17, # 58 Concordia
8.4019062534638870E-17, # 78 Diana
8.3548198507378080E-17, # 195 Eurykleia
8.3518243314079400E-17, # 79 Eurynome
8.2763495213043580E-17, # 322 Phaeo
8.0870792881032830E-17, # 362 Havnia
8.0635184290055780E-17, # 464 Megaira
7.9950510449165410E-17, # 72 Feronia
7.7695552096218000E-17, # 365 Corduba
7.6808847746999990E-17, # 454 Mathesis
7.6171409872842650E-17, # 1021 Flammario
7.5494816293144010E-17, # 49 Pales
7.3799853443752290E-17, # 424 Gratia
7.3526620138415910E-17, # 99 Dike
7.2753938533407130E-17, # 43 Ariadne
7.2538616434661140E-17, # 535 Montague
7.2398415223662110E-17, # 25 Phocaea
7.1382433975293480E-17, # 739 Mandeville
6.9841945302118660E-17, # 363 Padua
6.9600000000000000E-17, # 516 Amherstia
6.9118774732740100E-17, # 599 Luisa
6.8832710966054010E-17, # 304 Olga
6.6012607669307700E-17, # 82 Alkmene
6.5399257724402350E-17, # 465 Alekto
6.4677672132374610E-17, # 1036 Ganymed
6.4260399533186430E-17, # 604 Tekmessa
6.3407783495216970E-17, # 431 Nephele
6.2485687084637610E-17, # 1171 Rusthawelia
6.2392433107751650E-17, # 53 Kalypso
6.2251018387380330E-17, # 389 Industria
6.1530641540755240E-17, # 598 Octavia
6.1123219726683140E-17, # 569 Misa
5.8167921005732780E-17, # 760 Massinga
5.7960397015532350E-17, # 112 Iphigenia
5.7345815663368580E-17, # 369 Aeria
5.6430000000000000E-17, # 336 Lacadiera
5.6404017439762430E-17, # 63 Ausonia
5.5534656287469120E-17, # 442 Eichsfeldia
5.5258241903855270E-17, # 115 Thyra
5.4001214087434230E-17, # 415 Palatia
5.1407683428632200E-17, # 338 Budrosa
5.1143093942960570E-17, # 329 Svea
5.0911367830144640E-17, # 60 Echo
5.0188993452074920E-17, # 593 Titania
4.9982048125892480E-17, # 752 Sulamitis
4.9312955095007290E-17, # 77 Frigga
4.8033737856096100E-17, # 287 Nephthys
4.7743547738799830E-17, # 778 Theobalda
4.6886407201292200E-17, # 44 Nysa
4.5372886294250990E-17, # 177 Irma
4.3573746250771270E-17, # 75 Eurydike
4.3177157644920530E-17, # 172 Baucis
4.1172682683995420E-17, # 914 Palisana
3.9520590343452720E-17, # 224 Oceana
3.6737960799234180E-17, # 485 Genua
3.4515074511865960E-17, # 337 Devosa
3.3519192811280560E-17, # 111 Ate
2.8123577458657680E-17, # 547 Praxedis
2.7000759625981350E-17, # 118 Peitho
2.5580213924578190E-17, # 113 Amalthea
2.1949122177817370E-17, # 347 Pariana
2.1760595614159460E-17, # 584 Semiramis
2.0538576829083310E-17, # 198 Ampella
2.0038393509521630E-17, # 591 Irmgard
1.8780122958989280E-17, # 432 Pythia
1.4497976797693290E-17, # 623 Chimaera
1.3196141221701560E-17, # 132 Aethra
1.1716991540679400E-17, # 585 Bilkis
9.9000011897959020E-19  # 433 Eros
]

const de430_343ast_ids =
    [1, 4, 2, 10, 31, 704, 511, 15, 3, 16,
    65, 88, 48, 52, 451, 87, 7, 423, 29, 24,
    13, 790, 372, 107, 354, 96, 386, 39, 324, 11,
    94, 22, 120, 185, 14, 536, 9, 165, 19, 130,
    128, 187, 41, 532, 85, 444, 702, 137, 471, 45,
    6, 76, 154, 409, 145, 308, 121, 349, 144, 216,
    59, 259, 566, 747, 8, 93, 54, 405, 47, 489,
    106, 344, 168, 238, 268, 69, 712, 420, 104, 690,
    129, 618, 375, 814, 150, 196, 117, 139, 164, 225,
    147, 92, 146, 173, 27, 212, 596, 895, 141, 5,
    105, 56, 57, 419, 127, 490, 410, 654, 381, 74,
    388, 68, 505, 508, 89, 360, 909, 134, 481, 46,
    334, 469, 140, 776, 211, 241, 40, 34, 514, 28,
    705, 328, 175, 303, 159, 1093, 70, 42, 554, 194,
    95, 247, 466, 356, 156, 209, 51, 21, 772, 192,
    98, 91, 476, 162, 35, 426, 455, 804, 276, 595,
    346, 37, 788, 86, 602, 762, 17, 283, 18, 506,
    769, 233, 250, 773, 545, 12, 488, 23, 326, 275,
    203, 266, 221, 521, 751, 357, 230, 200, 191, 114,
    635, 36, 181, 176, 626, 148, 26, 50, 784, 675,
    786, 407, 393, 709, 980, 171, 694, 416, 491, 780,
    30, 335, 412, 404, 674, 713, 71, 377, 350, 110,
    223, 373, 100, 449, 38, 210, 498, 102, 1015, 84,
    90, 345, 791, 358, 663, 227, 32, 740, 206, 80,
    313, 366, 445, 236, 143, 503, 696, 1467, 240, 385,
    83, 62, 517, 109, 683, 849, 97, 160, 81, 201,
    387, 103, 680, 1107, 213, 135, 691, 20, 205, 667,
    124, 163, 568, 735, 58, 78, 195, 79, 322, 362,
    464, 72, 365, 454, 1021, 49, 424, 99, 43, 535,
    25, 739, 363, 516, 599, 304, 82, 465, 1036, 604,
    431, 1171, 53, 389, 598, 569, 760, 112, 369, 336,
    63, 442, 115, 415, 338, 329, 60, 593, 752, 77,
    287, 778, 44, 177, 75, 172, 914, 224, 485, 337,
    111, 547, 118, 113, 347, 584, 198, 591, 432, 623,
    132, 585, 433
]

# Matrix of J2 interactions included in DE430 ephemeris, according to Folkner et al., 2014
const UJ_interaction = fill(false, length(μ), length(μ))
UJ_interaction[2:11, su] .= true # per Folkner et al. (2014), Section III, 1st paragraph: Sun's J2 only interacts with the Moon and planets
UJ_interaction[union(1:ea-1,ea+1:7), ea] .= true # Earth's grav potential interacts with Sun, Mercury, Venus, Moon, Mars and Jupiter
UJ_interaction[union(1:mo-1,mo+1:7), mo] .= true # Moon's grav potential interacts with Sun, Mercury, Venus, Earth, Mars and Jupiter

const au = 1.495978707E8 # astronomical unit value in km
const yr = 365.25 # days in a Julian year
const daysec = 86_400 # number of seconds in a day
const clightkms = 2.99792458E5 # speed of light, km/sec
const c_au_per_day = daysec*(clightkms/au) # speed of light in au per day
const c_au_per_sec = clightkms/au # speed of light in au per sec
const c_cm_per_sec = 100_000*clightkms # speed of light in cm per sec

const sundofs = nbodyind(length(μ), su)
const earthdofs = nbodyind(length(μ), ea)

const J2000 = 2451545.0

const R_sun = 696000.0/au # Solar radius in au, value taken from DE430 docs

const α_p_sun = 286.13 # Sun's rotation pole right ascension (degrees)
const δ_p_sun = 63.87 # Sun's rotation pole declination (degrees)

const RSUN  = 6.9600000000000000E+05
const J2SUN = 2.1106088532726840E-07
const JS = [J2SUN*(RSUN/au)^2]

const RE     =  6.378136300E+03
const J2E    =  1.082625450E-03
const J3E    = -2.532410000E-06
const J4E    = -1.619898000E-06
const J5E    = -2.277345000E-07
const J2EDOT = -2.600000000E-11
const JE = [J2E*(RE/au)^2, J3E*(RE/au)^3, J4E*(RE/au)^4, J5E*(RE/au)^5]

const RM = 1.7380000000000000E+03
const J2M  =  2.0321568464952570E-04
const J3M  =  8.4597026974594570E-06
const J4M  = -9.7044138365700000E-06
const J5M  =  7.4221608384052890E-07
const J6M  = -1.3767531350969900E-05
const JM = [J2M*(RM/au)^2, J3M*(RM/au)^3, J4M*(RM/au)^4, J5M*(RM/au)^5, J6M*(RM/au)^6]

const JSEM = zeros(5, 6)
JSEM[su,2:2] = JS
JSEM[ea,2:5] = JE
JSEM[mo,2:6] = JM

const n1SEM = [2, 0, 0, 5, 6]
const n2M = 6

const C22M =  2.2382740590560020E-05

const C31M =  2.8480741195592860E-05
const S31M =  5.8915551555318640E-06
const C32M =  4.8449420619770600E-06
const S32M =  1.6844743962783900E-06
const C33M =  1.6756178134114570E-06
const S33M = -2.4742714379805760E-07
const C3M = [C31M, C32M, C33M]*(RM/au)^3
const S3M = [S31M, S32M, S33M]*(RM/au)^3

const C41M = -5.7048697319733210E-06
const S41M =  1.5789202789245720E-06
const C42M = -1.5912271792977430E-06
const S42M = -1.5153915796731720E-06
const C43M = -8.0678881596778210E-08
const S43M = -8.0349266627431070E-07
const C44M = -1.2692158612216040E-07
const S44M =  8.2964257754075220E-08
const C4M = [C41M, C42M, C43M, C44M]*(RM/au)^4
const S4M = [S41M, S42M, S43M, S44M]*(RM/au)^4

const C51M = -8.6629769308983560E-07
const S51M = -3.5272289393243820E-06
const C52M =  7.1199537967353330E-07
const S52M =  1.7107886673430380E-07
const C53M =  1.5399750424904520E-08
const S53M =  2.8736257616334340E-07
const C54M =  2.1444704319218450E-08
const S54M =  5.2652110720146800E-10
const C55M =  7.6596153884006140E-09
const S55M = -6.7824035473995330E-09
const C5M = [C51M, C52M, C53M, C54M, C55M]*(RM/au)^5
const S5M = [S51M, S52M, S53M, S54M, S55M]*(RM/au)^5

const C61M =  1.2024363601545920E-06
const S61M = -2.0453507141252220E-06
const C62M = -5.4703897324156850E-07
const S62M = -2.6966834353574270E-07
const C63M = -6.8785612757292010E-08
const S63M = -7.1063745295915780E-08
const C64M =  1.2915580402925160E-09
const S64M = -1.5361616966632300E-08
const C65M =  1.1737698784460500E-09
const S65M = -8.3465073195142520E-09
const C66M = -1.0913395178881540E-09
const S66M =  1.6844213702632920E-09
const C6M = [C61M, C62M, C63M, C64M, C65M, C66M]*(RM/au)^6
const S6M = [S61M, S62M, S63M, S64M, S65M, S66M]*(RM/au)^6

const CM = zeros(6, 6)
CM[3,1:3] = C3M
CM[4,1:4] = C4M
CM[5,1:5] = C5M
CM[6,1:6] = C6M

const SM = zeros(6, 6)
SM[3,1:3] = S3M
SM[4,1:4] = S4M
SM[5,1:5] = S5M
SM[6,1:6] = S6M

const R_moon = RM/au # lunar radius in au, value taken from DE430 docs
const β_L = 6.3102131934887270E-04 # Lunar moment parameter, β_L = (C_T-A_T)/B_T
const γ_L = 2.2773171480091860E-04 # Lunar moment parameter, γ_L = (B_T-A_T)/C_T
const k_2M = 0.024059 # potential Love number
const τ_M = 9.5830547273306690E-02 # time-lag for the lunar solid-body tide (days)
const α_c = 0.0007 # Ratio of polar moment of inertia of core to mean total polar moment of inertia
const f_c = 2.4623904789198150E-04 #oblateness of core
const k_ν_div_C_T = 1.6365616533709530E-08 # Friction coefficient between core and mantle, radian/day

const J2_M_und = 2.0321568464952570E-04 # undistorted lunar 2nd zonal harmonic coefficient

#diagonal elements of undistorted lunar mantle moment of inertia
const A_T = (2(1-β_L*γ_L)/(2β_L-γ_L+β_L*γ_L))*J2_M_und
const B_T = (2(1+γ_L)/(2β_L-γ_L+β_L*γ_L))*J2_M_und
const C_T = (2(1+β_L)/(2β_L-γ_L+β_L*γ_L))*J2_M_und

const k_ν = k_ν_div_C_T*(C_T*μ[mo]*(R_moon^2))

#diagonal elements of lunar core moment of inertia
const A_c = α_c*C_T*(1-f_c)
const B_c = α_c*C_T*(1-f_c)
const C_c = α_c*C_T
const C_c_m_A_c = (C_c - A_c)*(μ[mo]*R_moon^2)

#lunar undistorted total moment of inertia
const ITM_und = diagm([A_T, B_T, C_T]*μ[mo]*R_moon^2)

#lunar core moment of inertia
const I_c = diagm([A_c, B_c, C_c]*μ[mo]*R_moon^2)

const ld = 384402.0 # Lunar distance (km)
# 384399.014 km  mean semi-major axis value was retrieved from:
# J. G. Williams, D. H. Boggs, and W. M. Folkner (2013). “DE430 Lunar Orbit, Physical Librations and Surface Coordinates,” JPL IOM 335-JW,DB,WF-20130722-016 (internal document)
const n_moon = sqrt((μ[ea]+μ[mo])/((384399.014/au)^3)) # lunar mean motion (rad/day)

const k_20E = 0.335 # Potential Love number for long-period deformation
const k_21E = 0.32 # Potential Love number for diurnal deformation
const k_22E = 0.32 # Potential Love number for semi-diurnal deformation

const τ_0p = 0.0640 # Orbital time-lag for long-period deformation, days
const τ_1p = -0.044 # Orbital time-lag for diurnal deformation, days --> DE430 report states -0.0044, but tech comments and Williams et al. (2006) state -0.044
const τ_2p = -0.1000 # Orbital time-lag for semi-diurnal deformation, days

const τ_0 = 0.0 # Rotational time-lag for long-period deformation, days
const τ_1 = 7.3632190228041890E-03 # Rotational time-lag for diurnal deformation, days
const τ_2 = 2.5352978633388720E-03 # Rotational time-lag for semi-diurnal deformation, days

# standard value of nominal mean angular velocity of Earth (rad/day), ESAA 2014 Sec 7.4.3.3 p. 296: 7.2921151467e-5 rad/second
const ω_E = daysec*7.2921151467e-5 # (7.2921151467e-5 rad/sec)*daysec -> rad/day

# Earth/Moon mass ratio
const EMRAT = 8.1300569074190620E+01


text:
################################################################################
# airfoilprep.py WRAPPER CLASS
################################################################################
"""
  `Polar(Re, alpha, cl, cd, cm, x, y)`

Defines section lift, drag, and pitching moment coefficients as a
function of angle of attack at a particular Reynolds number. This object acts
as an interface for the Python object `Polar` from `airfoilprep.py`, with
the following properties:

  * `self.PyPolar.Re`     : (float) Reynolds number.
  * `self.PyPolar.alpha`  : (array) Angles of attacks.
  * `self.PyPolar.cl`     : (array) Lift coefficient at each AOA.
  * `self.PyPolar.cd`     : (array) Drag coefficient at each AOA.
  * `self.PyPolar.cm`     : (array) Moment coefficient at each AOA.

  # Arguments
  * Re::Int64               : Reynolds number of this polar
  * alpha::Array{Float64,1} : Angles of attack (deg)
  * cl::Array{Float64,1}    : Lift coefficient at each AOA.
  * cd::Array{Float64,1}    : Drag coefficient at each AOA.
  * cm::Array{Float64,1}    : Moment coefficient at each AOA.
  * x::Array{Float64,1}     : x-coordinates of airfoil geometry
  * y::Array{Float64,1}     : y-coordinates of airfoil geometry

NOTE: `alpha` input and output is always in degrees.
**NOTE 2: Airfoil points x,y must go from trailing edge around the top, then the
bottom and end back at the trailing edge.**
"""
mutable struct Polar
  # Initialization variables (USER INPUT)
  init_Re::Float64                        # Reynolds number of this polar
  init_alpha::Array{Float64,1}            # Angles of attack (deg)
  init_cl::Array{Float64,1}               # Lift coefficient at each AOA.
  init_cd::Array{Float64,1}               # Drag coefficient at each AOA.
  init_cm::Array{Float64,1}               # Moment coefficient at each AOA.
  Ma::Float64                             # Mach number
  npanels::Int                            # Number of panels
  ncrit::Float64                          # XFOIL's turbulence transition parameter
  x::Array{Float64,1}                     # x-coordinates of airfoil geometry
  y::Array{Float64,1}                     # y-coordinates of airfoil geometry
  xsepup::Array{Float64,1}                # Separation point over upper surface
  xseplo::Array{Float64,1}                # Separation point over lower surface

  # Internal variables
  pyPolar::PyCall.PyObject

  Polar(init_Re, init_alpha, init_cl, init_cd, init_cm;
            Ma=0.0, npanels=0, ncrit=0,
            x=Float64[], y=Float64[],
            xsepup=Float64[], xseplo=Float64[],
            pyPolar=prepy.Polar(init_Re, init_alpha, init_cl, init_cd, init_cm)
        ) = new(
            init_Re, init_alpha, init_cl, init_cd, init_cm,
            Ma, npanels, ncrit,
            x, y,
            xsepup, xseplo,
            pyPolar)
end

"Returns Re of this Polar"
function get_Re(self::Polar)
    return self.pyPolar.Re
end
"Returns Ma of this Polar"
function get_Ma(self::Polar)
    return self.Ma
end
"Returns number of panels of this Polar"
function get_npanels(self::Polar)
    return self.npanels
end
"Returns ncrit of this Polar"
function get_ncrit(self::Polar)
    return self.ncrit
end
"Returns (alpha, cl) points of this Polar"
function get_cl(self::Polar)
    return (self.pyPolar.alpha, self.pyPolar.cl)
end
"Returns (alpha, cd) points of this Polar"
function get_cd(self::Polar)
    return (self.pyPolar.alpha, self.pyPolar.cd)
end
"Returns (alpha, cm) points of this Polar"
function get_cm(self::Polar)
    return (self.pyPolar.alpha, self.pyPolar.cm)
end
"Returns (alpha, xsepup) points of this Polar"
function get_xsepup(self::Polar)
    return (self.init_alpha, self.xsepup)
end
"Returns (alpha, xseplo) points of this Polar"
function get_xseplo(self::Polar)
    return (self.init_alpha, self.xseplo)
end
"Returns (x, y) points of the airfoil geometry of this Polar"
function get_geometry(self::Polar)
    return (self.x, self.y)
end
"Returns a dummy polar object"
function dummy_polar()
    return Polar(-1, Float64[], Float64[], Float64[], Float64[])
end
"Reads a polar as downloaded from Airfoiltools.com"
function read_polar(file_name::String; path::String="", optargs...)
    header = ["Alpha","Cl","Cd","Cdp","Cm","Top_Xtr","Bot_Xtr"]
    data = DataFrames.DataFrame(CSV.File(joinpath(path,file_name), datarow=12, header=header))
    if data[:,1] == []
        error("Unusable data file. Possibly insufficient data in $file_name.  Data must start on row 12 (XFoil polar output formatting is assumed).")
    end
    polar = Polar(-1, data[:,1], data[:,2], data[:,3], data[:,5]; optargs...)
    return polar
end
"Reads a polar as saved from a Polar object"
function read_polar2(file_name::String; path::String="", optargs...)
    header = ["Alpha","Cl","Cd","Cm"]
    data = DataFrames.DataFrame(CSV.File(joinpath(path,file_name), datarow=2, header=header))
    polar = Polar(-1, data[:,1], data[:,2], data[:,3], data[:,4]; optargs...)
    return polar
end
"Saves a polar in Polar format"
function save_polar2(self::Polar, file_name::String; path::String="")
    _file_name = file_name*(occursin(".", file_name) ? "" : ".csv")
    f = open(joinpath(path,_file_name),"w")

    # Header
    header = ["Alpha","Cl","Cd","Cm"]
    for (i,h) in enumerate(header)
        write(f, h)
        write(f, i!=size(header)[1] ? "," : "\n")
    end

    # Data
    alpha, cl = get_cl(self)
    _, cd = get_cd(self)
    _, cm = get_cm(self)
    for (i,a) in enumerate(alpha)
        write(f, "$a,$(cl[i]),$(cd[i]),$(cm[i])\n")
    end

    close(f)
end

"""
  `correction3D(self::Polar, r_over_R::Float64, chord_over_r::Float64,
                        tsr::Float64; alpha_max_corr=30, alpha_linear_min=-5,
                        alpha_linear_max=5)`

Applies 3-D corrections for rotating sections from the 2-D data.

  Parameters
  ----------
  r_over_R : float
      local radial position / rotor radius
  chord_over_r : float
      local chord length / local radial location
  tsr : float
      tip-speed ratio
  alpha_max_corr : float, optional (deg)
      maximum angle of attack to apply full correction
  alpha_linear_min : float, optional (deg)
      angle of attack where linear portion of lift curve slope begins
  alpha_linear_max : float, optional (deg)
      angle of attack where linear portion of lift curve slope ends

  Returns
  -------
  polar : Polar
      A new Polar object corrected for 3-D effects

  Notes
  -----
  The Du-Selig method :cite:`Du1998A-3-D-stall-del` is used to correct lift, and
  the Eggers method :cite:`Eggers-Jr2003An-assessment-o` is used to correct drag.
"""
function correction3D(self::Polar, r_over_R::Float64, chord_over_r::Float64,
                      tsr; alpha_max_corr=30, alpha_linear_min=-5,
                      alpha_linear_max=5)

    new_pyPolar = self.pyPolar.correction3D(r_over_R, chord_over_r, tsr,
                  alpha_linear_min=alpha_linear_min,
                  alpha_linear_max=alpha_linear_max,
                  alpha_max_corr=alpha_max_corr)

    new_polar = _pyPolar2Polar(new_pyPolar; _get_nonpypolar_args(self)...)

    return new_polar
end

"""
  `extrapolate(self::Polar, cdmax::Float64; AR=nothing, cdmin=0.001,
                    nalpha=15)`
Extrapolates force coefficients up to +/- 180 degrees using Viterna's method
:cite:`Viterna1982Theoretical-and`.

  Parameters
  ----------
  cdmax : float
      maximum drag coefficient
  AR : float, optional
      aspect ratio = (rotor radius / chord_75% radius)
      if provided, cdmax is computed from AR
  cdmin: float, optional
      minimum drag coefficient.  used to prevent negative values that can
      sometimes occur with this extrapolation method
  nalpha: int, optional
      number of points to add in each segment of Viterna method

  Returns
  -------
  polar : Polar
      a new Polar object

  Notes
  -----
  If the current polar already supplies data beyond 90 degrees then
  this method cannot be used in its current form and will just return itself.

  If AR is provided, then the maximum drag coefficient is estimated as

  cdmax = 1.11 + 0.018*AR

"""
function extrapolate(self::Polar, cdmax::Float64; AR=nothing, cdmin=0.001,
                      nalpha=15)

    new_pyPolar = self.pyPolar.extrapolate(cdmax, AR=AR, cdmin=cdmin,
                                              nalpha=nalpha)

    new_polar = _pyPolar2Polar(new_pyPolar; _get_nonpypolar_args(self)...)

    return new_polar
end

"Plots this Polar. Give it `geometry=(x,y,1)` for ploting the airfoil geometry,
where x and y are the points of the airfoil (the third number gives the size
of the plot)"
function plot(self::Polar; geometry::Bool=true, label="", style=".-",
                              cdpolar=true, rfl_style="-", fig_id="polar_curves",
                              figsize=[7, 5], rfl_figfactor=2/3,
                              title_str="automatic", polar_optargs=[],
                              legend_optargs=[(:loc, "best")],
                              to_plot=["Cl", "Cd", "Cm"],
                              fig=nothing, axs=nothing)

    # Geometry
    if geometry
        x,y = get_geometry(self)
        plot_airfoil(x, y; label=label, style=rfl_style, figfactor=rfl_figfactor, fig_id=fig_id*"_rfl")
    end

    _fig = fig!=nothing ? fig : figure(fig_id, figsize=figsize.*[length(to_plot), 1])
    _axs = axs!=nothing ? axs : _fig.subplots(1, length(to_plot))

    alpha, cl = get_cl(self)
    _, cd = get_cd(self)
    _, cm = get_cm(self)

    if title_str !== "automatic"
        _fig.suptitle(title_str)
    end

    ttl = title_str=="automatic" ? "$(self.pyPolar.Re)" : title_str

    for (pi, this_plot) in enumerate(to_plot)

        ax = _axs[pi]

        if this_plot=="Cl"
            if title_str=="automatic"; ax.title.set_text("Lift curve at Re=$(self.pyPolar.Re)"); end;
            ax.plot(alpha, cl, style, label=label; polar_optargs...)
            ax.set_xlabel(L"Angle of attack $\alpha (^\circ)$")
            ax.set_ylabel(L"C_l")
            ax.grid(true, color="0.8", linestyle="--")

        elseif this_plot=="Cd"
            if title_str=="automatic"; ax.title.set_text("Drag polar at Re=$(self.pyPolar.Re)"); end;
            ax.plot( cdpolar ? cl : alpha, cd, style, label=label; polar_optargs...)
            if cdpolar
              ax.set_xlabel(L"C_l")
            else
              ax.set_xlabel(L"Angle of attack $\alpha (^\circ)$")
            end
            ax.set_ylabel(L"C_d")
            ax.grid(true, color="0.8", linestyle="--")

        elseif this_plot=="Cm"
            if title_str=="automatic"; ax.title.set_text("Moment curve at Re=$(self.pyPolar.Re)"); end;
            ax.plot(alpha, cm, style, label=label; polar_optargs...)
            ax.set_xlabel(L"Angle of attack $\alpha (^\circ)$")
            ax.set_ylabel(L"C_m")
            ax.grid(true, color="0.8", linestyle="--")
        end

        if pi==length(to_plot) && legend_optargs!=nothing; ax.legend(; legend_optargs...); end;
    end
end

"Compares two polars. Returns [err1, err2, err3] with an error between 0 and 1
of Cl, Cd, and Cm curves, respectively."
function compare(polar1::Polar, polar2::Polar; verbose::Bool=false)
    out = []
    curves = [:cl, :cd, :cm]

    # Determines the range of alphas to compare
    min_alpha1 = minimum(polar1.pyPolar.alpha)
    min_alpha2 = minimum(polar2.pyPolar.alpha)
    max_alpha1 = maximum(polar1.pyPolar.alpha)
    max_alpha2 = maximum(polar2.pyPolar.alpha)
    min_alpha = max(min_alpha1, min_alpha2)
    max_alpha = min(max_alpha1, max_alpha2)
    if max_alpha < min_alpha # Case that ranges don't intercept
        @warn("Requested comparison of polar that don't coincide. Returning -1.")
        return [-1, -1, -1]
    end

    if verbose; println("min_alpha:$min_alpha\tmax_alpha:$max_alpha"); end;

    # Creates splines of each curve on each polar
    splines = Dict()
    for polar in [polar1, polar2] # each polar
        this_splines = []
        for curve in curves # each curve
            this_spline = Dierckx.Spline1D(polar.pyPolar.alpha, getproperty(polar.pyPolar, curve))
            push!(this_splines, this_spline)
        end
        splines[polar] = this_splines
    end

    for (i, curve) in enumerate(curves)
        f1(x) = splines[polar1][i](x)
        f2(x) = splines[polar2][i](x)
        absf1(x) = abs(f1(x))
        absf2(x) = abs(f2(x))
        dif(x) = abs(f2(x) - f1(x))

        # Area of the difference
        num = quadgk(dif, min_alpha, max_alpha)[1]
        # Areas of both
        den = max(quadgk(absf1, min_alpha, max_alpha)[1],
                      quadgk(absf2, min_alpha, max_alpha)[1])
        push!(out, num/den)

        if verbose; println("num=$num\tden=$den"); end;
    end

    return out
end

"Returns a polar that is injective in alpha (angles are unique)"
function injective(polar::Polar; start_i::Int64=2)
    alpha, cl = get_cl(polar)
    _, cd = get_cd(polar)
    _, cm = get_cm(polar)
    na = length(alpha)

    # New polar data
    out_a, out_cl, out_cd, out_cm = zeros(na), zeros(na), zeros(na), zeros(na)
    # Iterates over data averaging repeated angles
    ave_a, ave_cl, ave_cd, ave_cm, cum_i = 0.0, 0.0, 0.0, 0.0, 0
    count = 0
    for (i,this_a) in enumerate(alpha)

        # Accumulates repeated values
        ave_a += this_a
        ave_cl += cl[i]
        ave_cd += cd[i]
        ave_cm += cm[i]
        cum_i += 1

        # Averages repeated values
        next_a = i!=size(alpha)[1] ? alpha[i+1] : nothing
        if i<start_i || this_a!=next_a
            count += 1
            out_a[count] = ave_a/cum_i
            out_cl[count] = ave_cl/cum_i
            out_cd[count] = ave_cd/cum_i
            out_cm[count] = ave_cm/cum_i
            ave_a, ave_cl, ave_cd, ave_cm, cum_i = 0.0, 0.0, 0.0, 0.0, 0
        end
    end

    # Injective polar
    new_polar = Polar(get_Re(polar), out_a[1:count], out_cl[1:count],
                                    out_cd[1:count], out_cm[1:count];
                                          _get_nonpypolar_args(polar)...)

    return new_polar
end

### INTERNAL FUNCTIONS #########################################################
"Given a Python Polar object, it returns a julia Polar object"
function _pyPolar2Polar(pyPolar::PyCall.PyObject; nonpypolar_args...)

    polar = Polar(pyPolar.Re, pyPolar.alpha, pyPolar.cl, pyPolar.cd,
                  pyPolar.cm; nonpypolar_args...)

    return polar
end

"Return arguments unrelated to airfoilpreppy"
function _get_nonpypolar_args(polar::Polar)

    ma = get_Ma(polar)
    npanels = get_npanels(polar)
    ncrit = get_ncrit(polar)
    x, y = get_geometry(polar)
    xsepup = get_xsepup(polar)[2]
    xseplo = get_xseplo(polar)[2]

    return ((:Ma, ma), (:npanels, npanels), (:ncrit, ncrit), (:x, x), (:y, y),
                                        (:xsepup, xsepup), (:xseplo, xseplo))
end
### END OF airfoilprep.py CLASS ################################################


text:
#=
"""
    ig = image_geom(...)

OBSOLETE
Constructor for `ImageGeom`, where `dx,dy,dz` and `fov` and `fovz` may have units

# Arguments
- `nx::Int = 128`
- `ny::Int = nx`
- `dx::RealU = ?` (must specify one of `dx` or `fov`)
- `dy::RealU = -dx`
- `offset_x::Real = 0` (unitless)
- `offset_y::Real = 0` (unitless)
- `fov::RealU = ?` (if specified, then `nx*dx=ny*dy`)
- `nz::Int = 0`
- `dz::RealU = ?` (need one of `dz` or `zfov` if `nz > 0`)
- `zfov::RealU = ?` (if specified, then `nz*dz`)
- `offset_z::Real = 0` (unitless)
- `offsets::Symbol = :none` or :dsp
- `mask::Union{Symbol,AbstractArray{Bool}} = :all` | `:circ` | `:all_but_edge_xy`
"""
function image_geom( ;
    nx::Int = 128,
    ny::Int = nx,
    dx::RealU = NaN,
    dy::RealU = NaN,
    offset_x::Real = 0,
    offset_y::Real = 0,
    fov::RealU = NaN,
    nz::Int = 0,
    dz::RealU = NaN,
    zfov::RealU = NaN,
    offset_z::Real = 0,
    offsets::Symbol = :none,
    mask::Union{Symbol,AbstractArray{Bool}}    = :all,
)

    # handle optional arguments

    is3 = !isnan(dz) || !isnan(zfov) || (nz > 0)

    # offsets
    if offsets === :dsp
        # check that no nonzero offset was specified
        (offset_x != 0 || offset_y != 0) && throw("offsets usage incorrect")
        offset_x = 0.5
        offset_y = 0.5
        if is3
            offset_z != 0 && throw("offset_z usage incorrect")
            offset_z = 0.5
        end
    elseif offsets != :none
        throw("offsets $offsets")
    end

    # transverse plane (2d) distances
    if true
        if isnan(fov)
            isnan(dx) && throw("dx or fov required")
        #    fov = nx * dx
        else
            !isnan(dx) && throw("dx and fov?")
            dx = fov / nx
            dy = fov / ny
        end

        dy = isnan(dy) ? -dx : dy # default dy
    end

    # 3D geometry
    if is3
        if isnan(zfov)
            isnan(dz) && throw("dz or zfov required")
        #    zfov = nz * dz
        else
            !isnan(dz) && throw("dz and zfov?")
            dz = zfov / nz
        end
    end

    # mask
    if mask === :all
        mask = is3 ? Trues(nx, ny, nz) : Trues(nx, ny)
    elseif mask === :circ
        mask = image_geom_circle(nx, ny, dx, dy)
        if nz > 0
            mask = repeat(mask, 1, 1, nz)
        end
    elseif mask === :all_but_edge_xy
        mask = trues(nx,ny,max(nz,1))
        mask[  1,   :, :] .= false
        mask[end,   :, :] .= false
        mask[  :,   1, :] .= false
        mask[  :, end, :] .= false
        if !is3
            mask = mask[:,:,1]
        end
    elseif isa(mask, Symbol)
        throw("mask symbol $mask")
    elseif size(mask,1) != nx || size(mask,2) != ny || (is3 && size(mask,3) != nz)
        throw("mask size $(size(mask)), nx=$nx ny=$ny nz=$nz")
    end

    # return the object (type unstable, but used rarely so ok)
    return is3 ?
        ImageGeom((nx, ny, nz), (dx, dy, dz),
            (offset_x, offset_y, offset_z), mask) :
        ImageGeom((nx, ny), (dx, dy),
            (offset_x, offset_y), mask)
end
=#


text:
struct BlockParams
  repeat::Int
  kernel::Tuple{Int, Int}
  stride::Int
  expansion_ratio::Int
  in_channels::Int
  out_channels::Int
end

struct GlobalParams
  width_coef::Real
  depth_coef::Real
  image_size::Tuple{Int, Int}

  depth_divisor::Int
  min_depth::Union{Nothing, Int}
end

# (width_coefficient, depth_coefficient, resolution)
get_efficientnet_coefficients(model_name::String) =
  Dict(
    "b0" => (1.0, 1.0, 224),
    "b1" => (1.0, 1.1, 240),
    "b2" => (1.1, 1.2, 260),
    "b3" => (1.2, 1.4, 300),
    "b4" => (1.4, 1.8, 380),
    "b5" => (1.6, 2.2, 456),
    "b6" => (1.8, 2.6, 528),
    "b7" => (2.0, 3.1, 600),
    "b8" => (2.2, 3.6, 672))[model_name]

function get_efficientnet_params(model_name)
  block_params = [
    BlockParams(1, (3, 3), 1, 1,  32,  16),
    BlockParams(2, (3, 3), 2, 6,  16,  24),
    BlockParams(2, (5, 5), 2, 6,  24,  40),
    BlockParams(3, (3, 3), 2, 6,  40,  80),
    BlockParams(3, (5, 5), 1, 6,  80, 112),
    BlockParams(4, (5, 5), 2, 6, 112, 192),
    BlockParams(1, (3, 3), 1, 6, 192, 320)]

  width_coef, depth_coef, resolution = get_efficientnet_coefficients(model_name)
  global_params = GlobalParams(
    width_coef, depth_coef, (resolution, resolution), 8, nothing)
  block_params, global_params
end

function round_filter(filters, global_params::GlobalParams)
  global_params.width_coef ≈ 1 && return filters

  depth_divisor = global_params.depth_divisor
  filters *= global_params.width_coef
  min_depth = global_params.min_depth
  min_depth = min_depth ≡ nothing ? depth_divisor : min_depth

  new_filters = max(min_depth, (floor(Int, filters + depth_divisor / 2) ÷ depth_divisor) * depth_divisor)
  new_filters < 0.9 * filters && (new_filters += global_params.depth_divisor)
  new_filters
end


text:
using TestPackage

function main()
    return 0
end

TestPackage.@run_if_deployed main()