text:
module Particles
export AbstractParticles

    abstract type AbstractParticles end
    using Juriba.Constants
    # const J = Juriba.Constants
    # using GSL
    # const GSL_CONST = GSL.Const

    # const cee  =           GSL_CONST.CGSM_SPEED_OF_LIGHT
    # const emgm=            GSL_CONST.CGSM_MASS_ELECTRON
    # const pmgm  =          GSL_CONST.CGSM_MASS_PROTON
    # const kboltz    =      GSL_CONST.CGSM_BOLTZMANN
    # const kboltz_kev2erg = 1.6022e-9
    # const gr_to_kev   =    5.6095883571872e+29
    # const me_kev    =      511.
    # const emerg   =        (GSL_CONST.CGSM_MASS_ELECTRON*GSL_CONST.CGSM_SPEED_OF_LIGHT^2.0)
    # const charg   =        4.8e-10
    # const sigtom =         GSL_CONST.CGSM_THOMSON_CROSS_SECTION
    # const herg   =         GSL_CONST.CGSM_PLANCKS_CONSTANT_H
    # const hkev   =         (GSL_CONST.CGSM_PLANCKS_CONSTANT_H*6.2415e8)
    # const mjy   =          1.e-26
    # const re0    =         2.81794e-13
    # const gconst  =        GSL_CONST.CGSM_GRAVITATIONAL_CONSTANT
    # const sbconst  =       GSL_CONST.CGSM_STEFAN_BOLTZMANN_CONSTANT
    # const aconst   =       7.56e-15

    struct pl_params
        s::Float64
        n::Float64
    end

    struct bkn_params
        s1::Float64
        s2::Float64
        brk::Float64
        max::Float64
        m::Float64
    end

    struct th_params
        t::Float64
        n::Float64
        m::Float64
    end

    struct k_params
        t::Float64
        k::Float64
    end

    struct injection_mixed_params
        s::Float64
        t::Float64
        nth::Float64
        npl::Float64
        m::Float64
        min::Float64
        max::Float64
        cutoff::Float64
    end

    struct injection_kappa_params
        t::Float64
        k::Float64
        n::Float64
        m::Float64
    end

    struct injection_pl_params
        s::Float64
        n::Float64
        m::Float64
        max::Float64
    end

    struct injection_bkn_params
        s1::Float64
        s2::Float64
        brk::Float64
        max::Float64
        m::Float64
        n::Float64
    end

    function count_particles(p::AbstractParticles)
        temp = 0.0
        for i ∈ 1:p.size-1
            temp = temp + 0.5*(p.p[i+1] - p.p[i])*(p.ndens[i+1] + p.ndens[i])
        end
        return temp
    end

    function count_particles_energy(p::AbstractParticles)
        temp = 0.0
        for i ∈ 1:p.size-1
            temp = temp + 0.5*(p.p[i+1] - p.p[i])*(p.gdens[i+1] + p.gdens[i])
        end
        return temp
    end

    function av_p(p::AbstractParticles)
        temp = 0.0
        for i ∈ 1:p.size-1
            temp = temp + 0.5*(p.p[i+1] - p.p[i])*
                            (p.p[i+1]*p.ndens[i+1] + p.p[i]*p.ndens[i])
        end
        return temp/count_particles(p)
    end

    function av_gamma(p::AbstractParticles)
        return sqrt((av_p(p)/(p.mass_gr*cee))^2 + 1)
    end

    function av_psq(p::AbstractParticles)
        temp = 0.0
        for i ∈ 1:p.size - 1
            temp = temp + 0.5*(p.p[i+1] - p.p[i])*p.p[i+1]^2*
                        p.ndens[i+1]*p[i]^2*p.ndens[i]
        end

        return temp/count_particles(p)
    end

    function av_gammasq(p::AbstractParticles)

        return sqrt(av_psq(p)/((p.mass_gr*cee)^2) + 1)

    end

    function initialize_gdens(p::AbstractParticles)
        for i ∈ 1:p.size
            p.gdens[i] = p.ndens[i]*p.gamma[i]*p.mass_gr*cee/sqrt(p.gamma[i]^2 - 1)
        end
    end

    function initialize_pdens(p::AbstractParticles)
        for i ∈ 1:p.size
            ndens[i] = p.gdens[i]*p.p[i]/((p.mass_gr*cee)^2*
                    sqrt(p.p[i]/(p.mass_gr*cee)^2 + 1))
        end
    end

    function gdens_differentiate(p::AbstractParticles)
        temp = Vector{Float64}(undef, p.size)

        for i in eachindex(temp)
            temp[i] = p.gdens[i]/p.gamma[i]
        end

        for i ∈ 1:p.size - 1
            p.gdens_diff[i] = (temp[i+1] - temp[i])/
                            (p.mass_gr*cee^2*p.gamma[i+1] - p.gamma[i])
        end

        p.gdens_diff[end] = p.gdens_diff[end-1]
    end

    function set_mass(p::AbstractParticles, m)
        p.mass_gr = m
        p.mass_kev = m*gr_to_kev
    end

end


text:
# SPDX-License-Identifier: MIT

module BayesianModels

import Distributions
import ForwardDiff
using LinearAlgebra
using SpecialFunctions
using Zygote

# Custom type compressed lower triangular matrix.
include("utils/clowertri.jl")

export vech
export CompressedLowerTriangular
export CompressedSymmetric

include("distributions/expfamily.jl")

export η
export μ
export ξ
export A
export kldiv
export sample
export stdparams
export unpack

include("distributions/normal.jl")
include("distributions/gamma.jl")
include("distributions/wishart.jl")
include("distributions/normalwishart.jl")
include("distributions/jointnormal.jl")

export JointNormalFixedCov
export Normal
export NormalDiag
export NormalIso
export Gamma
export Wishart
export NormalWishart

include("models/model.jl")

export loglikelihood
export newposterior
export elbo
export getstats

include("models/normal.jl")
include("models/gsm.jl")
include("models/gsmphonetic.jl")

export NormalModel
export GSM
export GSMPhonetic

end # module

ext: jl

max_stars_repo_name: BUTSpeechFIT/PPCA

max_issues_repo_name: BUTSpeechFIT/PPCA

max_forks_repo_name: BUTSpeechFIT/PPCA


text:
################################################################################
#                        BFV Scheme definition
################################################################################

struct BFVParams <: SHEShemeParams
    # The Cypertext ring over which operations are performed
    ℛ
    # The big ring used during multiplication
    ℛbig
    # The plaintext ring.
    ℛplain
    relin_window
    σ
    Δ
end
scheme_name(p::Type{BFVParams}) = "BFV"

ℛ_plain(p::BFVParams) = p.ℛplain
ℛ_cipher(p::BFVParams) = p.ℛ

function π⁻¹(params::BFVParams, plaintext)
    @fields_as_locals params::BFVParams
    params.Δ * oftype(zero(params.ℛ), ℛplain(plaintext))
end

function π(params::BFVParams, b)
    @fields_as_locals params::BFVParams
    ℛplain(map(x->coefftype(ℛplain)(convert(Integer, mod(SignedMod(divround(x, Δ)), modulus(base_ring(ℛplain))))), NTT.coeffs_primal(b)))
end

𝒩(params::BFVParams) = RingSampler(params.ℛ, DiscreteNormal(0, params.σ))
𝒢(params::BFVParams) = RingSampler(params.ℛ, DiscreteNormal(0, params.σ))

mul_expand(params::BFVParams, c::CipherText) = map(c->switch(params.ℛbig, c), c.cs)
function mul_contract(params::BFVParams, c)
    @fields_as_locals params::BFVParams
    map(c) do e
        switch(ℛ, multround(e, modulus(base_ring(ℛplain)), modulus(coefftype(ℛ))))
    end
end

################################################################################
#                 BFV Noise Modeling / Parameter Generation
################################################################################

# Matches parameter generation in PALISADE
function BFVParams(p, σ=8/√(2pi), α=9, r=1; eval_mult_count = 0, security = CryptParameters.HEStd_128_classic, relin_window=1)
    @assert r >= 1
    Berr = σ*√(α)
    Bkey = Berr
    δ(n) = 2*√(n)
    Vnorm(n) = Berr * (1 + 2*δ(n)*Bkey)

    function nRLWE(q)
        if isa(security, CryptParameters.StdSecurity)
            CryptParameters.std_ring_dim(CryptParameters.HEStd_error, security, ceil(log2(q)))
        else
            # The security parameter is interpreted as the hermite factor as
            # in PALISADE.
            log2(q / σ) / (4 * log2(security));
        end
    end

    n = 512
    q = 0

    if eval_mult_count > 0
        w = 2^r
        ϵ₁(n) = 4 / δ(n)*Bkey
        C₁(n) = (1 + ϵ₁(n))*δ(n)^2*p*Bkey
        C₂(n, qPrev) =
            δ(n)^2*p*Bkey*(Bkey + p^2) +
            δ(n)*(floor(log2(qPrev) / r) + 1)*w*Berr
        qBFV(n, qPrev) =
            p^2 + 2p*(
                C₁(n)^eval_mult_count * Vnorm(n) +
                eval_mult_count*C₁(n)^(eval_mult_count-1)*C₂(n, qPrev))

        qPrev = 1e6
        q = qBFV(n, qPrev)
        qPrev = q

        while nRLWE(q) > n
            while nRLWE(q) > n
                n *= 2
                # TODO: So in original, but is this right?
                # Shouldn't we set qPrev = q first.
                q = qBFV(n, qPrev)
                qPrev = q
            end

            q = qBFV(n, qPrev)

            while abs(q - qPrev) > 0.001q
                qPrev = q
                q = qBFV(n, qPrev)
            end
        end
    end

    bits = ceil(Integer, log2(q))+1
    T = bits > 256 ? Int512 : bits > 128 ? Int256 : Int128
    qPrime = nextprime(T(2)^(ceil(Int, log2(q))+1) + 1, 1; interval=2n)
    largebits = 2*ceil(Int, log2(q)) + ceil(Int, log2(p)) + 3
    Tlarge = largebits > 256 ? Int512 : largebits > 128 ? Int256 : Int128
    qLargeBig = nextprime(big(2)^largebits + 1, 1; interval=2n)
    qPrimeLarge = Tlarge(qLargeBig)

    Δ = div(qPrime, p)

    𝔽 = GaloisField(qPrime)
    ℛ = NegacyclicRing{𝔽, n}(GaloisFields.minimal_primitive_root(𝔽, 2n))
    𝔽big = GaloisField(qPrimeLarge)
    r = GaloisFields.minimal_primitive_root(𝔽big, 2n)
    ℛbig = NegacyclicRing{𝔽big, n}(r)

    BFVParams(ℛ, ℛbig, plaintext_space(ℛ, p), relin_window, σ, Δ)
end

"""
Compute the *invariant noise budget*, defined by:

        -log2(2‖v‖) = log2(q) - log2(q‖v‖) - 1.

If this quantity is >0, the ciphertext is expected to decrypt correctly with
high probability.

This notion of noise was first introduced by the SEAL HE library. See [CLP19]
for details.

[CLP19] Anamaria Costache, Kim Laine, and Rachel Player
        "Homomorphic noise growth in practice: comparing BGV and FV"
        https://eprint.iacr.org/2019/493.pdf
"""
function invariant_noise_budget(pk::PrivKey{BFVParams}, c::CipherText{<:Any, BFVParams})
    @fields_as_locals pk::PrivKey
    @fields_as_locals params::BFVParams

    b = c[1]
    spow = secret

    for i = 2:length(c)
        b += spow*c[i]
        spow *= secret
    end

    ℛplain = plaintext_space(params)

    function birem(x)
        r = rem(x, Δ)
        if r > div(Δ, 2)
            return Δ - r
        else
            return r
        end
    end

    # -log2(2‖v‖) = log(q) - log(t) - 1 - max_i log2(Δ |v_i|)
    log2(modulus(coefftype(ℛ))) - log2(modulus(coefftype(ℛplain))) - 1 -
        maximum(log2(birem(c.n)) for c in NTT.coeffs_primal(b))
end
invariant_noise_budget(kp::KeyPair, c::CipherText) =
    invariant_noise_budget(kp.priv, c)

export invariant_noise_budget

################################################################################
#                 BFV Computational utilities
################################################################################

function multround(e::SignedMod, a::Integer, b::Integer)
    div(e * a, b, RoundNearestTiesAway)
end
function multround(e::BigInt, a::Integer, b::Integer)
    div(e * a, b, RoundNearestTiesAway)
end
multround(e, a::Integer, b::fmpz) = multround(e, a, BigInt(b))
multround(e::fmpz, a::Integer, b::Integer) = multround(BigInt(e), a, b)
multround(e::fmpz, a::Integer, b::fmpz) = multround(BigInt(e), a, BigInt(b))

function multround(e, a::Integer, b)
    oftype(e, broadcast(NTT.coeffs_primal(e)) do x
        if isa(x, AbstractAlgebra.Generic.Res{fmpz})
            multround(BigInt(Nemo.lift(x)), a, b)
        else
            multround(SignedMod(x), a, b).x
        end
    end)
end

Nemo.modulus(e::GaloisFields.PrimeField) = GaloisFields.char(e)
Nemo.lift(e::GaloisFields.PrimeField) = e.n
Nemo.lift(e::Nemo.nmod) = lift(Nemo.ZZ, e)

divround(e::Integer, q::Integer) = div(e, q, RoundNearestTiesAway)
divround(e::fmpz, q::Integer) = divround(BigInt(e), q)
function divround(e, d::Integer)
    div(SignedMod(e), d, RoundNearestTiesAway).x
end

function switchel(T, e)
    q = modulus(e)
    halfq = q >> 1
    diff = modulus(T) > q ? modulus(T) - q : q - modulus(T)
    en = convert(Integer, e)
    if (q < modulus(T))
        if en > halfq
            return T(en + diff)
        else
            return T(en)
        end
    else
        if en > halfq
            return T(en - diff)
        else
            return T(en)
        end
    end
end

function switch(ℛ, e)
    ℛ(broadcast(NTT.coeffs_primal(e)) do x
        switchel(coefftype(ℛ), x)
    end)
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "aiger"
version = v"1.9.9"

# Collection of sources required to complete build
sources = [
    ArchiveSource("http://fmv.jku.at/aiger/aiger-1.9.9.tar.gz", "1e50d3db36f5dc5ed0e57aa4c448b9bcf82865f01736dde1f32f390b780350c7"),
    DirectorySource("./bundled")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
for f in ${WORKSPACE}/srcdir/patches/*.patch; do
    atomic_patch -p1 ${f}
done
cd aiger-1.9.9/
./configure.sh
make -j
mkdir $prefix/bin
PREFIX=$prefix make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:i686, libc=:glibc),
    Linux(:x86_64, libc=:glibc),
    Linux(:aarch64, libc=:glibc),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf),
    Linux(:powerpc64le, libc=:glibc),
    Linux(:i686, libc=:musl),
    Linux(:x86_64, libc=:musl),
    Linux(:aarch64, libc=:musl),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf),
    FreeBSD(:x86_64)
]

# The products that we will ensure are always built
products = [
    ExecutableProduct("aigdd", :aigdd),
    ExecutableProduct("aigreset", :aigreset),
    ExecutableProduct("aigmiter", :aigmiter),
    ExecutableProduct("aignm", :aignm),
    ExecutableProduct("soltostim", :soltostim),
    ExecutableProduct("aigmove", :aigmove),
    ExecutableProduct("aigflip", :aigflip),
    ExecutableProduct("aigunconstraint", :aigunconstraint),
    ExecutableProduct("aigstrip", :aigstrip),
    ExecutableProduct("aigsim", :aigsim),
    ExecutableProduct("aigtodot", :aigtodot),
    ExecutableProduct("aigtoblif", :aigtoblif),
    ExecutableProduct("aigunroll", :aigunroll),
    ExecutableProduct("wrapstim", :wrapstim),
    ExecutableProduct("aigor", :aigor),
    ExecutableProduct("aigtosmv", :aigtosmv),
    ExecutableProduct("aigtocnf", :aigtocnf),
    ExecutableProduct("andtoaig", :andtoaig),
    ExecutableProduct("aigjoin", :aigjoin),
    ExecutableProduct("smvtoaig", :smvtoaig),
    ExecutableProduct("aigfuzz", :aigfuzz),
    ExecutableProduct("aiginfo", :aiginfo),
    ExecutableProduct("bliftoaig", :bliftoaig),
    ExecutableProduct("aigsplit", :aigsplit),
    ExecutableProduct("aigtoaig", :aigtoaig),
    ExecutableProduct("aigand", :aigand)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
export cascade

"""
A `cascade` consists of the eight piles of cards on the main
part of a FreeCell tableau.
+ `cascade()` creates an empty cascade
+ `cascade(deck())` creates a full cascade from a deck of cards
"""
struct cascade
    piles::Vector{Vector{Card}}
    function cascade()
        p = Vector{Vector{Card}}(undef, 8)
        for k = 1:8
            p[k] = Card[]
        end
        new(p)
    end
end

function cascade(d::Vector{Card})
    C = cascade()
    for i = 1:4
        a = 7 * (i - 1) + 1
        b = 7i
        C.piles[i] = d[a:b]
    end

    for k = 5:8
        i = k - 4
        a = 28 + 6 * (i - 1) + 1
        b = 28 + 6i
        C.piles[k] = d[a:b]
    end

    return C
end

function string(C::cascade)::String
    result = ""
    for k = 1:8
        result *= row_string(C.piles[k])
        if k < 8
            result *= "\n"
        end
    end
    return result
end

function show(io::IO, C::cascade)
    print(io, string(C))
end

"""
`add_check(pile::Vector{Card}, PC::Card)` checks to see if `PC`
may be added to the bottom of the cascade pile `pile`.
"""
function add_check(pile::Vector{Card}, PC::Card)::Bool
    if length(pile) == 0
        return true
    end
    last = pile[end]
    return rank(last) == rank(PC) + 1 && color(last) != color(PC)
end

"""
`add_check(C::cascade, k::Int, PC::Card)::Bool`
checks to see if card `PC` may be added to pile `k`
of the cascade `C`.
"""
function add_check(C::cascade, k::Int, PC::Card)::Bool
    if k < 0 || k > 8      # no such pile!
        return false
    end

    return add_check(C.piles[k], PC)
end

"""
`add_card(C::cascade, k::Int, PC::Card)::Bool` adds the card `PC`
to pile `k` of `C`.
"""
function add_card(C::cascade, k::Int, PC::Card)::Bool
    if !add_check(C, k, PC)
        return false
    end
    push!(C.piles[k], PC)
    return true
end

"""
`get_card(C::cascade,k)` returns the bottom card in pile `k` of
cascade `C`.
"""
function get_card(C::cascade, k::Int)::Card
    p = C.piles[k]
    if length(p) == 0
        error("Attempt to get card from empty cascade pile $k")
    end
    return p[end]
end

"""
`rm_card(C::cascade, k::Int)::Card` removes the last card from pile `k`
of cascade `C`.
"""
function rm_card(C::cascade, k::Int)::Card
    PC = get_card(C, k)
    C.piles[k] = C.piles[k][1:end-1]
    return PC
end

"""
`score(pile::Vector{Card})` subtacts 1 for each card on the pile,
but adds 1 back in for builds at the bottom.
"""
function score(pile::Vector{Card})
    n = length(pile)

    if n == 0
        return 20

    end

    if n < 2
        return 5
    end

    rks = rank.(pile)
    drks = diff(rks)

    if all(drks .<= 0)
        return 0
    end

    i = findfirst(t -> t > 0, drks)  # find first ascent in pile (it's bad)
    return -n + i

    # result = -Int(round(n^1.5))

    # return result

    # if n < 2
    #     return result
    # end

    # run = 0
    # for k = n:-1:2
    #     if color(pile[k]) != color(pile[k-1]) && rank(pile[k]) == rank(pile[k-1]) - 1
    #         run += 1
    #     else
    #         break
    #     end
    # end
    # return result + run

end

"""
`score(C::cascade)` is the sum of the scores of the piles in `C`
"""
score(C::cascade) = sum(score(p) for p in C.piles)


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

using BAT
using Test

using LinearAlgebra, Random
using Distributions, PDMats
using ArraysOfArrays
using ValueShapes

struct test_mvdist <: Distribution{Multivariate, Continuous}
    d::Distribution{Multivariate, Continuous}
end

test_sampler(tmv::test_mvdist) = tmv.d

@testset "proposaldist" begin
    @testset "BAT.GenericProposalDist" begin
        gpd = BAT.GenericProposalDist{test_mvdist, typeof(test_sampler)}(
             test_mvdist(MvNormal([0.0, 1.0], [1. 0.5; 0.5 2])),
             test_sampler)

        @test typeof(gpd.d) <: Distribution{Multivariate, Continuous}
        @test typeof(gpd.s) <: MvNormal
        @test gpd.s.μ ≈ [0.0, 1.0]
        @test gpd.sampler_f == test_sampler

        gpd2 = BAT.GenericProposalDist(
             test_mvdist(MvNormal([0.0, 1.0], [1. 0.5; 0.5 2])),
             test_sampler)

        @test typeof(gpd.d) <: Distribution{Multivariate, Continuous}
        @test typeof(gpd.s) <: MvNormal
        @test gpd.s.μ ≈ [0.0, 1.0]
        @test gpd.sampler_f == test_sampler

        mvN_gpd = @inferred BAT.GenericProposalDist(MvNormal([0.0, 1.0], [1. 0.5; 0.5 2]))
        @test mvN_gpd.sampler_f == BAT.bat_sampler
        @test typeof(mvN_gpd.d) <: MvNormal

        mv_gpd = @inferred BAT.GenericProposalDist(MvTDist, 2, 5)
        @test typeof(mv_gpd.d) <: Distributions.GenericMvTDist
        @test length(mv_gpd.d) == 2
        @test mv_gpd.d.df == 5

        gpd3 = similar(gpd2, test_mvdist(MvNormal([0.0, 1.0], [1. 0.5; 0.5 2])))
        @test gpd3.sampler_f == test_sampler
        @test typeof(gpd3.d) <: test_mvdist

        gpd3 = @inferred similar(mv_gpd, MvNormal([0.0, 1.0], [1. 0.5; 0.5 2]))
        @test gpd3.sampler_f == BAT.bat_sampler
        @test typeof(gpd3.d) <: MvNormal

        apd = @inferred convert(BAT.AbstractProposalDist, gpd3, Float64, 2)
        @test typeof(apd) <: BAT.AbstractProposalDist
        @test_throws ArgumentError convert(BAT.AbstractProposalDist,
            gpd3, Float64, 3)

        d = MvTDist(1.5, zeros(2), PDMat([1. 0.5; 0.5 2]))
        gpd = @inferred BAT.GenericProposalDist(d)
        @test Matrix(BAT.get_cov(gpd)) ≈ [1. 0.5; 0.5 2]
        gpd = @inferred BAT.set_cov(gpd, Matrix(ScalMat(2, 1)))

        @test Matrix(BAT.get_cov(gpd)) ≈ Matrix{Float64}(I, 2, 2)

        @test totalndof(gpd) == 2
        @test issymmetric(gpd)

        d = MvTDist(1.5, ones(2), PDMat([1. 0.5; 0.5 2]))
        gpd = @inferred BAT.GenericProposalDist(d)
        @test !issymmetric(gpd)

        gpd = @inferred BAT.GenericProposalDist(MvTDist, Float64, 3, 2.0)
        @test totalndof(gpd) == 3
        @test typeof(gpd) <: BAT.GenericProposalDist{<:Distributions.GenericMvTDist, typeof(BAT.bat_sampler)}
        @test gpd.d.df ≈ 2.0

    end

    @testset "BAT.proposaldist_logpdf" begin
        d = MvTDist(1.5, zeros(2), PDMat([1. 0.5; 0.5 2]))
        gpd = @inferred BAT.GenericProposalDist(d)

        p = zeros(3)

        @test BAT.proposaldist_logpdf(
            gpd, [0.0,0.0], [0.1,-0.1]) ≈ -2.1441505
    end

    @testset "proposal_rand" begin
        d = MvTDist(1.5, zeros(2), PDMat([1. 0.5; 0.5 2]))
        gpd = @inferred BAT.GenericProposalDist(d)

        BAT.proposal_rand!(MersenneTwister(7002), gpd, zeros(2), [2.5, -1.0]) ≈
            [1.942373680136, 3.8454891831]
    end

    @testset "BAT.ProposalDistSpec" begin
        mvTps = @inferred BAT.MvTDistProposal(4.0)
        @test mvTps.df ≈ 4.0
        @test BAT.MvTDistProposal().df ≈ 1.0

        gpd = @inferred mvTps(Float64, 3)
        @test totalndof(gpd) == 3
        @test typeof(gpd) <: BAT.GenericProposalDist{<:Distributions.GenericMvTDist, typeof(BAT.bat_sampler)}
        @test gpd.d.df ≈ 4.0
    end
end


text:
using PyPlot;

function time_avg(ts, ys)
  avgs = zeros(length(ts));
  avgs[1] = ys[1];
  sum = 0.0;
  for i=2:length(ts)
    dt = ts[i] - ts[i-1];
    sum += dt * (ys[i] + ys[i-1]) / 2.0;
    avgs[i] = sum / ts[i];
  end
  return avgs;
end

#=
for datafile in filter(x->contains(x, "linear-spring_beta"), readdir("."))
  basename = split(datafile, ',')[1];
  data, headers = readdlm(datafile, ','; header=true);
  ts = vec(data[:,1]);
  Us = vec(data[:,2]);
  xs = vec(data[:,3]);
  xsqs = vec(data[:,4]);

  plot(ts, time_avg(ts, Us));
  xlabel("Number of steps");
  ylabel("\$ \\langle U \\rangle \$");
  savefig(basename * "_U.png");
  clf();

  plot(ts, time_avg(ts, xs));
  xlabel("Number of steps");
  ylabel("\$ \\langle x \\rangle \$");
  savefig(basename * "_x.png");
  clf();

  plot(ts, time_avg(ts, xsqs));
  xlabel("Number of steps");
  ylabel("\$ \\langle x^2 \\rangle \$");
  savefig(basename * "_xsq.png");
  clf();
end
=#

clf();
for (y, ystr, loc) in zip(2:4, ["U", "x", "x^2"], [4, 1, 4])
  for (beta, bstr) in zip([0.1; 1.0; 5.0; 10.0], ["00100"; "01000"; "05000"; "10000"])
    for (delta, dstr) in zip([0.05; 0.25; 1.25; 4.0], ["00050", "00250", "01250", "04000"])
      datafile = "linear-spring_beta-$(bstr)_delta-$(dstr).csv";
      data, headers = readdlm(datafile, ','; header=true);
      ts = vec(data[:,1]);
      ys = vec(data[:,y]);

      plot(ts, time_avg(ts, ys), label="\$\\delta = $delta\$");
    end
    xlabel("Number of steps");
    ylabel("\$ \\langle $ystr \\rangle \$");
    legend(; loc=loc);
    savefig("linear-spring_beta-$(bstr)_$ystr.png");
    clf();
  end
end

max_stars_repo_name: grasingerm/port-authority

max_issues_repo_name: grasingerm/port-authority

max_forks_repo_name: grasingerm/port-authority


text:
export ray_class_group, stable_subgroups

add_verbose_scope(:RayFacElem)
add_assert_scope(:RayFacElem)

###############################################################################
#
#  Map Type
#
###############################################################################

mutable struct MapRayClassGrp{T} <: Map{T, FacElemMon{Hecke.NfOrdIdlSet}, HeckeMap, MapRayClassGrp}
  header::Hecke.MapHeader
  modulus_fin::NfOrdIdl #The finite part of the modulus
  modulus_inf::Array{InfPlc,1} #The infinite part of the modulus
  fact_mod::Dict{NfOrdIdl, Int} #The factorization of the finite part of the modulus

  #Dictionaries to cache preimages. Used in the action on the ray class group
  prime_ideal_preimage_cache::Dict{NfOrdIdl, GrpAbFinGenElem}
  prime_ideal_cache::Array{NfOrdIdl, 1}

  evals::Array{NfOrdQuoRingElem,1}# Evaluations of the units and class group generators.
  quots::Array  #Quotients of the ring by p^n for p dividing the modulus
  idemps::Array{Tuple{NfOrdElem, NfOrdElem},1} #Idempotents for discrete logarithm
  coprime_elems::Array{nf_elem,1}

  tame_mult_grp::Dict{NfOrdIdl,Tuple{NfOrdElem,fmpz,Function}} #The multiplicative group, tame part
  wild_mult_grp::Dict{NfOrdIdl,Tuple{Array{NfOrdElem,1},Array{fmpz,1},Function}} #Multiplicative group, wild part

  function MapRayClassGrp{T}() where {T}
    z = new{T}()
    z.prime_ideal_preimage_cache = Dict{NfOrdIdl, GrpAbFinGenElem}()
    return z
  end
end

###############################################################################
#
#  Ray Class Group interface
#
###############################################################################

doc"""
***
    ray_class_group(m::NfOrdIdl, inf_plc::Array{InfPlc,1}=InfPlc[]; p_part,n_quo)

> Given a modulus with finite part $m$ and infinite part inf_plc, it returns the Ray Class Group Cl_m. If p_part is given, the function will return the largest quotient of the Ray Class Group of p-power order. If n_quo is given, it will return the quotient of the Ray Class Group by n

"""

function ray_class_group(m::NfOrdIdl, inf_plc::Array{InfPlc,1}=InfPlc[]; p_part=0, n_quo=0)

  if p_part!=0
    @assert isprime(fmpz(p_part))
    return ray_class_group_p_part(p_part, m, inf_plc)
  elseif n_quo!=0
    return ray_class_group_quo(n_quo,m,inf_plc)
  else
    return ray_class_group_fac_elem(m,inf_plc)
  end

end

###############################################################################
#
#  Functions for the evaluation of factored elements
#
###############################################################################

#
#  Multiple elements evaluation
#

function fac_elems_eval(O::NfOrd, Q::NfOrdQuoRing, elems::Array{FacElem{nf_elem, AnticNumberField},1}, lp::Dict{NfOrdIdl, Int}, exponent::fmpz)

  newelems=_preproc(O,elems,exponent)
  quots=[]
  idemps=Tuple{NfOrdElem, NfOrdElem}[]
  el=[Q(1) for i=1:length(newelems)]
  I=ideal(O,1)
  for (p,vp) in lp
    q=p^vp
    y, Qn=_eval_quo(O, newelems, p, q, anti_uniformizer(p), vp, exponent)
    push!(quots,Qn)
    a,b=idempotents(I,q)
    push!(idemps,(a,b))
    for i=1:length(el)
      el[i]=Q(y[i])*Q(a)+el[i]*Q(b)
    end
    I=I*q
  end
  return el, quots, idemps

end

#
#  Reduces the elements modulo the exponents and computes a representation as a product of elements in O
#

function _preproc(O::NfOrd, elems::Array{FacElem{nf_elem, AnticNumberField},1}, exponent::fmpz)

  assure_has_basis_mat_inv(O)
  M = O.tcontain
  newelems=FacElem{NfOrdElem, NfOrd}[]
  for el in elems
    x=Dict{NfOrdElem, fmpz}()
    for (f,k) in el.fac
      l=mod(k,exponent)
      if !iszero(l)
        elem_to_mat_row!(M.num, 1, M.den, f)
        M = mul!(M, M, O.basis_mat_inv)
        if M.den==1
          el=O(vec(Array(M.num)))
          if haskey(x,el)
            x[el]+= l
          else
            x[el]=l
          end
        else
          d=O(M.den)
          n=O(vec(Array(M.num)))
          if haskey(x,n)
            x[n]=mod(x[n]+l,exponent)
          else
            x[n]=l
          end
          if haskey(x,d)
            x[d]=mod(x[d]+exponent-l,exponent)
          else
            x[d]=exponent-l
          end
        end
      end
    end
    if !isempty(x)
      push!(newelems, FacElem(x))
    else
      push!(newelems,FacElem(Dict(O(1)=> 1)))
    end
  end
  return newelems

end

function _eval_quo(O::NfOrd, elems::Array{FacElem{NfOrdElem, NfOrd},1}, p::NfOrdIdl, q::NfOrdIdl, anti_uni::nf_elem, mult::Int, exp::fmpz)

  if mult==1
    if nbits(p.minimum)<64
      @vtime :RayFacElem 2 Q,mQ=ResidueFieldSmall(O,p)
      el=[Q(1) for i=1:length(elems)]
      for i=1:length(elems)
        J=elems[i]
        for (f,k) in J.fac
          act_el=f
          if mQ(act_el)!=0
            el[i]*=mQ(act_el)^k
            continue
          end
          val=valuation(act_el,p)
          act_el=O(act_el*(anti_uni^val),false)
          el[i]*= mQ(act_el)^k
        end
      end
    else
      @vtime :RayFacElem 2 Q,mQ=ResidueField(O,p)
      el=[Q(1) for i=1:length(elems)]
      for i=1:length(elems)
        J=elems[i]
        for (f,k) in J.fac
          if mQ(f)!=0
            el[i]*=mQ(f)^k
            continue
          end
          val=valuation(f,p)
          act_el=O(f*(anti_uni^val),false)
          el[i]*= mQ(act_el)^k
        end
      end
    end
    return [mQ\el[i] for i=1:length(el)], (Q,mQ)
  else
    @vtime :RayFacElem 2 Q,mQ=quo(O,q)
    el=[Q(1) for i=1:length(elems)]
    for i=1:length(elems)
      J=elems[i]
      for (f,k) in J.fac
        act_el=f
        if mod(act_el, p)!=0
          el[i]*=mQ(act_el)^k
          continue
        end
        val=valuation(act_el,p)
        act_el=O(act_el*(anti_uni^val),false)
        el[i]*=Q(act_el)^k
      end
    end
    return [el[i].elem for i=1:length(el)], Q
  end

end

#
#  Single element evaluation (for the disclog)
#

function _fac_elem_evaluation(O::NfOrd, Q::NfOrdQuoRing, quots::Array, idemps::Array, J::FacElem{nf_elem,AnticNumberField}, primes::Dict{NfOrdIdl, Int}, exponent::fmpz)

  assure_has_basis_mat_inv(O)
  M=O.tcontain
  element=Q(1)
  i=0
  #Reduce the exponents and reduce to elements in O
  x=Dict{NfOrdElem, fmpz}()
  for (f,k) in J.fac
    l=mod(k,exponent)
    if !iszero(l)
      elem_to_mat_row!(M.num, 1, M.den, f)
      M = mul!(M, M, O.basis_mat_inv)
      if M.den==1
        el=O(vec(Array(M.num)))
        if haskey(x,el)
          x[el]=mod(x[el]+l,exponent)
        else
          x[el]=l
        end
      else
        d=O(M.den)
        n=O(vec(Array(M.num)))
        if haskey(x,n)
          x[n]=mod(x[n]+l,exponent)
        else
          x[n]=l
        end
        if haskey(x,d)
          x[d]=mod(x[d]-l,exponent)
        else
          x[d]=exponent-l
        end
      end
    end
  end
  if isempty(x)
    return element.elem
  end
  tobeeval=FacElem(x)
  for (p,vp) in primes
    i+=1
    y=_eval_quo(O, quots[i], tobeeval, p, anti_uniformizer(p), vp)
    a,b=idemps[i]
    element=Q(Q(y)*Q(a)+element*Q(b))
  end
  return element.elem

end

function _eval_quo(O::NfOrd, Q1, J::FacElem{NfOrdElem, NfOrd}, p::NfOrdIdl, anti_uni::nf_elem,  mult::Int)
  if mult==1
    Q=Q1[1]
    mQ=Q1[2]
    el=Q(1)
    for (f,k) in J.fac
      act_el=f
      if mQ(act_el)!=0
        el*=mQ(act_el)^k
        continue
      end
      val=valuation(act_el,p)
      act_el=O(act_el*(anti_uni^val),false)
      el*= mQ(act_el)^k
    end
    return mQ\el
  else
    Q=Q1
    el=Q(1)
    for (f,k) in J.fac
      act_el=f
      if mod(act_el, p)!=0
        el*=Q(act_el)^k
        continue
      end
      val=valuation(act_el,p)
      act_el=O(act_el*(anti_uni^val),false)
      el*= Q(act_el)^k
    end
    return el.elem
  end

end

###############################################################################
#
#  Ray Class Group - Auxiliary functions
#
###############################################################################

#
# Function that finds the generators of the infinite part
#
function carlos_units(O::NfOrd)
  try
    c = _get_carlos_units_of_order(O)
    return c
  catch
    K=O.nf
    p=real_places(K)
    S=DiagonalGroup([2 for i=1:length(p)])

    function logS(x::Array{Int, 1})
      return S([x[i] > 0 ? 0 : 1 for i=1:length(x)])
    end

    s = typeof(S[1])[]
    g = elem_type(O)[]
    u, mu = sub(S, s, false)
    b = 10
    cnt = 0
    while b > 0
      a = rand(O, b)
      if a==0
        continue
      end
      emb=signs(K(a),p)
      t=S([emb[x]==1 ? 0 : 1 for x in collect(keys(emb))])
      if !Hecke.haspreimage(mu, t)[1]
        push!(s, t)
        push!(g, a)
        u, mu = sub(S, s, false)
        if order(u) == order(S)
          break
        end
      else
        cnt += 1
        if cnt > 1000
          b *= 2
          cnt = 0
        end
      end
    end
    if b <= 0
      b = 10
      cnt = 0
      bas = lll_basis(O)
      while true
        @assert b>0
        a = rand(bas, 1:b)
        if a==0
          continue
        end
        emb=signs(a,p)
        t=S([emb[x]==1 ? 0 : 1 for x in collect(keys(emb))])
        if !Hecke.haspreimage(mu, t)[1]
          push!(s, t)
          push!(g, O(a,false))
          u, mu = sub(S, s, false)
          if order(u) == order(S)
            break
          end
        else
          cnt += 1
          if cnt > 1000
            b *= 2
            cnt = 0
          end
        end
      end
    end
    hS = Hecke.GrpAbFinGenMap(S, S, vcat([x.coeff for x in s]))   # Change of coordinates so that the canonical basis elements are mapped to the elements found above
    r = elem_type(O)[]
    for i=1:length(p)
      y = haspreimage(hS,S[i])[2]
      push!(r, prod([g[i]^Int(y[i]) for i=1:length(p)]))
    end

    function exp(A::GrpAbFinGenElem)

      s=O(1)
      if iszero(A.coeff)
        return s
      end
      for i=1:length(p)
        if Int(A.coeff[1,i]) == 1
          s=s*r[i]
        end
      end
      return s
    end

    function log(B::nf_elem)
      emb=Hecke.signs(B,p)
      return S([emb[x]==1 ? 0 : 1 for x in collect(keys(emb))])
    end

    function log(B::FacElem{nf_elem})
      emb=Hecke.signs(B,p)
      return S([emb[x]==1 ? 0 : 1 for x in collect(keys(emb))])
    end

    _set_carlos_units_of_order(O, (S,exp,log))
    return (S,exp,log)
  end
end

function _infinite_primes(O::NfOrd, p::Array{InfPlc,1}, m::NfOrdIdl)

    K=O.nf
    if p==real_places(K)
      S,exp1,log1= carlos_units(O)
      function exp2(a::GrpAbFinGenElem)
        return m.gen_one*exp1(a)
      end
      return S, exp2, log1
    end

    S=DiagonalGroup([2 for i=1:length(p)])

    function logS(x::Array{Int, 1})
      return S([x[i] > 0 ? 0 : 1 for i=1:length(x)])
    end

    s = typeof(S[1])[]
    g = elem_type(O)[]
    u, mu = sub(S, s, false)
    b = 10
    cnt = 0
    while true
      @assert b > 0
      a = rand(m, b)
      if a==0
        continue
      end
      emb=signs(K(a),p)
      t=S([emb[x]==1 ? 0 : 1 for x in collect(keys(emb))])
      if !Hecke.haspreimage(mu, t)[1]
        push!(s, t)
        push!(g, a)
        u, mu = sub(S, s, false)
        if order(u) == order(S)
          break
        end
      else
        cnt += 1
        if cnt > 1000
          b *= 2
          cnt = 0
        end
        if b <= 0
          b = 10
          cnt = 0
          bas = lll_basis(O)
          while true
            @assert b>0
            a = rand(bas, 1:b)
            if a==0
              continue
            end
            emb=signs(a,p)
            t=S([emb[x]==1 ? 0 : 1 for x in collect(keys(emb))])
            if !Hecke.haspreimage(mu, t)[1]
              push!(s, t)
              push!(g, O(a, false))
              u, mu = sub(S, s, false)
              if order(u) == order(S)
                break
              end
            else
              cnt += 1
              if cnt > 1000
                b *= 2
                cnt = 0
              end
            end
          end
        end
      end
    end
    hS = Hecke.GrpAbFinGenMap(S, S, vcat([x.coeff for x in s]))   # Change of coordinates so that the canonical basis elements are mapped to the elements found above
    r = elem_type(O)[]
    for i=1:length(p)
      y = haspreimage(hS,S[i])[2]
      push!(r, prod([g[i]^Int(y[i]) for i=1:length(p)]))
    end

    function exp(A::GrpAbFinGenElem)

      s=O(m.gen_one)
      if iszero(A.coeff)
        return s
      end
      for i=1:length(p)
        if Int(A.coeff[1,i]) == 1
          s=s*r[i]
        end
      end
      return s
    end

    function log(B::nf_elem)
      emb=Hecke.signs(B,p)
      return S([emb[x]==1 ? 0 : 1 for x in collect(keys(emb))])
    end

    function log(B::FacElem{nf_elem})
      emb=Hecke.signs(B,p)
      return S([emb[x]==1 ? 0 : 1 for x in collect(keys(emb))])
    end
  return S, exp, log

end

#
#  Function that stores the principal generators element of the powers of the generators
#  in the class group map
#

function _assure_princ_gen(mC::MapClassGrp)

  if isdefined(mC, :princ_gens)
    return true
  end
  C=domain(mC)
  mC.princ_gens=Array{Tuple{FacElem{NfOrdIdl,NfOrdIdlSet}, FacElem{nf_elem, AnticNumberField}},1}(ngens(C))
  for i=1:ngens(C)
    I=FacElem(Dict(mC(C[i])=> fmpz(1)))
    pr=principal_gen_fac_elem(I^C.snf[i])
    mC.princ_gens[i]=(I,pr)
  end
  return true

end

#
#  Changes the exponential map of the class group so that the chosen representatives are coprime to the modulus
#

function _coprime_ideal(C::GrpAbFinGen, mC::Map, m::NfOrdIdl)

  O=parent(m).order
  K=nf(O)
  L=NfOrdIdl[]

  for i=1:ngens(C)
    a=mC(C[i])
    if iscoprime(a,m)
      push!(L,a)
    else
      J=inv(a)
      s=K(rand(J.num,5))//J.den  # Is the bound acceptable?
      I=s*a
      simplify(I)
      I = numerator(I)
      while !iscoprime(I,m)
        s=K(rand(J.num,5))//J.den
        I=s*a
        simplify(I)
        I = numerator(I)
      end
      push!(L,I)
    end
  end

  function exp(a::GrpAbFinGenElem)
    e=FacElem(Dict{NfOrdIdl,fmpz}(ideal(O,1) => fmpz(1)))
    for i=1:ngens(C)
      if Int(a.coeff[1,i])!= 0
        e*=FacElem(Dict(L[i] => a.coeff[1,i]))
      end
    end
    return e
  end

  return exp

end

function _elements_to_coprime_ideal(C::GrpAbFinGen, mC::Map, m::NfOrdIdl)

  O=parent(m).order
  K=nf(O)
  L=Array{NfOrdIdl,1}(ngens(C))
  el=Array{nf_elem,1}(ngens(C))

  for i=1:ngens(C)
    a=first(keys(mC.princ_gens[i][1].fac))
    if iscoprime(a,m)
      L[i]=a
      el[i]=K(1)
    else
      J=inv(a)
      s=K(rand(J.num,5))//J.den  # Is the bound acceptable?
      I=s*a
      simplify(I)
      I = numerator(I)
      while !iscoprime(I,m)
        s=K(rand(J.num,5))//J.den
        I=s*a
        simplify(I)
        I = numerator(I)
      end
      L[i]=I
      el[i]=s
    end
  end
  for i=1:ngens(C)
    @assert iscoprime(L[i],m)
  end

  function exp(a::GrpAbFinGenElem)
    e=Dict{NfOrdIdl,fmpz}()
    for i=1:ngens(C)
      if Int(a.coeff[1,i])!= 0
        e[L[i]]= a.coeff[1,i]
      end
    end
    if isempty(e)
      e[ideal(O,1)]=1
    end
    return FacElem(e)
  end

  return exp, el

end
function empty_ray_class(m::NfOrdIdl)
  O=order(parent(m))
  X=DiagonalGroup(Int[])
  function exp(a::GrpAbFinGenElem)
    return FacElem(Dict(ideal(O,1) => fmpz(1)))
  end

  function disclog(J::Union{NfOrdIdl, FacElem{NfOrdIdl}})
    return X(Int[])
  end

  mp=Hecke.MapRayClassGrp{typeof(X)}()
  mp.header = Hecke.MapHeader(X, FacElemMon(parent(m)) , exp, disclog)
  mp.modulus_fin=ideal(O,1)
  mp.modulus_inf=InfPlc[]

  return X,mp

end

function class_as_ray_class(C::GrpAbFinGen, mC::MapClassGrp, exp_class::Function,  m::NfOrdIdl, n::Integer)

  O=order(m)
  X,_=quo(C, n,false)
  function exp(a::GrpAbFinGenElem)
    return exp_class(a)
  end

  function disclog(J::NfOrdIdl)
    return X((mC\J).coeff)
  end

  function disclog(J::FacElem)
    a= X([0 for i=1:ngens(X)])
    for (f,k) in J.fac
      a+=k*disclog(f)
    end
    return a
  end

  mp=Hecke.MapRayClassGrp{typeof(X)}()
  mp.header = Hecke.MapHeader(X, FacElemMon(parent(m)) , exp, disclog)
  mp.modulus_fin=ideal(O,1)
  mp.modulus_inf=InfPlc[]
  mp.fact_mod=Dict{NfOrdIdl, Int}()

  return X,mp
end

function class_as_ray_class(C::GrpAbFinGen, mC::MapClassGrp, exp_class::Function,  m::NfOrdIdl)

    O=order(m)
    X=deepcopy(C)
    function exp(a::GrpAbFinGenElem)
      return exp_class(a)
    end

    function disclog(J::NfOrdIdl)
      return X((mC\J).coeff)
    end

    function disclog(J::FacElem)
      a= X([0 for i=1:ngens(X)])
      for (f,k) in J.fac
        a+=k*disclog(f)
      end
      return a
    end

    mp=Hecke.MapRayClassGrp{typeof(X)}()
    mp.header = Hecke.MapHeader(X, FacElemMon(parent(m)) , exp, disclog)
    mp.modulus_fin=ideal(O,1)
    mp.modulus_inf=InfPlc[]
    mp.fact_mod=Dict{NfOrdIdl, Int}()
    return X,mp

end

###################################################################################
#
#  Ray Class Group
#
###################################################################################

function ray_class_group_fac_elem(m::NfOrdIdl, inf_plc::Array{InfPlc,1}=InfPlc[])

#
# We compute the group using the sequence U -> (O/m)^* _> Cl^m -> Cl -> 1
# First of all, we compute all these groups with their own maps
#

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

  C, mC = class_group(O)

  exp_class=Hecke._coprime_ideal(C,mC,m)
  U, mU = unit_group_fac_elem(O)
  Q, pi= quo(O,m)
  G, mG=_multgrp_ray(Q)

  lp=Q.factor

  p = [ x for x in inf_plc if isreal(x) ]
  if !isempty(p)
    H,eH,lH=Hecke._infinite_primes(O,p,m)
    T=G
    G=direct_product(G,H)
  end

  @vprint :RayFacElem 1 "The multiplicative group is $G \n"
  @vprint :RayFacElem 1 "The class group is $C \n"
  @vprint :RayFacElem 1 "The units are $U \n"

  expon=exponent(G)

#
# We construct the relation matrix and evaluate units and relations with the class group in the quotient by m
# Then we compute the discrete logarithms
#

  R=zero_matrix(FlintZZ,ngens(C)+ngens(U)+ngens(G), ngens(C)+ngens(G))
  for i=1:ngens(C)
    R[i,i]=C.snf[i]
  end
  if issnf(G)
    for i=1:ngens(G)
      R[i+ngens(C),i+ngens(C)]=G.snf[i]
    end
  else
    for i=1:ngens(G)
      R[i+ngens(C),i+ngens(C)]=G.rels[i,i]
    end
  end

  @vprint :RayFacElem 1 "Collecting elements to be evaluated; first, units \n"
  evals=[]
  tobeeval=FacElem{nf_elem, AnticNumberField}[]
  if U.snf[1]==2
    push!(evals,O(-1))
  else
    push!(tobeeval, mU(U[1]))
  end
  append!(tobeeval,[mU(U[i]) for i=2:ngens(U)])

  @vprint :RayFacElem 1 "then principal ideal generators \n"
  princ_gens=[]
  for i=1:ngens(C)
    @vtime :RayFacElem 1 push!(princ_gens, Hecke.principal_gen_fac_elem((exp_class(C[i]))^(Int(order(C[i])))))
  end
  append!(tobeeval,princ_gens)

  @vprint :RayFacElem 1 "Time for elements evaluation: "
  @vtime :RayFacElem 1 ev,quots,idemps=fac_elems_eval(O,Q,tobeeval,lp,expon)
  append!(evals,ev)
  @vprint :RayFacElem 1 "\n"

  for i=1:ngens(U)
    @vprint :RayFacElem 1 "Disclog of unit $i \n"
    a=(mG\Q(evals[i])).coeff
    if !isempty(p)
      if i==1
        a=hcat(a, matrix(FlintZZ,1,length(p), [1 for i in p]))
      else
        b=lH(mU(U[i]))
        a=hcat(a, b.coeff)
      end
    end
    for j=1:ngens(G)
      R[i+ngens(G)+ngens(C),ngens(C)+j]=a[1,j]
    end
  end

#
# We compute the relation between generators of Cl and (O/m)^* in Cl^m
#

  for i=1: ngens(C)
    @vprint :RayFacElem 1 "Disclog of class group element $i \n"
    a=((mG\Q(evals[i+ngens(U)]))).coeff
    if !isempty(p)
      b=lH(princ_gens[i])
      a=hcat(a, b.coeff)
    end
    for j=1: ngens(G)
      R[i,ngens(C)+j]=-a[1,j]
    end
  end

  X=AbelianGroup(R)

#
# Discrete logarithm
#

  function disclog(J::FacElem)

    @vprint :RayFacElem 1 "Disc log of element $J \n"
    a= X([0 for i=1:ngens(X)])
    for (f,k) in J.fac
      a+=k*disclog(f)
    end
    return a
  end

  function disclog(J::NfOrdIdl)

    if isone(J)
    @vprint :RayFacElem 1 "J is one \n"
      return X([0 for i=1:ngens(X)])
    else
      L=mC\J
      @vprint :RayFacElem 1 "Disc log of element J in the Class Group: $(L.coeff) \n"
      s=exp_class(L)
      I=J* inv(s)
      @vprint :RayFacElem 1 "This ideal is principal: $I \n"
      z=principal_gen_fac_elem(I)
      el=_fac_elem_evaluation(O,Q,quots,idemps,z,lp,expon)
      @vprint :RayFacElem 1 "and 'generated' by $el \n"
      y=(mG\Q(el)).coeff
      @vprint :RayFacElem 1 "in the unit group, $y \n"
      if !isempty(p)
        b=lH(z)
        @vprint :RayFacElem 1 "the signs are $b \n"
        y=hcat(y, b.coeff)
      end
      return X(hcat(L.coeff,y))
    end
  end

#
# Exp map
#

  function expo(a::GrpAbFinGenElem)
    b=C([a.coeff[1,i] for i=1:ngens(C)])
    if isempty(p)
      c=G([a.coeff[1,i] for i=ngens(C)+1:ngens(X)])
      return exp_class(b)*ideal(O,pi\(mG(c)))
    else
      c=T([a.coeff[1,i] for i=ngens(C)+1:ngens(T)+ngens(C)])
      d=H([a.coeff[1,i] for i=ngens(T)+ngens(C)+1:ngens(X)])
      el=pi\(mG(c))
      @vprint :RayFacElem 1 "I have the element $el \n"
      @vprint :RayFacElem 1 "I want $(d.coeff) \n"
      # I need to modify $el$ so that it has the correct sign at the embeddings contained in primes
      vect=(lH(K(el))).coeff
      if vect==d.coeff
        return exp_class(b)*ideal(O,el)
      else
        correction=eH(d)
        while vect!=d.coeff
          el=el+correction
          vect=(lH(K(el))).coeff
        end
        return exp_class(b)*ideal(O,el)
      end
    end
  end

  mp=MapRayClassGrp{typeof(X)}()
  mp.header = Hecke.MapHeader(X, FacElemMon(parent(m)), expo, disclog)
  mp.modulus_fin=m
  mp.modulus_inf=p
  mp.fact_mod=Q.factor
  mp.tame_mult_grp=mG.tame
  mp.wild_mult_grp=mG.wild
  return X, mp

end

########################################################
#
#  Ray Class Group - p-part
#
########################################################

function ray_class_group_p_part(p::Integer, m::NfOrdIdl, inf_plc::Array{InfPlc,1}=InfPlc[])

  O=parent(m).order
  K=nf(O)
  Q,pi=quo(O,m)
  G, mG, lp = _mult_grp(Q,p)
  C, mC = class_group(O)
  _assure_princ_gen(mC)

  if p==2
    pr = [ x for x in inf_plc if isreal(x) ]
    if !isempty(pr)
      H,eH,lH=Hecke._infinite_primes(O,pr,m)
      T=G
      G =Hecke.direct_product(G,H)
    end
  end

  valclassp=Int(p^(valuation(order(C[ngens(C)]),p)))
  nonppartclass=div(order(C[ngens(C)]),valclassp)

  if valclassp==1 && order(G)==1
    return empty_ray_class(m)
  end

  C, mC, vect = _class_group_mod_n(C,mC,valclassp)
  U, mU = unit_group_fac_elem(O)
  exp_class, Kel = Hecke._elements_to_coprime_ideal(C,mC,m)

  if order(G)==1
    return class_as_ray_class(C,mC,exp_class,m)
  end

  n=ideal(O,1)
  for (q,vq) in lp
    n*=q^vq
  end

#
# We start to construct the relation matrix
#

  expo=exponent(G)
  inverse_d=gcdx(nonppartclass,expo)[2]

  R=zero_matrix(FlintZZ, ngens(C)+ngens(U)+ngens(G), ngens(C)+ngens(G))
  for i=1:ngens(C)
    R[i,i]=C.snf[i]
  end
  if issnf(G)
    for i=1:ngens(G)
      R[i+ngens(C),i+ngens(C)]=G.snf[i]
    end
  else
    for i=1:ngens(G)
      R[i+ngens(C),i+ngens(C)]=G.rels[i,i]
    end
  end

#
# We compute the relation matrix given by the image of the map U -> (O/m)^*
#

  @assert issnf(U)
  @vprint :RayFacElem 1 "Collecting elements to be evaluated; first, units \n"
  evals=NfOrdQuoRingElem[]
  tobeeval=FacElem{nf_elem, AnticNumberField}[]
  if gcd(U.snf[1],p)!=1
    if U.snf[1]==2
      push!(evals,Q(-1))
    else
      push!(tobeeval, mU(U[1]))
    end
  else
    push!(evals, Q(1))
  end
  append!(tobeeval,[mU(U[i]) for i=2:ngens(U)])

  @vprint :RayFacElem 1 "then principal ideal generators \n"
  for i=1:ngens(C)
    @vtime :RayFacElem 1 push!(tobeeval, mC.princ_gens[i][2]*Kel[i])
  end

  @vprint :RayFacElem 1 "Time for elements evaluation: "
  ev,quots,idemps=fac_elems_eval(O,Q,tobeeval,lp,expo)
  @vtime :RayFacElem 1 append!(evals, ev)
  @vprint :RayFacElem 1 "\n"

  for i=1:ngens(U)
    @vprint :RayFacElem 1 "Disclog of unit $i \n"
    a=(mG\(evals[i].elem)).coeff
    if p==2 && !isempty(pr)
      if i==1
        a=hcat(a, matrix(FlintZZ,1,length(pr), [1 for i in pr]))
      else
        b=lH(tobeeval[length(tobeeval)-ngens(C)-ngens(U)+i])
        a=hcat(a, b.coeff)
      end
    end
    for j=1:ngens(G)
      R[i+ngens(G)+ngens(C),ngens(C)+j]=a[1,j]
    end
  end

#
# We compute the relation between generators of Cl and (O/m)^* in Cl^m
#

  for i=1: ngens(C)
    @vprint :RayFacElem 1 "Disclog of class group element $i \n"
    invn=gcdx(vect[i], C.snf[i])[2]
    a=((mG\(evals[i+ngens(U)].elem))*invn).coeff
    if p==2 && !isempty(pr)
      b=lH(tobeeval[length(tobeeval)-ngens(C)+i])
      a=hcat(a, b.coeff)
    end
    for j=1: ngens(G)
      R[i,ngens(C)+j]=-a[1,j]
    end
  end

  X=AbelianGroup(R)

#
# Discrete logarithm
#

  function disclog(J::FacElem)
    a= X([0 for i=1:ngens(X)])
    for (f,k) in J.fac
      a+=k*disclog(f)
    end
    return a
  end

  function disclog(J::NfOrdIdl)
    if isone(J)
      return X([0 for i=1:ngens(X)])
    else
      L=mC\J
      s=exp_class(L)
      I=J* inv(s)
      I=I^Int(nonppartclass)
      z=principal_gen_fac_elem(I)
      el=Hecke._fac_elem_evaluation(O, Q, quots, idemps, z, lp, expo)
      y=((mG\el)*inverse_d).coeff
      if p==2 && !isempty(pr)
        b=lH(z)
        y=hcat(y, b.coeff)
      end
      return X(hcat(L.coeff,y))
    end
  end

#
# Exp map
#

  function expon(a::GrpAbFinGenElem)
    b=C([a.coeff[1,i] for i=1:ngens(C)])
    if p!=2 || isempty(pr)
      c=G([a.coeff[1,i] for i=ngens(C)+1:ngens(X)])
      return exp_class(b)*ideal(O,pi\(mG(c)))
    else
      c=T([a.coeff[1,i] for i=ngens(C)+1:ngens(T)+ngens(C)])
      d=H([a.coeff[1,i] for i=ngens(T)+ngens(C)+1: ngens(X)])
      el=pi\(mG(c))
      # I need to modify $el$ so that it has the correct sign at the embeddings contained in primes
      vect=(lH(K(el))).coeff
      if vect==d.coeff
        return exp_class(b)*ideal(O,el)
      else
        correction=eH(d)
        while vect!=d.coeff
          el=el+correction
          vect=(lH(K(el))).coeff
        end
        return exp_class(b)*ideal(O,el)
      end
    end
  end

  mp=Hecke.MapRayClassGrp{typeof(X)}()
  mp.header = Hecke.MapHeader(X, FacElemMon(parent(m)) , expon, disclog)
  mp.modulus_fin=n
  mp.modulus_inf=inf_plc
  mp.fact_mod=lp
  mp.tame_mult_grp=mG.tame
  mp.wild_mult_grp=mG.wild

  return X,mp
end

#####################################################################################################
#
#  Quotient by n of the Ray Class Group
#
#####################################################################################################

function _class_group_mod_n(C::GrpAbFinGen, mC::Hecke.MapClassGrp, n::Integer)

  @assert issnf(C)
  O=parent(mC(C[1])).order
  K=nf(O)
  if gcd(C.snf[ngens(C)],n)==1
    G=DiagonalGroup(Int[])
    function exp1(a::GrpAbFinGenElem)
      return ideal(O, O(1))
    end
    function disclog1(I::NfOrdIdl)
      return G(Int[])
    end
    mp=Hecke.MapClassGrp{typeof(G)}()
    mp.header=Hecke.MapHeader(G, mC.header.codomain,exp1,disclog1)
    mp.princ_gens=[(FacElem(Dict(ideal(O,1)=> fmpz(1))), FacElem(Dict(K(1)=> 1)))]
    return G,mp, fmpz[]

  else

    ind=1
    while gcd(order(C[ind]),n)==1
      ind+=1
    end

    vect=[gcd(C.snf[ind+j],n) for j=0:ngens(C)-ind]
    G=DiagonalGroup(vect)
    G.issnf=true
    G.snf=vect

    function exp2(a::GrpAbFinGenElem)
      x=ideal(O,1)
      for i=1:ngens(G)
        if a[i]!=0
          x*=numerator(evaluate(mC.princ_gens[ind+i-1][1]))^(Int(a[i]))
        end
      end
      return x
    end

    function disclog2(I::NfOrdIdl)
      x=mC\I
      y=G([0 for j=1:ngens(G)])
      for i=ind:ngens(C)
          y.coeff[1,i-ind+1]=x.coeff[1,i]
      end
      return y
    end

    mp=Hecke.MapClassGrp{typeof(G)}()
    mp.header=Hecke.MapHeader(G, mC.header.codomain, exp2, disclog2)
    mp.princ_gens=mC.princ_gens[ind:end]

    return G,mp, [divexact(C.snf[ind+j],gcd(C.snf[ind+j],n)) for j=0:ngens(C)-ind]
  end
end

function ray_class_group_quo(n::Integer, m::NfOrdIdl, inf_plc::Array{InfPlc,1}=InfPlc[])

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

  #
  #  Take the relevant part of the modulus
  #

  fac=factor(m)
  y1=Dict{NfOrdIdl,Int}()
  y2=Dict{NfOrdIdl,Int}()
  for (q,e) in fac
    if gcd(norm(q)-1,n)!=1
      y1[q]=Int(1)
      if gcd(norm(q),n)!=1 && e>=2
        y2[q]=Int(e)
      end
    elseif gcd(norm(q),n)!=1 && e>=2
      y2[q]=Int(e)
    end
  end
  return ray_class_group_quo(n, m, y1, y2, inf_plc)

end

function ray_class_group_quo(O::NfOrd, n_quo::Int, m::Int, wprimes::Dict{NfOrdIdl,Int}=Dict{NfOrdIdl, Int}(), inf_plc::Array{InfPlc,1}=InfPlc[])

  K=nf(O)
  d1=Dict{NfOrdIdl, Int}()
  lp=factor(m)
  for q in keys(lp.fac)
    lq=prime_decomposition(O,q)
    for (P,e) in lq
      d1[P]=1
    end
  end
  return ray_class_group_quo(n_quo, ideal(O,1), d1, wprimes, inf_plc, check_expo=true)

end

function ray_class_group_quo(n::Integer, m::NfOrdIdl, y1::Dict{NfOrdIdl,Int}, y2::Dict{NfOrdIdl,Int}, inf_plc::Array{InfPlc,1}=InfPlc[]; check_expo=false)
  # check_expo checks, before the computation of the units, if the exponent of the group can be n.
  # if it is lower for sure, it returns the trivial group.
  # I HAVE TO FIND A BETTER METHOD.
  O=parent(m).order
  K=nf(O)

  # Compute the modulus of the quotient
  I=ideal(O,1)
  for (q,vq) in y1
    I*=q^vq
  end
  for (q,vq) in y2
    I*=q^vq
  end
  lp=merge(max,y1,y2)

  Q,pi=quo(O,I)
  Q.factor=lp
  C, mC = class_group(O)
  _assure_princ_gen(mC)
  @vtime :RayFacElem 1 G, mG, tame, wild= _mult_grp_mod_n(Q,y1,y2,n)

  if mod(n,2)==0
    pr = [ x for x in inf_plc if isreal(x) ]
    if !isempty(pr)
      @vtime :RayFacElem 1 H,eH,lH=Hecke._infinite_primes(O,pr,I)
      T=G
      G =Hecke.direct_product(G,H)
    end
  end

  if gcd(C.snf[end],n)==1 && order(G)==1
    return empty_ray_class(m)
  end

  f=collect(keys(factor(fmpz(n)).fac))
  val=Array{Int,1}(length(f))
  for i=1:length(f)
    val[i]=valuation(C.snf[end],f[i])
  end
  valclass=1
  for i=1:length(f)
    if val[i]!=0
      valclass*=f[i]^(val[i])
    end
  end
  nonnclass=divexact(C.snf[end], valclass)

  C, mC, vect = _class_group_mod_n(C,mC,Int(valclass))

  if check_expo && exponent(C)*exponent(G)<n
    return empty_ray_class(m)
  end

  U, mU = unit_group_fac_elem(O)
  exp_class, Kel = Hecke._elements_to_coprime_ideal(C,mC,m)
  for i=1:ngens(C)
    @hassert :RayFacElem 1 iscoprime(numerator(evaluate(exp_class(C[i]))), m)
  end

  if order(G)==1
    return class_as_ray_class(C,mC,exp_class,m,n)
  end

#
# We start to construct the relation matrix
#

  expo=exponent(G)

  R=zero_matrix(FlintZZ, 2*ngens(C)+ngens(U)+2*ngens(G), ngens(C)+ngens(G))
  for i=1:ngens(C)
    R[i,i]=C.snf[i]
  end
  if issnf(G)
    for i=1:ngens(G)
      R[i+ngens(C),i+ngens(C)]=G.snf[i]
    end
  else
    for i=1:ngens(G)
      R[i+ngens(C),i+ngens(C)]=G.rels[i,i]
    end
  end
  for i=1:cols(R)
    R[ngens(C)+ngens(U)+ngens(G)+i,i]=n
  end

#
# We compute the relation matrix given by the image of the map U -> (O/m)^*
#

  @hassert :RayFacElem 1 issnf(U)
  @vprint :RayFacElem 1 "Collecting elements to be evaluated; first, units \n"
  evals=NfOrdQuoRingElem[]
  tobeeval=FacElem{nf_elem, AnticNumberField}[]
  if gcd(U.snf[1],n)!=1
    if U.snf[1]==2
      push!(evals,Q(-1))
    else
      push!(tobeeval, mU(U[1]))
    end
  else
    push!(evals,Q(1))
  end
  append!(tobeeval,[mU(U[i]) for i=2:ngens(U)])

  @vprint :RayFacElem 1 "then principal ideal generators \n"
  for i=1:ngens(C)
    push!(tobeeval, mC.princ_gens[i][2]*(Kel[i]^(C.snf[i]*vect[i])))
  end

  @vprint :RayFacElem 1 "Time for elements evaluation: "
  @vtime :RayFacElem 1 ev,quots,idemps=fac_elems_eval(O,Q,tobeeval,lp,gcd(expo,n))
  append!(evals, ev)
  @vprint :RayFacElem 1 "\n"

  for i=1:ngens(U)
    @vprint :RayFacElem 1 "Disclog of unit $i \n"
    a=(mG\(evals[i].elem)).coeff
    if mod(n,2)==0 && !isempty(pr)
      if i==1
        a=hcat(a, matrix(FlintZZ,1,length(pr), [1 for i in pr]))
      else
        b=lH(mU(U[i]))
        a=hcat(a, b.coeff)
      end
    end
    for j=1:ngens(G)
      R[i+ngens(G)+ngens(C),ngens(C)+j]=a[1,j]
    end
  end

#
# We compute the relation between generators of Cl and (O/m)^* in Cl^m
#

  for i=1:ngens(C)
    @vprint :RayFacElem 1 "Disclog of class group element $i \n"
    invn=invmod(vect[i],fmpz(expo))
    investigated=evaluate(mC.princ_gens[i][2]*(Kel[i]^(C.snf[i]*vect[i])))
    a=((mG\(evals[i+ngens(U)].elem))*invn).coeff
    if mod(n,2)==0 && !isempty(pr)
      b=lH(mC.princ_gens[i][2]*Kel[i])
      a=hcat(a, b.coeff)
    end
    for j=1: ngens(G)
      R[i,ngens(C)+j]=-a[1,j]
    end
  end

  X=AbelianGroup(R)

#
# Discrete logarithm
#
  inverse_d=invmod(fmpz(nonnclass),fmpz(expo))
  @assert gcd(fmpz(nonnclass),fmpz(expo))==1

  function disclog(J::FacElem{NfOrdIdl, NfOrdIdlSet})

    a= C([0 for i=1:ngens(C)])
    for (f,k) in J.fac
      a+=k*(mC\f)
    end
    Id=J* inv(exp_class(a))
    Id=Id^Int(nonnclass)
    z=principal_gen_fac_elem(Id)
    el=Hecke._fac_elem_evaluation(O, Q, quots, idemps, z, lp, gcd(expo,n))
    y=((mG\(el))*inverse_d).coeff
    if mod(n,2)==0 && !isempty(pr)
      b=lH(z)
      y=hcat(y, b.coeff)
    end
    return X(hcat(a.coeff,y))
  end

  function disclog(J::NfOrdIdl)

    @hassert :RayFacElem 1 iscoprime(J,I)
    if J.is_principal==1
      if isdefined(J,:princ_gen)
        el=J.princ_gen
        y=(mG\(el)).coeff
        if mod(n,2)==0 && !isempty(pr)
          b=lH(K(el))
          y=hcat(y, b.coeff)
        end
        return X(hcat(C([0 for i=1:ngens(C)]).coeff,y))
      elseif isdefined(J,:princ_gen_special)
        el=O(J.princ_gen_special[2])+O(J.princ_gen_special[3])
        y=(mG\(el)).coeff
        if mod(n,2)==0 && !isempty(pr)
          b=lH(K(el))
          y=hcat(y, b.coeff)
        end
        return X(hcat(C([0 for i=1:ngens(C)]).coeff,y))
      else
        z=principal_gen_fac_elem(J)
        el=Hecke._fac_elem_evaluation(O, Q, quots, idemps, z, lp, gcd(expo,n))
        y=(mG\(el)).coeff
        if mod(n,2)==0 && !isempty(pr)
          b=lH(z)
          y=hcat(y, b.coeff)
        end
        return X(hcat(C([0 for i=1:ngens(C)]).coeff,y))
      end
    else
      W=mC\J
      s=exp_class(W)
      for (el,v) in s.fac
        s.fac[el]=-nonnclass*v
      end
      if haskey(s.fac, J)
        s.fac[J]+=nonnclass
      else
        s.fac[J]=nonnclass
      end
      z=principal_gen_fac_elem(s)
      el=Hecke._fac_elem_evaluation(O, Q, quots, idemps, z, lp, gcd(expo,n))
      y=((mG\(el))*inverse_d).coeff
      if mod(n,2)==0 && !isempty(pr)
        b=lH(z)
        y=hcat(y, b.coeff)
      end
      return X(hcat(W.coeff,y))
    end

  end

#
# Exp map
#

  function expon(a::GrpAbFinGenElem)
    b=C([a.coeff[1,i] for i=1:ngens(C)])
    if mod(n,2)!=0  || isempty(pr)
      c=G([a.coeff[1,i] for i=ngens(C)+1:ngens(X)])
      return exp_class(b)*ideal(O,mG(c).elem)
    else
      c=T([a.coeff[1,i] for i=ngens(C)+1:ngens(T)+ngens(C)])
      d=H([a.coeff[1,i] for i=ngens(T)+ngens(C)+1: ngens(X)])
      el=mG(c).elem
      # I need to modify $el$ so that it has the correct sign at the embeddings contained in primes
      vect=(lH(K(el))).coeff
      if vect==d.coeff
        return exp_class(b)*ideal(O,el)
      else
        correction=eH(d)
        while vect!=d.coeff
          el=el+correction
          vect=(lH(K(el))).coeff
        end
        return exp_class(b)*ideal(O,el)
      end
    end
  end

  mp=Hecke.MapRayClassGrp{typeof(X)}()
  mp.header = Hecke.MapHeader(X, FacElemMon(parent(m)) , expon, disclog)
  mp.modulus_fin=I
  mp.evals=evals
  mp.quots=quots
  mp.idemps=idemps
  mp.coprime_elems=Kel
  mp.fact_mod=lp
  mp.tame_mult_grp=tame
  mp.wild_mult_grp=wild

  if mod(n,2)==0
    mp.modulus_inf=pr
  else
    mp.modulus_inf=inf_plc
  end
  return X,mp

end

##################################################################################
#
#  Ray Class Group over QQ
#
##################################################################################

function ray_class_groupQQ(O::NfOrd, modulus::Int, inf_plc::Bool, n_quo::Int)

  R=ResidueRing(FlintZZ, modulus, cached=false)
  U,mU=unit_group_mod(R, n_quo)
  if inf_plc
    function disc_log1(I::NfOrdIdl)
      @assert gcd(minimum(I),modulus)==1
      i=Int(I.minimum)
      return mU\(R(i))
    end

    function expon1(a::GrpAbFinGenElem)
      x=mU(a)
      return FacElem(Dict{NfOrdIdl, fmpz}(ideal(O,lift(x)) => 1))
    end

    mp=Hecke.MapRayClassGrp{typeof(U)}()
    mp.header = Hecke.MapHeader(U, FacElemMon(parent(ideal(O,1))) , expon1, disc_log1)
    mp.modulus_fin=ideal(O,modulus)
    mp.modulus_inf=real_places(nf(O))

    return U,mp

  elseif isodd(n_quo)

    function disc_log2(I::NfOrdIdl)
      @assert gcd(minimum(I),modulus)==1
      i=Int(I.minimum)
      return mU\(R(i))
    end

    function expon2(a::GrpAbFinGenElem)
      x=mU(a)
      return FacElem(Dict{NfOrdIdl, fmpz}(ideal(O,lift(x)) => 1))
    end

    mp=Hecke.MapRayClassGrp{typeof(U)}()
    mp.header = Hecke.MapHeader(U, FacElemMon(parent(ideal(O,1))) , expon2, disc_log2)
    mp.modulus_fin=ideal(O,modulus)
    mp.modulus_inf=[]

    return U,mp

  else

    Q,mQ=quo(U, [mU\(R(-1))])

    function disc_log(I::NfOrdIdl)
      i=Int(minimum(I))
      return mQ(mU\(R(i)))
    end

    function expon(a::GrpAbFinGenElem)
      x=mU(mQ\a)
      return FacElem(Dict{NfOrdIdl, fmpz}(ideal(O,x) => 1))
    end

    mp=Hecke.MapRayClassGrp{typeof(Q)}()
    mp.header = Hecke.MapHeader(Q, FacElemMon(parent(ideal(O,1))) , expon, disc_log)
    mp.modulus_fin=ideal(O,modulus)
    mp.modulus_inf=[]

    return Q,mp

  end

end

##################################################################################
#
#  Action of the Galois Group on the Ray Class Group
#
##################################################################################

function _aut_on_id(O::NfOrd, phi::Hecke.NfToNfMor, I::NfOrdIdl)

  K=nf(O)
  if I.is_principal==1
    if isdefined(I, :princ_gen)
      y=K(I.princ_gen)
      y=O(phi(y), false)
      return ideal(O,y)
    else
      y=K(I.gen_two)
      y=O(phi(y), false)
      J=ideal(O,I.gen_one,y)
      J.is_principal=1
      return J
    end
  else
    y=K(I.gen_two)
    y=O(phi(y), false)
    return ideal(O,I.gen_one,y)
  end
end

#
#  Find small primes that generate the ray class group (or a quotient)
#  It needs a map GrpAbFinGen -> NfOrdIdlSet
#
function find_gens(mR::MapRayClassGrp)

  O = order(codomain(mR))
  R = domain(mR)
  m = mR.modulus_fin
  mm = minimum(m)

  sR = GrpAbFinGenElem[]
  lp = NfOrdIdl[]
  q, mq = quo(R, sR,false)

  #
  #  First, generators of the multiplicative group.
  #  If the class group is trivial, they are enough
  #

  if !isempty(mR.fact_mod)
    totally_positive_generators(mR, Int(mm), true)
    tmg=mR.tame_mult_grp
    wld=mR.wild_mult_grp
    for (p,v) in tmg
      I=ideal(O,v[1])
      f=mR\I
      if iszero(mq(f))
        continue
      end
      push!(sR, f)
      push!(lp, I)
      q, mq = quo(R, sR, false)
      if order(q) == 1
        return lp,sR
      end
    end

    for (p,v) in wld
      for i=1:length(v[1])
        I=ideal(O,v[1][i])
        f=mR\I
        if iszero(mq(f))
          continue
        end
        push!(sR, f)
        push!(lp, I)
        q, mq = quo(R, sR, false)
        if order(q) == 1
          return lp, sR
        end
      end
    end
  end

  if !isempty(mR.modulus_inf)
    S, ex, lo=carlos_units(O)
    for i=1:length(mR.modulus_inf)
      pl=mR.modulus_inf[i]
      @assert isreal(pl)
      delta=mm*ex(S[i])
      el=1+delta
      con=abs_upper_bound(1/real(conjugates_arb(delta))[i], fmpz)
      el+=con*delta
      I=ideal(O,el)
      f=mR\I
      if iszero(mq(f))
        continue
      end
      push!(sR, f)
      push!(lp, I)
      q, mq = quo(R, sR, false)
      if order(q)==1
        return lp, sR
      end
    end

  end

  if isdefined(mR, :prime_ideal_cache)
    S = mR.prime_ideal_cache
  else
    S = prime_ideals_up_to(O, max(1000,100*clog(discriminant(O),10)^2))
    mR.prime_ideal_cache = S
  end
  q, mq = quo(R, sR, false)
  for P in S
    if !iscoprime(P,m)
      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
    if iszero(mq(f))
      continue
    end
    push!(sR, f)
    push!(lp, P)
    q, mq = quo(R, sR, false)
    if order(q) == 1
      break
    end
  end
  @assert order(q)==1
  return lp, sR
end

function _act_on_ray_class(mR::MapRayClassGrp, Aut::Array{Hecke.NfToNfMor,1}=Array{Hecke.NfToNfMor,1}())

  R=mR.header.domain
  O=mR.header.codomain.base_ring.order
  K=nf(O)
  if isempty(Aut)
    Aut=automorphisms(K)
  end
  if ngens(R)==0
    return GrpAbFinGenMap[]
  end

  lgens,subs=find_gens(mR)
  if isempty(lgens)
    return GrpAbFinGenMap[]
  end

  G=Array{GrpAbFinGenMap,1}(length(Aut))
  #
  #  Instead of applying the automorphisms to the elements given by mR, I choose small primes
  #  generating the group and study the action on them. In this way, I take advantage of the cache of the
  #  class group map
  #

  lgens,subs=find_gens(mR)

  if isempty(lgens)
    push!(G, GrpAbFinGenMap(R))
    return G
  end

  #
  #  Write the matrices for the change of basis
  #
  auxmat=zero_matrix(FlintZZ, ngens(R), length(lgens)+nrels(R))
  for i=1:length(lgens)
    for j=1:ngens(R)
      auxmat[j,i]=subs[i][j]
    end
  end
  if issnf(R)
    for i=1:ngens(R)
      auxmat[i,length(lgens)+i]=R.snf[i]
    end
  else
    for i=1:ngens(R)
      for j=1:nrels(R)
        auxmat[i,length(lgens)+j]=R.rels[j,i]
      end
    end
  end

  Ml=solve(auxmat,eye(auxmat,ngens(R)))'
  #
  #  Now, we compute the action on the group
  #

  for k=1:length(Aut)
    M=zero_matrix(FlintZZ,length(lgens), ngens(R))
    for i=1:length(lgens)
      @vtime :RayFacElem 3 J=_aut_on_id(O,Aut[k],lgens[i])
      @vtime :RayFacElem 3 elem=mR\J
      for j=1:ngens(R)
        M[i,j]=elem[j]
      end
    end
    G[k]= hom(R, R, sub(Ml,1:rows(Ml), 1:length(lgens))*M)
    @hassert :RayFacElem 1 isbijective(G[k])
  end

  return G

end

##################################################################################
#
#  Stable Subgroups function
#
##################################################################################

doc"""
***
    stable_subgroups(R::GrpAbFinGen, quotype::Array{Int,1}, act::Array{T, 1}; op=sub)

> Given a group R, an array of endomorphisms of the group and the type of the quotient, it returns all the stable
> subgroups of R such that the corresponding quotient has the required type.
"""

function stable_subgroups(R::GrpAbFinGen, quotype::Array{Int,1}, act::Array{T, 1}; op=sub) where T <: Map{GrpAbFinGen, GrpAbFinGen}

  c=lcm(quotype)
  Q,mQ=quo(R,c, false)
  lf=factor(order(Q)).fac
  list=[]
  for p in keys(lf)

    x=valuation(c,p)
    if x==0
      continue
    end
    G,mG=psylow_subgroup(Q, p, false)
    S,mS=snf(G)
    quotype_p=Int[]
    for i=1:length(quotype)
      v=valuation(quotype[i],p)
      if v>0
        push!(quotype_p,v)
      end
    end
    if Int(p)*quotype_p==S.snf
      plist1=GrpAbFinGenElem[c*R[i] for i=1:ngens(R)]
      push!(list, ([plist1]))
      continue
    end
    #
    #  Action on the group: we need to distinguish between FqGModule and ZpnGModule (in the first case the algorithm is more efficient)
    #

    if x==1

      F, _ = Nemo.FiniteField(Int(p), 1, "_", cached=false)
      act_mat=Array{fq_nmod_mat, 1}(length(act))
      for w=1:length(act)
        act_mat[w]=zero_matrix(F,ngens(S), ngens(S))
      end
      for w=1:ngens(S)
        el=mG(mS(S[w]))
        for z=1:length(act)
          elz=mS\(haspreimage(mG,act[z](el))[2])
          for l=1:ngens(S)
            act_mat[z][w,l]=elz[l]
          end
        end
      end
      M=FqGModule(act_mat)

      #
      #  Searching for submodules
      #

      ind=length(quotype_p)
      plist=submodules(M,ind)
      push!(list, (_lift_and_construct(x, mQ,mG,mS,c) for x in plist))

    else

      RR=ResidueRing(FlintZZ, Int(p^x), cached=false)
      act_mat=Array{nmod_mat,1}(length(act))
      auxmat1=hcat(mG.map', rels(Q)')
      auxmat2=mS.map*mG.map
      W=MatrixSpace(RR,ngens(S), ngens(S), false)
      for z=1:length(act)
        y=transpose(solve(auxmat1, (auxmat2*act[z].map)'))
        y=sub(y,1:ngens(S), 1:ngens(G))*mS.imap
        act_mat[z]=W(y)
      end

      #
      #  Searching for submodules
      #

      M=Hecke.ZpnGModule(S,act_mat)
      plist=submodules(M,typequo=quotype_p)
      push!(list, (_lift_and_construct(x, mQ,mG,mS,c) for x in plist))

    end
  end
  if isempty(list)
    return ([])
  end

  return ( op(R,vcat(c...)) for c in Iterators.product(list...))

end

function _lift_and_construct(A::nmod_mat, mQ::GrpAbFinGenMap, mG::GrpAbFinGenMap, mS::GrpAbFinGenMap, c::Int)

  R=mQ.header.domain
  newsub=GrpAbFinGenElem[c*R[i] for i=1:ngens(R)]
  for i=1:rows(A)
    y=sub(A,i:i,1:cols(A))
    if !iszero(y)
      push!(newsub,mQ\(mG(mS(mS.header.domain(lift(y))))))
    end
  end
  return newsub

end

function _lift_and_construct(A::fq_nmod_mat, mQ::GrpAbFinGenMap, mG::GrpAbFinGenMap, mS::GrpAbFinGenMap, c ::Int)

  R=mQ.header.domain
  newsub=GrpAbFinGenElem[c*R[i] for i=1:ngens(R)]
  for i=1:rows(A)
    z=zero_matrix(FlintZZ,1,cols(A))
    for j=1:cols(z)
      z[1,j]=FlintZZ(coeff(A[i,j],0))
    end
    push!(newsub,mQ\(mG(mS(mS.header.domain(z)))))
  end
  return newsub

end


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    CopyMapper

A mapping strategy in which data points are copied directly to the
domain at specified locations.
"""
struct CopyMapper{V1,V2} <: AbstractMapper
  orig::V1
  dest::V2
end

CopyMapper(dest) = CopyMapper(nothing, dest)
CopyMapper() = CopyMapper(nothing, nothing)

function Base.map(spatialdata::AbstractData{T,N},
                  domain::AbstractDomain{T,N},
                  targetvars::NTuple{K,Symbol},
                  mapper::CopyMapper) where {N,T,K}
  @assert targetvars ⊆ keys(variables(spatialdata)) "target variables must be present in spatial data"

  # dictionary with mappings
  mappings = Dict(var => Dict{Int,Int}() for var in targetvars)

  # retrieve origin and destination indices
  orig = mapper.orig ≠ nothing ? mapper.orig : 1:npoints(spatialdata)
  dest = mapper.dest ≠ nothing ? mapper.dest : 1:npoints(spatialdata)

  @assert length(orig) == length(dest) "invalid mapping specification"

  for i in eachindex(orig, dest)
    # save pair if there is data for variable
    for var in targetvars
      if isvalid(spatialdata, orig[i], var)
        push!(mappings[var], dest[i] => orig[i])
      end
    end
  end

  mappings
end


text:
#=

@author : Victoria Gregory, John Stachurski

=#

struct LakeModel{TF <: AbstractFloat}
    λ::TF
    α::TF
    b::TF
    d::TF
    g::TF
    A::Matrix{TF}
    A_hat::Matrix{TF}
end

"""
Constructor with default values for `LakeModel`

##### Fields of `LakeModel`

  - λ : job finding rate
  - α : dismissal rate
  - b : entry rate into labor force
  - d : exit rate from labor force
  - g : net entry rate
  - A : updates stock
  - A_hat : updates rate

"""
function LakeModel(;λ::AbstractFloat=0.283,
                    α::AbstractFloat=0.013,
                    b::AbstractFloat=0.0124,
                    d::AbstractFloat=0.00822)

    g = b - d
    A = [(1-λ) * (1-d) + b  (1-d) * α + b;
         (1-d) * λ          (1-d) * (1-α)]
    A_hat = A ./ (1 + g)

    return LakeModel(λ, α, b, d, g, A, A_hat)
end

"""
Finds the steady state of the system :math:`x_{t+1} = \hat A x_{t}`

##### Arguments

  - lm : instance of `LakeModel`
  - tol: convergence tolerance

##### Returns

  - x : steady state vector of employment and unemployment rates

"""

function rate_steady_state(lm::LakeModel, tol::AbstractFloat=1e-6)
    x = 0.5 * ones(2)
    error = tol + 1
    while (error > tol)
        new_x = lm.A_hat * x
        error = maximum(abs, new_x - x)
        x = new_x
    end
    return x
end

"""
Simulates the the sequence of Employment and Unemployent stocks

##### Arguments

  - X0 : contains initial values (E0, U0)
  - T : number of periods to simulate

##### Returns

  - X_path : contains sequence of employment and unemployment stocks

"""

function simulate_stock_path{TF<:AbstractFloat}(
                    lm::LakeModel, X0::AbstractVector{TF}, T::Integer)
    X_path = Array{TF}(2, T)
    X = copy(X0)
    for t in 1:T
        X_path[:, t] = X
        X = lm.A * X
    end
    return X_path
end

"""
Simulates the the sequence of employment and unemployent rates.

##### Arguments

  - X0 : contains initial values (E0, U0)
  - T : number of periods to simulate

##### Returns

  - X_path : contains sequence of employment and unemployment rates

"""
function simulate_rate_path{TF<:AbstractFloat}(lm::LakeModel,
                                               x0::Vector{TF}, T::Integer)
    x_path = Array{TF}(2, T)
    x = copy(x0)
    for t in 1:T
        x_path[:, t] = x
        x = lm.A_hat * x
    end
    return x_path
end


text:
# Main GaussianProcess type

abstract type GPBase end

"""
    predict_f(gp::GPBase, X::Matrix{Float64}[; full_cov::Bool = false])

Return posterior mean and variance of the Gaussian Process `gp` at specfic points which are
given as columns of matrix `X`. If `full_cov` is `true`, the full covariance matrix is
returned instead of only variances.
"""
function predict_f(gp::GPBase, x::MatF64; full_cov::Bool=false)
    size(x,1) == gp.dim || throw(ArgumentError("Gaussian Process object and input observations do not have consistent dimensions"))
    if full_cov
        return _predict(gp, x)
    else
        ## Calculate prediction for each point independently
            μ = Array{Float64}(undef, size(x,2))
            σ2 = similar(μ)
        for k in 1:size(x,2)
            m, sig = _predict(gp, x[:,k:k])
            μ[k] = m[1]
            σ2[k] = max(diag(sig)[1], 0.0)
        end
        return μ, σ2
    end
end

# 1D Case for prediction
predict_f(gp::GPBase, x::VecF64; full_cov::Bool=false) = predict_f(gp, x'; full_cov=full_cov)

function make_posdef!(m::MatF64)
    n = size(m, 1)
    size(m, 2) == n || throw(ArgumentError("Covariance matrix must be square"))
    for _ in 1:10 # 10 chances
        try
            return m, cholesky(m)
        catch
            # that wasn't (numerically) positive definite,
            # so let's add some weight to the diagonal
            ϵ = 1e-6 * tr(m) / n
            @inbounds for i in 1:n
                m[i, i] += ϵ
            end
        end
    end
    m, cholesky(m)
end

"""
    tolerant_PDMat(Σ::Matrix{Float64})

Try to encode covariance matrix `Σ` as positive definite matrix of type `PDMat`.

Sometimes covariance matrices of Gaussian processes are positive definite mathematically
but have negative eigenvalues numerically. To resolve this issue, small weights are added
to the diagonal (and hereby all eigenvalues are raised by that amount mathematically)
until all eigenvalues are positive numerically.
"""
function tolerant_PDMat(Sigma_raw::MatF64)
    Sigma_raw, chol = make_posdef!(Sigma_raw)
    PDMat(Sigma_raw, chol)
end


text:
# BBoxArray

export BBoxArray
struct BBoxArray{D} <: BBoxSet{D}
    bboxes::Vector{BBox{D}}
    BBoxArray{D}() where {D} = new(BBox{D}[])
    BBoxArray{D}(bboxes::Vector{BBox{D}}) where {D} = new(bboxes)
end

BBoxArray(bboxes::Vector{BBox{D}}) where {D} = BBoxArray{D}(bboxes)

function invariant(bbarray::BBoxArray)
    any(isempty, bbarray.bboxes) && return false
    for i in 1:length(bbarray.bboxes), j in (i + 1):length(bbarray.bboxes)
        intersects(bbarray.bboxes[i], bbarray.bboxes[j]) && return false
    end
    return true
end

Base.ndims(bbarray::BBoxArray{D}) where {D} = D

Base.eltype(bbarray::BBoxArray) = Int

Base.isempty(bbarray::BBoxArray) = isempty(bbarray.bboxes)

Base.length(bbarray::BBoxArray) = sum(length, bbarray.bboxes)

function shift!(bbarray::BBoxArray{D}, offset::SVector{D,Int}) where {D}
    for i in 1:length(bbarray.bboxes)
        bbarray.bboxes[i] = shift(bbarray.bboxes[i], offset)
    end
    return bbarray
end

function shift(bbarray::BBoxArray{D}, offset::SVector{D,Int}) where {D}
    return shift!(BBoxArray(copy(bbarray.bboxes)))
end

function Base.in(pos::SVector{D,Int}, bbarray::BBoxArray{D}) where {D}
    return any(bbox -> pos in bbox, bbarray.bboxes)
end

# in, intersect, union, setdiff, all with BBox and BBoxArray


text:
module Calc

using Test
using DebuggerFramework
using AbstractTrees

include("utilities.jl")

# This file implementes a simple calculator. It is meant to illustrate
# and test line-wise stepping, display of local variables and
# language-specic REPL prompts.

@enum(TokenKind,
    # Punctuation
    LPAREN,
    RPAREN,
    # Operators
    PLUS,
    MINUS,
    TIMES,
    DIV,
    # Literals
    INT,
    # Backreferences
    BACKREF,
    # Meta
    WHITESPACE,
    NEWLINE,
    ENDMARKER,
    ERROR,
)
is_op_kind(kind::TokenKind) = kind in (PLUS, MINUS, TIMES, DIV)

# Lexer
struct Token
    kind::TokenKind
    startbyte::UInt
    endbyte::UInt
    val::Union{Int, Nothing}
end

mutable struct Lexer
    io::IO
    current_char::Char
    charstore::IOBuffer
    token_start_pos::UInt
    current_row::Int
    current_col::Int
end
Lexer(io::IO) = Lexer(io, '\0', IOBuffer(), 0, 0, 0)
peekchar(l::Lexer) = peekchar(l.io)
dpeekchar(l::Lexer) = dpeekchar(l.io)

function emit(l::Lexer, kind::TokenKind)
    text = String(take!(l.charstore))
    val = kind == INT ? parse(Int, text) :
          kind == BACKREF ? parse(Int, text[2:end]) :
          nothing
    Token(kind, l.token_start_pos, position(l.io), val)
end

function readchar(l::Lexer)
    l.current_char = readchar(l.io)
    write(l.charstore, l.current_char)
    if l.current_char == '\n'
        l.current_row += 1
        l.current_col = 1
    elseif !eof(l.current_char)
        l.current_col += 1
    end
    return l.current_char
end

const operators = Dict(
    '+' => PLUS,
    '-' => MINUS,
    '*' => TIMES,
    '/' => DIV,
)

function accept(l::Lexer, f::Union{Function, Char, Vector{Char}, String})
    c = peekchar(l)
    if isa(f, Function)
        ok = f(c)
    elseif isa(f, Char)
        ok = c == f
    else
        ok = c in f
    end
    ok && readchar(l)
    return ok
end

function accept_number(l::Lexer, f::F) where {F}
    while true
        pc, ppc = dpeekchar(l)
        if pc == '_' && !f(ppc)
            return
        elseif f(pc) || pc == '_'
            readchar(l)
        else
            return
        end
    end
end

"""
accept_batch(l::Lexer, f)
Consumes all following characters until `accept(l, f)` is `false`.
"""
function accept_batch(l::Lexer, f)
    ok = false
    while accept(l, f)
        ok = true
    end
    return ok
end

function lex_digit(l::Lexer)
    accept_number(l, isdigit)
    emit(l, INT)
end

is_non_nl_whitespace(c) = c != '\n' && iswhitespace(c)

function lex_whitespace(l::Lexer)
    accept_batch(l, is_non_nl_whitespace)
    return emit(l, WHITESPACE)
end

eof(io::IO) = Base.eof(io)
eof(c::Char) = c === EOF_CHAR

function start_token!(l::Lexer)
    l.token_start_pos = position(l.io)
end

function Base.popfirst!(l::Lexer)
    start_token!(l)
    c = readchar(l)
    if eof(c)
        return emit(l, ENDMARKER)
    elseif c in keys(operators)
        return emit(l, operators[c])
    elseif c == '('
        return emit(l, LPAREN)
    elseif c == ')'
        return emit(l, RPAREN)
    elseif isdigit(c)
        return lex_digit(l)
    elseif is_non_nl_whitespace(c)
        return lex_whitespace(l)
    elseif c == '_'
        accept_number(l, isdigit)
        return emit(l, BACKREF)
    elseif c == '\n'
        return emit(l, NEWLINE)
    else
        return emit(l, ERROR)
    end
end
Base.iterate(l::Lexer, state=nothing) = eof(l.io) ? nothing : (popfirst!(l),state)
Base.IteratorSize(::Type{Lexer}) = Base.SizeUnknown()

# CST Parser
struct CSTNode
    fullspan::UInt
    span::UnitRange{UInt}
end

abstract type CalcExpr; end

# Leaf nodes
struct Punctuation
    n::CSTNode
    kind::TokenKind
end
function Base.show(io::IO, punct::Punctuation)
    print(io, Punctuation, " ", punct.kind)
    print(io, "  ", punct.n.fullspan, " (", convert(UnitRange{Int}, punct.n.span), ")")
end
AbstractTrees.printnode(io::IO, lit::Punctuation) = show(io, lit)

struct Literal <: CalcExpr
    n::CSTNode
    value::Int64
end
function Base.show(io::IO, lit::Literal)
    print(io, Literal, "[", lit.value, "]")
    print(io, "  ", lit.n.fullspan, " (", convert(UnitRange{Int}, lit.n.span), ")")
end
AbstractTrees.printnode(io::IO, lit::Literal) = show(io, lit)

struct Backref <: CalcExpr
    n::CSTNode
    value::Int64
end
function Base.show(io::IO, br::Backref)
    print(io, Backref, "[_", br.value, "]")
    print(io, "  ", br.n.fullspan, " (", convert(UnitRange{Int}, br.n.span), ")")
end
AbstractTrees.printnode(io::IO, lit::Backref) = show(io, Backref)

function spansum(args...)
    fullspan = sum(arg->arg.n.fullspan, args)
    span = first(first(args).n.span):(
        fullspan - (last(args).n.fullspan -
                    last(last(args).n.span)))
    CSTNode(fullspan, span)
end

# Higher level nodes
struct Operator <: CalcExpr
    n::CSTNode
    lhs::CalcExpr
    kind::Punctuation
    rhs::CalcExpr
end
function Operator(lhs, kind, rhs)
    Operator(spansum(lhs, kind, rhs), lhs, kind, rhs)
end
AbstractTrees.children(o::Operator) = (o.lhs, o.kind, o.rhs)

struct ParenthesizedExpr <: CalcExpr
    n::CSTNode
    lparen::Punctuation
    expr::CalcExpr
    rparen::Punctuation
end
function ParenthesizedExpr(lparen, expr, rparen)
    ParenthesizedExpr(spansum(lparen, expr, rparen), lparen, expr, rparen)
end
AbstractTrees.children(o::ParenthesizedExpr) = (o.lparen, o.expr, o.rparen)

function AbstractTrees.printnode(io::IO, x::T) where T <: CalcExpr
    print(io, T, "  ", x.n.fullspan, " (", convert(UnitRange{Int}, x.n.span), ")")
end
function Base.show(io::IO, x::CalcExpr)
    if get(io, :compact, false)
        AbstractTrees.printnode(io, x)
    else
        AbstractTrees.print_tree(io, x, 3)
    end
end

mutable struct Parser
    l::Lexer
    first_leading_ws::Union{Token, Nothing}
    saw_newline::Bool
    current_token::Token
    had_error::Bool
end

function Parser(io::IO)
    l = Lexer(io)
    ws = nothing
    t = popfirst!(l)
    if t.kind == WHITESPACE || t.kind == NEWLINE
        ws = t
        while t.kind == WHITESPACE || t.kind == NEWLINE
            t = popfirst!(l)
        end
    end
    return Parser(l, ws, false, t, false)
end
tok(p::Parser) = p.current_token

macro propagate_error(p, x)
    quote
        x = $(esc(x))
        if $(esc(p)).had_error
            return nothing
        end
        x
    end
end

function emit_error(p::Parser, text)
    p.had_error = true
    return nothing
end

function INSTANCE(p::Parser)
    t = tok(p)
    fullspan = 0
    if p.first_leading_ws != nothing
        fullspan += t.startbyte - p.leading_ws.startbyte
    end
    p.first_leading_ws = nothing
    p.saw_newline = false
    span_start = fullspan
    # Gobble up whitespace
    nt = popfirst!(p.l)
    while nt.kind == WHITESPACE
        nt = popfirst!(p.l)
    end
    fullspan = nt.startbyte - t.endbyte
    # If the next token is a newline, add it to our trailing whitespace
    if nt.kind == NEWLINE
        p.saw_newline = true
        fullspan += 1
        nt = popfirst!(p.l)
        # Accumulate any leading whitespace now and get us to the first non-ws token
        if nt.kind == NEWLINE || nt.kind == WHITESPACE
            p.first_leading_ws = nt
            while nt.kind == NEWLINE || nt.kind == WHITESPACE
                nt = popfirst!(p.l)
            end
        end
    end
    p.current_token = nt
    t_span = t.endbyte - t.startbyte
    node = CSTNode(fullspan + t_span, span_start .+ (1:t_span))
    if t.kind == INT
        return Literal(node, t.val)
    elseif t.kind == BACKREF
        return Backref(node, t.val)
    elseif t.kind in (LPAREN, RPAREN) || is_op_kind(t.kind)
        return Punctuation(node, t.kind)
    else
        error("Unrecognized leaf token")
    end
end

function parse_operand(p::Parser)
    t = tok(p)
    if t.kind == LPAREN
        lparen = INSTANCE(p)
        expr = @propagate_error p parse_expr(p)
        rparen = INSTANCE(p)
        return ParenthesizedExpr(lparen, expr, rparen)
    elseif t.kind == INT || t.kind == BACKREF
        return INSTANCE(p)
    else
        return emit_error(p, "Any operand must always be a parenthesized expression, a backref or a literal")
    end
end

function parse_expr(p::Parser)
    lhs = @propagate_error p parse_operand(p)
    t = tok(p)
    if p.saw_newline || tok(p).kind == ENDMARKER
        return lhs
    end
    if !is_op_kind(t.kind)
        return emit_error(p, "Invalid non-operator after LHS expression")
    end
    op = INSTANCE(p)
    rhs = @propagate_error p parse_operand(p)
    return Operator(lhs, op, rhs)
end

function parse_toplevel(p::Parser)
    expr = @propagate_error p parse_expr(p)
    if !(p.saw_newline || tok(p) == ENDMARKER)
        return emit_error(p, "Extra token after the end of an expression")
    end
    expr
end

function parseall(data::String)
    p = Parser(IOBuffer(data))
    code = CalcExpr[]
    while !p.had_error && tok(p).kind != ENDMARKER
        expr = @propagate_error p parse_toplevel(p)
        push!(code, expr)
    end
    ParsedSourceCode(data, code)
end

# Interpreter
struct ParsedSourceCode
   text::String
   parsed::Vector{CalcExpr}
end

mutable struct InterpreterState
    backrefs::Vector{Any}
    code::ParsedSourceCode
    valstack::Vector{Any}
    it::Any
    pc::Int
end

function interesting_node_iterator(ast)
    Iterators.filter(x->!(x isa Calc.Punctuation || x isa Calc.ParenthesizedExpr), PostOrderDFS(ast))
end

function enter(code::ParsedSourceCode)
    it = interesting_node_iterator(code.parsed[1])
    InterpreterState(
        Any[],
        code,
        Any[],
        (it,),
        1
    )
end

const op_to_func_map = Dict(
    PLUS => +,
    MINUS => -,
    TIMES => *,
    DIV => /,
)

iterate_copy(it) = iterate(it)
iterate_copy(it, state) = iterate(it, deepcopy(state))

function _step!(s::InterpreterState)
    _step!(s, iterate_copy(s.it...)...)
end

function _step!(s::InterpreterState, expr, ns)
    if expr isa Literal
        push!(s.valstack, expr.value)
    elseif expr isa Backref
        push!(s.valstack, s.backrefs[expr.value])
    elseif expr isa Operator
        op2 = pop!(s.valstack)
        op1 = pop!(s.valstack)
        push!(s.valstack, op_to_func_map[expr.kind.kind](op1, op2))
    else
        error("Unrecognized expression")
    end
    s.it = (s.it[1], ns)
end

function step!(s::InterpreterState)
    y = iterate_copy(s.it...)
    if y !== nothing
        _step!(s)
    end
    it, state = s.it
    if iterate_copy(s.it...) === nothing
        if s.pc == length(s.code.parsed)
            return false
        end
        s.pc = s.pc + 1
        it = interesting_node_iterator(s.code.parsed[s.pc])
        s.it = (it,)
        push!(s.backrefs, s.valstack[end])
        empty!(s.valstack)
    end
    return true
end

# Debugger
function skip_literals!(s::InterpreterState)
    while true
        expr, ns = iterate_copy(s.it...)
        if !(expr isa Literal || expr isa Backref)
            break
        end
        step!(s)
    end
end

function DebuggerFramework.debug(code::ParsedSourceCode, args...)
    s = enter(code)
    skip_literals!(s)
    DebuggerFramework.RunDebugger([s], args...)
end

DebuggerFramework.locdesc(s::InterpreterState) = "statement $(s.pc)"
function DebuggerFramework.locinfo(s::InterpreterState)
    it, state = s.it
    pc_range = 1:(iterate_copy(s.it...) === nothing ? s.pc : s.pc - 1)
    offset = isempty(pc_range) ? 0 : sum(i->s.code.parsed[i].n.fullspan, pc_range)
    DebuggerFramework.BufferLocInfo(s.code.text, DebuggerFramework.compute_line(
        DebuggerFramework.SourceFile(s.code.text), offset), 0, 1)
end

function DebuggerFramework.print_next_state(io::IO, state, s::InterpreterState)
    y = iterate_copy(s.it...)
    if y === nothing
        if s.pc == length(s.code.parsed)
            return false
        end
        expr = first(interesting_node_iterator(s.code.parsed[s.pc + 1]))
    else
        expr, _ = y
    end
    print(io, "About to run: ")
    if expr isa Literal
        print(io, expr.value)
    elseif expr isa Backref
        print(io, "_", expr.value)
    elseif expr isa Operator
        print(io, s.valstack[end-1], " ", op_to_func_map[expr.kind.kind], " ", s.valstack[end])
    end
    println(io)
    return true
end

function DebuggerFramework.execute_command(state, s::InterpreterState, cmd::Union{Val{:s},Val{:si}}, command)
    if cmd == Val{:si}()
        if !step!(s)
            popfirst!(state.stack)
            return false
        end
        return true
    elseif cmd == Val{:s}()
        if !step!(s)
            popfirst!(state.stack)
            return false
        end
        skip_literals!(s)
        return true
    end
    return false
end

function DebuggerFramework.eval_code(state, s::InterpreterState, command)
    p = Parser(IOBuffer(command))
    expr = parse_expr(p)
    ns = enter(ParsedSourceCode(command, CalcExpr[expr]))
    ns.backrefs = copy(s.backrefs)
    while step!(ns); end
    ns.valstack[end]
end

using REPL
using REPL.LineEdit
function DebuggerFramework.language_specific_prompt(state, frame::InterpreterState)
    if haskey(state.language_modes, :calc)
        return state.language_modes[:calc]
    end
    calc_prompt = LineEdit.Prompt(DebuggerFramework.promptname(state.level, "calc");
        # Copy colors from the prompt object
        prompt_prefix = state.repl.prompt_color,
        prompt_suffix = (state.repl.envcolors ? Base.input_color : state.repl.input_color),
        on_enter = REPL.return_callback)
    calc_prompt.hist = state.main_mode.hist
    calc_prompt.hist.mode_mapping[:calc] = calc_prompt

    calc_prompt.on_done = (s,buf,ok)->begin
        if !ok
            LineEdit.transition(s, :abort)
            return false
        end
        xbuf = copy(buf)
        command = String(take!(buf))
        ok, result = DebuggerFramework.eval_code(state, command)
        REPL.print_response(state.repl, ok ? result : result[1], ok ? nothing : result[2], true, true)
        println(state.repl.t)
        LineEdit.reset_state(s)
    end
    calc_prompt.keymap_dict = LineEdit.keymap([REPL.mode_keymap(state.main_mode);state.standard_keymap])
    state.language_modes[:calc] = calc_prompt
    return calc_prompt
end

end


text:
#
# bilateral.jl --
#
# Implements the bilateral filter.
#
#------------------------------------------------------------------------------
#
# This file is part of the `LocalFilters.jl` package licensed under the MIT
# "Expat" License.
#
# Copyright (C) 2018, Éric Thiébaut.
#

module BilateralFilter

export bilateralfilter!, bilateralfilter!

using ..LocalFilters
using ..LocalFilters: Neighborhood, RectangularBox, Kernel, axes, _store!

bilateralfilter(A::AbstractArray{T}, args...) where {T<:Real} =
    # Provide type for computations and result.
    bilateralfilter(float(T), A, args...)

function bilateralfilter(::Type{T}, A::AbstractArray,
                         args...) where {T<:AbstractFloat}
    return bilateralfilter!(T, similar(Array{T}, axes(A)), A, args...)
end

"""
```julia
bilateralfilter([T,] A, Fr, Gs, ...)
```

yields the result of applying the bilateral filter on array `A`.

Argument `Fr` specifies the range kernel for smoothing differences in
intensities, it is a function which takes two values from `A` as arguments and
returns a nonnegative value.

Arguments `Gs, ...` specify the spatial kernel for smoothing differences in
coordinates.

Optional argument `T` can be used to force the floating-point type used for
(most) computations.

The in-place version is:

```julia
bilateralfilter!([T,] dst, A, Fr, Gs, ...)
```

which stores in `dst` the result of applying the bilateral filter on array `A`.

See [wikipedia](https://en.wikipedia.org/wiki/Bilateral_filter).

"""
function bilateralfilter!(dst::AbstractArray{Td,N},
                          A::AbstractArray{Ta,N},
                          args...) where {Td<:Real, Ta, N}
    # Provide type for computations.
    return bilateralfilter!(float(Td), dst, A, args...)
end

@doc @doc(bilateralfilter!) bilateralfilter

function bilateralfilter!(::Type{T},
                          dst::AbstractArray{Td,N},
                          A::AbstractArray{Ta,N},
                          Fr::Function,
                          Gs::Kernel{Tg,N}) where {T<:AbstractFloat, Td, Ta,
                                                   Tg<:Real, N}
    # The state is the tuple: (central_value, numerator, denominator).
    return localfilter!(dst, A, Kernel{T}(Gs),
                        (val) -> (val, zero(T), zero(T)),
                        (v, val, ker) -> _update(v, val, ker,
                                                 convert(T, Fr(val, v[1]))),
                        _final!)
end

function bilateralfilter!(::Type{T},
                          dst::AbstractArray{Td,N},
                          A::AbstractArray{Ta,N},
                          Fr::Function,
                          Gs::Function,
                          B::RectangularBox{N}
                          ) where {T<:AbstractFloat, Td, Ta, N}
    return bilateralfilter!(T, dst, A, Fr, Kernel(T, Gs, B))
end

function bilateralfilter!(::Type{T},
                          dst::AbstractArray{Td,N},
                          A::AbstractArray{Ta,N},
                          Fr::Function,
                          Gs::Function,
                          width::Integer
                          ) where {T<:AbstractFloat, Td<:Real, Ta<:Real, N}
    @assert width > 0 && isodd(width) "width or region of interest must be at least one and odd"
    h = Int(width) >> 1
    I = CartesianIndex(ntuple(i -> h, N))
    B = RectangularBox{N}(-I, I)
    return bilateralfilter!(T, dst, A, Fr, Gs, B)
end

# Range filter specifed by its standard deviation.
function bilateralfilter!(::Type{T},
                          dst::AbstractArray{Td,N},
                          A::AbstractArray{Ta,N},
                          σr::Real,
                          args...) where {T<:AbstractFloat,
                                          Td<:Real, Ta<:Real, N}
    @assert isfinite(σr) && σr > 0
    qr = _gaussfactor(T, σr)
    return bilateralfilter!(T, dst, A,
                            (v, v0) -> _gausswindow(qr, v - v0), args...)
end

# Distance filter specifed by its standard deviation and the size of the ROI.
function bilateralfilter!(::Type{T},
                          dst::AbstractArray{Td,N},
                          A::AbstractArray{Ta,N},
                          Fr::Function,
                          σs::Real, B) where {T<:AbstractFloat,
                                              Td<:Real, Ta<:Real, N}
    @assert isfinite(σs) && σs > 0
    qs = _gaussfactor(T, σs)
    return bilateralfilter!(T, dst, A, Fr, (i) -> _gausswindow(qs, i), B)
end

_gaussfactor(::Type{T}, σ::Real) where {T<:AbstractFloat} =
    -1/(2*convert(T, σ)^2)

# η = -1/2σ²
_gausswindow(η::T, I::CartesianIndex{N}) where {T<:AbstractFloat,N} =
    exp(sum((k) -> I[k]^2, 1:N)*η)

_gausswindow(η::T, x::Real, x0::Real) where {T<:AbstractFloat} =
    _gausswindow(η, x - x0)

_gausswindow(η::T, x::Real) where {T<:AbstractFloat} =
    exp(convert(T, x*x)*η)

function _update(v::Tuple{V,T,T}, val::V,
                 ws::T, wr::T) where {V, T<:AbstractFloat}
    w = wr*ws
    return (v[1], v[2] + convert(T, val)*w, v[3] + w)
end

_final!(dst, i, v::Tuple{V,T,T}) where {T<:AbstractFloat,V} =
    _store!(dst, i, (v[3] > zero(T) ? v[2]/v[3] : zero(T)))

end # module


text:
using CompScienceMeshes
using Test
using SauterSchwabQuadrature

pI = point(1,5,3)
pII = point(2,5,3)
pIII = point(7,1,0)
pVI = point(10,11,12)
pVII = point(10,11,13)
pVIII = point(11,11,12)

Sourcechart = simplex(pI,pII,pIII)
Testchart = simplex(pVI,pVII,pVIII)

Accuracy = 12
pd = PositiveDistance(SauterSchwabQuadrature._legendre(Accuracy,0.0,1.0))

function integrand(x,y)
			return(((x-pI)'*(y-pVII))*exp(-im*1*norm(x-y))/(4pi*norm(x-y)))
end

function INTEGRAND(û,v̂)
   n1 = neighborhood(Testchart, û)
   n2 = neighborhood(Sourcechart, v̂)
   x = cartesian(n1)
   y = cartesian(n2)
   output = integrand(x,y)*jacobian(n1)*jacobian(n2)
   return(output)
end

result = sauterschwab_parameterized(INTEGRAND, pd)-
             verifintegral2(Sourcechart, Testchart, integrand, Accuracy)

@test norm(result) < 1.e-3