text:
"""

        composite_and_radialaverage_withmask(azimuthal_average,composite_addition,buf,radiusbins,surface_pressure,variable_of_interest)
    """

function composite_and_radialaverage!(azimuthal_average,composite_addition,buf1,buf2,centers_labels_and_cyclones,radiusbins,variable_of_interest;maskcyclones = true)
    totalcyclonecount = 0
    nd = ndims(variable_of_interest)
    for timeindex in 1:size(variable_of_interest,nd)
        centers_and_labels,cyclones = centers_labels_and_cyclones[timeindex]
        if !isnothing(centers_and_labels[1])
            count = add_allcyclones!(composite_addition,buf1,buf2,selectdim(variable_of_interest,nd,timeindex),cyclones,centers_and_labels;maskcyclones)
            totalcyclonecount += count
        end
    end
    if !iszero(totalcyclonecount)
        composite_addition ./= totalcyclonecount
        if nd == 4
            for idx in 1:(length(radiusbins) - 1)
                azimuthal_average[idx,:] .= averageallindistance((radiusbins[idx],radiusbins[idx+1]),composite_addition,(256,256), 2000 )
            end
        elseif nd == 3
            for idx in 1:(length(radiusbins) - 1)
                azimuthal_average[idx] = averageallindistance((radiusbins[idx],radiusbins[idx+1]),composite_addition,(256,256), 2000 )
            end
        end

    end
end

function gather_diagnostic_3d(data_dir,exp_name,desired_var)
    diagnostic = zeros(Float32,512,512,80,600)
    for iter in 1:12
        file_name = string(exp_name,"_nosmoothing",iter+12,".jld")
        file_path = joinpath(data_dir,file_name)
        file = h5dopen(file_path, "r")
        index_initial = (iter-1)*(50)+1
        index_end = iter*50
        var_from_file = read(file, desired_var)
        diagnostic[:,:,:,index_initial:index_end] .=  var_from_file[:,:,:,:]
        close(file)
    end
    return diagnostic
end

"""
        getcomposites(file2d,file3d,diagnostics_dir,exp_,initial_timeindex,final_timeindex,partnumber,outputInterval)

    """
function get_composites(file2d,file3d,exp_name)

    @info "Starting composite routine"
    diagnostics_data_dir = "/global/cscratch1/sd/aramreye/for_postprocessing/ApeBudgetOutputs/"
    composites_data_dir = "/global/cscratch1/sd/aramreye/for_postprocessing/CompositeOutputs/"
    azimuthal_averages_data_dir = "/global/cscratch1/sd/aramreye/for_postprocessing/AzimuthualAverageOutputs/"

    variables_diag = ["convec_heating_anomaly","rad_heating_anomaly","buoyancy_anomaly"]
    variables_orig_3d = ["U","V","QV","TABS","QRAD","PP","W"]
    variables_orig_2d = ["USFC","VSFC","PW","Prec","LHF","SHF"]

    initial_timeindex = 841
    #    initial_timeindex = 1195
    final_timeindex = 1200

    iterator_time_2d    = initial_timeindex*2-1:2:final_timeindex*2
    iterator_time_3d    = initial_timeindex:1:final_timeindex

    radiusbins = 1000:2000:512000;
    ds3d                = Dataset(file3d)
    ds2d                = Dataset(file2d)
    z                   = variable(ds3d,"z")[:]                       :: Array{Float32,1}
    t                   = variable(ds3d,"time")[iterator_time_3d]     :: Array{Float32,1}
    surface_pressure    = variable(ds2d,"PSFC")[:,:,iterator_time_2d] :: Array{Float32,3}
    pressure_anomaly    = surface_pressure .- mean(surface_pressure,dims=(1,2))

    ##Create buffers to reuse on allocating functions

    size_var_of_interest_2d = size(surface_pressure)[1:end-1]
    size_var_of_interest_3d = (size_var_of_interest_2d...,length(z))
    ## Create 3D buffers
    #var_of_interest_3d      = zeros(size_var_of_interest_3d)
    composite_addition_3d   = zeros(Float32,size_var_of_interest_3d)
    buf_3d_2                = similar(composite_addition_3d)
    buf_3d_1                = similar(composite_addition_3d)
    azimuthal_average_3d    = zeros(Float32,length(radiusbins) - 1,length(z))
    ## Create 2D buffers
    #var_of_interest_2d      = zeros(size_var_of_interest_2d)
    composite_addition_2d   = zeros(Float32,size_var_of_interest_2d)
    buf_2d_1                = similar(composite_addition_2d)
    buf_2d_2                = similar(composite_addition_2d)
    azimuthal_average_2d    = zeros(Float32,length(radiusbins) - 1)
    @info size_var_of_interest_2d size_var_of_interest_3d

    centers_labels_and_cyclones = @views [
        detect_cyclones!(buf_2d_1, buf_2d_2, pressure_anomaly[:,:,timeindex],-5,2000) for
        timeindex in 1:size(surface_pressure,3)]; pressure_anomaly = [];

    composite_filename_withmask = string(composites_data_dir,exp_name        ,"_withmask.jld")
    composite_filename_nomask   = string(composites_data_dir,exp_name        ,"_nomask.jld")
    azimuthal_filename_withmask = string(azimuthal_averages_data_dir,exp_name,"_withmask.jld")
    azimuthal_filename_nomask   = string(azimuthal_averages_data_dir,exp_name,"_nomask.jld")

    ### *** add allcyclones ***
    @info "Creating composite for surface pressure "
    flush(stdout)
    composite_addition_2d .= 0.0
    buf_2d_1              .= 0.0
    buf_2d_2              .= 0.0
    azimuthal_average_2d  .= 0.0
    composite_and_radialaverage!(azimuthal_average_2d,composite_addition_2d,buf_2d_1,
    buf_2d_2,centers_labels_and_cyclones,radiusbins,surface_pressure;maskcyclones = true)

    @info "Writing"
    flush(stdout)
    jldopen(composite_filename_withmask, "w") do file
        write(file,"surface_pressure",composite_addition_2d)
    end

    jldopen(azimuthal_filename_withmask, "w") do file
        write(file,"surface_pressure",azimuthal_average_2d)
    end
    GC.gc()
    GC.gc()
    ## Create products without mask
    @info "Creating composite without mask for: surface_pressure"
    flush(stdout)
    composite_addition_2d .= 0.0
    buf_2d_1              .= 0.0
    buf_2d_2              .= 0.0
    azimuthal_average_2d  .= 0.0
    composite_and_radialaverage!(azimuthal_average_2d,composite_addition_2d,buf_2d_1,buf_2d_2
    ,centers_labels_and_cyclones,radiusbins,surface_pressure;maskcyclones = false)
    @info "Writing"
    flush(stdout)
    jldopen(composite_filename_nomask, "w") do file
        write(file,"surface_pressure",composite_addition_2d)
    end
    jldopen(azimuthal_filename_nomask, "w") do file
        write(file,"surface_pressure",azimuthal_average_2d)
    end
    GC.gc()
    GC.gc()

    for current_variable in variables_orig_3d
        ##NCDatasets.load!(ds["PSFC"].var,a,:,:,:)

        var_of_interest = variable(ds3d,current_variable)[:,:,:,iterator_time_3d]    :: Array{Float32,4}
        composite_addition_3d .= 0.0
        buf_3d_1              .= 0.0
        buf_3d_2              .= 0.0
        azimuthal_average_3d  .= 0.0

        ## Create products with mask
        @info "Creating composite for: " current_variable
        flush(stdout)
        composite_and_radialaverage!(azimuthal_average_3d,composite_addition_3d,buf_3d_1,buf_3d_2,centers_labels_and_cyclones,radiusbins,var_of_interest;maskcyclones = true)
        @info "Writing"
        flush(stdout)
        jldopen(composite_filename_withmask, "r+",mmaparrays=true) do file
            write(file,current_variable,composite_addition_3d)
        end
        jldopen(azimuthal_filename_withmask, "r+",mmaparrays=true) do file
            write(file,current_variable,azimuthal_average_3d)
        end
        GC.gc()
        GC.gc()
        ## Create products without mask
        @info "Creating composite without mask for: " current_variable
        flush(stdout)
        composite_addition_3d .= 0.0
        buf_3d_1              .= 0.0
        buf_3d_2              .= 0.0
        azimuthal_average_3d  .= 0.0
        composite_and_radialaverage!(azimuthal_average_3d,composite_addition_3d,buf_3d_1,buf_3d_2,centers_labels_and_cyclones,radiusbins,var_of_interest;maskcyclones = false)
        @info "Writing"
        flush(stdout)
        jldopen(composite_filename_nomask, "r+",mmaparrays=true) do file
            write(file,current_variable,composite_addition_3d)
        end
        jldopen(azimuthal_filename_nomask, "r+",mmaparrays=true) do file
            write(file,current_variable,azimuthal_average_3d)
        end
        GC.gc()
    end

    for current_variable in variables_orig_2d
        composite_addition_2d .= 0.0
        buf_2d_1              .= 0.0
        buf_2d_2              .= 0.0
        azimuthal_average_2d  .= 0.0

        var_of_interest = variable(ds2d,current_variable)[:,:,iterator_time_2d]    :: Array{Float32,3}

        ## Create products with mask
        @info "Creating composite for: " current_variable
        flush(stdout)
        composite_and_radialaverage!(azimuthal_average_2d,composite_addition_2d,buf_2d_1,buf_2d_2,centers_labels_and_cyclones,radiusbins,var_of_interest;maskcyclones = true)
        @info "Writing"
        flush(stdout)
        jldopen(composite_filename_withmask, "r+",mmaparrays=true) do file
            write(file,current_variable,composite_addition_2d)
        end
        jldopen(azimuthal_filename_withmask, "r+",mmaparrays=true) do file
            write(file,current_variable,azimuthal_average_2d)
        end
        GC.gc()
        GC.gc()
        ## Create products without mask
        @info "Creating composite without mask for: " current_variable
        flush(stdout)
        composite_addition_2d .= 0.0
        buf_2d_1              .= 0.0
        buf_2d_2              .= 0.0
        azimuthal_average_2d  .= 0.0
        composite_and_radialaverage!(azimuthal_average_2d,composite_addition_2d,buf_2d_1,buf_2d_2,centers_labels_and_cyclones,radiusbins,var_of_interest;maskcyclones = false)
        @info "Writing"
        flush(stdout)
        jldopen(composite_filename_nomask, "r+",mmaparrays=true) do file
            write(file,current_variable,composite_addition_2d)
        end
        jldopen(azimuthal_filename_nomask, "r+",mmaparrays=true) do file
            write(file,current_variable,azimuthal_average_2d)
        end
        GC.gc()
    end

    for current_variable in variables_diag
        var_of_interest = gather_diagnostic_3d(diagnostics_data_dir,exp_name,current_variable)   :: Array{Float32,4}
        composite_addition_3d .= 0.0
        buf_3d_1              .= 0.0
        buf_3d_2              .= 0.0
        azimuthal_average_3d  .= 0.0
        ## Create products with mask
        @info "Creating composite for: " current_variable
        flush(stdout)
        composite_and_radialaverage!(azimuthal_average_3d,composite_addition_3d,buf_3d_1,buf_3d_2,centers_labels_and_cyclones,radiusbins,var_of_interest;maskcyclones = true)
        @info "Writing"
        flush(stdout)
        jldopen(composite_filename_withmask, "r+",mmaparrays=true) do file
            write(file,current_variable,composite_addition_3d)
        end
        jldopen(azimuthal_filename_withmask, "r+",mmaparrays=true) do file
            write(file,current_variable,azimuthal_average_3d)
        end
        GC.gc()
        GC.gc()
        ## Create products without mask
        @info "Creating composite without mask for: " current_variable
        flush(stdout)
        composite_addition_3d .= 0.0
        buf_3d_1              .= 0.0
        buf_3d_2              .= 0.0
        azimuthal_average_3d  .= 0.0
        composite_and_radialaverage!(azimuthal_average_3d,composite_addition_3d,buf_3d_1,buf_3d_2,centers_labels_and_cyclones,radiusbins,var_of_interest;maskcyclones = false)
        @info "Writing"
        flush(stdout)
        jldopen(composite_filename_nomask, "r+",mmaparrays=true) do file
            write(file,current_variable,composite_addition_3d)
        end
        jldopen(azimuthal_filename_nomask, "r+",mmaparrays=true) do file
            write(file,current_variable,azimuthal_average_3d)
        end
        GC.gc()
        GC.gc()
    end

    close(ds3d)
    close(ds2d)

    return nothing

end


text:
using TransformUtils
using Test

@testset "test in-place and efficient inverse SE3..." begin

dest = zeros(3)
src = [1.0;0;0]
wTb = SE3([1.0;0;0],SO3(Matrix{Float64}(LinearAlgebra.I,3,3)))

doinversetransform!(dest, wTb, src)

@test norm(dest[1:3]) < 1e-10

# using BenchmarkTools
# @btime for i in 1:1000 doinversetransform!(dest, wTb, src); end
# Profile.clear()
# @profiler for i in 1:1000000 doinversetransform!(dest, wTb, src); end

end

ext: jl


text:
export marginalrule

@marginalrule typeof(*)(:A_in) (m_out::NormalDistributionsFamily, m_A::PointMass, m_in::F) where { F <: NormalDistributionsFamily } = begin
    A = mean(m_A)
    W = A' * precision(m_out) * A
    b_in = convert(promote_variate_type(F, NormalWeightedMeanPrecision), A' * weightedmean(m_out), W)
    q_in = prod(ProdAnalytical(), b_in, m_in)
    return (A = m_A, in = q_in)
end


text:
module LIBMP

const SUCCESS = 0
const FAILURE = 1

const libmp = ""

macro apicall(fun, argtypes, args...)
    isa(fun, QuoteNode) || error("first argument to @apicall should be a symbol")

    return quote
        status = ccall(($fun, libcuda), Cint,
                       $(esc(argtypes)), $(map(esc, args)...))

        if status != SUCCESS
            error()
        end
    end
end

const PARAM_VERSION = Cint(0)
const NUM_PARAMS = Cint(1)

function query_param(param)
    val = Ref{Cint}()
    @apicall(:mp_query_param, (Cint, Ref{Cint}), param, val)
    return val[]
end

# Opaque handles
struct mp_reg end
struct mp_request end
struct mp_window end

const Reg = Ptr{mp_reg};
const Request = Ptr{mp_request};
const Window = Ptr{mp_request};

const INIT_DEFAULT       = Cint(0)
const INIT_WQ_ON_GPU     = Cint(1)
const INIT_RX_CQ_ON_GPU  = Cint(2)
const INIT_TX_CQ_ON_GPU  = Cint(3)
const INIT_DBREC_ON_GPU  = Cint(4)

"""
    init(comm, peers, flags, gpu_id)

Initialize the MP library

# Arguments:
- comm: MPI communicator to use to bootstrap connection establishing
- peers: array of MPI ranks with which to establish a connection
- flags: combination of mp_init_flags
"""
function init(comm, peers, flags, gpu_id)
    @apicall(:mp_init, (MPI.CComm, Ptr{Cint}, Cint, Cint, Cint), comm, peers, length(peers), flags, gpu_id)
end

function finalize()
    @apicall(:mp_finalize, ())
end

function register(addr, length, reg)
    @apicall(:mp_register, (Ptr{Cvoid}, Csize_t, Ptr{Reg}), addr, length, reg)
end

function deregister(reg)
    @apicall(:mp_deregister, (Ptr{Reg},), reg)
end

function irecv(buf, size, peer, mp_req, req)
    @apicall(:mp_irecv, (Ptr{Cvoid}, Cint, Cint, Ptr{Req}, Ptr{Request}), buf, size, peer, mp_reg, req)
end

function isend(buf, size, peer, mp_req, req)
    @apicall(:mp_irecv, (Ptr{Cvoid}, Cint, Cint, Ptr{Reg}, Ptr{Request}), buf, size, peer, mp_reg, req)
end

function wait(reg)
    @apicall(:mp_wait, (Ptr{Reg},), reg)
end

function wait_all(regs::Vector{Reg})
    @apicall(:mp_wait_all, (UInt32, Ptr{Reg},), length(regs), regs)
end

function progress_all(regs::Vector{Reg})
    @apicall(:mp_progress_all, (UInt32, Ptr{Reg},), length(regs), regs)
end

# /*
# * CUDA stream synchronous primitives
# */
# int mp_send_on_stream  (void *buf, int size, int peer, mp_reg_t *mp_reg,
#                   mp_request_t *req, cudaStream_t stream);
# int mp_isend_on_stream (void *buf, int size, int peer, mp_reg_t *mp_reg,
#                   mp_request_t *req, cudaStream_t stream);
# //int mp_irecv_on_stream (void *buf, int size, int peer, mp_reg_t *mp_reg,
# //                        mp_request_t *req, cudaStream_t stream);

# /* vector sends/recvs
# * caveats: all blocks are within same registration
# */
# int mp_isendv(struct iovec *v, int nblocks, int peer, mp_reg_t *mp_reg, mp_request_t *req);
# int mp_irecvv(struct iovec *v, int nblocks, int peer, mp_reg_t *mp_reg, mp_request_t *req);

# int mp_isendv_on_stream (struct iovec *v, int nblocks, int peer, mp_reg_t *mp_reg,
#        mp_request_t *req, cudaStream_t stream);

# /*
# * GPU synchronous functions
# */
# int mp_wait_on_stream (mp_request_t *req, cudaStream_t stream);
# int mp_wait_all_on_stream (uint32_t count, mp_request_t *req, cudaStream_t stream);

# /* Split API to allow for batching of operations issued to the GPU
# */

function send_prepare(buf, size, peer, reg, req)
    @apicall(:mp_send_prepare, (Ptr{Cvoid}, Cint, Cint, Ptr{Reg}, Ptr{Request}), buf, size, peer, reg, req)
end
# int mp_sendv_prepare (struct iovec *v, int nblocks, int peer, mp_reg_t *mp_reg,
#                    mp_request_t *req);

# int mp_send_post_on_stream (mp_request_t *req, cudaStream_t stream);
# int mp_isend_post_on_stream (mp_request_t *req, cudaStream_t stream);
# int mp_send_post_all_on_stream (uint32_t count, mp_request_t *req, cudaStream_t stream);
# int mp_isend_post_all_on_stream (uint32_t count, mp_request_t *req, cudaStream_t stream);

# /*
# * One-sided communication primitives
# */

# /* window creation */
# int mp_window_create(void *addr, size_t size, mp_window_t *window_t);
# int mp_window_destroy(mp_window_t *window_t);

# enum mp_put_flags {
# MP_PUT_INLINE  = 1<<0,
# MP_PUT_NOWAIT  = 1<<1, // don't generate a CQE, req cannot be waited for
# };

# int mp_iput (void *src, int size, mp_reg_t *src_reg, int peer, size_t displ, mp_window_t *dst_window_t, mp_request_t *req, int flags);
# int mp_iget (void *dst, int size, mp_reg_t *dst_reg, int peer, size_t displ, mp_window_t *src_window_t, mp_request_t *req);

# int mp_iput_on_stream (void *src, int size, mp_reg_t *src_reg, int peer, size_t displ, mp_window_t *dst_window_t, mp_request_t *req, int flags, cudaStream_t stream);

# int mp_put_prepare (void *src, int size, mp_reg_t *src_reg, int peer, size_t displ, mp_window_t *dst_window_t, mp_request_t *req, int flags);

# int mp_iput_post_on_stream (mp_request_t *req, cudaStream_t stream);

# int mp_iput_post_all_on_stream (uint32_t count, mp_request_t *req, cudaStream_t stream);

# /*
# * Memory related primitives
# */

# enum mp_wait_flags {
# MP_WAIT_GEQ = 0,
# MP_WAIT_EQ,
# MP_WAIT_AND,
# };

# int mp_wait32(uint32_t *ptr, uint32_t value, int flags);
# int mp_wait32_on_stream(uint32_t *ptr, uint32_t value, int flags, cudaStream_t stream);

# static inline int mp_wait_dword_geq_on_stream(uint32_t *ptr, uint32_t value, cudaStream_t stream)
# {
# return mp_wait32_on_stream(ptr, value, MP_WAIT_GEQ, stream);
# }

# static inline int mp_wait_dword_eq_on_stream(uint32_t *ptr, uint32_t value, cudaStream_t stream)
# {
# return mp_wait32_on_stream(ptr, value, MP_WAIT_EQ, stream);
# }

# /*
# *
# */

# typedef struct mp_desc_queue *mp_desc_queue_t;

# int mp_desc_queue_alloc(mp_desc_queue_t *dq);
# int mp_desc_queue_free(mp_desc_queue_t *dq);
# int mp_desc_queue_add_send(mp_desc_queue_t *dq, mp_request_t *req);
# int mp_desc_queue_add_wait_send(mp_desc_queue_t *dq, mp_request_t *req);
# int mp_desc_queue_add_wait_recv(mp_desc_queue_t *dq, mp_request_t *req);
# int mp_desc_queue_add_wait_value32(mp_desc_queue_t *dq, uint32_t *ptr, uint32_t value, int flags);
# int mp_desc_queue_add_write_value32(mp_desc_queue_t *dq, uint32_t *ptr, uint32_t value);
# int mp_desc_queue_post_on_stream(cudaStream_t stream, mp_desc_queue_t *dq, int flags);

end


text:
using simulate_quantum_memory
code = RMCode();

# File I/O
timestamp  = convert(Int64, floor(time()*1000));
salt = getpid();
f = open("ffwd-$salt-$timestamp.csv", "w");

# Settings
NUM_TRIALS_PER_RATE = 100;
RATES = logspace(-2, -6, 10);
MODELGRAINS = ceil(5 ./ RATES);
MAX_TS = RATES .^ (-4);

write(f, "# NUM_TRIALS_PER_RATE = $NUM_TRIALS_PER_RATE\n");

# Trials
write(f, "rate, modelgrain, max_t, time, theta_hat, theta_true, exit_status, p_uncorrectable\n")
close(f);
for j=1:length(RATES)
    rate = RATES[j];
    modelgrain = MODELGRAINS[j];
    max_t = MAX_TS[j];
    print("Simulating $NUM_TRIALS_PER_RATE trials of rate = $rate\n")
    for k=1:NUM_TRIALS_PER_RATE
        model = SingleAngleModel(rate, modelgrain);
        EXPORT = single_angle_fast_forward(model, max_t, report=false)
        time = EXPORT["time"];
        exit_status = EXPORT["exit_status"];
        p_uncorrectable = EXPORT["p_uncorrectable"];
        theta_hat = model.theta_hat;
        theta_true = model.theta_true;
        print("$rate, $modelgrain, $max_t, $time, $theta_hat, $theta_true, $exit_status, $p_uncorrectable\n")
        f = open("ffwd-$salt-$timestamp.csv", "a");
        write(f, "$rate, $modelgrain, $max_t, $time, $theta_hat, $theta_true, $exit_status, $p_uncorrectable\n");
        close(f);
    end
end
f = open("ffwd-$salt-$timestamp.csv", "a");
write(f, "# EOF");
close(f);


text:
A=[1. 2. 3.; 4. 5. 6.; 7. 8. 9.]

    function cost(X::ProductRepr)
        return cost([submanifold_components(X)...])
    end
    function cost(X::Array{Matrix{Float64},1})
        return -0.5 * norm(transpose(X[1]) * A * X[2])^2
    end

    function egrad(X::Array{Matrix{Float64},1})
        U = X[1]
        V = X[2]
        AV = A*V
        AtU = transpose(A)*U
        return [ -AV*(transpose(AV)*U), -AtU*(transpose(AtU)*V) ];
    end

    function rgrad(M::ProductManifold, X::ProductRepr)
        eG = egrad([submanifold_components(M,X)...])
        x = [submanifold_components(M,X)...]
        return Manifolds.ProductRepr(project.(M.manifolds, x, eG)...)
    end

    function e2rHess(M::Grassmann, x, ξ, eGrad::Matrix{T},Hess::Matrix{T}) where T<:Union{U, Complex{U}} where U<:AbstractFloat
        pxHess = project(M,x,Hess)
        xtGrad = x'*eGrad
        ξxtGrad = ξ*xtGrad
        return pxHess - ξxtGrad
    end

    function eHess(X::Array{Matrix{Float64},1}, H::Array{Matrix{Float64},1})
        U = X[1]
        V = X[2]
        Udot = H[1]
        Vdot = H[2]
        AV = A*V
        AtU = transpose(A)*U
        AVdot = A*Vdot
        AtUdot = transpose(A)*Udot
        return [ -(AVdot*transpose(AV)*U + AV*transpose(AVdot)*U + AV*transpose(AV)*Udot),
                 -(AtUdot*transpose(AtU)*V + AtU*transpose(AtUdot)*V + AtU*transpose(AtU)*Vdot)
            ]
    end

    function rhess(M::ProductManifold, X::ProductRepr, H::ProductRepr)
        x = [submanifold_components(M,X)...]
        h = [submanifold_components(M,H)...]
        eG = egrad(x)
        eH = eHess(x,h)
        return Manifolds.ProductRepr(e2rHess.(M.manifolds, x, h, eG, eH)...)
    end

@testset "Manopt Trust-Region" begin
    import Random: seed!
    seed!(42);

    N = Grassmann(3,2)
    M = N × N

    x = random_point(M)

    @test_throws ErrorException trust_regions(M, cost, rgrad, x, rhess; ρ_prime = 0.3)
    @test_throws ErrorException trust_regions(M, cost, rgrad, x, rhess; Δ_bar = -0.1)
    @test_throws ErrorException trust_regions(M, cost, rgrad, x, rhess; Δ = -0.1)
    @test_throws ErrorException trust_regions(M, cost, rgrad, x, rhess; Δ_bar = 0.1, Δ = 0.11)

    X = trust_regions(M, cost, rgrad, x, rhess; Δ_bar=4*sqrt(2*2) )
    opt = trust_regions(M, cost, rgrad, x, rhess; Δ_bar=4*sqrt(2*2), return_options=true )
    @test isapprox(M,X,get_solver_result(opt))

    @test cost(X) + 142.5 ≈ 0 atol=10.0^(-13)

    XuR = trust_regions(M, cost, rgrad, x, rhess;
        Δ_bar=4*sqrt(2*2),
        useRandom = true
    )

    @test cost(XuR) + 142.5 ≈ 0 atol=10.0^(-12)

    XaH = trust_regions(
        M,
        cost,
        rgrad,
        x,
        (p,x,ξ) -> approxHessianFD(
            p,
            x,
            x -> rgrad(p,x),
            ξ;
            stepsize=2^(-9),
            transport=ProductVectorTransport(ProjectionTransport(),ProjectionTransport())
        );
        stopping_criterion = StopWhenAny(
            StopAfterIteration(2000),
            StopWhenGradientNormLess(10^(-6))
        ),
        Δ_bar=4*sqrt(2*2),
    )
    @test cost(XaH) + 142.5 ≈ 0 atol=10.0^(-10)

    ξ = random_tangent(M,x)
    @test_throws ErrorException getHessian(SubGradientProblem(M,cost,rgrad),x, ξ)

    # Test the random step trust region
    p = HessianProblem(M, cost, rgrad, rhess, (M,x,ξ) -> ξ)
    o = TrustRegionsOptions(x, StopAfterIteration(2000), 10.0^(-8),
        sqrt(manifold_dimension(M)), retract, true, 0.1, 1000.)
    @test step_solver!(p,o,0) == nothing

    η = truncated_conjugate_gradient_descent(M, cost, rgrad, x, ξ, rhess, 0.5)
    ηOpt = truncated_conjugate_gradient_descent(M, cost, rgrad, x, ξ, rhess, 0.5; return_options=true)
    @test submanifold_components(get_solver_result(ηOpt)) == submanifold_components(η)
end


text:
import .Plots: Plots

function Base.convert(t::Type{PNG}, plot::Plots.Plot)
    pb = PipeBuffer()
    show(pb, MIME("image/png"), plot)
    return PNG(pb)
end

preprocess(name, plot::Plots.Plot, data) = preprocess(name, convert(PNG, plot), data)
preprocess(name, plots::AbstractArray{<:Plots.Plot}, data) = begin
    for (i, plot)=enumerate(plots)
        preprocess(name*"/$i", plot, data)
    end
    return data
end

"""
    log_image(logger, name, plot::Plots.Plot; [step=current_step])

Renders the Plots' and sends it to TensorBoard as an image with tag `name`.
"""
log_image(lg::TBLogger, name::AbstractString, img::Plots.Plot; step=nothing) =
    log_keyval(lg, name, img, step)

log_image(lg::TBLogger, name::AbstractString, img::AbstractArray{<:Plots.Plot}; step=nothing) =
    log_keyval(lg, name, img, step)


text:
# ---
# title: 453. Minimum Moves to Equal Array Elements
# id: problem453
# author: Indigo
# date: 2021-02-02
# difficulty: Easy
# categories: Math
# link: <https://leetcode.com/problems/minimum-moves-to-equal-array-elements/description/>
# hidden: true
# ---
#
# Given a **non-empty** integer array of size _n_ , find the minimum number of
# moves required to make all array elements equal, where a move is incrementing
# _n_ \- 1 elements by 1.
#
# **Example:**
#
#
#
#     Input:
#     [1,2,3]
#
#     Output:
#     3
#
#     Explanation:
#     Only three moves are needed (remember each move increments two elements):
#
#     [1,2,3]  =>  [2,3,3]  =>  [3,4,3]  =>  [4,4,4]
#
#
#
## @lc code=start
using LeetCode

function min_moves_453(nums::Vector{Int})
    sum(nums) - minimum(nums) * length(nums)
end
## @lc code=end


text:
__precompile__()

module ConditionalJuMP

using JuMP
using JuMP: GenericAffExpr, GenericRangeConstraint, Variable, AbstractJuMPScalar
using IntervalArithmetic: Interval
import Base: <=, ==, >=, !, &, hash, show
using Printf
using Nullables

export @disjunction,
    @implies,
    @?,
    @switch,
    @ifelse,
    warmstart!,
    interval,
    upperbound,
    lowerbound

include("macros.jl")

const NArg{N} = NTuple{N, Any}

"""
Simplification function that chooses the appropriate algorithm based on the number
of variables.
"""
function simplify!(e::JuMP.GenericAffExpr{T, Variable}) where T
    n = length(e.vars)
    ((n < 100) ? simplify_inplace! : simplify_dict!)(e)
end

"""
Naive O(N^2) simplification. Slower for very large expressions, but allocates
no memory and is solidly faster for expressions with < 100 variables.
"""
function simplify_inplace!(e::JuMP.GenericAffExpr{T, Variable}) where T
    i1 = 1
    iend = length(e.vars)
    while i1 < iend
        i2 = i1 + 1
        while i2 <= iend
            if e.vars[i1].col == e.vars[i2].col
                e.coeffs[i1] += e.coeffs[i2]
                e.vars[i2] = e.vars[iend]
                e.coeffs[i2] = e.coeffs[iend]
                iend -= 1
            else
                i2 += 1
            end
        end
        i1 += 1
    end
    resize!(e.vars, iend)
    resize!(e.coeffs, iend)
    if iszero(e.constant)
        # Ensure that we always use the canonical (e.g. positive) zero
        e.constant = zero(e.constant)
    end
    e
end

"""
O(N) simplification, but with a substantially larger constant cost due to the
need to construct a Dict.
"""
function simplify_dict!(e::JuMP.GenericAffExpr{T, Variable}) where T
    vars = Variable[]
    coeffs = T[]

    var_map = Dict{Int, Int}()
    for i in eachindex(e.vars)
        v, c = e.vars[i], e.coeffs[i]
        idx = v.col
        if c != 0
            if haskey(var_map, idx)
                coeffs[var_map[idx]] += c
            else
                push!(vars, v)
                push!(coeffs, c)
                var_map[idx] = length(vars)
            end
        end
    end
    constant = e.constant
    if iszero(constant)
        constant = zero(constant)
    end
    e.vars = vars
    e.coeffs = coeffs
    e.constant = constant
    AffExpr(vars, coeffs, constant)
end

function normalize!(c::GenericRangeConstraint)
    @assert c.ub != Inf || c.lb != -Inf "Can't construct a constraint with infinite lower and upper bounds"

    # terms = c.terms
    # lb = c.lb
    # ub = c.ub
    if c.ub == Inf
        # then flip the constraint over so its upper bound is finite
        c.lb, c.ub = -c.ub, -c.lb
        c.terms.coeffs .*= -1
        c.terms.constant *= -1
        # terms = -terms
    end
    if c.ub != 0
        x = c.ub
        c.ub = 0
        c.terms.constant -= x
        c.lb -= x
    end
    if iszero(c.lb)
        # Ensure that we don't have negative zeros so that hashing is reliable
        c.lb = zero(c.lb)
    end
    if iszero(c.ub)
        c.ub = zero(c.ub)
    end
    simplify!(c.terms)
    c
end

struct Constraint{T}
    c::GenericRangeConstraint{GenericAffExpr{T, Variable}}
    hash::UInt
    function Constraint{T}(c::GenericAffExpr{T}, lb, ub) where T
        constraint = GenericRangeConstraint{GenericAffExpr{T, Variable}}(copy(c), lb, ub)
        normalize!(constraint)
        h = _hash(constraint, UInt(0))
        new(constraint, h)
    end
end

function Constraint(c::GenericAffExpr{T1}, lb::T2, ub::T3) where {T1, T2, T3}
    Constraint{promote_type(T1, T2, T3)}(c, lb, ub)
end

# work-around because JuMP doesn't properly define hash()
function _hash(x::GenericRangeConstraint, h::UInt)
    h = hash(x.lb, h)
    h = hash(x.ub, h)
    h = hash(x.terms.constant, h)
    for v in x.terms.vars
        h = hash(v, h)
    end
    for c in x.terms.coeffs
        h = hash(c, h)
    end
    h
end

Base.hash(c::Constraint, h::UInt) = hash(c.hash, h)

function (==)(e1::Constraint, e2::Constraint)
    (e1.c.lb == e2.c.lb) || return false
    (e1.c.ub == e2.c.ub) || return false
    e1.c.terms.constant == e2.c.terms.constant || return false
    length(e1.c.terms.vars) == length(e2.c.terms.vars) || return false
    for i in eachindex(e1.c.terms.vars)
        e1.c.terms.vars[i] == e2.c.terms.vars[i] || return false
        e1.c.terms.coeffs[i] == e2.c.terms.coeffs[i] || return false
    end
    return true
end

show(io::IO, h::Constraint) = print(io, "$(h.c.lb) <= $(h.c.terms) <= $(h.c.ub)")

const Conditional = Vector{Constraint{Float64}}

show(io::IO, c::Conditional) = print(io, join(c, " & "))

Conditional(::typeof(<=), x, y) = Conditional([Constraint(x - y, -Inf, 0)])
Conditional(::typeof(>=), x, y) = Conditional([Constraint(y - x, -Inf, 0)])
Conditional(::typeof(==), x, y) = Conditional([Constraint(x - y, 0, 0)])
(&)(c1::Conditional, c2::Conditional) = union(c1, c2)

struct UnhandledComplementException <: Exception
    c::Conditional
end

function Base.showerror(io::IO, e::UnhandledComplementException)
    print(io, "The complement of condition $(e.c) cannot be automatically determined. You will need to manually specify a disjunction covering this condition and all of its alternatives")
end

# lb <= x <= ub
# -> x <= lb || x >= ub
hascomplement(c::Conditional) = (length(c) == 1) && (first(c).c.lb == -Inf || first(c).c.ub == Inf)
function complement(c::Conditional)
    hascomplement(c) || throw(UnhandledComplementException(c))
    constr = first(c).c
    if constr.lb == -Inf
        new_constr = Constraint(constr.terms, constr.ub, Inf)
    else
        @assert constr.ub == Inf
        new_constr = Constraint(constr.terms, -Inf, constr.lb)
    end
    Conditional([new_constr])
end
(!)(c::Conditional) = complement(c)

lowerbound(x::Number) = x
upperbound(x::Number) = x
lowerbound(x::Variable) = JuMP.getlowerbound(x)
upperbound(x::Variable) = JuMP.getupperbound(x)

interval(x::Number, simplify=true) = Interval(x, x)
interval(x::Variable, simplify=true) = Interval(JuMP.getlowerbound(x), JuMP.getupperbound(x))

function interval(e::JuMP.GenericAffExpr, needs_simplification=true)
    if needs_simplification
        simplify!(e)
    end
    if isempty(e.coeffs)
        return Interval(e.constant, e.constant)
    else
        result = Interval(e.constant, e.constant)
        for i in eachindex(e.coeffs)
            var = e.vars[i]
            coef = e.coeffs[i]
            result += Interval(coef, coef) * Interval(getlowerbound(var), getupperbound(var))
        end
        return result
    end
end

upperbound(e::GenericAffExpr, simplify=true) = upperbound(interval(e, simplify))
lowerbound(e::GenericAffExpr, simplify=true) = lowerbound(interval(e, simplify))
lowerbound(i::Interval) = i.lo
upperbound(i::Interval) = i.hi

const Implication = Pair{Conditional, Conditional}

struct Indicator
    var::Variable
    sense::Bool

    function Indicator(var::Variable, sense::Bool)
        @assert JuMP.getcategory(var) == :Bin
        new(var, sense)
    end
end

Indicator(v::Variable) = Indicator(v, true)

(!)(i::Indicator) = Indicator(i.var, !i.sense)

variable(i::Indicator) = i.var

struct IndicatorMap
    model::Model
    indicators::Dict{Conditional, Indicator}
    disjunctions::Vector{Vector{Implication}}
    y_idx::Base.RefValue{Int}
    z_idx::Base.RefValue{Int}
end

IndicatorMap(m::Model) = IndicatorMap(m,
    Dict{Conditional, AffExpr}(),
    Vector{Vector{Implication}}(),
    Ref(1),
    Ref(1))

function getindmap!(m::Model)::IndicatorMap
    if haskey(m.ext, :indicator_map)
        return m.ext[:indicator_map]
    else
        i = IndicatorMap(m)
        m.ext[:indicator_map] = i
        return i
    end
end

function newcontinuousvar(m::IndicatorMap)
    varname = @sprintf("y_{anon%d}", m.y_idx[])
    var = @variable(m.model, basename=varname)
    m.y_idx[] = m.y_idx[] + 1
    var
end

function newcontinuousvar(m::IndicatorMap, size::Dims)
    varname = @sprintf("{y_{anon%d}}", m.y_idx[])
    var = reshape(@variable(m.model, [1:prod(size)], basename=varname), size)
    m.y_idx[] = m.y_idx[] + 1
    var
end

newcontinuousvar(m::Model, args...) = newcontinuousvar(getindmap!(m), args...)

function newbinaryvar(m::IndicatorMap)
    var = @variable(m.model, basename=@sprintf("z_{anon%d}", m.z_idx[]), category=:Bin)
    m.z_idx[] = m.z_idx[] + 1
    var
end

newbinaryvar(m::Model, args...) = newbinaryvar(getindmap!(m), args...)

function getindicator!(m::IndicatorMap, c::Conditional, can_create=true)
    if haskey(m.indicators, c)
        return m.indicators[c]
    else
        if !can_create
            @show c
            error("Not allowed to create a new variable here. Something has gone wrong")
        end
        z = Indicator(newbinaryvar(m))
        # z = convert(AffExpr, newbinaryvar(m))
        implies!(m.model, z, c)
        m.indicators[c] = z

        if hascomplement(c)
            m.indicators[!c] = !z
        end
        return z
    end
end

getindicator!(m::Model, c::Conditional) = getindicator!(getindmap!(m), c)

function implies!(m::Model, z::Indicator, c::Conditional)
    for constraint in c
        implies!(m, z, constraint)
    end
end

struct UnboundedVariableException <: Exception
    terms::AffExpr
end

function Base.showerror(io::IO, e::UnboundedVariableException)
    print(io, "Cannot create an implication involving the expression: $(e.terms) because not all of the variables in it are fully bounded. Please use `JuMP.setlowerbound()` and `JuMP.setupperbound()` to set finite bounds for all variables appearing in this expression.")
end

function append(v::AbstractVector{T}, x::T) where T
    result = similar(v, length(v) + 1)
    for i in 1:length(v)
        result[i] = v[i]
    end
    result[end] = x
    result
end

function implies!(m::Model, z::Indicator, constraint::Constraint)
    expr_interval = interval(constraint.c.terms, false)
    if constraint.c.ub != Inf
        M = upperbound(expr_interval) - constraint.c.ub
        isfinite(M) || throw(UnboundedVariableException(constraint.c.terms))
        if z.sense
            @constraint m constraint.c.terms <= constraint.c.ub + M * (1 - z.var)
            constant = constraint.c.terms.constant - M
            expr = AffExpr(append(constraint.c.terms.vars, z.var), append(constraint.c.terms.coeffs, M), 0.0)
        else
            # @constraint m constraint.c.terms <= constraint.c.ub + M * (z.var)
            constant = constraint.c.terms.constant
            expr = AffExpr(append(constraint.c.terms.vars, z.var), append(constraint.c.terms.coeffs, -M), 0.0)
        end
        JuMP.addconstraint(m, LinearConstraint(expr, -Inf, constraint.c.ub - constant))
    end
    if constraint.c.lb != -Inf
        M = constraint.c.lb - lowerbound(expr_interval)
        isfinite(M) || throw(UnboundedVariableException(constraint.c.terms))
        if z.sense
            # @constraint m constraint.c.terms >= constraint.c.lb - M * (1 - z.var)
            constant = constraint.c.terms.constant + M
            expr = AffExpr(append(constraint.c.terms.vars, z.var), append(constraint.c.terms.coeffs, -M), 0.0)
        else
            # @constraint m constraint.c.terms >= constraint.c.lb - M * (z.var)
            constant = constraint.c.terms.constant
            expr = AffExpr(append(constraint.c.terms.vars, z.var), append(constraint.c.terms.coeffs, M), 0.0)
        end
        JuMP.addconstraint(m, LinearConstraint(expr, constraint.c.lb - constant, Inf))
    end
end

disjunction!(m::Model, args...) = disjunction!(getindmap!(m), args...)

function disjunction!(indmap::IndicatorMap, imps::NTuple{1, Implication})
    z = Indicator(newbinaryvar(indmap))
    lhs, rhs = imps[1]
    implies!(indmap.model, z, lhs)
    implies!(indmap.model, z, rhs)
    implies!(indmap.model, !z, !lhs)
    _fix(z, true)
end

function disjunction!(indmap::IndicatorMap, imps::NTuple{2, Implication})
    z = getindicator!(indmap, first(imps[1]))
    implies!(indmap.model, z, last(imps[1]))
    indmap.indicators[first(imps[2])] = !z
    implies!(indmap.model, !z, first(imps[2]))
    implies!(indmap.model, !z, last(imps[2]))
    push!(indmap.disjunctions, collect(imps))
end

function disjunction!(indmap::IndicatorMap, imps::Union{Tuple, AbstractArray})
    zs = getindicator!.(indmap, first.(imps))
    implies!.(indmap.model, zs, last.(imps))
    JuMP.addconstraint(indmap.model,
        LinearConstraint(AffExpr(collect(variable.(zs)), ones(length(zs)), 0.0), 1.0, 1.0))
    # @constraint(indmap.model, sum(variable.(zs)) == 1)
    push!(indmap.disjunctions, collect(imps))
end

function  disjunction!(indmap::IndicatorMap,
                       conditions::Union{Tuple{Vararg{Conditional}},
                                         <:AbstractArray{Conditional}})
    disjunction!(indmap, (c -> (c => Conditional())).(conditions))
end

Broadcast.broadcastable(i::IndicatorMap) = Ref(i)
Broadcast.broadcastable(m::Model) = Ref(m)

implies!(m::Model, imps::Implication...) = disjunction!(m, imps)

implies!(m::Model, imps::Pair...) = implies!(m, convert.(Implication, imps)...)

Base.convert(::Type{Conditional}, ::Nothing) = Conditional()

function implies!(m::Model, imp::Implication)
    c1, c2 = imp
    comp1 = !c1
    disjunction!(m, (imp, (comp1 => Conditional())))
end

getmodel(x::JuMP.Variable) = x.m
getmodel(x::JuMP.GenericAffExpr) = first(x.vars).m
getmodel(c::Constraint) = getmodel(c.c.terms)

function getmodel(c::Conditional)::Model
    for arg in c
        if !isempty(arg.c.terms.vars)
            return getmodel(arg)
        end
    end
    error("Could not find JuMP Model in conditional $c")
end

nan_to_null(x) = isnan(x) ? Nullable{typeof(x)}() : Nullable{typeof(x)}(x)

"""
Like JuMP.getvalue, but returns a Nullable{T}() for unset variables instead
of throwing a warning
"""
_getvalue(x::Variable) = nan_to_null(JuMP._getValue(x))
_getvalue(x::Number) = nan_to_null(x)

function _getvalue(x::JuMP.GenericAffExpr{T, Variable}) where {T}
    result::Nullable{T} = x.constant
    for i in eachindex(x.coeffs)
        result = result .+ x.coeffs[i] .* _getvalue(x.vars[i])
    end
    result
end

function _getvalue(h::Constraint)
    v = _getvalue(h.c.terms)
    max.(v .- h.c.ub, h.c.lb .- v)
end

_getvalue(c::Conditional) = maximum(_getvalue, c)

_setvalue(v::Variable, x) = JuMP.setvalue(v, x)
function _setvalue(s::JuMP.GenericAffExpr, x)
    @assert length(s.vars) == 1
    _setvalue(s.vars[1], (x - s.constant) / s.coeffs[1])
end

function _getvalue(i::Indicator)
    v = _getvalue(i.var)
    if i.sense
        return v
    else
        return 1 .- v
    end
end

Base.show(io::IO, cv::Implication) = print(io, "(", first(cv), ") implies (", last(cv), ")")

Base.@pure isjump(x) = false
# Base.@pure isjump(args::Tuple) = any(isjump, args)
Base.@pure isjump(x, args...) = isjump(x) || isjump(args...)
Base.@pure isjump(arg::Pair) = isjump(arg.first)
Base.@pure isjump(c::Constraint) = true
Base.@pure isjump(c::Conditional) = true
Base.@pure isjump(x::AbstractJuMPScalar) = true
Base.@pure isjump(x::AbstractArray{<:AbstractJuMPScalar}) = true

function _conditional(op::Op, args...) where {Op <: Union{typeof(<=), typeof(>=), typeof(==), typeof(in)}}
    if isjump(args...)
        Conditional(op, args...)
    else
        op(args...)
    end
end

_conditional(op, args...) = op(args...)

function _all_conditional(args)
    (&)(args...)
end

function switch!(m::Model, args::Pair{<:Conditional}...)
    y = newcontinuousvar(m)
    conditions = first.(args)
    values = last.(args)
    setlowerbound(y, minimum(lowerbound, values))
    setupperbound(y, maximum(upperbound, values))
    disjunction!(m, map(cv -> (cv[1] => @?(y == cv[2])), args))
    y
end

function switch!(m::Model, args::Pair{<:Conditional, <:AbstractArray}...)
    y = newcontinuousvar(m, size(last(args[1])))
    conditions = first.(args)
    values = last.(args)
    for I in eachindex(y)
        setlowerbound(y[I], minimum(v -> lowerbound(v[I]), values))
        setupperbound(y[I], maximum(v -> upperbound(v[I]), values))
    end
    disjunction!(m, map(cv -> (cv[1] => @?(y .== cv[2])), args))
    y
end

function ifelse_conditional(c::Conditional, v1, v2)
    @assert size(v1) == size(v2)
    m = getmodel(c)
    y = newcontinuousvar(m)
    setlowerbound(y, min(lowerbound(v1), lowerbound(v2)))
    setupperbound(y, max(upperbound(v1), upperbound(v2)))
    disjunction!(m, (c => @?(y == v1), !c => @?(y == v2)))
    y
end

function ifelse_conditional(c::Conditional, v1::AbstractArray, v2::AbstractArray)
    @assert size(v1) == size(v2)
    m = getmodel(c)
    y = newcontinuousvar(m, size(v1))
    setlowerbound.(y, min.(lowerbound.(v1), lowerbound.(v2)))
    setupperbound.(y, max.(upperbound.(v1), upperbound.(v2)))
    @disjunction(m, (c => y .== v1), (!c => (y .== v2)))
    # disjunction!(m, (c => @?(y .== v1), !c => @?(y .== v2)))
    y
end

isfixed(v::Variable) = JuMP.getcategory(v) == :Fixed
isfixed(s::JuMP.GenericAffExpr) = all(isfixed, s.vars)
isfixed(i::Indicator) = isfixed(i.var)

_fix(v::Variable, x) = JuMP.fix(v, x)
_fix(i::Indicator, x::Bool) = _fix(i.var, x ⊻ (!i.sense))

function _unfix(v::Variable)
    JuMP.setcategory(v, :Bin)
    setlowerbound(v, 0)
    setupperbound(v, 1)
end

_unfix(i::Indicator) = _unfix(i.var)

function _setvalue(m::Model, c::Conditional)
    if length(c) != 1
        return nothing
    end
    constraint = first(c).c
    if isfinite(constraint.lb) && constraint.lb == constraint.ub
        _setvalue(constraint.terms, constraint.ub)
    end
    return nothing
end

_setvalue(i::Indicator, x::Bool) = _setvalue(i.var, x ⊻ !(i.sense))

function warmstart!(m::Model, fix=false)
    indmap = getindmap!(m)
    while true
        modified = false
        for implications in indmap.disjunctions
            violations = Nullable{Float64}[_getvalue(first(i)) for i in implications]
            if !any(isnull, violations)
                best_match = argmin(get.(violations))
                for i in eachindex(violations)
                    imp = implications[i]
                    lhs, rhs = imp
                    z = getindicator!(indmap, lhs, false)
                    satisfied = i == best_match
                    if fix
                        if !isfixed(z)
                            modified = true
                        end
                        _fix(z, satisfied)
                    else
                        if isnull(_getvalue(z))
                            modified = true
                        end
                        _unfix(z)
                        _setvalue(z, satisfied)
                    end
                    if satisfied
                        _setvalue(m, rhs)
                    end
                end
            end
        end
        if !modified
            break
        end
    end
end

function switch(args::Pair...)
    if isjump(args...)
        switch!(getmodel(args[1].first), args...)
    else
        for arg in args
            if arg.first
                return last(arg)
            end
        end
        error("One of the cases in the switch statement must always match")
    end
end

"""
Gurobi doesn't understand constant terms in the objective, which can
affect the way the relative MIP gap is interpreted. To work around this,
we'll replace the constant term k with an affine term 1*x for a new
variable x fixed to equal k.
"""
function handle_constant_objective!(m::Model)
    x = @variable(m, basename="x_{objective_constant}")
    obj = getobjective(m)
    k = obj.aff.constant
    JuMP.fix(x, k)
    JuMP.setobjective(m, getobjectivesense(m), obj - k + x)
    @assert getobjective(m).aff.constant == 0
end

end


text:
function index_problems(n)

  reps = parse(Int64, ENV["LAMP_REPS"])

  m = Int(floor(n / 10))
  A = randn(m, m)
  B = randn(m, m)
  C = Vector(randn(m))

  Benchmarker.add_data(csv, "loop_nai", Benchmarker.measure(reps, naive_loop, A, B, C))
  Benchmarker.add_data(csv, "loop_rec", Benchmarker.measure(reps, recommended_loop, A, B, C))

  @test(isapprox(naive_loop(A,B,C)[2], recommended_loop(A,B,C)[2]))
end

function naive_loop( A::DenseArray{Float64, 2},  B::DenseArray{Float64, 2},  C::Vector{Float64})
  result = copy(C)
  start::Float64 = 0.0
  finish::Float64 = 0.0
  Benchmarker.cachescrub()
  Benchmarker.gcscrub()
  start = time_ns()

  temp::DenseArray{Float64, 2} = zeros(size(A, 1), size(B, 2))
  for i in range(1, size(result, 1), step=1)
    temp = A*B
    @inbounds result[i] = temp[i, i]
  end

  finish = time_ns()
  return ((finish-start)*1e-9, result)
end

function recommended_loop( A::DenseArray{Float64, 2},  B::DenseArray{Float64, 2},  C::Vector{Float64})
  result = copy(C)
  start::Float64 = 0.0
  finish::Float64 = 0.0
  Benchmarker.cachescrub()
  Benchmarker.gcscrub()
  start = time_ns()

  temp::DenseArray{Float64, 2} = A*B
  for i in range(1, size(result, 1), step=1)
    @inbounds result[i] = temp[i, i]
  end

  finish = time_ns()
  return ((finish-start)*1e-9, result)
end


text:
using Documenter
using Literate

using Pkg
if haskey(ENV, "GITLAB_CI")
    Pkg.add([PackageSpec(name = x; rev = "master")
             for x in ["CUDAapi", "GPUArrays", "CUDAnative", "NNlib", "CUDAdrv"]])
end

using CuArrays

# generate tutorials
OUTPUT = joinpath(@__DIR__, "src/tutorials/generated")
Literate.markdown(joinpath(@__DIR__, "src/tutorials/intro.jl"), OUTPUT)

makedocs(
    modules = [CuArrays],
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    sitename = "CuArrays.jl",
    pages = [
        "Home" => "index.md",
        "Tutorials"  => [
            "tutorials/generated/intro.md"
        ],
    ],
    doctest = true
)


text:
# Single Layer
function (SL::SingleLayerKernel{T,S})(x,y)::T  where {T,S<:Elastostatic}
    N = ambient_dimension(S)
    x==y && return zero(T)
    μ = SL.op.μ
    λ = SL.op.λ
    ν = λ/(2*(μ+λ))
    r = x .- y
    d = norm(r)
    RRT = r*transpose(r) # r ⊗ rᵗ
    if N==2
        ID = SMatrix{2,2,Float64,4}(1,0,0,1)
        return 1/(8π*μ*(1-ν))*(-(3-4*ν)*log(d)*ID + RRT/d^2)
        # return (λ + 3μ)/(4*π*(N-1)*μ*(λ+2μ))* (-log(d)*ID + (λ+μ)/(λ+3μ)*RRT/d^2)
    elseif N==3
        ID = SMatrix{3,3,Float64,9}(1,0,0,0,1,0,0,0,1)
        return 1/(16π*μ*(1-ν)*d)*((3-4*ν)*ID + RRT/d^2)
    end
end

# Double Layer Kernel
function (DL::DoubleLayerKernel{T,S})(x,y,ny)::T where {T,S<:Elastostatic}
    N = ambient_dimension(S)
    x==y && return zero(T)
    μ = DL.op.μ
    λ = DL.op.λ
    ν = λ/(2*(μ+λ))
    r = x .- y
    d = norm(r)
    RRT = r*transpose(r) # r ⊗ rᵗ
    drdn = -dot(r,ny)/d
    dγdny = 1/d*drdn
    if N==2
        ID = SMatrix{2,2,Float64,4}(1,0,0,1)
        return -1/(4π*(1-ν)*d)*(drdn*((1-2ν)*ID + 2*RRT/d^2) + (1-2ν)/d*(r*transpose(ny) - ny*transpose(r)))
        # return μ/(2*(N-1)*π*(λ+2μ))*( (ID + N*(λ+μ)/(μ*d^2)*RRT)*dγdny - 1/d^N*(r*transpose(ny)-ny*transpose(r)) )
    elseif N==3
        ID = SMatrix{3,3,Float64,9}(1,0,0,0,1,0,0,0,1)
        return -1/(8π*(1-ν)*d^2)*(drdn * ((1-2*ν)*ID + 3*RRT/d^2) + (1-2*ν)/d*(r*transpose(ny) - ny*transpose(r)))
    end
end

# Adjoint Double Layer Kernel
function (ADL::AdjointDoubleLayerKernel{T,S})(x,y,nx)::T where {T,S<:Elastostatic}
    DL = DoubleLayerKernel{T}(ADL.op)
    return -DL(x,y,nx) |> transpose
    # N = ambient_dim(S)
    # x==y && return zero(T)
    # μ = ADL.op.μ
    # λ = ADL.op.λ
    # ν = λ/(2*(μ+λ))
    # r = x .- y
    # d = norm(r)
    # RRT = r*transpose(r) # r ⊗ rᵗ
    # if N==2
    #     ID = Mat{2,2,Float64,4}(1,0,0,1)
    #     return 1/(4π*(1-ν)*d)*((1-2ν)/d*(nx*transpose(r) + dot(r,nx)*ID - r*transpose(nx)) + 2/d^3*dot(r,nx)*RRT)
    # elseif N==3
    #     ID = Mat{3,3,Float64,9}(1,0,0,0,1,0,0,0,1)
    #     return 1/(8π*(1-ν)*d^2)*((1-2ν)/d*(-nx*transpose(r) + dot(r,nx)*ID + r*transpose(nx)) + 3/d^3*dot(r,nx)*RRT)
    # end
end

# Hypersingular kernel
function (HS::HyperSingularKernel{T,S})(x,y,nx,ny)::T where {T,S<:Elastostatic}
    N = ambient_dimension(S)
    x==y && return zero(T)
    μ = HS.op.μ
    λ = HS.op.λ
    ν = λ/(2*(μ+λ))
    r = x .- y
    d = norm(r)
    RRT = r*transpose(r) # r ⊗ rᵗ
    drdn    = dot(r,ny)/d
    if N==2
        ID = SMatrix{2,2,Float64,4}(1,0,0,1)
        return μ/(2π*(1-ν)*d^2)* (2*drdn/d*( (1-2ν)*nx*transpose(r) + ν*(dot(r,nx)*ID + r*transpose(nx)) - 4*dot(r,nx)*RRT/d^2 ) +
                                  2*ν/d^2*(dot(r,nx)*ny*transpose(r) + dot(nx,ny)*RRT) +
                                  (1-2*ν)*(2/d^2*dot(r,nx)*r*transpose(ny) + dot(nx,ny)*ID + ny*transpose(nx)) -
                                  (1-4ν)*nx*transpose(ny)
                                  )
    elseif N==3
        ID = SMatrix{3,3,Float64,9}(1,0,0,0,1,0,0,0,1)
        return μ/(4π*(1-ν)*d^3)* (3*drdn/d*( (1-2ν)*nx*transpose(r) + ν*(dot(r,nx)*ID + r*transpose(nx)) - 5*dot(r,nx)*RRT/d^2 ) +
                                  3*ν/d^2*(dot(r,nx)*ny*transpose(r) + dot(nx,ny)*RRT) +
                                  (1-2*ν)*(3/d^2*dot(r,nx)*r*transpose(ny) + dot(nx,ny)*ID + ny*transpose(nx)) -
                                  (1-4ν)*nx*transpose(ny)
                                  )
    end
end

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
################################################################################
#
#  Basic field access
#
################################################################################

@doc Markdown.doc"""
    order(a::NfRelOrdIdl) -> NfRelOrd

Returns the order of $a$.
"""
order(a::NfRelOrdIdl) = a.order

@doc Markdown.doc"""
    nf(a::NfRelOrdIdl) -> NumField

Returns the number field, of which $a$ is an integral ideal.
"""
nf(a::NfRelOrdIdl) = nf(order(a))

################################################################################
#
#  Parent
#
################################################################################

parent(a::NfRelOrdIdl) = a.parent

function check_parent(x::NfRelOrdIdl, y::NfRelOrdIdl)
   if order(x) !== order(y)
     error("Ideals do not have the same order.")
   end
end

################################################################################
#
#  "Assure" functions for fields
#
################################################################################

function assure_has_basis_pmatrix(a::Union{NfRelOrdIdl, NfRelOrdFracIdl})
  if isdefined(a, :basis_pmatrix)
    return nothing
  end
  if !isdefined(a, :pseudo_basis)
    error("No pseudo_basis and no basis_pmatrix defined.")
  end
  pb = pseudo_basis(a, copy = false)
  L = nf(order(a))
  M = zero_matrix(base_field(L), degree(L), degree(L))
  C = Vector{frac_ideal_type(order_type(base_field(L)))}()
  for i = 1:degree(L)
    elem_to_mat_row!(M, i, pb[i][1])
    push!(C, deepcopy(pb[i][2]))
  end
  M = M*basis_mat_inv(order(a), copy = false)
  a.basis_pmatrix = pseudo_hnf(PseudoMatrix(M, C), :lowerleft, true)
  return nothing
end

function assure_has_pseudo_basis(a::Union{NfRelOrdIdl, NfRelOrdFracIdl})
  if isdefined(a, :pseudo_basis)
    return nothing
  end
  if !isdefined(a, :basis_pmatrix)
    error("No pseudo_basis and no basis_pmatrix defined.")
  end
  P = basis_pmatrix(a, copy = false)
  B = basis_nf(order(a), copy = false)
  L = nf(order(a))
  K = base_field(L)
  pseudo_basis = Vector{Tuple{elem_type(L), frac_ideal_type(order_type(K))}}()
  for i = 1:degree(L)
    t = L()
    for j = 1:degree(L)
      t += P.matrix[i, j]*B[j]
    end
    push!(pseudo_basis, (t, deepcopy(P.coeffs[i])))
  end
  a.pseudo_basis = pseudo_basis
  return nothing
end

function assure_has_basis_matrix(a::Union{NfRelOrdIdl, NfRelOrdFracIdl})
  if isdefined(a, :basis_matrix)
    return nothing
  end
  a.basis_matrix = basis_pmatrix(a).matrix
  return nothing
end

function assure_has_basis_mat_inv(a::Union{NfRelOrdIdl, NfRelOrdFracIdl})
  if isdefined(a, :basis_mat_inv)
    return nothing
  end
  a.basis_mat_inv = inv(basis_matrix(a, copy = false))
  return nothing
end

################################################################################
#
#  Pseudo basis / basis pseudo-matrix
#
################################################################################

@doc Markdown.doc"""
      pseudo_basis(a::NfRelOrdIdl{T, S}) -> Vector{Tuple{NumFieldElem{T}, S}}
      pseudo_basis(a::NfRelOrdFracIdl{T, S}) -> Vector{Tuple{NumFieldElem{T}, S}}

Returns the pseudo-basis of $a$.
"""
function pseudo_basis(a::Union{NfRelOrdIdl, NfRelOrdFracIdl}; copy::Bool = true)
  assure_has_pseudo_basis(a)
  if copy
    return deepcopy(a.pseudo_basis)
  else
    return a.pseudo_basis
  end
end

@doc Markdown.doc"""
      basis_pmatrix(a::NfRelOrdIdl) -> PMat
      basis_pmatrix(a::NfRelOrdFracIdl) -> PMat

Returns the basis pseudo-matrix of $a$.
"""
function basis_pmatrix(a::Union{NfRelOrdIdl, NfRelOrdFracIdl}; copy::Bool = true)
  assure_has_basis_pmatrix(a)
  if copy
    return deepcopy(a.basis_pmatrix)
  else
    return a.basis_pmatrix
  end
end

################################################################################
#
#  (Inverse) basis matrix
#
################################################################################

@doc Markdown.doc"""
      basis_matrix(a::NfRelOrdIdl{T, S}) -> Generic.Mat{T}
      basis_matrix(a::NfRelOrdFracIdl{T, S}) -> Generic.Mat{T}

Returns the basis matrix of $a$.
"""
function basis_matrix(a::Union{NfRelOrdIdl, NfRelOrdFracIdl}; copy::Bool = true)
  assure_has_basis_matrix(a)
  if copy
    return deepcopy(a.basis_matrix)
  else
    return a.basis_matrix
  end
end

@doc Markdown.doc"""
      basis_mat_inv(a::NfRelOrdIdl{T, S}) -> Generic.Mat{T}
      basis_mat_inv(a::NfRelOrdFracIdl{T, S}) -> Generic.Mat{T}

Returns the inverse of the basis matrix of $a$.
"""
function basis_mat_inv(a::Union{NfRelOrdIdl, NfRelOrdFracIdl}; copy::Bool = true)
  assure_has_basis_mat_inv(a)
  if copy
    return deepcopy(a.basis_mat_inv)
  else
    return a.basis_mat_inv
  end
end

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

function show(io::IO, s::NfRelOrdIdlSet)
  print(io, "Set of ideals of ")
  print(io, s.order)
end

function show(io::IO, a::NfRelOrdIdl)
  compact = get(io, :compact, false)
  if compact
    print(io, "Ideal with basis pseudo-matrix\n")
    show(IOContext(io, :compact => true), basis_pmatrix(a, copy = false))
  else
    print(io, "Ideal of\n")
    show(IOContext(io, :compact => true), order(a))
    print(io, "\nwith basis pseudo-matrix\n")
    show(IOContext(io, :compact => true), basis_pmatrix(a, copy = false))
  end
end

################################################################################
#
#  Parent object overloading and user friendly constructors
#
################################################################################

function defines_ideal(O::NfRelOrd{T, S}, M::PMat{T, S}) where {T, S}
  K = base_field(nf(O))
  coeffs = basis_pmatrix(O, copy = false).coeffs
  I = PseudoMatrix(identity_matrix(K, degree(O)), deepcopy(coeffs))
  return _spans_subset_of_pseudohnf(M, I, :lowerleft)
end

@doc Markdown.doc"""
    ideal(O::NfRelOrd, M::PMat, check::Bool = true, M_in_hnf::Bool = false) -> NfRelOrdIdl

Creates the ideal of $\mathcal O$ with basis pseudo-matrix $M$. If check is set,
then it is checked whether $M$ defines an ideal. If M_in_hnf is set, then it is
assumed that $M$ is already in lower left pseudo HNF.
"""
function ideal(O::NfRelOrd{T, S}, M::PMat{T, S}, check::Bool = true, M_in_hnf::Bool = false) where {T, S}
  if check
    !defines_ideal(O, M) && error("The pseudo-matrix does not define an ideal.")
  end
  !M_in_hnf ? M = pseudo_hnf(M, :lowerleft, true) : nothing
  return NfRelOrdIdl{T, S}(O, M)
end

@doc Markdown.doc"""
    ideal(O::NfRelOrd, M::Generic.Mat, check::Bool = true) -> NfRelOrdIdl

Creates the ideal of $\mathcal O$ with basis matrix $M$. If check is set,
then it is checked whether $M$ defines an ideal.
"""
function ideal(O::NfRelOrd{T, S}, M::Generic.Mat{T}, check::Bool = true) where {T, S}
  coeffs = deepcopy(basis_pmatrix(O, copy = false).coeffs)
  return ideal(O, PseudoMatrix(M, coeffs), check)
end

@doc Markdown.doc"""
    ideal(O::NfRelOrd{T, S}, x::NfRelElem{T}, y::NfRelElem{T}, a::S, b::S, check::Bool = true) -> NfRelOrdIdl{T, S}

Creates the ideal $x\cdot a + y\cdot b$ of $\mathcal O$. If check is set,
then it is checked whether these elements define an ideal.
"""
function ideal(O::NfRelOrd{T, S}, x::NumFieldElem{T}, y::NumFieldElem{T}, a::S, b::S, check::Bool = true) where {T, S}
  d = degree(O)
  pb = pseudo_basis(O, copy = false)
  M = zero_matrix(base_field(nf(O)), 2*d, d)
  C = Array{S}(undef, 2*d)
  for i = 1:d
    elem_to_mat_row!(M, i, pb[i][1]*x)
    C[i] = pb[i][2]*a
  end
  for i = (d + 1):2*d
    elem_to_mat_row!(M, i, pb[i - d][1]*y)
    C[i] = pb[i - d][2]*b
  end
  M = M*basis_mat_inv(O, copy = false)
  PM = PseudoMatrix(M, C)
  if check
    !defines_ideal(O, PM) && error("The elements do not define an ideal.")
  end
  PM = sub(pseudo_hnf(PM, :lowerleft), (d + 1):2*d, 1:d)
  return NfRelOrdIdl{T, S}(O, PM)
end

function ideal(O::NfRelOrd{T, S}, x::NumFieldElem{T}, y::NumFieldElem{T}, a::NfOrdIdl, b::NfOrdIdl, check::Bool = true) where {T, S}
  aa = frac_ideal(order(a), a, fmpz(1))
  bb = frac_ideal(order(b), b, fmpz(1))
  return ideal(O, x, y, aa, bb, check)
end

function ideal(O::NfRelOrd{T, S}, x::NumFieldElem{T}, y::NumFieldElem{T}, a::NfRelOrdIdl, b::NfRelOrdIdl, check::Bool = true) where {T, S}
  aa = frac_ideal(order(a), basis_pmatrix(a), true)
  bb = frac_ideal(order(b), basis_pmatrix(b), true)
  return ideal(O, x, y, aa, bb, check)
end

@doc Markdown.doc"""
    ideal(O::NfRelOrd{T, S}, x::NfRelOrdElem{T}) -> NfRelOrdIdl{T, S}
    *(O::NfRelOrd{T, S}, x::NfRelOrdElem{T}) -> NfRelOrdIdl{T, S}
    *(x::NfRelOrdElem{T}, O::NfRelOrd{T, S}) -> NfRelOrdIdl{T, S}

Creates the ideal $x\cdot \mathcal O$ of $\mathcal O$.
"""
function ideal(O::NfRelOrd{T, S}, x::NfRelOrdElem) where {T, S}
  x = O(x)
  d = degree(O)
  pb = pseudo_basis(O, copy = false)
  M = zero_matrix(base_field(nf(O)), d, d)
  if iszero(x)
    return NfRelOrdIdl{T, S}(O, PseudoMatrix(M, [ deepcopy(pb[i][2]) for i = 1:d ]))
  end
  for i = 1:d
    elem_to_mat_row!(M, i, pb[i][1]*nf(O)(x))
  end
  M = M*basis_mat_inv(O, copy = false)
  PM = PseudoMatrix(M, [ deepcopy(pb[i][2]) for i = 1:d ])
  PM = pseudo_hnf(PM, :lowerleft)
  return NfRelOrdIdl{T, S}(O, PM)
end

function ideal(O::NfRelOrd, x::Union{ Int, fmpz, NfOrdElem })
  return ideal(O, O(x))
end

*(O::NfRelOrd, x::T) where { T <: Union{ Int, fmpz, NfOrdElem, NfRelOrdElem } } = ideal(O, x)

*(x::T, O::NfRelOrd) where { T <: Union{ Int, fmpz, NfOrdElem, NfRelOrdElem } } = ideal(O, x)

@doc Markdown.doc"""
    ideal(O::NfRelOrd{T, S}, a::S, check::Bool = true) -> NfRelOrdIdl{T, S}

Creates the ideal $a \cdot \mathcal O$ of $\mathcal O$. If check is set,
then it is checked whether $a$ defines an (integral) ideal.
"""
function ideal(O::NfRelOrd{T, S}, a::S, check::Bool = true) where {T, S}
  d = degree(O)
  pb = pseudo_basis(O, copy = false)
  M = identity_matrix(base_field(nf(O)), d)
  PM = PseudoMatrix(M, [ a*pb[i][2] for i = 1:d ])
  if check
    !defines_ideal(O, PM) && error("The coefficient ideal does not define an ideal.")
  end
  PM = pseudo_hnf(PM, :lowerleft)
  return NfRelOrdIdl{T, S}(O, PM)
end

function ideal(O::NfRelOrd{nf_elem, NfOrdFracIdl}, a::NfOrdIdl, check::Bool = true)
  aa = frac_ideal(order(a), a, fmpz(1))
  return ideal(O, aa, check)
end

function ideal(O::NfRelOrd, a::NfRelOrdIdl, check::Bool = true)
  @assert order(a) == order(pseudo_basis(O, copy = false)[1][2])

  aa = frac_ideal(order(a), basis_pmatrix(a), true)
  return ideal(O, aa, check)
end

@doc Markdown.doc"""
    *(O::NfRelOrd{T, S}, a::S) -> NfRelOrdIdl{T, S}
    *(a::S, O::NfRelOrd{T, S}) -> NfRelOrdIdl{T, S}

Creates the ideal $a \cdot \mathcal O$ of $\mathcal O$.
"""
*(O::NfRelOrd{T, S}, a::S) where {T, S} = ideal(O, a)

*(a::S, O::NfRelOrd{T, S}) where {T, S} = ideal(O, a)

*(O::NfRelOrd, a::Union{NfOrdIdl, NfRelOrdIdl}) = ideal(O, a)

*(a::Union{NfOrdIdl, NfRelOrdIdl}, O::NfRelOrd) = ideal(O, a)

################################################################################
#
#  Deepcopy
#
################################################################################

function Base.deepcopy_internal(a::NfRelOrdIdl{T, S}, dict::IdDict) where {T, S}
  z = NfRelOrdIdl{T, S}(a.order)
  for x in fieldnames(typeof(a))
    if x != :order && x != :parent && isdefined(a, x)
      setfield!(z, x, Base.deepcopy_internal(getfield(a, x), dict))
    end
  end
  z.order = a.order
  z.parent = a.parent
  return z
end

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

function copy(a::NfRelOrdIdl)
  return a
end

################################################################################
#
#  Equality
#
################################################################################

@doc Markdown.doc"""
    ==(a::NfRelOrdIdl, b::NfRelOrdIdl) -> Bool

Returns whether $a$ and $b$ are equal.
"""
function ==(a::NfRelOrdIdl, b::NfRelOrdIdl)
  order(a) !== order(b) && return false
  return basis_pmatrix(a, copy = false) == basis_pmatrix(b, copy = false)
end

################################################################################
#
#  iszero/isone
#
################################################################################

iszero(a::NfRelOrdIdl) = iszero(basis_matrix(a, copy = false)[1, 1])

isone(a::NfRelOrdIdl) = isone(minimum(a))

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

# Assumes, that det(basis_matrix(a)) == 1
function assure_has_norm(a::NfRelOrdIdl{T, S}) where {T, S}
  if a.has_norm
    return nothing
  end
  if iszero(a)
    O = order(basis_pmatrix(a, copy = false).coeffs[1])
    a.norm = O()*O
    a.has_norm = true
    return nothing
  end
  c = basis_pmatrix(a, copy = false).coeffs
  d = inv_coeff_ideals(order(a), copy = false)
  n = c[1]*d[1]
  for i = 2:degree(order(a))
    n *= c[i]*d[i]
  end
  if T == nf_elem
    simplify(n)
    @assert n.den == 1
    a.norm = n.num
  else
    @assert denominator(n) == 1
    a.norm = ideal_type(order(n))(order(n), basis_pmatrix(n, copy = false))
  end
  a.has_norm = true
  return nothing
end

@doc Markdown.doc"""
    norm(a::NfRelOrdIdl) -> NfOrdIdl

Returns the norm of $a$.
"""
function norm(a::NfRelOrdIdl; copy::Bool = true)
  assure_has_norm(a)
  if copy
    return deepcopy(a.norm)
  else
    return a.norm
  end
end

function norm(a::NfRelOrdIdl, k::Union{ NfRel, AnticNumberField, NfRel_ns })
  n = norm(a)
  while nf(order(n)) != k
    n = norm(n)
  end
  return n
end

function norm(a::NfRelOrdIdl, k::FlintRationalField)
  n = norm(a)
  while !(n isa fmpz)
    n = norm(n)
  end
  return n
end

function absolute_norm(a::NfRelOrdIdl)
  return norm(a, FlintQQ)
end

################################################################################
#
#  Ideal addition / GCD
#
################################################################################

@doc Markdown.doc"""
    +(a::NfRelOrdIdl, b::NfRelOrdIdl) -> NfRelOrdIdl

Returns $a + b$.
"""
function +(a::NfRelOrdIdl{T, S}, b::NfRelOrdIdl{T, S}) where {T, S}
  check_parent(a, b)
  d = degree(order(a))
  H = vcat(basis_pmatrix(a), basis_pmatrix(b))
  if T == nf_elem
    m = norm(a) + norm(b)
    H = sub(pseudo_hnf_full_rank_with_modulus(H, m, :lowerleft), (d + 1):2*d, 1:d)
  else
    H = sub(pseudo_hnf(H, :lowerleft), (d + 1):2*d, 1:d)
  end
  return ideal(order(a), H, false, true)
end

################################################################################
#
#  Ideal multiplication
#
################################################################################

@doc Markdown.doc"""
    *(a::NfRelOrdIdl, b::NfRelOrdIdl) -> NfRelOrdIdl

Returns $a \cdot b$.
"""
function *(a::NfRelOrdIdl{T, S}, b::NfRelOrdIdl{T, S}) where {T, S}
  check_parent(a, b)
  if iszero(a) || iszero(b)
    return order(a)()*order(a)
  end
  pba = pseudo_basis(a, copy = false)
  pbb = pseudo_basis(b, copy = false)
  ma = basis_matrix(a, copy = false)
  mb = basis_matrix(b, copy = false)
  L = nf(order(a))
  K = base_field(L)
  d = degree(order(a))
  M = zero_matrix(K, d^2, d)
  C = Array{frac_ideal_type(order_type(K)), 1}(undef, d^2)
  t = L()
  for i = 1:d
    for j = 1:d
      t = mul!(t, pba[i][1], pbb[j][1])
      elem_to_mat_row!(M, (i - 1)*d + j, t)
      C[(i - 1)*d + j] = simplify(pba[i][2]*pbb[j][2])
    end
  end
  PM = PseudoMatrix(M, C)
  PM.matrix = PM.matrix*basis_mat_inv(order(a), copy = false)
  if T == nf_elem
    m = norm(a)*norm(b)
    H = sub(pseudo_hnf_full_rank_with_modulus(PM, m, :lowerleft), (d*(d - 1) + 1):d^2, 1:d)
  else
    H = sub(pseudo_hnf(PM, :lowerleft), (d*(d - 1) + 1):d^2, 1:d)
  end
  return ideal(order(a), H, false, true)
end

Base.:(^)(A::NfRelOrdIdl, e::Int) = Base.power_by_squaring(A, e)

################################################################################
#
#  Ad hoc multiplication
#
################################################################################

@doc Markdown.doc"""
    *(a:NfRelOrdIdl{T, S}, x::T) -> NfRelOrdIdl{T, S}

Returns the ideal $x\cdot a$.
"""
function *(a::NfRelOrdIdl{T, S}, x::T) where {T, S}
  if iszero(x)
    return order(a)()*order(a)
  end

  return ideal(order(a), x*basis_pmatrix(a), true, true)
end

*(x::T, a::NfRelOrdIdl{T, S}) where {T, S} = a*x

function *(a::Union{NfRelOrdIdl, NfRelOrdFracIdl}, x::Union{ Int, fmpz })
  if iszero(x)
    return order(a)()*order(a)
  end

  return typeof(a)(order(a), x*basis_pmatrix(a))
end

*(x::Union{ Int, fmpz}, a::Union{NfRelOrdIdl, NfRelOrdFracIdl}) = a*x

################################################################################
#
#  Intersection / LCM
#
################################################################################

@doc Markdown.doc"""
    intersect(a::NfRelOrdIdl, b::NfRelOrdIdl) -> NfRelOrdIdl

Returns $a \cap b$.
"""
function intersect(a::NfRelOrdIdl{T, S}, b::NfRelOrdIdl{T, S}) where {T, S}
  check_parent(a, b)
  d = degree(order(a))
  Ma = basis_pmatrix(a)
  Mb = basis_pmatrix(b)
  M1 = hcat(Ma, deepcopy(Ma))
  z = zero_matrix(base_ring(Ma.matrix), d, d)
  M2 = hcat(PseudoMatrix(z, Mb.coeffs), Mb)
  M = vcat(M1, M2)
  if T == nf_elem
    m = intersect(norm(a), norm(b))
    H = sub(pseudo_hnf_full_rank_with_modulus(M, m, :lowerleft), 1:d, 1:d)
  else
    H = sub(pseudo_hnf(M, :lowerleft), 1:d, 1:d)
  end
  return ideal(order(a), H, false, true)
end

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

@doc Markdown.doc"""
      inv(a::NfRelOrdIdl) -> NfRelOrdFracIdl
      inv(a::NfRelOrdFracIdl) -> NfRelOrdFracIdl

Computes the inverse of $a$, that is, the fractional ideal $b$ such that
$ab = O$, where $O$ is the ambient order of $a$. $O$ must be maximal.
"""
function inv(a::Union{NfRelOrdIdl{T, S}, NfRelOrdFracIdl{T, S}}) where {T, S}
  if !ismaximal(order(a))
    error("Not implemented (yet).")
  end
  @assert !iszero(a)
  O = order(a)
  d = degree(O)
  pb = pseudo_basis(a, copy = false)
  bmO = basis_matrix(O, copy = false)
  bmOinv = basis_mat_inv(O, copy = false)
  M = bmO*representation_matrix(pb[1][1])*bmOinv
  for i = 2:d
    M = hcat(M, bmO*representation_matrix(pb[i][1])*bmOinv)
  end
  invcoeffs = inv_coeff_ideals(O, copy = false)
  C = Array{S}(undef, d^2)
  for i = 1:d
    for j = 1:d
      C[(i - 1)*d + j] = simplify(pb[i][2]*invcoeffs[j])
    end
  end
  PM = PseudoMatrix(transpose(M), C)
  PM = sub(pseudo_hnf(PM, :upperright, true), 1:d, 1:d)
  N = inv(transpose(PM.matrix))
  PN = PseudoMatrix(N, [ simplify(inv(I)) for I in PM.coeffs ])
  PN = pseudo_hnf(PN, :lowerleft, true)
  return frac_ideal(O, PN, true)
end

################################################################################
#
#  Division
#
################################################################################

@doc Markdown.doc"""
      divexact(a::NfRelOrdIdl, b::NfRelOrdIdl) -> NfRelOrdFracIdl

Returns $ab^{-1}$.
"""
function divexact(a::NfRelOrdIdl{T, S}, b::NfRelOrdIdl{T, S}) where {T, S}
  O = order(a)
  return frac_ideal(O, basis_pmatrix(a, copy = false), true)*inv(b)
end

//(a::NfRelOrdIdl{T, S}, b::NfRelOrdIdl{T, S}) where {T, S} = divexact(a, b)

################################################################################
#
#  P-radical
#
################################################################################

# Returns an element x with v_p(x) = v_p(a) for all p in primes.
function element_with_valuation(a::T, primes::Vector{T}) where {T <: Union{NfOrdIdl, NfRelOrdIdl}}
  products = Vector{T}()
  for p in primes
    push!(products, a*p)
  end
  foundOne = false
  x = order(a)()
  while !foundOne
    x = rand(a, 2^61) # magic number
    foundOne = true
    for p in products
      if x in p
        foundOne = false
        break
      end
    end
  end
  return x
end

# Algorithm V.8. and VI.8. in "Berechnung relativer Ganzheitsbasen mit dem
# Round-2-Algorithmus" by C. Friedrichs.
@doc Markdown.doc"""
      pradical(O::NfRelOrd, P::NfOrdIdl) -> NfRelOrdIdl

Given a prime ideal $P$, this function returns the $P$-radical
$\sqrt{P\mathcal O}$ of $\mathcal O$, which is
just $\{ x \in \mathcal O \mid \exists k \in \mathbf Z_{\geq 0} \colon x^k
\in P\mathcal O \}$. It is not checked that $P$ is prime.
"""
function pradical(O::NfRelOrd, P::Union{NfOrdIdl, NfRelOrdIdl})
  d = degree(O)
  L = nf(O)
  K = base_field(L)
  OK = maximal_order(K)
  pb = pseudo_basis(O, copy = false)

  is_absolute = (typeof(K) == AnticNumberField)

  # Compute a pseudo basis of O with integral ideals:
  basis_mat_int = zero_matrix(K, d, d)
  pbint = Vector{Tuple{elem_type(L), typeof(P)}}()
  for i = 1:d
    t = divexact(pb[i][1], denominator(pb[i][2]))
    if is_absolute
      push!(pbint, (t, deepcopy(numerator(pb[i][2]))))
    else
      push!(pbint, (t, numerator(pb[i][2])))
    end
    elem_to_mat_row!(basis_mat_int, i, t)
  end
  if is_absolute
    Oint = typeof(O)(L, PseudoMatrix(basis_mat_int, [ frac_ideal(OK, pbint[i][2], fmpz(1)) for i = 1:d ]))
  else
    Oint = typeof(O)(L, PseudoMatrix(basis_mat_int, [ frac_ideal(OK, basis_pmatrix(pbint[i][2], copy = false)) for i = 1:d ]))
  end

  if is_absolute
    pOK = ideal(OK, OK(minimum(P)))
  else
    pOK = minimum(P, copy = false)*OK
  end
  prime_ideals = factor(pOK)

  elts_with_val = Vector{elem_type(OK)}(undef, d)
  for i = 1:d
    elts_with_val[i] = element_with_valuation(pbint[i][2], [ p for (p, e) in prime_ideals ])
  end
  F, mF = ResidueField(OK, P)
  mmF = extend(mF, K)
  A = zero_matrix(F, d, d)

  # If the prime number in P is too small one can't use the trace.
  if is_absolute
    p = minimum(P)
  else
    p = prime_number(P)
  end
  if p <= d
    q = order(F)
    k = clog(fmpz(degree(Oint)), q)
    for i = 1:d
      t = Oint((L(K(elts_with_val[i]))*pbint[i][1])^(q^k))
      ar = coordinates(t)
      for j = 1:d
        A[j, i] = mmF(divexact(ar[j], K(elts_with_val[j])))
      end
    end
  else
    for i = 1:d
      for j = i:d
        t = L(K(elts_with_val[i]))*pbint[i][1]*L(K(elts_with_val[j]))*pbint[j][1]
        A[i, j] = mF(OK(tr(t)))
        A[j, i] = deepcopy(A[i, j])
      end
    end
  end

  B = nullspace(A)[2]
  M1 = zero_matrix(K, d, d)
  imF = pseudo_inv(mF)
  # Write a basis of the kernel of A in the rows of M1.
  for i = 1:ncols(B)
    for j = 1:nrows(B)
      M1[i, j] = K(imF(B[j, i])*elts_with_val[j])
    end
  end

  M1 = zero_matrix(K, d, d)
  imF = pseudo_inv(mF)
  # Write a basis of the kernel of A in the rows of M1.
  for j = 1:nrows(B)
    t = K(denominator(pb[j][2], copy = false))
    for i = 1:ncols(B)
      M1[i, j] = divexact(K(imF(B[j, i])*elts_with_val[j]), t)
    end
  end
  PM1 = PseudoMatrix(M1)
  PM2 = PseudoMatrix(identity_matrix(K, d), [ pb[i][2]*deepcopy(P) for i = 1:d ])
  m = det(PM2)
  PM = sub(pseudo_hnf_full_rank_with_modulus(vcat(PM1, PM2), numerator(m, copy = false), :lowerleft), (d + 1):2*d, 1:d)

  return ideal(O, PM, false, true)
end

################################################################################
#
#  Ring of multipliers
#
################################################################################

# Algorithm VII.1. in "Berechnung relativer Ganzheitsbasen mit dem
# Round-2-Algorithmus" by C. Friedrichs.
@doc Markdown.doc"""
    ring_of_multipliers(a::NfRelOrdIdl) -> NfRelOrd

Computes the order $(a : a)$, which is the set of all $x \in K$
with $xa \subseteq a$, where $K$ is the ambient number field
of $a$.
"""
function ring_of_multipliers(a::NfRelOrdIdl{T1, T2}) where {T1, T2}
  O = order(a)
  K = base_field(nf(O))
  d = degree(O)
  pb = pseudo_basis(a, copy = false)
  S = basis_mat_inv(O, copy = false)*basis_mat_inv(a, copy = false)
  M = basis_matrix(O, copy = false)*representation_matrix(pb[1][1])*S
  for i = 2:d
    M = hcat(M, basis_matrix(O, copy = false)*representation_matrix(pb[i][1])*S)
  end
  invcoeffs = [ simplify(inv(pb[i][2])) for i = 1:d ]
  C = Array{T2}(undef, d^2)
  for i = 1:d
    for j = 1:d
      if i == j
        C[(i - 1)*d + j] = K(1)*order(pb[i][2])
      else
        C[(i - 1)*d + j] = simplify(pb[i][2]*invcoeffs[j])
      end
    end
  end
  PM = PseudoMatrix(transpose(M), C)
  PM = sub(pseudo_hnf(PM, :upperright, true), 1:d, 1:d)
  N = inv(transpose(PM.matrix))*basis_matrix(O, copy = false)
  PN = PseudoMatrix(N, [ simplify(inv(I)) for I in PM.coeffs ])
  return NfRelOrd{T1, T2}(nf(O), PN)
end

################################################################################
#
#  Absolute to relative
#
################################################################################

function relative_ideal(a::NfOrdIdl, m::NfToNfRel)
  L = codomain(m)
  Labs = domain(m)
  @assert nf(order(a)) == Labs
  K = base_field(L)
  O = relative_order(order(a), m)
  B = basis(a, copy = false)
  d = degree(L)
  dabs = degree(Labs)
  M = zero_matrix(K, dabs, d)
  for i = 1:dabs
    elem_to_mat_row!(M, i, m(Labs(B[i])))
  end
  M = M*basis_mat_inv(O, copy = false)
  PM = sub(pseudo_hnf(PseudoMatrix(M), :lowerleft, true), (dabs - d + 1):dabs, 1:d)
  return ideal(O, PM, false, true)
end

################################################################################
#
#  Index divisors
#
################################################################################

function isindex_divisor(O::NfRelOrd, p::Union{NfOrdIdl, NfRelOrdIdl})
  f = nf(O).pol
  return valuation(discriminant(f), p) != valuation(discriminant(O), p)
end

################################################################################
#
#  Prime decomposition
#
################################################################################

function prime_decomposition(O::NfRelOrd, p::Union{NfOrdIdl, NfRelOrdIdl})
  if isindex_divisor(O, p)
    return prime_dec_index(O, p)
  end

  return prime_dec_nonindex(O, p)
end

function prime_dec_nonindex(O::NfRelOrd, p::Union{NfOrdIdl, NfRelOrdIdl})
  L = nf(O)
  OK = order(p)
  @assert OK == O.basis_pmatrix.coeffs[1].order
  @assert OK.ismaximal == 1
  a = gen(L)
  K = base_field(L)
  f = L.pol

  Kx = parent(f)
  Fp, mF = ResidueField(OK, p)
  mmF = extend(mF, K)
  immF = pseudo_inv(mmF)
  Fy, y = PolynomialRing(Fp,"y", cached=false)
  fmodp = Hecke.nf_elem_poly_to_fq_poly(Fy, mmF, f)
  fac = factor(fmodp)
  result = Array{Tuple{ideal_type(O), Int}, 1}()
  for (q, e) in fac
    g = Hecke.fq_poly_to_nf_elem_poly(Kx, immF, q)
    ga = g(a)
    P = ideal(O, L(1), ga, p, OK(1)*OK)
    P.is_prime = 1
    P.splitting_type = (e, degree(q))
    P.minimum = deepcopy(p)
    P.non_index_div_poly = q
    Oga = O(ga)
    # TODO: Warum funktioniert das? Muss uniformizer(p) ein p-uniformizer sein?
    if iszero(Oga)
      @assert e == 1
      P.p_uniformizer = O(uniformizer(p))
    else
      if e != 1
        P.p_uniformizer = Oga
      else
        if valuation(Oga, P) == 1
          P.p_uniformizer = Oga
        else
          P.p_uniformizer = Oga + O(uniformizer(p))
        end
      end
    end
    push!(result, (P, e))
  end
  return result
end

function prime_dec_index(O::NfRelOrd, p::Union{NfOrdIdl, NfRelOrdIdl})
  L = nf(O)
  K = base_field(L)
  pbasisO = pseudo_basis(O, copy = false)
  pO = p*O

  Ip = pradical(O, p)
  A, OtoA = AlgAss(O, Ip, p)
  AtoO = pseudo_inv(OtoA)
  AA = decompose(A)

  result = Vector{Tuple{ideal_type(O), Int}}()
  m = PseudoMatrix(zero_matrix(K, 1, degree(O)))
  for (B, BtoA) in AA
    f = dim(B)
    idem = BtoA(B[1]) # Assumes that B == idem*A
    M = representation_matrix(idem)
    ker = left_kernel_basis(M)
    N = basis_pmatrix(Ip)
    for i = 1:length(ker)
      b = coordinates(AtoO(A(ker[i])))
      for j = 1:degree(O)
        m.matrix[1, j] = b[j]
      end
      N = vcat(N, deepcopy(m))
    end
    N = sub(pseudo_hnf(N, :lowerleft), nrows(N) - degree(O) + 1:nrows(N), 1:degree(O))
    P = ideal(O, N, false, true)
    P.is_prime = 1
    e = valuation(pO, P)
    P.splitting_type = (e, f)
    P.minimum = deepcopy(p)
    push!(result, (P, e))
  end

  return result
end

# Returns all prime ideals in O containing the prime number p
function prime_ideals_over(O::NfRelOrd, p::Union{ Int, fmpz })
  pdec = prime_ideals_over(base_ring(O), p)

  primes = Vector{ideal_type(O)}()
  for q in pdec
    qdec = prime_decomposition(O, q)
    append!(primes, [ qdec[i][1] for i = 1:length(qdec) ])
  end

  return primes
end

################################################################################
#
#  Reduction of element modulo ideal
#
################################################################################

function mod!(a::NfRelOrdElem, I::NfRelOrdIdl)
  O = order(I)
  b = coordinates(a, copy = false)
  PM = basis_pmatrix(I, copy = false) # PM is assumed to be in lower left pseudo hnf
  t = parent(b[1])()
  t1 = parent(b[1])()
  for i = degree(O):-1:1
    t = add!(t, mod(b[i], PM.coeffs[i]), -b[i])
    for j = 1:i
      t1 = mul!(t1, t, PM.matrix[i, j])
      b[j] = add!(b[j], b[j], t1)
    end
  end

  t = nf(O)()
  B = basis_nf(O, copy = false)
  zero!(a.elem_in_nf)
  for i = 1:degree(O)
    t = mul!(t, B[i], nf(O)(b[i]))
    a.elem_in_nf = add!(a.elem_in_nf, a.elem_in_nf, t)
  end

  return a
end

function mod(a::NfRelOrdElem, I::NfRelOrdIdl)
  return mod!(deepcopy(a), I)
end

function mod!(a::NfRelOrdElem, Q::RelOrdQuoRing)
  return mod!(a, ideal(Q))
end

function mod(a::NfRelOrdElem, Q::RelOrdQuoRing)
  return mod(a, ideal(Q))
end

################################################################################
#
#  Valuation
#
################################################################################

function valuation_naive(A::NfRelOrdIdl{T, S}, B::NfRelOrdIdl{T, S}) where {T, S}
  @assert order(A.basis_pmatrix.coeffs[1]) == order(B.basis_pmatrix.coeffs[1])
  @assert !iszero(A) && !iszero(B)
  O = order(A)
  Afrac = frac_ideal(O, basis_pmatrix(A), true)
  Bi = inv(B)
  i = 0
  C = Afrac*Bi
  @assert C != Afrac
  while isintegral(C)
    C = C*Bi
    i += 1
  end
  return i
end

valuation(A::NfRelOrdIdl{T, S}, B::NfRelOrdIdl{T, S}) where {T, S} = valuation_naive(A, B)

function valuation_naive(a::NfRelOrdElem{T}, B::NfRelOrdIdl{T, S}) where {T, S}
  @assert !iszero(a)
  @assert order(parent(a).basis_pmatrix.coeffs[1]) == order(B.basis_pmatrix.coeffs[1])
  @assert order((a * parent(a)).basis_pmatrix.coeffs[1]) == order(B.basis_pmatrix.coeffs[1])
  return valuation(a*parent(a), B)
end

valuation(a::NfRelOrdElem{T}, B::NfRelOrdIdl{T, S}) where {T, S} = valuation_naive(a, B)

valuation(a::fmpz, B::NfRelOrdIdl) = valuation(order(B)(a), B)

################################################################################
#
#  Factorization into prime ideals
#
################################################################################

function factor(A::NfRelOrdIdl{T, S}) where {T, S}
  n = norm(A)
  normFactors = factor(n)
  result = Dict{NfRelOrdIdl{T, S}, Int}()
  O = order(A)
  for p in keys(normFactors)
    prime_dec = prime_decomposition(O, p)
    for (P, e) in prime_dec
      v = valuation(A, P)
      if v != 0
        result[P] = v
        n = n//norm(P)^v
        simplify(n)
      end
      if isone(n)
        return result
      end
    end
  end
  return result
end

################################################################################
#
#  Minimum
#
################################################################################

@doc Markdown.doc"""
      minimum(A::NfRelOrdIdl) -> NfOrdIdl
      minimum(A::NfRelOrdIdl) -> NfRelOrdIdl

Returns the ideal $A \cap O$ where $O$ is the maximal order of the coefficient
ideals of $A$.
"""
function minimum(A::NfRelOrdIdl; copy::Bool = true)
  assure_has_minimum(A)
  if copy
    return deepcopy(A.minimum)
  else
    return A.minimum
  end
end

function assure_has_minimum(A::NfRelOrdIdl)
  if isdefined(A, :minimum)
    return nothing
  end
  @assert isone(basis_pmatrix(A, copy = false).matrix[1, 1])
  @assert isone(basis_pmatrix(order(A), copy = false).matrix[1, 1])

  M = deepcopy(basis_pmatrix(A, copy = false).coeffs[1])
  M = simplify(M)
  A.minimum = numerator(M)
  return nothing
end

################################################################################
#
#  Order modulo prime ideal
#
################################################################################

function ResidueField(O::NfRelOrd{T, S}, P::NfRelOrdIdl{T, S}) where {T, S}
  @assert order(P) == O
  @assert P.is_prime == 1
  mF = NfRelOrdToFqMor{T, S}(O, P)
  return codomain(mF), mF
end

################################################################################
#
#  Idempotents
#
################################################################################

@doc Markdown.doc"""
    idempotents(x::NfRelOrdIdl, y::NfRelOrdIdl) -> NfRelOrdElem, NfRelOrdElem

Returns a tuple `(e, f)` consisting of elements `e in x`, `f in y` such that
`1 = e + f`.
>
If the ideals are not coprime, an error is raised.
"""
function idempotents(x::NfRelOrdIdl{T, S}, y::NfRelOrdIdl{T, S}) where {T, S}
  check_parent(x, y)

  O = order(x)
  mx = minimum(x, copy = false)
  my = minimum(y, copy = false)
  g = mx + my
  if isone(g)
    u, v = idempotents(mx, my)
    return O(u), O(v)
  end

  d = degree(O)
  L = nf(O)
  K = base_field(L)
  OK = maximal_order(K)
  M = zero_matrix(K, 2*d + 1, 2*d + 1)

  M[1, 1] = K(1)
  z = coordinates(one(O))
  for i = 1:d
    M[1, i + 1] = z[i]
  end
  for i = 1:d
    for j = 1:d
      M[i + 1, j + 1] = deepcopy(basis_matrix(x, copy = false)[i, j])
      M[i + 1 + d, j + 1] = deepcopy(basis_matrix(y, copy = false)[i, j])
    end
    M[i + 1, i + d + 1] = K(1)
  end

  #=
    M is now
    ( 1 |  1  |  0  )
    ( 0 | M_x | I_d )
    ( 0 | M_y |  0  )
  =#

  coeffsx = deepcopy(basis_pmatrix(x, copy = false).coeffs)
  coeffsy = deepcopy(basis_pmatrix(y, copy = false).coeffs)
  C = [ K(1)*OK, coeffsx..., coeffsy... ]
  PM = PseudoMatrix(M, C)
  PM = pseudo_hnf(PM, :upperright)

  for i = 2:(d + 1)
    if !iszero(PM.matrix[1, i])
      error("Ideals are not coprime")
    end
  end

  pbx = pseudo_basis(x, copy = false)
  u = pbx[1][1]*PM.matrix[1, d + 2]
  for i = 2:d
    u += pbx[i][1]*PM.matrix[1, d + 1 + i]
  end

  @assert -u in x
  @assert u + 1 in y

  return O(-u), O(u + 1)
end

################################################################################
#
#  Inclusion of elements in ideals
#
################################################################################

@doc Markdown.doc"""
    in(x::NfRelOrdElem, y::NfRelOrdIdl)
    in(x::NumFieldElem, y::NfRelOrdIdl)
    in(x::fmpz, y::NfRelOrdIdl)

Returns whether $x$ is contained in $y$.
"""
function in(x::NfRelOrdElem, y::NfRelOrdIdl)
  parent(x) !== order(y) && error("Order of element and ideal must be equal")
  O = order(y)
  b_pmat = basis_pmatrix(y, copy = false)
  t = matrix(base_field(nf(O)), 1, degree(O), coordinates(x))
  t = t*basis_mat_inv(y, copy = false)
  for i = 1:degree(O)
    if !(t[1, i] in b_pmat.coeffs[i])
      return false
    end
  end
  return true
end

function in(x::NumFieldElem, y::NfRelOrdIdl)
  parent(x) !== nf(order(y)) && error("Number field of element and ideal must be equal")
  return in(order(y)(x),y)
end

in(x::fmpz, y::NfRelOrdIdl) = in(order(y)(x),y)

################################################################################
#
#  (Anti-)Uniformizer
#
################################################################################

@doc Markdown.doc"""
    uniformizer(P::NfRelOrdIdl) -> NfRelOrdElem

Returns an element $u \in P$ with valuation(u, P) == 1.
"""
function uniformizer(P::NfRelOrdIdl)
  @assert P.is_prime == 1

  if P.splitting_type[1] == 1
    return order(P)(uniformizer(minimum(P, copy = false)))
  end

  r = 500 # hopefully enough
  z = rand(P, r)
  while true
    if !iszero(z) && valuation(z, P) == 1
      break
    end
    z = rand(P, r)
  end
  return z
end

function _is_p_uniformizer(z::NfRelOrdElem, P::T, primes::Vector{T}) where {T <: NfRelOrdIdl}
  if iszero(z)
    return false
  end
  if valuation(z, P) != 1
    return false
  end
  for PP in primes
    if valuation(z, PP) != 0
      return false
    end
  end
  return true
end

@doc Markdown.doc"""
    p_uniformizer(P::NfRelOrdIdl) -> NfRelOrdElem

Returns an element $u \in P$ with valuation(u, P) == 1 and valuation 0 at all
other prime ideals lying over minimum(P).
"""
function p_uniformizer(P::NfRelOrdIdl)
  @assert P.is_prime == 1

  if isdefined(P, :p_uniformizer)
    return P.p_uniformizer
  end

  p = minimum(P, copy = false)
  prime_dec = prime_decomposition(order(P), p)
  primes = Vector{typeof(P)}()
  for (PP, e) in prime_dec
    if PP != P
      push!(primes, PP)
    end
  end
  r = 500
  z = rand(P, r)
  while !_is_p_uniformizer(z, P, primes)
    z = rand(P, r)
  end
  P.p_uniformizer = z
  return z
end

@doc Markdown.doc"""
    anti_uniformizer(P::NfRelOrdIdl) -> NumFieldElem

Returns an element $a$ in the number field containing $P$ with valuation(a, P) == -1
and non-negative valuation at all other prime ideals.
"""
function anti_uniformizer(P::NfRelOrdIdl{T, S}) where {T, S}
  @assert P.is_prime == 1

  if isdefined(P, :anti_uniformizer)
    return P.anti_uniformizer
  end

  p = minimum(P, copy = false)
  # We need a pseudo basis of O, where the coefficient ideals have valuation
  # 0 at p.
  O = order(P)
  N = basis_matrix(O)
  NN = basis_mat_inv(O)
  d = Vector{T}(undef, degree(O))
  a = elem_in_nf(uniformizer(p))
  for i = 1:degree(O)
    v = valuation(pseudo_basis(O, copy = false)[i][2], p)
    if !iszero(v)
      d[i] = a^v
      mul_row!(N, i, d[i])
      mul_col!(NN, i, inv(d[i]))
    else
      d[i] = one(base_field(nf(O)))
    end
  end

  u = elem_in_nf(p_uniformizer(P))
  M = representation_matrix(u)
  M = N*M*NN

  F, mF = ResidueField(order(p), p)
  mmF = extend(mF, nf(order(p)))
  immF = pseudo_inv(mmF)
  Mp = zero_matrix(F, nrows(M), ncols(M))
  for i = 1:nrows(M)
    for j = 1:ncols(M)
      Mp[i, j] = mmF(M[i, j])
    end
  end
  K = left_kernel_basis(Mp)
  @assert length(K) > 0
  x = nf(O)()
  for i = 1:degree(O)
    x += immF(K[1][i])*pseudo_basis(O, copy = false)[i][1]*d[i]
  end
  P.anti_uniformizer = x*anti_uniformizer(p)
  return P.anti_uniformizer
end

################################################################################
#
#  Random elements
#
################################################################################

function rand(a::NfRelOrdIdl, B::Int)
  pb = pseudo_basis(a, copy = false)
  z = nf(order(a))()
  for i = 1:degree(order(a))
    t = rand(pb[i][2], B)
    z += t*pb[i][1]
  end
  return order(a)(z)
end

################################################################################
#
#  Prime number in a prime ideal
#
################################################################################

function prime_number(p::NfRelOrdIdl)
  @assert p.is_prime == 1
  m = minimum(p, copy = false)
  if typeof(m) == NfOrdIdl
    return minimum(m)
  else
    return prime_number(m)
  end
end

################################################################################
#
#  Something with coprime
#
################################################################################

function coprime_to(I::NfRelOrdFracIdl, p::NfRelOrdIdl)
  pi = anti_uniformizer(p)
  a = rand(I, 500)
  l = valuation(a, p)
  @assert l >= 0
  if l > 0
    a = pi^l*a
  end
  @assert valuation(a, p) == 0
  return a
end

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

function Base.hash(A::NfRelOrdIdl, h::UInt)
  return Base.hash(basis_pmatrix(A, copy = false), h)
end

################################################################################
#
#  Approximation
#
################################################################################

# See also approximate_nonnegative and approximate_simple in NfOrd/Ideal/Prime.jl

# Returns x in K such that v_p(x) = v[i] for p = primes[i] and v_p(x) \geq 0 for all other primes p.
# Algorithm 1.7.8 in Hoppe: Normal forms over Dedekind domains
function approximate(v::Vector{Int}, primes::Vector{ <: NfRelOrdIdl })
  @assert length(v) == length(primes)
  @assert length(primes) > 0

  O = order(primes[1])

  # Make the set primes complete: add all prime ideals lying over the same prime numbers
  prime_numbers = Set{fmpz}()
  for p in primes
    push!(prime_numbers, prime_number(p))
  end

  primes2 = Vector{ideal_type(O)}()
  for p in prime_numbers
    pdec = prime_ideals_over(O, p)
    append!(primes2, pdec)
  end

  v2 = zeros(Int, length(primes2))

  D = Dict([ (primes[i], v[i]) for i = 1:length(primes) ])

  for i = 1:length(primes2)
    if haskey(D, primes2[i])
      v2[i] = D[primes2[i]]
    end
  end

  a_pos, a_neg = _approximate_simple(v2, primes2)

  # Take care of the additional negative valuations coming from a_neg^(-1)
  c = fmpq(absolute_norm(a_neg))
  for i = 1:length(primes)
    if v[i] >= 0
      continue
    end

    c *= fmpq(absolute_norm(primes[i]))^v[i]
  end

  return divexact(c*elem_in_nf(a_pos), elem_in_nf(a_neg))
end


text:
#=
Copyright (c) 2015, Intel Corporation

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

    * Redistributions of source code must retain the above copyright notice,
      this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of Intel Corporation nor the names of its contributors
      may be used to endorse or promote products derived from this software
      without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=#

export TransformLayer
using HDF5

type TransformEnsemble <: JuliaEnsemble
    net           :: Net
    name          :: Symbol
    neurons       :: Array{DataNeuron}
    mean          :: Array{Float32}
    random_mirror :: Bool
    scale         :: Float32
    connections   :: Vector{Connection}
    phase         :: Phase
    net_subgroup  :: Cint
end

function TransformLayer(net::Net, name::Symbol,
                        input_ensemble::AbstractEnsemble; crop=false,
                        random_mirror=false, scale=1.0f0, mean_file=nothing)
    @assert ndims(input_ensemble) == 3
    if crop != false
        neurons = Array(DataNeuron, crop..., size(input_ensemble, 3))
    else
        neurons = Array(DataNeuron, size(input_ensemble)...)
    end
    mean = zeros(Float32, size(input_ensemble)...)
    if mean_file != nothing
        h5open(mean_file, "r") do h5
            mean = h5["mean"]
            if ndims(mean) == 4
                mean = mean[:,:,:,1] .* scale
                mean = reshape(mean, size(mean)[1:3])
            else
                mean = mean[:,:,:] .* scale
            end
            @assert size(mean) == size(input_ensemble)
        end
    end
    for i = 1:length(neurons)
        neurons[i] = DataNeuron(0.0f0)
    end
    ens = TransformEnsemble(net, name, neurons, mean, random_mirror, scale, [], TrainTest, convert(Cint, 1))
    add_ensemble(net, ens)
    add_connections(net, input_ensemble, ens, (i, j, k) -> (i:i, j:j, k:k))
    ens
end

function forward(ens::TransformEnsemble, input::Array, value::Array, phase::Phase)
    if phase == Train
        w_off = abs(rand(Int)) % (size(input,1)-size(value, 1)+1)
        h_off = abs(rand(Int)) % (size(input,2)-size(value, 2)+1)
    else
        w_off = div(size(input,1)-size(input, 1), 2)
        h_off = div(size(input,2)-size(input, 2), 2)
    end
    for i in 1:size(input, 4)
        input[:,:,:,i] -= ens.mean
    end
    crop_w, crop_h, channels, num = size(value)
    if rand(UInt) % 2 == 0 || phase == Test
        for n = 1:num
            for c = 1:channels
                for h = 1:crop_h
                    @simd for w = 1:crop_w
                        @inbounds value[w,h,c,n] = input[w+w_off,h+h_off,c,n] * ens.scale
                    end
                end
            end
        end
    else
        for n = 1:num
            for c = 1:channels
                for h = 1:crop_h
                    @simd for w = 1:crop_w
                        @inbounds value[crop_w-w+1,h,c,n] = input[w+w_off,h+h_off,c,n] * ens.scale
                    end
                end
            end
        end
    end
end

function backward(ens::TransformEnsemble, args...)
end