text:
# Variable selection using power likelihood
module VarSel
using Distributions

# w ~ Beta(aw,bw)
a_w(p) = 1.0
b_w(p) = 2*p
# prior on beta's given w: N(0,1/L0) w.p. w, and 0 w.p. 1-w.
sigma = 1.0
L0 = 1/sigma^2

# prior on lambda: Gamma(a,b)
a,b = 1.0,1.0

# induced prior on # of nonzero coefs k: BetaBinomial(p,aw,bw)
log_betabern(k,n,a,b) = lgamma(a+k) + lgamma(b+n-k) - lgamma(a+b+n) - lgamma(a) - lgamma(b) + lgamma(a+b)
log_betabin(k,n,a,b) = lgamma(n+1) - lgamma(k+1) - lgamma(n-k+1) + log_betabern(k,n,a,b)
prior_k(k,p) = exp(log_betabin(k,p,a_w(p),b_w(p)))

# run sampler
function run(x,y,n_total,n_keep,zeta)
    # INPUTS
    # x = n-by-p matrix of covariates
    # y = n-by-1 array of responses
    # n_total = number of MCMC samples
    # n_keep = number of MCMC samples to record
    # zeta = power to raise likelihood to

    n,p = size(x)
    @assert(length(y)==n)
    aw = a_w(p)
    bw = b_w(p)

    # initialize
    if true
        beta = zeros(p)
    else # sample from prior
        w = rand(Beta(aw,bw))
        z = (rand(p).<w)
        beta = z.*randn(p)/sqrt(L0)
    end
    lambda = 1.0
    delta = y - x*beta
    k = sum(beta.!=0)
    sx = vec(sum(x.^2,1))

    # initialize record-keeping
    @assert(n_total>=n_keep)
    keepers = int(linspace(1,n_total,n_keep))
    beta_r = zeros(p,n_keep)
    lambda_r = zeros(n_keep)
    i_keep = 1

    for rep = 1:n_total
        # update lambda
        A = a + 0.5*n*zeta
        B = b + 0.5*sum(delta.^2)*zeta
        lambda = rand(Gamma(A,1/B))

        # update beta
        for j = 1:p
            delta += beta[j]*x[:,j]
            k -= int(beta[j]!=0)
            L = L0 + zeta*lambda*sx[j]
            M = (zeta*lambda/L)*dot(delta,x[:,j])
            log_p0 = log(bw + p-1 - k)  # unnormalized prob of 0
            log_p1 = log(aw + k) + 0.5*log(L0) - 0.5*log(L) + 0.5*L*M*M  # unnormalized prob of nonzero
            p0 = 1/(1 + exp(log_p1 - log_p0))
            if rand() < p0
                beta[j] = 0.0
            else
                beta[j] = randn()/sqrt(L) + M
                delta -= beta[j]*x[:,j]
                k += 1
            end
        end

        # record
        if rep==keepers[i_keep]
            lambda_r[i_keep] = lambda
            beta_r[:,i_keep] = beta
            i_keep += 1
        end
    end

    return beta_r,lambda_r,keepers
end

end # module

max_stars_repo_name: jwmi/CoarsenedPosterior

max_issues_repo_name: jwmi/CoarsenedPosterior

max_forks_repo_name: jwmi/CoarsenedPosterior


text:
@testset "Population" begin

    function test_interface(p, members)
        test_iterator(p, members)

        ep = evolve(NoOpEvolution(), p)
        test_iterator(ep, p)

        @test generation(ep) == generation(p) + 1
    end

    function test_iterator(p1, p2)
        @test length(p1) == length(p2)
        @test size(p1) == size(p2)
        @test collect(p1) == collect(p2)
    end

    test_interface(Population(1:10), 1:10)

    @testset "Fitness calc" begin
        mutable struct PopulationMockFitness <: AbstractFitness
            cnt::Int
        end
        struct PopulationMockCand <: AbstractCandidate
            nr::Int
        end

        function NaiveGAflux._fitness(f::PopulationMockFitness, c::PopulationMockCand)
            f.cnt +=1
            return f.cnt + c.nr
        end
        fpop = fitness(Population(PopulationMockCand.(10:10:100)), PopulationMockFitness(0))
        for (i, c) in enumerate(fpop)
            @test fitness(c) == 10i + i
        end

    end

    @testset "Persistent pop" begin
        testdir = "testPersistentPopulation"
        mems = PersistentArray(testdir, 4, identity)
        pop1 = Population(mems)
        @test generation(pop1) == 1
        pop2 = evolve(NoOpEvolution(), pop1)
        @test generation(pop2) == 2

        try

            persist(pop2)
            pop2 = Population(mems)

            @test generation(pop2) == 2

            pop3 = evolve(NoOpEvolution(), pop2)
            @test generation(pop3) == 3

            persist(pop3)
            pop3 = Population(mems)

            @test generation(pop3) == 3
        finally
            rm(testdir, force=true, recursive=true)
        end
    end
end

ext: jl


text:
function mean(x::Array)
    sum(x) / length(x)
end

function distance(arr::Array, x::Float64)
    mean([(y-x)^2 for y in arr])
end

function argsort(arr::Array)
    sortperm(arr, rev=false)
end

function ordinal_patterns(x::Array, m=3::Int, t=1::Int)
    partition = [x[i:i+(m-1)*t] for i in 1:t:length(x)-(m-1)*t]

    permutation = map(argsort, partition)

    counts = Dict{Array, Float64}()
    sum_w = 0.0
    for perm in permutation
        counts[perm] = get(counts, perm, 0) + 1.0
        sum_w += 1.0
    end

    p = [w / sum_w for w in values(counts)]
end

function weight_patterns(x::Array, m=3::Int, t=1::Int)
    partition = [x[i:i+(m-1)*t] for i in 1:t:length(x)-(m-1)*t]

    X = map(mean, partition)
    weight = map(distance, partition, X)

    permutation = map(argsort, partition)

    counts = Dict{Array, Float64}()
    sum_w = 0.0
    for (perm, w) in zip(permutation, weight)
        counts[perm] = get(counts, perm, 0) + w
        sum_w += w
    end

    p = [w / sum_w for w in values(counts)]
end

function weight_entropy(x::Array, m=3::Int, t=1::Int)
    y = weight_patterns(x, m, t)
    - sum(x * log(x) for x in y)
end

function bandt_and_pompe(x::Array, m=3::Int, t=1::Int)
    y = ordinal_patterns(x, m, t)
    - sum(x * log2(x) for x in y)
end

function bandt_and_pompe_normal(x::Array, m=3::Int, t=1::Int)
    bandt_and_pompe(x, m, t) / log2(factorial(m))
end

function renyi_entropy(x::Array, q=1.::Float, m=3::Int, t=1::Int)
    if q == 1
        return bandt_and_pompe(x, m, t)
    end
    y = ordinal_patterns(x, m, t)
    log2(sum(x^q for x in y)) / (1-q)
end

function complexity_entropy(x::Array, m=3::Int, t=1::Int)
    n = factorial(m)
    op = ordinal_patterns(x, m, t)
    h1 = -sum(x * log(x) for x in op)

    Q = -1 / (
        (0.5+0.5/n)*log(0.5+(0.5/n)) +
        (0.5/n)*log(0.5/n)*(n-1) +
        0.5*log(n)
    )
    op2 = vcat(0.5.*op .+ 0.5/n, [0.5/n for i in 1:(n - length(op))])
    h2 = -sum(x * log(x) for x in op2)

    JSD = h2 .- 0.5.*h1 .- 0.5*log(n)
    comp_JS = Q * JSD * h1/log(n)

    a = [h1/log(n), JSD, comp_JS]
end


text:
const __bodyfunction__ = Dict{Method,Any}()

# Find keyword "body functions" (the function that contains the body
# as written by the developer, called after all missing keyword-arguments
# have been assigned values), in a manner that doesn't depend on
# gensymmed names.
# `mnokw` is the method that gets called when you invoke it without
# supplying any keywords.
function __lookup_kwbody__(mnokw::Method)
    function getsym(arg)
        isa(arg, Symbol) && return arg
        @assert isa(arg, GlobalRef)
        return arg.name
    end

    f = get(__bodyfunction__, mnokw, nothing)
    if f === nothing
        fmod = mnokw.module
        # The lowered code for `mnokw` should look like
        #   %1 = mkw(kwvalues..., #self#, args...)
        #        return %1
        # where `mkw` is the name of the "active" keyword body-function.
        ast = Base.uncompressed_ast(mnokw)
        if isa(ast, Core.CodeInfo) && length(ast.code) >= 2
            callexpr = ast.code[end-1]
            if isa(callexpr, Expr) && callexpr.head == :call
                fsym = callexpr.args[1]
                if isa(fsym, Symbol)
                    f = getfield(fmod, fsym)
                elseif isa(fsym, GlobalRef)
                    if fsym.mod === Core && fsym.name === :_apply
                        f = getfield(mnokw.module, getsym(callexpr.args[2]))
                    elseif fsym.mod === Core && fsym.name === :_apply_iterate
                        f = getfield(mnokw.module, getsym(callexpr.args[3]))
                    else
                        f = getfield(fsym.mod, fsym.name)
                    end
                else
                    f = missing
                end
            else
                f = missing
            end
        else
            f = missing
        end
        __bodyfunction__[mnokw] = f
    end
    return f
end

function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    Base.precompile(Tuple{typeof(get_template!),String,Dict{String, String},Dict{Symbol, Any}})   # time: 0.21933882
    isdefined(GLAbstraction, Symbol("#63#68")) && Base.precompile(Tuple{getfield(GLAbstraction, Symbol("#63#68"))})   # time: 0.086304404
    Base.precompile(Tuple{typeof(gluniform),Int32,Observable{Vec{2, Float32}}})   # time: 0.06477806
    Base.precompile(Tuple{Core.kwftype(typeof(Type)),NamedTuple{(:x_repeat,), Tuple{Symbol}},Type{Texture},Vector{Float16}})   # time: 0.03858548
    Base.precompile(Tuple{typeof(signal_convert),Type,Observable{Vec2{Int32}}})   # time: 0.038141333
    Base.precompile(Tuple{typeof(gl_convert),Type{GLBuffer},Observable{Vec{2, Float32}}})   # time: 0.031721786
    Base.precompile(Tuple{typeof(toglsltype_string),GLBuffer{Point{3, Float32}}})   # time: 0.03062245
    Base.precompile(Tuple{typeof(signal_convert),Type,Observable{Float32}})   # time: 0.022515636
    Base.precompile(Tuple{typeof(compile_program),Vector{Shader},Vector{Tuple{Int64, String}}})   # time: 0.021999607
    let fbody = try __lookup_kwbody__(which(TextureParameters, (Type,Int64,))) catch missing end
        if !ismissing(fbody)
            precompile(fbody, (Symbol,Symbol,Symbol,Symbol,Symbol,Float32,Type{TextureParameters},Type,Int64,))
        end
    end   # time: 0.019532876
    Base.precompile(Tuple{typeof(gl_convert),Observable{Any}})   # time: 0.01664884
    Base.precompile(Tuple{typeof(gl_convert),Observable{Bool}})   # time: 0.0149467
    Base.precompile(Tuple{typeof(gl_convert),Observable{RGBA{Float32}}})   # time: 0.014625993
    Base.precompile(Tuple{typeof(const_lift),typeof(length),Observable{Vector{Point{2, Float32}}}})   # time: 0.014371832
    Base.precompile(Tuple{typeof(gl_convert),Type{GLBuffer},Observable{Vector{Point{3, Float32}}}})   # time: 0.01388594
    Base.precompile(Tuple{typeof(gl_convert),Observable{Int64}})   # time: 0.01345424
    Base.precompile(Tuple{typeof(const_lift),typeof(length),Observable{Vector{Point{3, Float32}}}})   # time: 0.011757969
    Base.precompile(Tuple{typeof(gl_convert),Type{GLBuffer},Observable{Vector{Vec{2, Float32}}}})   # time: 0.011719398
    isdefined(GLAbstraction, Symbol("#LazyShader#49#50")) && Base.precompile(Tuple{getfield(GLAbstraction, Symbol("#LazyShader#49#50")),Base.Iterators.Pairs{Symbol, Dict{String, String}, Tuple{Symbol}, NamedTuple{(:view,), Tuple{Dict{String, String}}}},Type{LazyShader},String,Vararg{String, N} where N})   # time: 0.010064982
    Base.precompile(Tuple{typeof(gl_convert),Type{GLBuffer},Observable{Vector{Vec{4, Float32}}}})   # time: 0.008807838
    isdefined(GLAbstraction, Symbol("#63#68")) && Base.precompile(Tuple{getfield(GLAbstraction, Symbol("#63#68"))})   # time: 0.008238689
    Base.precompile(Tuple{typeof(gl_convert),Type{GLBuffer},Observable{Vector{Point{2, Float32}}}})   # time: 0.008220657
    Base.precompile(Tuple{typeof(gl_convert),Type{GLBuffer},Observable{Vector{RGBA{Float32}}}})   # time: 0.007779626
    Base.precompile(Tuple{typeof(toglsltype_string),GLBuffer{Vec{4, Float32}}})   # time: 0.007245759
    Base.precompile(Tuple{typeof(isa_gl_struct),Observable{Int64}})   # time: 0.00689424
    Base.precompile(Tuple{typeof(toglsltype_string),Texture{Float16, 2}})   # time: 0.006558549
    Base.precompile(Tuple{typeof(gl_convert),Type{GLBuffer},Observable{Vector{Float32}}})   # time: 0.006534443
    Base.precompile(Tuple{typeof(gl_convert),SMatrix{4, 4, Float32, 16}})   # time: 0.0058618
    Base.precompile(Tuple{typeof(isa_gl_struct),Observable{Float32}})   # time: 0.005198346
    Base.precompile(Tuple{typeof(isa_gl_struct),Observable{Any}})   # time: 0.004741108
    Base.precompile(Tuple{typeof(toglsltype_string),Texture{Float16, 1}})   # time: 0.004597425
    Base.precompile(Tuple{typeof(toglsltype_string),GLBuffer{RGBA{Float32}}})   # time: 0.00437181
    Base.precompile(Tuple{typeof(isa_gl_struct),Observable{Vec2{Int32}}})   # time: 0.004354946
    Base.precompile(Tuple{typeof(isa_gl_struct),Observable{Int32}})   # time: 0.004301191
    Base.precompile(Tuple{typeof(isa_gl_struct),Observable{SMatrix{4, 4, Float32, 16}}})   # time: 0.004184706
    Base.precompile(Tuple{typeof(isa_gl_struct),Observable{Vec{3, Float32}}})   # time: 0.004105391
    Base.precompile(Tuple{typeof(isa_gl_struct),Observable{Vec{2, Float32}}})   # time: 0.004101697
    Base.precompile(Tuple{typeof(toglsltype_string),GLBuffer{Point{2, Float32}}})   # time: 0.003962705
    Base.precompile(Tuple{typeof(toglsltype_string),GLBuffer{Vec{2, Float32}}})   # time: 0.003796558
    Base.precompile(Tuple{typeof(isa_gl_struct),Observable{RGBA{Float32}}})   # time: 0.003778398
    Base.precompile(Tuple{typeof(isa_gl_struct),Observable{Bool}})   # time: 0.00368915
    Base.precompile(Tuple{typeof(isa_gl_struct),Dict{Symbol, Any}})   # time: 0.003686279
    Base.precompile(Tuple{typeof(isa_gl_struct),Nothing})   # time: 0.003181062
    Base.precompile(Tuple{typeof(isa_gl_struct),Symbol})   # time: 0.003076398
    Base.precompile(Tuple{typeof(isa_gl_struct),Bool})   # time: 0.003013663
    Base.precompile(Tuple{typeof(gl_convert),Vec{2, Float32}})   # time: 0.002978386
    Base.precompile(Tuple{typeof(mustache2replacement),String,Dict{String, String},Dict{Symbol, Any}})   # time: 0.002830558
    Base.precompile(Tuple{typeof(gluniform),Int32,Int64,Texture{RGBA{N0f8}, 2}})   # time: 0.002434047
    Base.precompile(Tuple{typeof(map_texture_paramers),Tuple{Symbol, Symbol}})   # time: 0.00121505
end


text:
using Plots

# Model
include_model("bimanipulation_block")

# add obstacle
function ϕ_func(model::BimanipulationBlock, q)
    p = view(q, 1:2)
	θ = q[3]
    R = rotation_matrix(θ)

	xl = 0.25
	xu = 0.75
	h = 0.05

	k1 = p + R * model.contact_corner_offset[1]
	k2 = p + R * model.contact_corner_offset[2]
	k3 = p + R * model.contact_corner_offset[3]
	k4 = p + R * model.contact_corner_offset[4]

	SVector{4}([k1[2] - ((k1[1] >= xl && k1[1] <= xu) ? h : 0.0),
			    k2[2] - ((k2[1] >= xl && k2[1] <= xu) ? h : 0.0),
			    k3[2] - ((k3[1] >= xl && k3[1] <= xu) ? h : 0.0),
			    k4[2] - ((k4[1] >= xl && k4[1] <= xu) ? h : 0.0)])
end

# Horizon
T = 11

# Time step
tf = 1.0
h = tf / (T-1)

# Bounds
_uu = Inf * ones(model.m)
_uu[model.idx_u] .= Inf
# _uu[model.idx_u[1:8]] .= 1000.0
# _uu[model.idx_u[3:4]] .= 0.0
# _uu[model.idx_u[7:8]] .= 0.0
# _uu[model.idx_λ] .= 1.0
_ul = zeros(model.m)
_ul[model.idx_u] .= -Inf
# _ul[model.idx_u[1:8]] .= -1000.0
# _ul[model.idx_u[3:4]] .= 0.0
# _ul[model.idx_u[7:8]] .= 0.0
# _ul[model.idx_λ] .= 1.0

_ul[model.idx_u] .= 0.0
_uu[model.idx_u] .= 0.0

_ul[model.idx_u[9:10]] = x0[t][nq .+ (1:2)] + model.control_input_offset[1]
_ul[model.idx_u[11:12]] = x0[t][nq .+ (1:2)] + model.control_input_offset[3]
_ul[model.idx_u[1:2]] = [-10.0; 0.5 * model.mass * model.g * h]
_ul[model.idx_u[5:6]] = [10.0; 0.5 * model.mass * model.g * h]

_uu[model.idx_u[9:10]] = x0[t][nq .+ (1:2)] + model.control_input_offset[1]
_uu[model.idx_u[11:12]] = x0[t][nq .+ (1:2)] + model.control_input_offset[3]
_uu[model.idx_u[1:2]] = [-10.0; 0.5 * model.mass * model.g * h]
_uu[model.idx_u[5:6]] = [10.0; 0.5 * model.mass * model.g * h]

ul, uu = control_bounds(model, T, _ul, _uu)

# Initial and final states
q1 = [0.0, 2.0 * r, 0.0]
x1 = [q1; q1]
qT = [0.0; 2.0 * r; 0.0]
xT = [qT; qT]
qM = [0.0, 2.0 * r, 0.0]

# Trajectory initialization
q0 = [q1, linear_interpolation(q1, qM, 13)..., linear_interpolation(qM, qT, 12)..., qT]
q0 = [q1 for t = 1:T+1]
x0 = configuration_to_state(q0)

u0 = [0.0 * rand(model.m) for t = 1:T-1] # random controls
for t = 1:T-1
	u0[t][model.idx_u] .= 0.0
	u0[t][model.idx_u[9:10]] = x0[t][nq .+ (1:2)] + model.control_input_offset[1]
	u0[t][model.idx_u[11:12]] = x0[t][nq .+ (1:2)] + model.control_input_offset[3]
	u0[t][model.idx_u[1:2]] = [-1.0; 0.5 * model.mass * model.g * h]
	u0[t][model.idx_u[5:6]] = [1.0; 0.5 * model.mass * model.g * h]
end

xl, xu = state_bounds(model, T, x1 = x1, xT = xT)

# qm1 = zeros(nq)
# vm1 = zeros(nq)
# qm2 = zeros(nq)
# vm2 = zeros(nq)
#
# D1L1, D2L1 = lagrangian_derivatives(model, qm1, vm1)
# D1L2, D2L2 = lagrangian_derivatives(model, qm2, vm2)
# 0.5 * h * D1L1 + D2L1 + 0.5 * h * D1L2 - D2L2
# B_func(model, q1)' * u0[1][model.idx_u[1:8]]
# k = control_kinematics_func(model, q1)

# Objective
include_objective(["velocity", "control_velocity", "nonlinear_stage"])

function l_stage(x, u, t)
	q2 = view(x, nq .+ (1:nq))
	k = control_kinematics_func(model, q2)
	k1 = k[1:2]
	p1 = u[model.idx_u[9:10]]

	k3 = k[5:6]
	p2 = u[model.idx_u[11:12]]

    J = 0.0
	u_ctrl_err = u[model.idx_u] - u0[t][model.idx_u]
	J += 1.0 * u_ctrl_err' * u_ctrl_err
	# J += 1000.0 * (k1 - p1)' * (k1 - p1)
	# J += 1000.0 * (k3 - p2)' * (k3 - p2)

    return J
end

l_terminal(x) = 0.0
obj_stage = nonlinear_stage_objective(l_stage, l_terminal)

# obj_velocity = velocity_objective(
#     [t > T / 2 ? Diagonal(1.0 * ones(model.nq)) : Diagonal(1.0 * ones(model.nq)) for t = 1:T-1],
#     model.nq,
#     h = h,
#     idx_angle = collect([3]))
#
obj_tracking = quadratic_tracking_objective(
    [Diagonal(1.0 * ones(model.n)) for t = 1:T],
    # [Diagonal(0.1 * ones(model.m)) for t = 1:T-1],
	[Diagonal([0.0 * ones(model.nu);
		zeros(model.nc);
		zeros(model.nb);
		zeros(model.m - model.nu - model.nc - model.nb)]) for t = 1:T-1],
    [x0[t] for t = 1:T],
    [zeros(model.m) for t = 1:T])

# obj_ctrl_vel = control_velocity_objective(Diagonal([1.0 * ones(model.nu); 0.0 * ones(model.nc + model.nb); zeros(model.m - model.nu - model.nc - model.nb)]))

obj_penalty = PenaltyObjective(1.0e5, model.m)
obj = MultiObjective([obj_penalty])#, obj_stage, obj_tracking])

# Constraints
include_constraints(["contact", "stage"])
t_idx = vcat([t for t = 1:T-1])

function control_comp_con!(c, x, u, t)
	q = x[nq .+ (1:nq)]
	s = u[model.idx_s]
    k = control_kinematics_func(model, q)
	k_input1 = u[model.idx_u][9:10]
	k_input2 = u[model.idx_u][11:12]

	u_ctrl = u[model.idx_u]

	e1 = k[1:2] - k_input1
	e2 = k[3:4] - k_input1
	e3 = k[5:6] - k_input1
	e4 = k[7:8] - k_input1

	d1 = e1' * e1
	d2 = e2' * e2
	d3 = e3' * e3
	d4 = e4' * e4

	# input 1
	c[1] = s[1] - u_ctrl[1] * d1
	c[2] = s[1] + u_ctrl[1] * d1

	c[3] = s[1] - u_ctrl[2] * d1
	c[4] = s[1] + u_ctrl[2] * d1

	c[5] = s[1] - u_ctrl[3] * d2
	c[6] = s[1] + u_ctrl[3] * d2

	c[7] = s[1] - u_ctrl[4] * d2
	c[8] = s[1] + u_ctrl[4] * d2

	c[9] = s[1] - u_ctrl[5] * d3
	c[10] = s[1] + u_ctrl[5] * d3

	c[11] = s[1] - u_ctrl[6] * d3
	c[12] = s[1] + u_ctrl[6] * d3

	c[13] = s[1] - u_ctrl[7] * d4
	c[14] = s[1] + u_ctrl[7] * d4

	c[15] = s[1] - u_ctrl[8] * d4
	c[16] = s[1] + u_ctrl[8] * d4

	# input 2
	f1 = k[1:2] - k_input2
	f2 = k[3:4] - k_input2
	f3 = k[5:6] - k_input2
	f4 = k[7:8] - k_input2

	g1 = f1' * f1
	g2 = f2' * f2
	g3 = f3' * f3
	g4 = f4' * f4

	c[17] = s[1] - u_ctrl[1] * g1
	c[18] = s[1] + u_ctrl[1] * g1

	c[19] = s[1] - u_ctrl[2] * g1
	c[20] = s[1] + u_ctrl[2] * g1

	c[21] = s[1] - u_ctrl[3] * g2
	c[22] = s[1] + u_ctrl[3] * g2

	c[23] = s[1] - u_ctrl[4] * g2
	c[24] = s[1] + u_ctrl[4] * g2

	c[25] = s[1] - u_ctrl[5] * g3
	c[26] = s[1] + u_ctrl[5] * g3

	c[27] = s[1] - u_ctrl[6] * g3
	c[28] = s[1] + u_ctrl[6] * g3

	c[29] = s[1] - u_ctrl[7] * g4
	c[30] = s[1] + u_ctrl[7] * g4

	c[31] = s[1] - u_ctrl[8] * g4
	c[32] = s[1] + u_ctrl[8] * g4

    nothing
end

n_ctrl_comp = 32
con_ctrl_comp = stage_constraints(control_comp_con!, n_ctrl_comp, (1:32), t_idx)

function control_limits_con!(c, x, u, t)
	q = x[nq .+ (1:nq)]
	θ = q[3]
	R = rotation_matrix(θ)
	u_ctrl = u[model.idx_u]

	c[1] = -1.0 * (R' * u_ctrl[1:2])[1]
	c[2] = (model.μ[end] * u_ctrl[1])^2.0 - u_ctrl[2]^2.0

	c[3] = -1.0 * (R' * u_ctrl[3:4])[2]
	c[4] = (model.μ[end] * u_ctrl[4])^2.0 - u_ctrl[3]^2.0

	c[5] = (R' * u_ctrl[5:6])[1]
	c[6] = (model.μ[end] * u_ctrl[5])^2.0 - u_ctrl[6]^2.0

	c[7] = (R' * u_ctrl[7:8])[2]
	c[8] = (model.μ[end] * u_ctrl[8])^2.0 - u_ctrl[7]^2.0
end

n_ctrl_lim = 8
con_ctrl_lim = stage_constraints(control_limits_con!, n_ctrl_lim, (1:8), t_idx)

con_contact = contact_constraints(model, T)
con = multiple_constraints([con_contact])#, con_ctrl_comp, con_ctrl_lim])

# Problem
prob = trajectory_optimization_problem(model,
               obj,
               T,
               h = h,
               xl = xl,
               xu = xu,
               ul = ul,
               uu = uu,
               con = con)

# Pack trajectories into vector
z0 = pack(x0, u0, prob)

#NOTE: may need to run examples multiple times to get good trajectories
# Solve nominal problem
@time z̄, info = solve(prob, copy(z0), tol = 1.0e-3, c_tol = 1.0e-3)

check_slack(z̄, prob)
x̄, ū = unpack(z̄, prob)
q = state_to_configuration(x̄)
u = [u[model.idx_u] for u in ū]
γ = [u[model.idx_λ] for u in ū]
b = [u[model.idx_b] for u in ū]
h̄ = h

include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
render(vis)
# open(vis)

visualize!(vis, model,
    # q0, u0,
	q, u,
	u_mag = 1.0,
    Δt = h)

# setobject!(vis["obstacle"],
# 	GeometryBasics.HyperRectangle(Vec(0.0, 0.0, 0.0), Vec(0.5, 0.5, 0.05)),
# 	MeshPhongMaterial(color = RGBA(0.7, 0.7, 0.7, 1.0)))
# settransform!(vis["obstacle"], Translation(0.25, -0.25, 0))

plot(hcat(q...)', color = :red, width = 1.0, labels = "")
plot(hcat([u..., u[end]]...)[model.idx_u[1:8], :]', linetype = :steppost, color = :black, width = 1.0, labels = "")

max_issues_repo_name: thowell/DirectMotionPlanning


text:
using Test
using GEMPIC

@testset "Maxwell 2D" begin
    function evaluate_spline_2d(nx1, nx2, degs, dofs)
        deg1, deg2 = degs
        vals = collect(reshape(dofs, nx1, nx2))
        for j in 1:nx2
            vals[:, j] .= GEMPIC.eval_uniform_periodic_spline_curve(deg1, vals[:, j])
        end
        for i in 1:nx1
            vals[i, :] .= GEMPIC.eval_uniform_periodic_spline_curve(deg2, vals[i, :])
        end
        return vec(vals)
    end

    x1min, x1max = 0.0, 2π
    nx1 = 16
    x2min, x2max = 0.0, 2π
    nx2 = 32

    mesh = TwoDGrid(x1min, x1max, nx1, x2min, x2max, nx2)

    deg = 3
    delta_t = 0.01
    nsteps = 300

    maxwell = TwoDMaxwell(mesh, deg)

    efield = [zeros(nx1 * nx2) for _ in 1:3]
    bfield = deepcopy(efield)
    efield_val = deepcopy(efield)
    bfield_val = deepcopy(efield)
    efield_ref = deepcopy(efield)
    bfield_ref = deepcopy(efield)

    x = LinRange(x1min, x1max, nx1 + 1)[1:(end - 1)] .* transpose(ones(nx2))
    y = ones(nx1) .* transpose(LinRange(x2min, x2max, nx2 + 1)[1:(end - 1)])

    w1 = sqrt(3)
    w2 = sqrt(3)

    rho = zeros(nx1 * nx2)
    rho_ref = zeros(nx1 * nx2)
    time = 0.0

    sin_k = (x, y) -> sin((x + y) - w1 * time)
    cos_k = (x, y) -> cos((x + y) - w1 * time)

    rho .= compute_rhs_from_function(maxwell, cos_k, 1, 0)

    compute_e_from_rho!(efield, maxwell, rho)

    efield_val1 = evaluate_spline_2d(nx1, nx2, (deg - 1, deg), efield[1])
    efield_val2 = evaluate_spline_2d(nx1, nx2, (deg, deg - 1), efield[2])
    efield_val3 = evaluate_spline_2d(nx1, nx2, (deg, deg), efield[3])

    efield_ref[1] .= vec(sin_k.(x, y) ./ 2)
    efield_ref[2] .= efield_ref[1]
    efield_ref[3] .= 0.0

    @test efield_val1 ≈ efield_ref[1] rtol = 1e-4
    @test efield_val2 ≈ efield_ref[2] rtol = 1e-4
    @test efield_val3 ≈ efield_ref[3]

    e1(x, y) = cos(x) * sin(y) * sin(sqrt(2) * time) / sqrt(2)
    e2(x, y) = -sin(x) * cos(y) * sin(sqrt(2) * time) / sqrt(2)
    b3(x, y) = -cos(x) * cos(y) * cos(sqrt(2) * time)

    time = -0.5 * delta_t

    bfield[1] .= l2projection(maxwell, e1, 1, 2)
    bfield[2] .= l2projection(maxwell, e2, 2, 2)
    bfield[3] .= l2projection(maxwell, b3, 3, 2)

    time = 0.0

    efield[1] .= l2projection(maxwell, e1, 1, 1)
    efield[2] .= l2projection(maxwell, e2, 2, 1)
    efield[3] .= l2projection(maxwell, b3, 3, 1)
    efield[3] .*= -1

    for istep in 1:nsteps
        compute_b_from_e!(bfield, maxwell, delta_t, efield)
        compute_e_from_b!(efield, maxwell, delta_t, bfield)
    end

    bfield_val[1] = evaluate_spline_2d(nx1, nx2, (deg, deg - 1), bfield[1])
    bfield_val[2] = evaluate_spline_2d(nx1, nx2, (deg - 1, deg), bfield[2])
    bfield_val[3] = evaluate_spline_2d(nx1, nx2, (deg - 1, deg - 1), bfield[3])

    time = (nsteps - 0.5) * delta_t
    bfield_ref[1] .= vec(e1.(x, y))
    bfield_ref[2] .= vec(e2.(x, y))
    bfield_ref[3] .= vec(b3.(x, y))

    @test bfield_ref[1] ≈ bfield_val[1] rtol = 1e-4
    @test bfield_ref[2] ≈ bfield_val[2] rtol = 1e-4
    @test bfield_ref[3] ≈ bfield_val[3] rtol = 1e-3

    efield_val[1] = evaluate_spline_2d(nx1, nx2, (deg - 1, deg), efield[1])
    efield_val[2] = evaluate_spline_2d(nx1, nx2, (deg, deg - 1), efield[2])
    efield_val[3] = evaluate_spline_2d(nx1, nx2, (deg, deg), efield[3])

    time = nsteps * delta_t

    efield_ref[1] .= vec(e1.(x, y))
    efield_ref[2] .= vec(e2.(x, y))
    efield_ref[3] .= -vec(b3.(x, y))

    @test efield_ref[1] ≈ efield_val[1] rtol = 1e-4
    @test efield_ref[2] ≈ efield_val[2] rtol = 1e-4
    @test efield_ref[3] ≈ efield_val[3] rtol = 1e-3

    efield[1] = l2projection(maxwell, cos_k, 1, 1)

    error2 = GEMPIC.inner_product(maxwell, efield[1], efield[1], 1, 1) - 2 * pi^2
    println(" Error in L2 norm squared: $error2")
    @test error2 ≈ 0 atol = 1e-5

    rho .= compute_rhs_from_function(maxwell, sin_k, 1, 1)

    compute_e_from_j!(efield[1], maxwell, rho, 1)

    efield_val1 = evaluate_spline_2d(nx1, nx2, [deg - 1, deg], efield[1])

    for i in eachindex(x, y)
        efield_ref[1][i] = cos_k(x[i], y[i]) - sin_k(x[i], y[i])
    end

    @show maximum(abs.(efield_val1 .- efield_ref[1]))
    @test efield_val1 ≈ efield_ref[1] atol = 1e-2

    time = 0.0
    rho_ref = compute_rhs_from_function(maxwell, cos_k, 1, 0)
    rho_ref .*= 2.0

    efield[1] = l2projection(maxwell, sin_k, 1, 1)
    efield[2] = l2projection(maxwell, sin_k, 2, 1)
    efield[3] = l2projection(maxwell, sin_k, 3, 1)

    compute_rho_from_e!(rho, maxwell, efield)

    @show maximum(abs.(rho .- rho_ref))

    @test rho ≈ rho_ref
end


text:
function printinrange(a, b)
    println(a)
    if a == b
        return
    else
        return printinrange(a + 1, b)
    end
end

# printinrange(2, 10)

function fibonacci(n)
    if n == 0
        return 0
    elseif n == 1
        return 1
    else
        return fibonacci(n - 1) + fibonacci(n - 2)
    end
end

# println(fibonacci(3))

const table = Dict(0=>0, 1=>1)
function fibmemo(n)
    if n in keys(table)
        return table[n]
    else
        res = table[n] = fibmemo(n - 1) + fibmemo(n - 2)
        return res
    end
end

# println(fibmemo(1000))


text:
"""
    extract_col_feats(df, cols)

find mean, std, minimum, maximum in df[!, col]
default value of columns are all numeric columns except date
"""
function extract_col_statvals(df::DataFrame, cols::Array{Symbol, 1})
    syms = []
    types = []
    vals = []
    for col in cols
        μ, σ = mean_and_std(skipmissing(df[!, col]))
        minval = minimum(skipmissing(df[!, col]))
        maxval = maximum(skipmissing(df[!, col]))

        push!(syms, String(col))
        push!(types, "μ")
        push!(vals, μ)

        push!(syms, String(col))
        push!(types, "σ")
        push!(vals, σ)

        push!(syms, String(col))
        push!(types, "minimum")
        push!(vals, minval)

        push!(syms, String(col))
        push!(types, "maximum")
        push!(vals, maxval)
    end

    ndsparse((
        symbol = syms,
        types = types),
        (value = vals,))
end

"""
    zscore!(df, col, new_col)

Apply zscore (normalization) to dataframe `df`
No need to implement `zscore`, just implement `zscore!``.
"""
function zscore!(df::DataFrame, col::Symbol, new_col::Symbol)
    to_be_normalized = df[!, col]
    df[!, new_col] = zscore(to_be_normalized)
end

function zscore!(df::DataFrame, col::Symbol, new_col::Symbol, μ::Real, σ::Real)
    to_be_normalized = df[!, col]
    df[!, new_col] = zscore(to_be_normalized, μ, σ)
end

function zscore!(df::DataFrame, cols::Array{Symbol, 1}, new_cols::Array{Symbol, 1})
    for (col, new_col) in zip(cols, new_cols)
        zscore!(df, col, new_col)
    end
end

function zscore!(df::DataFrame, cols::Array{Symbol, 1}, new_cols::Array{Symbol, 1}, μσs::Array{Real, 1})
    for (col, new_col) in zip(cols, new_cols)
        μ, σ = μσs[String(col), "μ"].value, μσs[String(col), "σ"].value
        zscore!(df, col, new_col, μ, σ)
    end
end

"""
    unzscore(A, μ, σ)

unzscore in Array using given μ and σ
"""
unzscore(A, μ::Real, σ::Real) = A .* σ .+ μ

"""
    unzscore!(df, zcol, ocol, μ, σ)

revert zscore normalization (single column, zcol -> ocol)
"""
function unzscore!(df::DataFrame, zcol::Symbol, ocol::Symbol, μ::Real, σ::Real)
    df[!, ocol] = unzscore(df[!, zcol], μ, σ)
end

"""
    unzscore!(df, zcols, ocols, μσs)

revert zscore normalization (single column, zcol -> ocol)
"""
function unzscore!(df::DataFrame, zcols::Array{Symbol, 1}, ocols::Array{Symbol, 1}, stattb::AbstractNDSparse)
    for (zcol, ocol) in zip(zcols, ocols)
        unzscore!(df, zcol, ocol,
            stattb[String(ocol), "μ"].value, stattb[String(ocol), "σ"].value)
    end
end

minmax_scaling(X::AbstractVector, a::F, b::F) where F<:AbstractFloat =
    a .+ (b - a) .* (X .- minimum(X)) ./ (maximum(X) - minimum(X))

minmax_scaling!(df::DataFrame, ocol::Symbol, mcol::Symbol,
    a::F, b::F) where F<:AbstractFloat =
        df[!, mcol] = minmax_scaling(df[!, ocol], a, b)

"""
    minmax_scaling!(df, ocols, mcols, a, b)

min-max normalization from df[!, ocols] to df[!, mcols], new range is (a, b)
https://en.wikipedia.org/wiki/Feature_scaling#Rescaling_(min-max_normalization)

Y = a + (X - min(X)) * (b - a) / (maximum(X) - minimum(X))
"""
function minmax_scaling!(df::DataFrame, ocols::Array{Symbol, 1}, mcols::Array{Symbol, 1},
    a::F, b::F) where F<:AbstractFloat

    for (ocol, mcol) in zip(ocols, mcols)
        minmax_scaling!(df, ocol, mcol, a, b)
    end
end

unminmax_scaling(Y::AbstractVector, minY::Real, maxY::Real, a::F, b::F) where F<:AbstractFloat =
    (Y .- a) .* ((maxY - minY) / (b - a)) .+ minY

"""
    unminmax_scaling!(df, mcol, ocol,  minY, maxY, a, b)

revert min-max normalization (single column)
"""
function unminmax_scaling!(df::DataFrame,
    mcol::Symbol, ocol::Symbol, minY::Real, maxY::Real, a::F, b::F) where F<:AbstractFloat
    df[!, ocol] = unminmax_scaling(df[!, mcol], minY, maxY, a, b)
end

"""
    unminmax_scaling!(df, mcols, ocols, minY, maxY, a, b)

revert min-max normalization (multiple column)
"""
function unminmax_scaling!(df::DataFrame, zcols::Array{Symbol, 1}, ocols::Array{Symbol, 1},
    stattb::AbstractNDSparse, a::F, b::F) where F<:AbstractFloat
    for (zcol, ocol) in zip(zcols, ocols)
        unminmax_scaling!(df, zcol, ocol,
        stattb[String(ocol), "minimum"].value, stattb[String(ocol), "maximum"].value, a, b)
    end
end

"""
    exclude_elem(cols, target_col)

Exclude element and return new splited array
"""
function exclude_elem(cols, target_col)
    new_cols = copy(cols)
    deleteat!(new_cols, new_cols .== target_col)

    new_cols
end

"""
    findrow(df, col, val)

Find fist row number in df[!, `col`] as `val` by brute-force
"""
function findrow(df::DataFrame, col::Symbol, val::Union{<:Real, DateTime, ZonedDateTime})

    idx = 0
    for row in eachrow(df)
        idx += 1
        if (row[col] == val)
            return idx
        end
    end

    idx = 0

    idx
end

"""
    WHO_PM10(val::Real)
    return WHO PM10 level

1: Good
2: Normal
3: Bad
4: Very Bad
"""
function WHO_PM10(val::Real)
    if 0 <= val < 31
        return 1
    elseif 31 <= val < 81
        return 2
    elseif 81 <= val < 151
        return 3
    else
        return 4
    end
end

"""
    WHO_PM25(val::Real)
    return WHO PM25 level

1: Good
2: Normal
3: Bad
4: Very Bad
"""
function WHO_PM25(val::Real)
    if 0 <= val < 16
        return 1
    elseif 16 <= val < 36
        return 2
    elseif 36 <= val < 76
        return 3
    else
        return 4
    end
end


text:
using Documenter, PCAs

makedocs(;
    modules=[PCAs],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/DarioSarra/PCAs.jl/blob/{commit}{path}#L{line}",
    sitename="PCAs.jl",
    authors="DarioSarra",
    assets=String[],
)

deploydocs(;
    repo="github.com/DarioSarra/PCAs.jl",
)


text:
using StructuralCausalModels

ProjDir = @__DIR__
cd(ProjDir) #do

d_string = "dag {A -> {E Z}; B -> {D Z}; Z -> {D E}; E -> D}"

dag = DAG("conf_triangles", d_string);
show(dag)

to_ggm(dag) |> display
println()

fname = ProjDir * "/conf_triangles.dot"
to_graphviz(dag, fname)
Sys.isapple() && run(`open -a GraphViz.app $(fname)`)

println("Basis set:")
bs = basis_set(dag)
display(bs)
#=
E ⊥ B | A, Z
D ⊥ A | B, E, Z
A ⊥ B
=#

f = :A; s = :B;

e = d_separation(dag, f, s)
println("d_separation($(dag.name), $f, $s) = $e\n")

f = [:B]; s = [:E];

e = d_separation(dag, f, s; c=[:A, :Z])
println("d_separation($(dag.name), $f, $s; cond=[:A, :Z]) = $e\n")

ap = all_paths(dag, :E, :D)
ap |> display
println()

bp = backdoor_paths(dag, ap, :E)
bp |> display
println()

adjustmentsets = adjustment_sets(dag, :E, :D)
println("Adjustment sets:")
adjustmentsets |> display

#end


text:
export left_order, right_order, normred, locally_free_basis, islocally_free

@doc Markdown.doc"""
    order(a::AlgAssRelOrdIdl) -> AlgAssRelOrd

Returns the order of which $a$ is a (possibly fractional) left or right ideal
if such an order is known.
"""
function order(a::AlgAssRelOrdIdl)
  if isdefined(a, :order)
    return a.order
  else
    error("No order defined")
  end
end

@doc Markdown.doc"""
    algebra(a::AlgAssRelOrdIdl) -> AbsAlgAss

Returns the algebra containing $a$.
"""
algebra(a::AlgAssRelOrdIdl) = a.algebra

_algebra(a::AlgAssRelOrdIdl) = algebra(a)

# The basis matrix is (should be) in lowerleft HNF, so if the upper left corner
# is not zero, then the matrix has full rank.
isfull_lattice(a::AlgAssRelOrdIdl) = !iszero(basis_matrix(a, copy = false)[1, 1])

################################################################################
#
#  Construction
#
################################################################################

@doc Markdown.doc"""
    ideal(A::AbsAlgAss, M::PMat, M_in_hnf::Bool = false) -> AlgAssRelOrdIdl

Returns the ideal in $A$ with basis pseudo-matrix $M$.
If `M_in_hnf == true`, it is assumed that $M$ is already in lower left pseudo HNF.
"""
function ideal(A::AbsAlgAss{S}, M::PMat{S, T}, M_in_hnf::Bool = false) where { S <: NumFieldElem, T }
  !M_in_hnf ? M = pseudo_hnf(M, :lowerleft) : nothing
  return AlgAssRelOrdIdl{S, T, typeof(A)}(A, M)
end

@doc Markdown.doc"""
    ideal(A::AbsAlgAss, O::AlgAssRelOrd, M::PMat, side::Symbol = :nothing,
          M_in_hnf::Bool = false)
      -> AlgAssRelOrdIdl

Returns the ideal of $O$ in $A$ with basis pseudo-matrix $M$ (in the basis of
$A$).
If the ideal is known to be a right/left/twosided ideal of $O$, `side` may be
set to `:right`/`:left`/`:twosided` respectively.
If `M_in_hnf == true`, it is assumed that $M$ is already in lower left pseudo HNF.
"""
function ideal(A::AbsAlgAss{S}, O::AlgAssRelOrd{S, T, U}, M::PMat{S, T}, side::Symbol = :nothing, M_in_hnf::Bool = false) where { S <: NumFieldElem, T, U }
  a = ideal(A, M, M_in_hnf)
  a.order = O
  _set_sidedness(a, side)
  if ismaximal_known(O) && ismaximal(O)
    if side == :left || side == :twosided
      a.left_order = O
    end
    if side == :right || side == :twosided
      a.right_order = O
    end
  end
  return a
end

@doc Markdown.doc"""
    ideal(O::AlgAssRelOrd, x::AbsAlgAssElem) -> AlgAssRelOrdIdl
    ideal(O::AlgAssRelOrd, x::AlgAssRelOrdElem) -> AlgAssRelOrdIdl

Returns the twosided principal ideal of $O$ generated by $x$.
"""
function ideal(O::AlgAssRelOrd{S, T, U}, x::AbsAlgAssElem{S}) where {S, T, U}
  A = algebra(O)
  @assert parent(x) === A
  if iszero(x)
    return _zero_ideal(A, O)
  end

  if iscommutative(O)
    a = ideal(O, x, :left)
    a.isright = 1
    return a
  end

  t1 = A()
  t2 = A()
  M = zero_matrix(base_ring(A), degree(O)^2, degree(O))
  pb = pseudo_basis(O, copy = false)
  C = Vector{T}(undef, degree(O)^2)
  for i = 1:degree(O)
    i1d = (i - 1)*degree(O)
    C[i1d + i] = pb[i][2]^2
    for j = (i + 1):degree(O)
      C[i1d + j] = simplify(pb[i][2]*pb[j][2])
      C[(j - 1)*degree(O) + i] = deepcopy(C[i1d + j])
    end
  end

  for i = 1:degree(O)
    t1 = mul!(t1, pb[i][1], x)
    ii = (i - 1)*degree(O)
    for j = 1:degree(O)
      t2 = mul!(t2, t1, pb[j][1])
      elem_to_mat_row!(M, ii + j, t2)
    end
  end
  M = sub(pseudo_hnf(PseudoMatrix(M, C), :lowerleft), nrows(M) - degree(O) + 1:nrows(M), 1:ncols(M))

  return ideal(A, O, M, :twosided, true)
end

ideal(O::AlgAssRelOrd{S, T, U}, x::AlgAssRelOrdElem{S, T, U}) where { S, T, U } = ideal(O, elem_in_algebra(x, copy = false))

@doc Markdown.doc"""
    ideal(O::AlgAssRelOrd, x::AbsAlgAssElem, side::Symbol) -> AlgAssRelOrdIdl
    ideal(O::AlgAssRelOrd, x::AlgAssRelOrdElem, side::Symbol) -> AlgAssRelOrdIdl

Returns the ideal $O \cdot x$ if `side == :left`, and $x \cdot O$ if
`side == :right`.
"""
function ideal(O::AlgAssRelOrd{S, T, U}, x::AbsAlgAssElem{S}, side::Symbol) where {S, T, U}
  A = algebra(O)
  @assert parent(x) === A
  if iszero(x)
    return _zero_ideal(A, O)
  end
  if side == :twosided
    return ideal(O, x)
  end
  @assert side == :left || side == :right "Option :$(side) for side not implemented"

  M = zero_matrix(base_ring(A), degree(O), degree(O))
  pb = pseudo_basis(O, copy = false)
  for i = 1:degree(O)
    if side == :left
      elem_to_mat_row!(M, i, pb[i][1]*x)
    elseif side == :right
      elem_to_mat_row!(M, i, x*pb[i][1])
    end
  end
  M = PseudoMatrix(M, deepcopy(basis_pmatrix(O, copy = false).coeffs))

  return ideal(A, O, M, side)
end

ideal(O::AlgAssRelOrd{S, T, U}, x::AlgAssRelOrdElem{S, T, U}, action::Symbol) where { S, T, U } = ideal(O, elem_in_algebra(x, copy = false), action)

@doc Markdown.doc"""
    *(O::AlgAssRelOrd, x::AbsAlgAssElem) -> AlgAssRelOrdIdl
    *(O::AlgAssRelOrd, x::AlgAssRelOrdElem) -> AlgAssRelOrdIdl
    *(O::AlgAssRelOrd, x::Int) -> AlgAssRelOrdIdl
    *(O::AlgAssRelOrd, x::fmpz) -> AlgAssRelOrdIdl
    *(x::AbsAlgAssElem, O::AlgAssRelOrd) -> AlgAssRelOrdIdl
    *(x::AlgAssRelOrdElem, O::AlgAssRelOrd) -> AlgAssRelOrdIdl
    *(x::Int, O::AlgAssRelOrd) -> AlgAssRelOrdIdl
    *(x::fmpz, O::AlgAssRelOrd) -> AlgAssRelOrdIdl

Returns the ideal $O \cdot x$ or $x \cdot O$ respectively.
"""
*(O::AlgAssRelOrd{S, T, U}, x::AbsAlgAssElem{S}) where {S, T, U} = ideal(O, x, :left)
*(x::AbsAlgAssElem{S}, O::AlgAssRelOrd{S, T, U}) where {S, T, U} = ideal(O, x, :right)
*(O::AlgAssRelOrd{S, T, U}, x::AlgAssRelOrdElem{S, T, U}) where {S, T, U} = ideal(O, x, :left)
*(x::AlgAssRelOrdElem{S, T, U}, O::AlgAssRelOrd{S, T, U}) where {S, T, U} = ideal(O, x, :right)
*(O::AlgAssRelOrd{S, T, U}, x::Union{ Int, fmpz }) where {S, T, U} = ideal(O, O(x), :left)
*(x::Union{ Int, fmpz }, O::AlgAssRelOrd{S, T, U}) where {S, T, U} = ideal(O, O(x), :right)

@doc Markdown.doc"""
    ideal(O::AlgAssRelOrd, a::NfOrdFracIdl) -> AlgAssRelOrdIdl
    ideal(O::AlgAssRelOrd, a::NfRelOrdFracIdl) -> AlgAssRelOrdIdl

Returns the ideal $a \cdot O$ where $a$ is a fractional ideal of `base_ring(O)`.
"""
function ideal(O::AlgAssRelOrd{S, T, U}, a::T) where {S, T, U}
  @assert order(a) === base_ring(O)
  d = degree(O)
  pb = pseudo_basis(O, copy = false)
  M = basis_matrix(O)
  PM = PseudoMatrix(M, [ a*pb[i][2] for i = 1:d ])
  return ideal(algebra(O), O, PM, :twosided)
end

@doc Markdown.doc"""
    ideal(O::AlgAssRelOrd, a::NfAbsOrdIdl) -> AlgAssRelOrdIdl
    ideal(O::AlgAssRelOrd, a::NfRelOrdIdl) -> AlgAssRelOrdIdl

Returns the ideal $a \cdot O$ where $a$ is an ideal of `base_ring(O)`.
"""
function ideal(O::AlgAssRelOrd{nf_elem, NfOrdFracIdl}, a::NfAbsOrdIdl)
  @assert order(a) === base_ring(O)
  aa = fractional_ideal(order(a), a, fmpz(1))
  return ideal(O, aa)
end

function ideal(O::AlgAssRelOrd, a::NfRelOrdIdl)
  @assert order(a) === base_ring(O)

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

@doc Markdown.doc"""
    *(O::AlgAssRelOrd, a::NfAbsOrdIdl) -> AlgAssRelOrdIdl
    *(O::AlgAssRelOrd, a::NfRelOrdIdl) -> AlgAssRelOrdIdl
    *(O::AlgAssRelOrd, a::NfOrdFracIdl) -> AlgAssRelOrdIdl
    *(O::AlgAssRelOrd, a::NfRelOrdFracIdl) -> AlgAssRelOrdIdl
    *(a::NfAbsOrdIdl, O::AlgAssRelOrd) -> AlgAssRelOrdIdl
    *(a::NfRelOrdIdl, O::AlgAssRelOrd) -> AlgAssRelOrdIdl
    *(a::NfOrdFracIdl, O::AlgAssRelOrd) -> AlgAssRelOrdIdl
    *(a::NfRelOrdFracIdl, O::AlgAssRelOrd) -> AlgAssRelOrdIdl

Returns the ideal $a \cdot O$ where $a$ is a (fractional) ideal of `base_ring(O)`.
"""
*(O::AlgAssRelOrd{S, T, U}, a::T) where {S, T, U} = ideal(O, a)

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

*(O::AlgAssRelOrd, a::Union{NfAbsOrdIdl, NfRelOrdIdl}) = ideal(O, a)

*(a::Union{NfAbsOrdIdl, NfRelOrdIdl}, O::AlgAssRelOrd) = ideal(O, a)

@doc Markdown.doc"""
    ideal_from_lattice_gens(A::AbsAlgAss, v::Vector{ <: AbsAlgAssElem })
      -> AlgAssRelOrdIdl

Returns the ideal in $A$ generated by the elements of `v` as a lattice over
`base_ring(maximal_order(base_ring(A)))`.
"""
function ideal_from_lattice_gens(A::AbsAlgAss{S}, v::Vector{ <: AbsAlgAssElem{S} }) where { S <: NumFieldElem }
  M = zero_matrix(base_ring(A), max(length(v), dim(A)), dim(A))
  for i = 1:length(v)
    elem_to_mat_row!(M, i, v[i])
  end
  PM = pseudo_hnf(PseudoMatrix(M), :lowerleft)
  if length(v) >= dim(A)
    PM = sub(PM, (nrows(PM) - dim(A) + 1):nrows(PM), 1:dim(A))
  end

  return ideal(A, PM, true)
end
@doc Markdown.doc"""
    ideal_from_lattice_gens(A::AbsAlgAss, O::AlgAssRelOrd,
                            v::Vector{ <: AbsAlgAssElem },
                            side::Symbol = :nothing)
      -> AlgAssRelOrdIdl

Returns the ideal of $O$ in $A$ generated by the elements of `v` as a lattice
over `base_ring(O)`.
If the ideal is known to be a right/left/twosided ideal of $O$, `side` may be
set to `:right`/`:left`/`:twosided` respectively.
"""
function ideal_from_lattice_gens(A::AbsAlgAss{S}, O::AlgAssRelOrd{S, T, U}, v::Vector{ <: AbsAlgAssElem{S} }, side::Symbol = :nothing) where { S <: NumFieldElem, T, U }
  a = ideal_from_lattice_gens(A, v)
  a.order = O
  _set_sidedness(a, side)
  return a
end

###############################################################################
#
#  Zero ideal
#
###############################################################################

function _zero_ideal(A::AbsAlgAss{S}) where { S <: NumFieldElem }
  a = ideal(A, PseudoMatrix(zero_matrix(base_ring(A), dim(A), dim(A))), true)
  a.iszero = 1
  return a
end

function _zero_ideal(A::AbsAlgAss{S}, O::AlgAssRelOrd{S, T, U}) where { S <: NumFieldElem, T, U }
  a = _zero_ideal(A)
  a.order = O
  _set_sidedness(a, :twosided)
  return a
end

function iszero(a::AlgAssRelOrdIdl)
  if a.iszero == 0
    if iszero_row(basis_matrix(a, copy = false), dim(algebra(a)))
      a.iszero = 1
    else
      a.iszero = 2
    end
  end
  return a.iszero == 1
end

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

function show(io::IO, a::AlgAssRelOrdIdl)
  print(io, "Ideal in ")
  show(IOContext(io, :compact => true), algebra(a))
  print(io, " with basis pseudo-matrix\n")
  show(IOContext(io, :compact => true), basis_pmatrix(a, copy = false))
end

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

function Base.deepcopy_internal(a::AlgAssRelOrdIdl, dict::IdDict)
  b = typeof(a)(algebra(a))
  for i in fieldnames(typeof(a))
    if isdefined(a, i)
      if i != :algebra && i != :order && i != :right_order && i != :left_order && i != :basis_pmatrix_wrt && i != :norm && i != normred
        setfield!(b, i, Base.deepcopy_internal(getfield(a, i), dict))
      end
    end
  end
  if isdefined(a, :order)
    b.order = order(a)
  end
  if isdefined(a, :right_order)
    b.right_order = right_order(a)
  end
  if isdefined(a, :left_order)
    b.left_order = left_order(a)
  end
  for (O, M) in a.basis_pmatrix_wrt
    b.basis_pmatrix_wrt[O] = deepcopy(M)
  end
  for (O, n) in a.norm
    b.norm[O] = deepcopy(n)
  end
  for (O, n) in a.normred
    b.normred[O] = deepcopy(n)
  end

  return b
end

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

function assure_has_pseudo_basis(a::AlgAssRelOrdIdl)
  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)
  A = algebra(a)
  K = base_ring(A)
  pb = Vector{Tuple{elem_type(A), fractional_ideal_type(order_type(K))}}()
  for i = 1:dim(A)
    t = elem_from_mat_row(A, P.matrix, i)
    push!(pb, (t, deepcopy(P.coeffs[i])))
  end
  a.pseudo_basis = pb
  return nothing
end

function assure_has_basis_pmatrix_wrt(a::AlgAssRelOrdIdl, O::AlgAssRelOrd)
  if haskey(a.basis_pmatrix_wrt, O)
    return nothing
  end

  PM = basis_pmatrix(a)
  PM.matrix = PM.matrix*basis_mat_inv(O, copy = false)
  PM = pseudo_hnf(PM, :lowerleft)
  a.basis_pmatrix_wrt[O] = PM
  return nothing
end

function assure_has_basis_matrix(a::AlgAssRelOrdIdl)
  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::AlgAssRelOrdIdl)
  if isdefined(a, :basis_mat_inv)
    return nothing
  end
  @assert isfull_lattice(a) "The ideal is not a full lattice"
  a.basis_mat_inv = inv(basis_matrix(a, copy = false))
  return nothing
end

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

@doc Markdown.doc"""
    pseudo_basis(a::AlgAssRelOrdIdl; copy::Bool = true)

Returns the pseudo basis of $a$, i. e. a vector $v$ of pairs $(e_i, a_i)$ such
that $a = \bigoplus_i a_i e_i$, where $e_i$ is an element of `algebra(a)`
and $a_i$ is a fractional ideal of `base_ring(left_order(a))`.
"""
function pseudo_basis(a::AlgAssRelOrdIdl; 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::AlgAssRelOrdIdl; copy::Bool = true) -> PMat

Returns the basis pseudo-matrix of $a$ with respect to the basis of the algebra.
"""
function basis_pmatrix(a::AlgAssRelOrdIdl; copy::Bool = true)
  if copy
    return deepcopy(a.basis_pmatrix)
  else
    return a.basis_pmatrix
  end
end

@doc Markdown.doc"""
    absolute_basis(a::AlgAssRelOrdIdl) -> Vector{ <: AbsAlgAssElem }

Returns a basis of $a$ as $\mathbb Z$-module.
"""
function absolute_basis(a::AlgAssRelOrdIdl)
  A = algebra(a)
  pb = pseudo_basis(a, copy = false)
  res = Vector{elem_type(A)}()
  for i = 1:dim(A)
    for b in basis(pb[i][2])
      push!(res, b*pb[i][1])
    end
  end
  return res
end

@doc Markdown.doc"""
    basis_pmatrix_wrt(a::AlgAssRelOrdIdl, O::AlgAssRelOrd; copy::Bool = true)
      -> PMat

Returns the basis pseudo-matrix of $a$ with respect to the basis of $O$.
"""
function basis_pmatrix_wrt(a::AlgAssRelOrdIdl, O::AlgAssRelOrd; copy::Bool = true)
  assure_has_basis_pmatrix_wrt(a, O)
  if copy
    return deepcopy(a.basis_pmatrix_wrt[O])
  else
    return a.basis_pmatrix_wrt[O]
  end
end

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

@doc Markdown.doc"""
    basis_matrix(a::AlgAssRelOrdIdl; copy::Bool = true) -> MatElem

Returns the basis matrix of $a$, that is the basis pseudo-matrix of $a$ without
the coefficient ideals.
"""
function basis_matrix(a::AlgAssRelOrdIdl; 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::AlgAssRelOrdIdl; copy::Bool = true) -> MatElem

Returns the inverse of the basis matrix of $a$.
"""
function basis_mat_inv(a::AlgAssRelOrdIdl; 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

################################################################################
#
#  Arithmetic
#
################################################################################

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

Returns $a + b$.
"""
function +(a::AlgAssRelOrdIdl{S, T, U}, b::AlgAssRelOrdIdl{S, T, U}) where {S, T, U}
  @assert algebra(a) === algebra(b)
  if iszero(a)
    return deepcopy(b)
  elseif iszero(b)
    return deepcopy(a)
  end

  d = dim(algebra(a))
  M = vcat(basis_pmatrix(a, copy = false), basis_pmatrix(b, copy = false))
  M = sub(pseudo_hnf(M, :lowerleft), (d + 1):2*d, 1:d)
  c = ideal(algebra(a), M, true)
  if isdefined(a, :order) && isdefined(b, :order) && order(a) === order(b)
    c.order = order(a)
  end
  return c
end

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

Returns $a \cdot b$.
"""
function *(a::AlgAssRelOrdIdl{S, T, U}, b::AlgAssRelOrdIdl{S, T, U}) where {S, T, U}
  @assert algebra(a) === algebra(b)

  if iszero(a)
    return deepcopy(a)
  elseif iszero(b)
    return deepcopy(b)
  end

  A = algebra(a)
  d = dim(A)
  d2 = d^2

  pba = pseudo_basis(a, copy = false)
  pbb = pseudo_basis(b, copy = false)
  M = zero_matrix(base_ring(A), d2, d)
  C = Array{fractional_ideal_type(order_type(base_ring(A))), 1}(undef, d2)
  t = one(A)
  for i = 1:d
    i1d = (i - 1)*d
    for j = 1:d
      t = mul!(t, pba[i][1], pbb[j][1])
      elem_to_mat_row!(M, i1d + j, t)
      C[i1d + j] = simplify(pba[i][2]*pbb[j][2])
    end
  end

  H = sub(pseudo_hnf(PseudoMatrix(M, C), :lowerleft), (d2 - d + 1):d2, 1:d)
  c = ideal(A, H, true)

  if _left_order_known_and_maximal(a)
    c.left_order = left_order(a)
  end
  if _right_order_known_and_maximal(b)
    c.right_order = right_order(b)
  end
  if isdefined(a, :order) && isdefined(b, :order) && order(a) === order(b)
    c.order = order(a)
  end

  return c
end

@doc Markdown.doc"""
    ^(a::AlgAssRelOrdIdl, e::Int) -> AlgAssRelOrdIdl
    ^(a::AlgAssRelOrdIdl, e::fmpz) -> AlgAssRelOrdIdl

Returns $a^e$.
"""
^(A::AlgAssRelOrdIdl, e::Int) = Base.power_by_squaring(A, e)
^(A::AlgAssRelOrdIdl, e::fmpz) = Base.power_by_squaring(A, BigInt(e))

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

Returns $a \cap b$.
"""
function intersect(a::AlgAssRelOrdIdl{S, T, U}, b::AlgAssRelOrdIdl{S, T, U}) where {S, T, U}
  @assert algebra(a) === algebra(b)
  d = dim(algebra(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)
  H = sub(pseudo_hnf(M, :lowerleft), 1:d, 1:d)
  c = ideal(algebra(a), H, true)

  if isdefined(a, :order) && isdefined(b, :order) && order(a) === order(b)
    c.order = order(a)
  end
  return c
end

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

@doc Markdown.doc"""
    *(a::AlgAssRelOrdIdl, x::NfAbsOrdElem) -> AlgAssRelOrdIdl
    *(x::NfAbsOrdElem, a::AlgAssRelOrdIdl) -> AlgAssRelOrdIdl
    *(a::AlgAssRelOrdIdl, x::NfRelOrdElem) -> AlgAssRelOrdIdl
    *(x::NfRelOrdElem, a::AlgAssRelOrdIdl) -> AlgAssRelOrdIdl
    *(a::AlgAssRelOrdIdl, x::Int) -> AlgAssRelOrdIdl
    *(x::Int, a::AlgAssRelOrdIdl) -> AlgAssRelOrdIdl
    *(a::AlgAssRelOrdIdl, x::fmpz) -> AlgAssRelOrdIdl
    *(x::fmpz, a::AlgAssRelOrdIdl) -> AlgAssRelOrdIdl
    *(a::AlgAssRelOrdIdl{S, T, U}, x::S) where { S, T, U } -> AlgAssRelOrdIdl{S, T, U}
    *(x::S, a::AlgAssRelOrdIdl{S, T, U}) where { S, T, U } -> AlgAssRelOrdIdl{S, T, U}
    *(a::AlgAssRelOrdIdl{S, T, U}, x::AbsAlgAssElem{S}) where { S, T, U }
      -> AlgAssRelOrdIdl{S, T, U}
    *(x::AbsAlgAssElem{S}, a::AlgAssRelOrdIdl{S, T, U}) where { S, T, U }
      -> AlgAssRelOrdIdl{S, T, U}
    *(a::AlgAssRelOrdIdl{S, T, U}, x::AlgAssRelOrdElem{S, T, U}) where { S, T, U }
      -> AlgAssRelOrdIdl{S, T, U}
    *(x::AlgAssRelOrdElem{S, T, U}, a::AlgAssRelOrdIdl{S, T}) where { S, T, U }
      -> AlgAssRelOrdIdl{S, T, U}

Returns the ideal $a*x$ respectively $x*a$.
"""
function *(a::AlgAssRelOrdIdl{S, T, U}, x::Union{ Int, fmpz, NfAbsOrdElem, NfRelOrdElem, S }) where { S, T, U }
  if iszero(x)
    return _zero_ideal(algebra(a))
  end
  b = ideal(algebra(a), x*basis_pmatrix(a, copy = false), true)
  if isdefined(a, :left_order)
    b.left_order = left_order(a)
  end
  if isdefined(a, :right_order)
    b.right_order = right_order(a)
  end
  if isdefined(a, :order)
    b.order = order(a)
  end
  return b
end

# Let's assume the base ring of algebra(a) is commutative (it should be a field)
*(x::Union{ Int, fmpz, NfAbsOrdElem, NfRelOrdElem, S }, a::AlgAssRelOrdIdl{S, T, U}) where { S, T, U } = a*x

function *(a::AlgAssRelOrdIdl{S, T, U}, x::AbsAlgAssElem{S}) where { S, T, U }
  if iszero(x)
    return _zero_ideal(algebra(a))
  end
  A = algebra(a)
  M = zero_matrix(base_ring(A), dim(A), dim(A))
  for i = 1:dim(A)
    t = pseudo_basis(a, copy = false)[i][1]*x
    elem_to_mat_row!(M, i, t)
  end
  PM = PseudoMatrix(M, [ deepcopy(pseudo_basis(a, copy = false)[i][2]) for i = 1:dim(A) ])
  b = ideal(A, PM)
  if isdefined(a, :left_order)
    b.left_order = left_order(a)
  end
  return b
end

function *(x::AbsAlgAssElem{S}, a::AlgAssRelOrdIdl{S, T, U}) where { S, T, U }
  if iszero(x)
    return _zero_ideal(algebra(a))
  end
  A = algebra(a)
  M = zero_matrix(base_ring(A), dim(A), dim(A))
  for i = 1:dim(A)
    t = x*pseudo_basis(a, copy = false)[i][1]
    elem_to_mat_row!(M, i, t)
  end
  PM = PseudoMatrix(M, [ deepcopy(pseudo_basis(a, copy = false)[i][2]) for i = 1:dim(A) ])
  b = ideal(A, PM)
  if isdefined(a, :right_order)
    b.right_order = right_order(a)
  end
  return b
end

function *(a::AlgAssRelOrdIdl{S, T, U}, x::AlgAssRelOrdElem{S, T, U}) where { S, T, U }
  b = a*elem_in_algebra(x, copy = false)
  if isdefined(a, :order) && parent(x) === order(a) && isright_ideal(a)
    b.order = order(a)
    b.isright = a.isright
    b.isleft = a.isleft
  end
  return b
end

function *(x::AlgAssRelOrdElem{S, T, U}, a::AlgAssRelOrdIdl{S, T, U}) where { S, T, U }
  b = elem_in_algebra(x, copy = false)*a
  if isdefined(a, :order) && parent(x) === order(a) && isleft_ideal(a)
    b.order = order(a)
    b.isright = a.isright
    b.isleft = a.isleft
  end
  return b
end

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

@doc Markdown.doc"""
    in(x::AbsAlgAssElem, a::AlgAssRelOrdIdl) -> Bool
    in(x::AlgAssRelOrdElem, a::AlgAssRelOrdIdl) -> Bool

Returns `true` if $x$ is in $a$ and `false` otherwise.
"""
function in(x::AbsAlgAssElem{S}, a::AlgAssRelOrdIdl{S, T, U}) where { S, T, U }
  parent(x) !== algebra(a) && error("Algebra of element and ideal must be equal")
  A = algebra(a)
  b_pmat = basis_pmatrix(a, copy = false)
  t = matrix(base_ring(A), 1, dim(A), coeffs(x, copy = false))
  t = t*basis_mat_inv(a, copy = false)
  for i = 1:dim(A)
    if !(t[1, i] in b_pmat.coeffs[i])
      return false
    end
  end
  return true
end

in(x::AlgAssRelOrdElem, a::AlgAssRelOrdIdl) = in(elem_in_algebra(x, copy = false), a)

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

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

Returns `true` if $a$ and $b$ are equal and `false` otherwise.
"""
function ==(a::AlgAssRelOrdIdl{S, T, U}, b::AlgAssRelOrdIdl{S, T, U}) where { S, T, U }
  algebra(a) !== algebra(b) && return false
  return basis_pmatrix(a, copy = false) == basis_pmatrix(b, copy = false)
end

################################################################################
#
#  isleft/isright
#
################################################################################

# functions isright_ideal and isleft_ideal are in AlgAss/Ideal.jl

# This only works if a.order is defined, otherwise order(a) throws an error.

function _test_ideal_sidedness(a::AlgAssRelOrdIdl, O::AlgAssRelOrd, side::Symbol)
  b = ideal(O, one(O))

  if side == :right
    c = a*b
  elseif side == :left
    c = b*a
  else
    error("side must be either :left or :right")
  end

  return _spans_subset_of_pseudohnf(basis_pmatrix(c, copy = false), basis_pmatrix(a, copy = false), :lowerleft)
end

function _test_ideal_sidedness(a::AlgAssRelOrdIdl, side::Symbol)
  O = order(a)
  b = ideal(O, one(O))

  if side == :right
    c = a*b
  elseif side == :left
    c = b*a
  else
    error("side must be either :left or :right")
  end

  return _spans_subset_of_pseudohnf(basis_pmatrix(c, copy = false), basis_pmatrix(a, copy = false), :lowerleft)
end

################################################################################
#
#  Ring of multipliers, left and right order
#
################################################################################

# This computes a basis pseudo-matrix for \{ x \in A | bx \subseteq a \} if
# side == :left or \{ x \in A | xb \subseteq a \} if side == :right.
function _colon_raw(a::AlgAssRelOrdIdl{S, T, U}, b::AlgAssRelOrdIdl{S, T, U}, side::Symbol) where { S, T, U }
  @assert isfull_lattice(a) && isfull_lattice(b)
  A = algebra(a)
  @assert A === algebra(b)
  K = base_ring(A)
  d = dim(A)
  pba = pseudo_basis(a, copy = false)
  pbb = pseudo_basis(b, copy = false)
  M = zero_matrix(K, d^2, d)
  for i = 1:d
    N = representation_matrix(pbb[i][1], side)*basis_mat_inv(a, copy = false)
    for s = 1:d
      for t = 1:d
        M[t + (i - 1)*d, s] = N[s, t]
      end
    end
  end
  invcoeffs = [ simplify(inv(pba[i][2])) for i = 1:d ]
  C = Array{T}(undef, d^2)
  for i = 1:d
    for j = 1:d
      C[(i - 1)*d + j] = simplify(pbb[i][2]*invcoeffs[j])
    end
  end
  PM = PseudoMatrix(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 ])
  return PN
end

function ring_of_multipliers(a::AlgAssRelOrdIdl, action::Symbol = :left)
  PM = _colon_raw(a, a, action)
  return Order(algebra(a), PM)
end

@doc Markdown.doc"""
    left_order(a::AlgAssRelOrdIdl) -> AlgAssRelOrd

Returns the largest order of which $a$ is a left ideal, that is
$\{ x \in A \mid xa \subseteq a\}$.
"""
function left_order(a::AlgAssRelOrdIdl)
  if isdefined(a, :left_order)
    return a.left_order
  end

  a.left_order = ring_of_multipliers(a, :right)
  return a.left_order
end

@doc Markdown.doc"""
    right_order(a::AlgAssRelOrdIdl) -> AlgAssRelOrd

Returns the largest order of which $a$ is a right ideal, that is
$\{ x \in A \mid ax \subseteq a\}$.
"""
function right_order(a::AlgAssRelOrdIdl)
  if isdefined(a, :right_order)
    return a.right_order
  end

  a.right_order = ring_of_multipliers(a, :left)
  return a.right_order
end

# Checks whether the left order is known and if it is known to be maximal
function _left_order_known_and_maximal(a::AlgAssRelOrdIdl)
  if !isdefined(a, :left_order)
    return false
  end
  O = left_order(a)
  if ismaximal_known(O) && ismaximal(O)
    return true
  end
  return false
end

# Checks whether the right order is known and if it is known to be maximal
function _right_order_known_and_maximal(a::AlgAssRelOrdIdl)
  if !isdefined(a, :right_order)
    return false
  end
  O = right_order(a)
  if ismaximal_known(O) && ismaximal(O)
    return true
  end
  return false
end

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

function mod!(x::AlgAssRelOrdElem, a::AlgAssRelOrdIdl)
  A = algebra(a)
  O = parent(x)
  @assert order(a) === O
  b = coordinates(x, copy = false)
  PM = basis_pmatrix_wrt(a, O, 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 = A()
  B = pseudo_basis(O, copy = false)
  zero!(x.elem_in_algebra)
  for i = 1:degree(O)
    t = mul!(t, B[i][1], A(b[i]))
    x.elem_in_algebra = add!(x.elem_in_algebra, x.elem_in_algebra, t)
  end

  return x
end

@doc Markdown.doc"""
    mod(x::AlgAssRelOrdElem, a::AlgAssRelOrdIdl) -> AlgAssRelOrdElem

Returns $y$ in `parent(x)` such that $x \equiv y \mod a$ and the coefficients
of $y$ are reduced modulo $a$.
Assumes `parent(x) == order(a)` and that $a$ is an integral ideal of `order(a)`.
"""
function mod(x::AlgAssRelOrdElem, a::AlgAssRelOrdIdl)
  return mod!(deepcopy(x), a)
end

function mod!(x::AlgAssRelOrdElem, Q::RelOrdQuoRing)
  return mod!(x, ideal(Q))
end

function mod(x::AlgAssRelOrdElem, Q::RelOrdQuoRing)
  return mod(x, ideal(Q))
end

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

# Assumes, that det(basis_matrix(a)) == 1
function assure_has_norm(a::AlgAssRelOrdIdl{S, T, U}, O::AlgAssRelOrd{S, T, U}) where { S, T, U }
  if haskey(a.norm, O)
    return nothing
  end
  if iszero(a)
    a.norm[O] = _algebra(O)()*O
    return nothing
  end

  # Now consider a as an ideal of O (but without computing a basis pseudo matrix
  # in the basis of O).
  c = basis_pmatrix(a, copy = false).coeffs
  d = inv_coeff_ideals(O, copy = false)
  n = c[1]*d[1]
  for i = 2:dim(algebra(a))
    n *= c[i]*d[i]
  end
  # The basis matrix of a in the basis of O would be
  # basis_matrix(a)*basis_mat_inv(O), so its determinant is
  # inv(det(basis_matrix(O))), since det(basis_matrix(a)) == 1.
  n *= inv(det(basis_matrix(O, copy = false)))
  simplify(n)
  a.norm[O] = n
  return nothing
end

@doc Markdown.doc"""
    norm(a::AlgAssRelOrdIdl{S, T, U}, O::AlgAssRelOrd{S, T, U}; copy::Bool = true)
      where { S, T, U } -> T

Returns the norm of $a$ considered as an (possibly fractional) ideal of $O$.
"""
function norm(a::AlgAssRelOrdIdl{S, T, U}, O::AlgAssRelOrd{S, T, U}; copy::Bool = true) where { S, T, U }
  assure_has_norm(a, O)
  if copy
    return deepcopy(a.norm[O])
  else
    return a.norm[O]
  end
end

@doc Markdown.doc"""
    norm(a::AlgAssRelOrdIdl{S, T, U}; copy::Bool = true) where { S, T, U } -> T

Returns the norm of $a$ considered as an (possibly fractional) ideal of
`order(a)`.
"""
norm(a::AlgAssRelOrdIdl; copy::Bool = true) = norm(a, order(a), copy = copy)

function assure_has_normred(a::AlgAssRelOrdIdl{S, T, U}, O::AlgAssRelOrd{S, T, U}) where { S, T, U }
  if haskey(a.normred, O)
    return nothing
  end
  if iszero(a)
    a.normred[O] = norm(a, O)
    return nothing
  end

  A = algebra(a)
  m = isqrt(dim(A))
  @assert m^2 == dim(A)
  N = norm(a, O, copy = false)
  b, I = ispower(N, m)
  @assert b "Cannot compute reduced norm. Maybe the algebra is not simple and central?"
  a.normred[O] = I
  return nothing
end

@doc Markdown.doc"""
    normred(a::AlgAssRelOrdIdl{S, T}, O::AlgAssRelOrd{S, T}; copy::Bool = true)
      where { S, T, U } -> T

Returns the reduced norm of $a$ considered as an (possibly fractional) ideal
of $O$.
It is assumed that the algebra containing $a$ is simple and central.
"""
function normred(a::AlgAssRelOrdIdl{S, T, U}, O::AlgAssRelOrd{S, T, U}; copy::Bool = true) where { S, T, U }
  @assert issimple(algebra(a)) && iscentral(algebra(a)) "Only implemented for simple and central algebras"
  assure_has_normred(a, O)
  if copy
    return deepcopy(a.normred[O])
  else
    return a.normred[O]
  end
end

@doc Markdown.doc"""
    normred(a::AlgAssRelOrdIdl{S, T, U}; copy::Bool = true) where { S, T, U } -> T

Returns the reduced norm of $a$ considered as an (possibly fractional) ideal
of `order(a)`.
"""
normred(a::AlgAssRelOrdIdl; copy::Bool = true) = normred(a, order(a), copy = copy)

################################################################################
#
#  Locally free basis
#
################################################################################

@doc Markdown.doc"""
    locally_free_basis(a::AlgAssRelOrdIdl, p::Union{ NfAbsOrdIdl, NfRelOrdIdl })
      -> AlgAssRelOrdElem

Returns an element $x$ of the order $O$ of $a$ such that $a_p = O_p \cdot x$
where $p$ is a prime ideal of `base_ring(O)`.
See also `islocally_free`.
"""
locally_free_basis(I::AlgAssRelOrdIdl, p::Union{ NfAbsOrdIdl, NfRelOrdIdl }) = locally_free_basis(order(I), I, p)

@doc Markdown.doc"""
    locally_free_basis(O::AlgAssRelOrd, a::AlgAssRelOrdIdl,
                       p::Union{ NfAbsOrdIdl, NfRelOrdIdl }) -> AlgAssRelOrdElem

Returns an element $x$ of $O$ such that $a_p = O_p \cdot x$ where $p$ is a
prime ideal of `base_ring(O)`.
It is assumed that $a$ is an ideal of $O$ and $a \subseteq O$.
See also `islocally_free`.
"""
function locally_free_basis(O::AlgAssRelOrd, I::AlgAssRelOrdIdl, p::Union{ NfAbsOrdIdl, NfRelOrdIdl })
  b, x = islocally_free(O, I, p)
  if !b
    error("The ideal is not locally free at the prime")
  end
  return x
end

@doc Markdown.doc"""
    islocally_free(a::AlgAssRelOrdIdl, p::Union{ NfAbsOrdIdl, NfRelOrdIdl })
      -> Bool, AlgAssRelOrdElem

Returns a tuple `(true, x)` with an element $x$ of the order $O$ of $a$ such
that $a_p = O_p x$ if $a$ is locally free at $p$, and `(false, 0)` otherwise.
$p$ is a prime ideal of `base_ring(O)`.
See also `locally_free_basis`.
"""
islocally_free(I::AlgAssRelOrdIdl, p::Union{ NfAbsOrdIdl, NfRelOrdIdl }) = islocally_free(order(I), I, p)

# See Bley, Wilson "Computations in relative algebraic K-groups", section 4.2
@doc Markdown.doc"""
    islocally_free(O::AlgAssRelOrd, a::AlgAssRelOrdIdl,
                   p::Union{ NfAbsOrdIdl, NfRelOrdIdl }) -> Bool, AlgAssRelOrdElem

Returns a tuple `(true, x)` with an element $x$ of $O$ such that $a_p = O_p x$
if $a$ is locally free at $p$, and `(false, 0)` otherwise.
$p$ is a prime ideal of `base_ring(O)`.
It is assumed that $a$ is an ideal of $O$ and $a \subseteq O$.
See also `locally_free_basis`.
"""
function islocally_free(O::AlgAssRelOrd, I::AlgAssRelOrdIdl, p::Union{ NfAbsOrdIdl, NfRelOrdIdl })
  pO = p*O
  OpO, toOpO = AlgAss(O, pO, p)
  J = radical(OpO)
  OJ, toOJ = quo(OpO, J)
  decOJ = decompose(OJ)

  pI = pO*I
  IpI, toIpI = AlgAss(I, pI, p)
  basisIpI = elem_type(O)[ toIpI\IpI[i] for i = 1:dim(IpI) ]
  gensJ = Vector{elem_type(IpI)}()
  for b in basis(J, copy = false)
    bb = toOpO\b
    for c in basisIpI
      push!(gensJ, toIpI(bb*c))
    end
  end
  JinIpI = ideal_from_gens(IpI, gensJ)
  IJ, toIJ = quo(IpI, JinIpI)

  a = O()
  for i = 1:length(decOJ)
    A, AtoOJ = decOJ[i]
    B, BtoA = _as_algebra_over_center(A)
    _T = AlgMat{elem_type(base_ring(B)), dense_matrix_type(base_ring(B))}
    C, BtoC = _as_matrix_algebra(B)::Tuple{_T, AbsAlgAssMor{AlgAss{elem_type(base_ring(B))}, _T, dense_matrix_type(base_ring(B))}}
    e = toOpO\(toOJ\(AtoOJ(BtoA(BtoC\C[1]))))
    basiseIJ = Vector{elem_type(IJ)}()
    for b in basisIpI
      bb = toIJ(toIpI(e*b))
      if !iszero(bb)
        push!(basiseIJ, bb)
      end
    end

    # Construct an Fq-basis for e*IJ where Fq \cong centre(A)
    Z, ZtoA = center(A)
    basisZ = elem_type(O)[ toOpO\(toOJ\(AtoOJ(ZtoA(Z[i])))) for i = 1:dim(Z) ]

    basiseIJoverZ = Vector{elem_type(O)}()
    M = zero_matrix(base_ring(IJ), dim(Z), dim(IJ))
    MM = zero_matrix(base_ring(IJ), 0, dim(IJ))
    r = 0
    for i = 1:length(basiseIJ)
      b = toIpI\(toIJ\basiseIJ[i])

      for j = 1:dim(Z)
        bb = toIJ(toIpI(basisZ[j]*b))
        elem_to_mat_row!(M, j, bb)
      end

      N = vcat(MM, M)
      s = rank(N)
      if s > r
        push!(basiseIJoverZ, b)
        MM = N
        r = s
      end
      if r == length(basiseIJ)
        break
      end
    end

    if length(basiseIJoverZ) != degree(C)
      # I is not locally free
      return false, O()
    end

    for i = 1:length(basiseIJoverZ)
      a += mod(toOpO\(toOJ\(AtoOJ(BtoA(BtoC\C[i]))))*basiseIJoverZ[i], pI)
    end
  end

  return true, mod(a, pI)
end

################################################################################
#
#  p-Radical
#
################################################################################

# See Friedrichs: "Berechnung von Maximalordnungen über Dedekindringen", Algorithmus 5.1
@doc Markdown.doc"""
    pradical(O::AlgAssRelOrd, p::Union{ NfAbsOrdIdl, NfRelOrdIdl })
      -> AlgAssRelOrdIdl

Returns the ideal $\sqrt{p \cdot O}$ where $p$ is a prime ideal of `base_ring(O)`.
"""
function pradical(O::AlgAssRelOrd, p::Union{ NfAbsOrdIdl, NfRelOrdIdl })
  K = base_ring(algebra(O))
  OK = maximal_order(K)
  pO = p*O
  OpO, OtoOpO = AlgAss(O, pO, p)
  J = radical(OpO)

  if isempty(basis(J, copy = false))
    return pO
  end

  N = basis_pmatrix(pO, copy = false)
  m = numerator(det(N), copy = false)
  t = PseudoMatrix(zero_matrix(K, 1, degree(O)))
  for b in basis(J, copy = false)
    bb = OtoOpO\b
    elem_to_mat_row!(t.matrix, 1, elem_in_algebra(bb, copy = false))
    N = vcat(N, deepcopy(t))
  end
  N = sub(pseudo_hnf_full_rank_with_modulus(N, m, :lowerleft), nrows(N) - degree(O) + 1:nrows(N), 1:degree(O))
  return ideal(algebra(O), O, N, :twosided, true)
end

################################################################################
#
#  Primes lying over a prime
#
################################################################################

# See Friedrichs: "Berechnung von Maximalordnungen über Dedekindringen", Algorithmus 5.23
@doc Markdown.doc"""
    prime_ideals_over(O::AlgAssRelOrd, p::Union{ NfAbsOrdIdl, NfRelOrdIdl })
      -> Vector{AlgAssRelOrdIdl}

Returns all prime ideals of $O$ lying over $p$ where $p$ is a prime ideal
of `base_ring(O)`.
"""
function prime_ideals_over(O::AlgAssRelOrd, p::Union{ NfAbsOrdIdl, NfRelOrdIdl })
  return _prime_ideals_over(O, pradical(O, p), p)
end

# prad is expected to be pradical(O, p). If strict_containment == true and prad
# is already prime, we return an empty array.
function _prime_ideals_over(O::AlgAssRelOrd, prad::AlgAssRelOrdIdl, p::Union{ NfAbsOrdIdl, NfRelOrdIdl }; strict_containment::Bool = false)
  K = base_ring(algebra(O))
  OK = order(p)

  A, OtoA = AlgAss(O, prad, p)
  decA = decompose(A)

  if length(decA) == 1
    if strict_containment
      return ideal_type(O)[]
    else
      return [ prad ]
    end
  end

  lifted_components = Vector{typeof(basis_pmatrix(prad, copy = false))}()
  for k = 1:length(decA)
    N = zero_matrix(K, dim(decA[k][1]), degree(O))
    for i = 1:dim(decA[k][1])
      b = OtoA\(decA[k][2](decA[k][1][i]))
      elem_to_mat_row!(N, i, elem_in_algebra(b, copy = false))
    end
    push!(lifted_components, PseudoMatrix(N))
  end

  primes = Vector{ideal_type(O)}()
  for i = 1:length(decA)
    N = basis_pmatrix(prad, copy = false)
    m = numerator(det(N), copy = false)
    for j = 1:length(decA)
      if i == j
        continue
      end

      N = vcat(N, lifted_components[j])
    end
    N = sub(pseudo_hnf_full_rank_with_modulus(N, m, :lowerleft), nrows(N) - degree(O) + 1:nrows(N), 1:degree(O))
    push!(primes, ideal(algebra(O), O, N, :twosided, true))
  end

  return primes
end

################################################################################
#
#  Maximal integral ideals
#
################################################################################

# Computes any maximal integral ideal with left order O (if side = :left) or
# right order O (if side = :right) which contains p.
# Assumes (so far?) that the algebra is simple and O is maximal.
function maximal_integral_ideal(O::AlgAssRelOrd, p::Union{ NfAbsOrdIdl, NfRelOrdIdl }, side::Symbol)
  A = algebra(O)
  @assert issimple(A)
  @assert ismaximal(O)

  K = base_ring(algebra(O))
  OK = base_ring(O)

  P = pradical(O, p) # if the algebra is simple, then the pradical is the unique prime lying over p

  # P is the Jacobson radical of O/pO, so O/P is a simple algebra
  B, OtoB = quo(O, P, p)
  C, BtoC, CtoB = _as_algebra_over_center(B)
  D, CtoD = _as_matrix_algebra(C)

  n = degree(D)
  if isone(n)
    return P
  end

  N = basis_pmatrix(P)
  m = numerator(det(N), copy = false)
  t = zero_matrix(K, 1, degree(O))
  # Now we only need to lift a basis for diag(1, ..., 1, 0)*D (side = :left) or
  # D*diag(1, ..., 1, 0) (side = :right) since these are maximal ideals of D.
  if side == :left
    jMax = n - 1
    iMax = n
  elseif side == :right
    jMax = n
    iMax = n - 1
  else
    error("Option :$(side) for side not implemented")
  end
  for j = 1:jMax
    jn = (j - 1)*n
    for i = 1:iMax
      b = (OtoB\(CtoB(CtoD\D[jn + i])))
      elem_to_mat_row!(t, 1, elem_in_algebra(b, copy = false))
      N = vcat(N, PseudoMatrix(deepcopy(t), [ K(1)*OK ]))
    end
  end
  N = sub(pseudo_hnf_full_rank_with_modulus(N, m, :lowerleft), nrows(N) - degree(O) + 1:nrows(N), 1:degree(O))

  M = ideal(algebra(O), O, N, side, true)
  if side == :left
    M.left_order = O # O is maximal
  else
    M.right_order = O
  end
  return M
end

# Constructs a maximal integral ideal M of O such that M\cap R = p and I\subseteq M,
# where O is the left order (if side = :left) or right order (if side = :right)
# of I. It is assumed that I \subseteq O.
# M is a left ideal of O if side = :left and a right ideal if side = :right.
# Assumes (so far?) that the algebra is simple and O is maximal.
function maximal_integral_ideal_containing(I::AlgAssRelOrdIdl, p::Union{ NfAbsOrdIdl, NfRelOrdIdl }, side::Symbol)
  if side == :left
    O = left_order(I)
  elseif side == :right
    O = right_order(I)
  else
    error("Option :$(side) for side not implemented")
  end

  @assert issimple(algebra(O))
  @assert ismaximal(O)

  n = normred(I, O)
  if valuation(n, p) == 0
    error("Cannot find a maximal ideal for the given prime")
  end
  if n == p
    # The ideal is maximal iff its reduced norm is prime
    return I
  end

  P = pradical(O, p) # if the algebra is simple, then the pradical is the unique prime lying over p
  J = I + P
  if normred(J, O) == p
    return J
  end

  OP, toOP = quo(O, P, p)
  B, OPtoB, BtoOP = _as_algebra_over_center(OP)
  C, toC = _as_matrix_algebra(B)

  JinC = ideal_from_gens(C, [ toC(OPtoB(toOP(O(b)))) for b in absolute_basis(J) ])
  y = left_principal_generator(JinC)
  m = matrix(y)
  r = rref!(m)
  k = degree(C)
  @assert r < k - 1 # Otherwise J would be maximal, which we have checked...

  # We need to "add" pivots to m to get rank k - 1
  row = 1
  col = 1
  while r < k - 1
    if !iszero(m[row, col])
      row += 1
      continue
    end

    r += 1
    m[r, col] = 1
    col += 1
  end
  @assert rank(m) == k - 1
  if side == :left
    c = C*C(m)
  else
    c = C(m)*C
  end
  basis_c = basis(c, copy = false)

  t = zero_matrix(base_ring(algebra(O)), length(basis_c), degree(O))
  for i = 1:length(basis_c)
    b = toOP\(BtoOP(toC\(basis_c[i])))
    elem_to_mat_row!(t, i, elem_in_algebra(b, copy = false))
  end
  PM = vcat(basis_pmatrix(P), PseudoMatrix(t))
  n = numerator(det(basis_pmatrix(P, copy = false)), copy = false)
  PM = sub(pseudo_hnf_full_rank_with_modulus(PM, n, :lowerleft), length(basis_c) + 1:nrows(PM), 1:ncols(PM))

  M = ideal(algebra(O), O, PM, side, true)
  @assert normred(M, O) == p
  if side == :left
    M.left_order = O # O is maximal
  else
    M.right_order = O
  end

  return M
end

################################################################################
#
#  Factorization into maximal ideals
#
################################################################################

# Computes maximal integral ideals M_1, ..., M_k such that I = M_1\cdots M_k,
# and such that the right order of any factor is the left order of the next
# factor.
# Always considers I as an ideal of its left order.
function factor(I::AlgAssRelOrdIdl)
  O = left_order(I)
  @assert ismaximal(O)

  factors = Vector{ideal_type(O)}()
  n = normred(J, O)
  fac_n = factor(n)
  primes = collect(keys(fac_n))
  sort!(primes, lt = (p, q) -> minimum(p, copy = false) < minimum(q, copy = false))
  fac_n[primes[end]] -= 1 # We don't need to find the "last" maximal ideal
  for p in primes
    for i = 1:fac_n[p]
      M = maximal_integral_ideal_containing(I, p, :left)
      push!(factors, M)
      I = divexact_left(I, M)
    end
  end
  push!(factors, J)

  return factors
end

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

@doc Markdown.doc"""
    rand(a::AlgAssRelOrdIdl, B::Int) -> AbsAlgAssElem

Returns a random element of $a$ whose coefficient size is controlled by $B$.
"""
function rand(a::AlgAssRelOrdIdl, B::Int)
  pb = pseudo_basis(a, copy = false)
  z = algebra(a)()
  for i = 1:degree(order(a))
    t = rand(pb[i][2], B)
    z += t*pb[i][1]
  end
  return z
end

################################################################################
#
#  Coprime representative
#
################################################################################

# Returns x \in A with Ix + a*O == O and Ix \subseteq O.
# a should be an ideal of base_ring(O).
function integral_coprime_representative(O::AlgAssRelOrd, I::AlgAssRelOrdIdl, a::Union{ NfAbsOrdIdl, NfRelOrdIdl })
  A = algebra(O)
  d = denominator(I, O)
  I = d*I

  if one(O) in I + a*O
    return one(A)*d
  end

  fac_a = factor(a)
  primes = collect(keys(fac_a))

  x = A()
  vals = ones(Int, length(primes))
  for i = 1:length(primes)
    p = primes[i]
    vals[i] = 0
    z = approximate_nonnegative(vals, primes)
    vals[i] = 1
    if one(O) in I + p*O
      # I is already coprime to this prime
      x += A(elem_in_nf(z, copy = false))
      continue
    end
    b, g = islocally_free(O, I, p)
    @assert b "No local generator found for $p"
    ig = inv(elem_in_algebra(g, copy = false))
    Ig = I*ig
    y = coprime_denominator(Ig, O, p)
    x += ig*elem_in_nf(y, copy = false)*elem_in_nf(z, copy = false)
  end
  return x*d
end

################################################################################
#
#  Inverses
#
################################################################################

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

Returns an ideal $b$ in `algebra(a)` such that `a*b == left_order(a)`
and `b*a == right_order(b)`.
"""
function inv(a::AlgAssRelOrdIdl)
  if isdefined(a, :right_order)
    O = right_order(a)
    PM = _colon_raw(O(1)*O, a, :right)
  else
    O = left_order(a)
    PM = _colon_raw(O(1)*O, a, :left)
  end
  b = ideal(algebra(a), PM)
  if isdefined(a, :left_order)
    b.right_order = left_order(a)
  end
  if isdefined(a, :right_order)
    b.left_order = right_order(a)
  end
  if isdefined(a, :order)
    b.order = order(a)
  end
  return b
end

################################################################################
#
#  Divexact
#
################################################################################

@doc Markdown.doc"""
    divexact_left(a::AlgAssRelOrdIdl, b::AlgAssRelOrdIdl) -> AlgAssRelOrdIdl

Returns an ideal $c$ such that $a = b \cdot c$.
"""
function divexact_left(a::AlgAssRelOrdIdl{S, T, U}, b::AlgAssRelOrdIdl{S, T, U}) where { S, T, U }
  @assert algebra(a) === algebra(b)
  PM = _colon_raw(a, b, :left)
  c = ideal(algebra(a), PM)

  if _right_order_known_and_maximal(a)
    c.right_order = right_order(a)
  end
  if _right_order_known_and_maximal(b)
    c.left_order = right_order(b)
  end

  return c
end

@doc Markdown.doc"""
    divexact_right(a::AlgAssRelOrdIdl, b::AlgAssRelOrdIdl) -> AlgAssRelOrdIdl

Returns an ideal $c$ such that $a = c \cdot b$.
"""
function divexact_right(a::AlgAssRelOrdIdl{S, T, U}, b::AlgAssRelOrdIdl{S, T, U}) where { S, T, U }
  @assert algebra(a) === algebra(b)
  PM = _colon_raw(a, b, :right)
  c = ideal(algebra(a), PM)

  if _left_order_known_and_maximal(a)
    c.left_order = left_order(a)
  end
  if _left_order_known_and_maximal(b)
    c.right_order = left_order(b)
  end

  return c
end

colon(a::AlgAssRelOrdIdl{S, T, U}, b::AlgAssRelOrdIdl{S, T, U}) where { S, T, U } = divexact_right(a, b)

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

@doc Markdown.doc"""
    denominator(a::AlgAssRelOrdIdl, O::AlgAssRelOrd) -> fmpz

Returns the smallest positive integer $d$ such that $da$ is contained in $O$.
"""
function denominator(a::AlgAssRelOrdIdl{S, T, U}, O::AlgAssRelOrd{S, T, U}) where { S, T, U }
  @assert algebra(a) === algebra(O)
  if iszero(a)
    return fmpz(1)
  end

  n = dim(algebra(a))
  PM = basis_pmatrix_wrt(a, O, copy = false)
  pb = pseudo_basis(O, copy = false)
  inv_coeffs = inv_coeff_ideals(O, copy = false)
  d = fmpz(1)
  for i = 1:n
    for j = 1:i
      d = lcm(d, denominator(simplify(PM.matrix[i, j]*PM.coeffs[i]*inv_coeffs[j])))
    end
  end
  return d
end

@doc Markdown.doc"""
    denominator(a::AlgAssRelOrdIdl) -> fmpz

Returns the smallest positive integer $d$ such that $da$ is contained in
`order(a)`.
"""
denominator(a::AlgAssRelOrdIdl) = denominator(a, order(a))

# Assumes that I is "locally integral at p", i. e. I_p \subseteq O_p.
# Returns x in R\setminus p such that Ix \subseteq O, where R = order(p)
function coprime_denominator(I::AlgAssRelOrdIdl{T1, T2}, O::AlgAssRelOrd{T1, T2}, p::Union{ NfAbsOrdIdl, NfRelOrdIdl }) where { T1, T2 }
  basis_O, basis_I, MO, MI = coprime_bases(O, I, p)
  OK = order(p)
  S = Dict{ideal_type(OK), Int}()
  for i = 1:degree(O)
    ai = inv(basis_O[i][2])
    for j = 1:degree(O)
      c = basis_I[j][2]*ai*MI[i, j]
      if iszero(norm(c))
        continue
      end
      d = denominator(c)
      facD = factor(d)
      for (q, e) in facD
        qdec = prime_decomposition(OK, q)
        for (Q, _) in qdec
          v = valuation(c, Q)
          if v >= 0
            continue
          end
          if haskey(S, Q)
            f = S[Q]
            S[Q] = max(f, -v)
          else
            S[Q] = -v
          end
        end
      end
    end
  end
  if haskey(S, p)
    error("The ideal is not locally integral at p")
  end
  primes = Vector{ideal_type(OK)}()
  vals = Vector{Int}()
  for (q, e) in S
    push!(primes, q)
    push!(vals, e)
  end
  push!(primes, p)
  push!(vals, 0)
  return approximate_nonnegative(vals, primes)
end

# Assumes N \subseteq M.
# Returns pseudo bases (a_i, alpha_i)_i of M and (b_i, beta_i)_i of N, a basis
# matrix of M_p and a basis matrix (b_{ij})_{i,j} of N_p in the basis of M_p,
# such that M_p = \bigoplus_i R_p alpha_i and N_p = \bigoplus_i R_p beta_i,
# where R = order(p) and M = \bigoplus a_i alpha_i and N = \bigoplus b_i beta_i,
# beta_i = sum_j b_{ij} alpha_j.
# Top level functions to avoid "type mix-ups" (like AlgAssRelOrd together with NfRelOrdIdl).
coprime_bases(O::NfRelOrd{S, T, U}, I::NfRelOrdIdl{S, T, U}, p::Union{ NfAbsOrdIdl, NfRelOrdIdl }) where { S, T, U } = _coprime_bases(O, I, p)

coprime_bases(I::NfRelOrdIdl{S, T, U}, J::NfRelOrdIdl{S, T, U}, p::Union{ NfAbsOrdIdl, NfRelOrdIdl }) where { S, T, U } = _coprime_bases(I, J, p)

coprime_bases(O::AlgAssRelOrd{S, T, U}, I::AlgAssRelOrdIdl{S, T, U}, p::Union{ NfAbsOrdIdl, NfRelOrdIdl }) where { S, T, U } = _coprime_bases(O, I, p)

coprime_bases(I::AlgAssRelOrdIdl{S, T, U}, J::AlgAssRelOrdIdl{S, T, U}, p::Union{ NfAbsOrdIdl, NfRelOrdIdl }) where { S, T, U } = _coprime_bases(I, J, p)

function _coprime_bases(M::Union{ NfRelOrd{S, T}, AlgAssRelOrd{S, T, U}, NfRelOrdIdl{S, T}, AlgAssRelOrdIdl{S, T, U} }, N::Union{ NfRelOrdIdl{S, T}, AlgAssRelOrdIdl{S, T, U} }, p::Union{ NfAbsOrdIdl, NfRelOrdIdl }) where { S, T, U }
  A = _algebra(M)
  @assert A === _algebra(N)
  OK = order(p)
  u = elem_in_nf(uniformizer(p), copy = false)
  iu = inv(u)
  basis_M = Vector{Tuple{elem_type(A), fractional_ideal_type(OK)}}()
  for (b, c) in pseudo_basis(M, copy = false)
    v = valuation(c, p)
    if v == 0
      push!(basis_M, (deepcopy(b), deepcopy(c)))
    elseif v < 0
      v = -v
      b = b*iu^v
      c = c*u^v
      push!(basis_M, (b, c))
    else
      b = b*u^v
      c = c*iu^v
      push!(basis_M, (b, c))
    end
  end

  basis_N = Vector{Tuple{elem_type(A), fractional_ideal_type(OK)}}()
  for (b, c) in pseudo_basis(N, copy = false)
    v = valuation(c, p)
    if v == 0
      push!(basis_N, (deepcopy(b), deepcopy(c)))
    elseif v < 0
      v = -v
      b = b*iu^v
      c = c*u^v
      push!(basis_N, (b, c))
    else
      b = b*u^v
      c = c*iu^v
      push!(basis_N, (b, c))
    end
  end

  mat_M = zero_matrix(_base_ring(A), dim(A), dim(A))
  for i = 1:dim(A)
    elem_to_mat_row!(mat_M, i, basis_M[i][1])
  end

  mat_N = zero_matrix(_base_ring(A), dim(A), dim(A))
  for i = 1:dim(A)
    elem_to_mat_row!(mat_N, i, basis_N[i][1])
  end
  mat_N = mat_N*inv(mat_M)

  return basis_M, basis_N, mat_M, mat_N
end

################################################################################
#
#  Preimage
#
################################################################################

function _as_ideal_of_smaller_algebra(m::AbsAlgAssMor, I::AlgAssRelOrdIdl)
  A = domain(m)
  B = codomain(m)
  K = base_ring(A)
  @assert K === base_ring(B)
  @assert dim(A) <= dim(B)
  @assert algebra(I) === B
  OA = maximal_order(A)
  # Transport OA to B
  M = zero_matrix(base_ring(B), dim(B), dim(B))
  c = NfOrdFracIdl[]
  PB = pseudo_basis(OA, copy = false)
  for i = 1:dim(A)
    t = m(PB[i][1])
    elem_to_mat_row!(M, i, t)
    push!(c, PB[i][2])
  end
  for i in (dim(A)+1):dim(B)
    push!(c, one(K) * maximal_order(K))
  end
  J = _intersect_modules(basis_pmatrix(I), pseudo_matrix(M, c), false)
  # Compute the intersection of M and I
  # Map the basis to A
  gen = elem_type(A)[]
  for i = 1:nrows(J)
    t = elem_from_mat_row(B, J.matrix, i)
    b, s = haspreimage(m, t)
    @assert b
    for el in absolute_basis(J.coeffs[i]) # replace with generators
      push!(gen, A(el) * s)
    end
  end
  J = ideal_from_lattice_gens(A, OA, gen)
  return J
end


text:
using Random: shuffle
using Distributed: @everywhere, @distributed
@everywhere include("gru_dynamic_struct.jl")

create_model_definitions(layer) =
begin
    eval(Meta.parse("@everywhere " * "struct Model\n" * *(["l$i::Layer\n" for i in 1:length(layers) +1]...) * "end"))
    eval(Meta.parse("@everywhere " * "(model::Model)(io) =\n" * "begin\n" * *(["io = model.l$i(io)\n" for i in 1:length(layers) +1]...) * "end"))
end

soften = arr -> (begin
    new_arr = []
    append!(new_arr, softmax(arr[1:Int(length(arr)/4)]))
    # soft = softmax(arr[1:Int(length(arr)/4)])
    # for i in 1:Int(length(arr)/4)
    #     push!(new_arr, soft[i])
    # end
    append!(new_arr, sigm.(arr[Int(length(arr)/4):Int(length(arr)*3/4)]))
    # for i in Int(length(arr)/4):Int(length(arr)*3/4)
    #     push!(new_arr, sigm(arr[i]))
    # end
new_arr
end)

prop(model, x) =
begin
    for mfield in fieldnames(Model)
        layer = getfield(model, mfield)
        layer.state = zeros(1, length(layer.bs))
    end
soften([model(t) for t in x][end])
end

cross_entropy(out, label) =
begin
    - sum(label .* log.(out))
end

train!(model, datas, lr) =
begin

    result = @distributed vcat for (x,y) in shuffle(datas)
        d = @diff cross_entropy(prop(model, x), y)
        grads = []
        for mfield in fieldnames(Model)
            layer = getfield(model, mfield)
            for lfield in fieldnames(Layer)
                push!(grads, grad(d, getfield(layer, lfield)))
            end
        end
        grads, value(d)
    end

    loss = 0.0
    for (g,l) in result
        loss += l
        i = 0
        for mfield in fieldnames(Model)
            layer = getfield(model, mfield)
            for lfield in fieldnames(Layer)
                i +=1
                setfield!(layer, lfield, Param(getfield(layer, lfield) - g[i] .* lr))
            end
        end
    end
    @show loss
end


text:
@testset "evaluation" begin
    cube = randn(rng, 10, 101, 101) .+ 100
    angles = 90 .* rand(rng, 10)
    res = Median()(cube, angles)
    res2 = Median()(cube, angles, cube)
    @test res == res2
    @test median(res) ≈ 0 atol = 1 / sqrt(length(res)) # within 1σ
end


text:
# this file is defined to test @changeprecision T include(...)
foo(x) = x/3


text:
using ProgressMeter, Dierckx, JLD

let
	splineknots = [fill(0.2, 4); 0.4:0.1:1.8; fill(2.0, 4)]

	splinecoeff(x, y) = Spline1D(x, y).c
	global interpolatespline(x, coeff) = Spline1D(splineknots, coeff, 3, 3, 0.0)(x)

	global function makecalibrationcache()
		lambdas = splineknots[1]:0.1:splineknots[end]
		years = 1800:10:2010
		ylen = length(years)
		states, forcingfactor, fertilization = precalchistory(lambdas, years)
		statemat = reshapestates(states)
		forccoeff = splinecoeff(lambdas, forcingfactor)
		fertcoeff = splinecoeff(lambdas, fertilization)
		statecoeff = Array{Float64,3}(undef, 75,length(forccoeff),ylen)
		for ivar = 1:75, iyear=1:ylen
			statecoeff[ivar,:,iyear] = splinecoeff(lambdas, statemat[ivar,:,iyear])
		end
		path = joinpath(@__DIR__, "..")
		JLD.save("$path/cachedclimatehistory.jld", "statecoeff", statecoeff, "forccoeff", forccoeff, "fertcoeff", fertcoeff)
	end

	function precalchistory(lambdas, years)
		len = length(lambdas)
		states = Vector{Vector{ClimateState}}(undef, len)
		forcingfactor = Vector{Float64}(undef, len)
		fertilization = Vector{Float64}(undef, len)
		progressmeter = Progress(len, 1)
		for i=1:len
			results, p = solveclimate(; lambda=lambdas[i], timestep=0.001, lastyear=2010, usecache=false)
			states[i] = results[iyear.(years)]
			forcingfactor[i] = p.aerosolforcingfactor
			fertilization[i] = p.fertilization
			next!(progressmeter)
		end
		return states, forcingfactor, fertilization
	end

	function reshapestates(states::Vector{Vector{ClimateState}})::Array{Float64,3}
		llen = length(states)
		ylen = length(states[1])
		out = Array{Float64,3}(undef, 75,llen,ylen)
		for ilambda=1:llen, iyear=1:ylen
			out[:,ilambda,iyear] = state2vector(states[ilambda][iyear])
		end
		out
	end

	# if !isfile("$path/cachedclimatehistory.jld2")
	# 	println("Precalculating historical model results...")
	# 	println("(This will take a minute or two but will make the model run much faster.)")
	# 	makecalibrationcache()
	# end

	path = joinpath(@__DIR__, "..")
	cache = JLD.load("$path/cachedclimatehistory.jld")
	global const cached_coeff_state = cache["statecoeff"]
	global const cached_coeff_forcing = cache["forccoeff"]
	global const cached_coeff_fertilization = cache["fertcoeff"]

end

max_stars_repo_name: niclasmattsson/ClimateCalculator

max_issues_repo_name: niclasmattsson/ClimateCalculator

max_forks_repo_name: niclasmattsson/ClimateCalculator


text:
#MAIN ENTRY POINT
#final columns for the decompressed file
prepcrsp(;refreshcrsp=true)::DataFrame = loadcrspdata(refreshcrsp=refreshcrsp)

function loadcrspdata(crspname=CRSP_NAME; datapath=DATA_PATH, refreshcrsp=true,
  incsvstream::Function = IN_CSV_STREAM,
  injlsstream::Function = IN_JLS_STREAM, outjlsstream::Function = OUT_JLS_STREAM)::DataFrame
  local p::String = "$datapath\\$crspname"
  local crsp::DataFrame

  #makes a serialization object so we don't have to keep parsing the CSV
  if refreshcrsp || (!isfile("$p.jls.lz4"))
    crsp = incsvstream("$p.csv.lz4") |> CSV.read |> DataFrame

    #println(describe(crsp))

    #preprocesscrsp!(crsp)

    crsp = preprocesscrsp(crsp)

    outjlsstream("$p.jls.lz4") do s
      serialize(s, crsp)
    end

  else
    injlsstream("$p.jls.lz4") do s
      crsp = deserialize(s)
    end

  end

  @info "CRSP data loaded and saved into file $crspname.jls.lz4."

  return crsp
end

function preprocesscrsp(crsp::DataFrame;
  crspdateformat::DateFormat = CRSP_DATE_FORMAT,
  crspcolumns::Vector{Symbol} = RETAINED_COLS_CRSP,
  crsptype::Symbol = CRSP_TYPE,
  testoutput::Bool = TEST_OUTPUT,
  usedcrspfields::Vector{Symbol} = USED_CRSP_FIELDS,
  time2stale::T = TIME_2_STALE)::DataFrame where T<:DatePeriod

  local scrsps::GroupedDataFrame

  #rename columns to lower case

  names!(crsp, (s::Symbol->Symbol(lowercase(string(s)))).(names(crsp)))

  #parse some columns
  crsp.ret = (s::MString->parsecrsp(Float64, s)).(crsp.ret)
  crsp.dlret= (s::MString->parsecrsp(Float64, s)).(crsp.dlret)
  crsp.siccd = (s::MString->parsecrsp(Int, s)).(crsp.siccd)

  rename!(crsp, :ret=>:retlisted)
  #make sure we capture delisting returns

  #treat -100% returns as missing values
  crsp.retlisted .= (r->
    !ismissing(r) && (r > -1.0) ? r : missing).(crsp.retlisted)
  crsp.dlret .= (r->
    !ismissing(r) && (r > -1.0) ? r : missing).(crsp.dlret)

  crsp.ret = ((retlisted,dlret)->
    ismissing(retlisted) && ismissing(dlret) ? missing :
      ismissing(dlret) ? retlisted :
      ismissing(retlisted) ? dlret :
      (1+retlisted)*(1+dlret) - 1.
    ).(crsp.retlisted,crsp.dlret)

  #drop rows with missing data
  #below stopped working
  #deleterows!(crsp, (!).(completecases(crsp[:,[:ret, :SHRCD, :EXCHCD, :SICCD]])))

  crsp = crsp[completecases(crsp[:,[:ret, :shrcd, :exchcd, :siccd]]),:]

  #NOTE: SHRCD==10,11 : US Common Shares
  #NOTE:  :EXCHCD ∈ [1,2,3] : NYSE, Amex, NASDAQ
  #NOTE: r[:SICCD] ∉ 6000:6999 : No financial companies
  #this improves performance
  local shrcdvalid::Set{Int} = Set{Int}((10,11))
  local exchcdvalid::Set{Int} = Set{Int}((1,2,3))
  local siccdinvalid::Set{Int} = Set{Int}(6000:6999)
  filter!(r::DataFrameRow->
    (r[:shrcd] ∈ shrcdvalid) &&
    (r[:exchcd] ∈ exchcdvalid ) &&
    (r[:siccd] ∉ siccdinvalid),
    crsp)

  #Now work on the price
  #in this case, use the mean of the two values
  #price cannot be zero
  crsp.prc .= (
    prc->coalesce(iszero(prc),true) ? missing : prc).(crsp.prc)
  crsp.dlprc .= (
    dlprc->coalesce(iszero(dlprc),true) ? missing : dlprc).(crsp.dlprc)

  crsp.price = ((prc,dlprc)->
    ismissing(prc) && ismissing(dlprc) ? missing :
      ismissing(dlprc) ? prc :
      ismissing(prc) ? dlprc :
      (prc+dlprc)/2
  ).(crsp.prc, crsp.dlprc)

  #add in the alternative price if nothing else is available
  if crsptype == :monthly
    crsp.altprc .= (
      altprc->coalesce(iszero(altprc),true) ? missing : altprc).(crsp.altprc)
    crsp.altprc .= (abs).(crsp.altprc)
    crsp.price .= ((price,altprc)->
      ismissing(price) ? altprc : price).(crsp.price,crsp.altprc)
  end

  #negative values computed via bid-ask spread method
  crsp.price .= (abs).(crsp.price)

  #adjust the units and compute the market cap
  crsp.mktcap = Vector{MInt}((round).(crsp.price .* crsp.shrout .* 1000))
  crsp.mktcap .= (i::MInt->coalesce(iszero(i),true) ? missing : i).(crsp.mktcap)

  #println(describe(crsp))
  crsp.date = (s->parsecrsp(Date, s)).(crsp.date)
  #println(describe(crsp))

  #get log returns
  crsp.lret = (log).(1.0 .+ crsp.ret)

  #now dedup crsp
  dedupctr::Int = 0
  crsp.keep = trues(size(crsp,1))

  scrsps = groupby(crsp, [:permno, :date])
  @mpar for i ∈ 1:length(scrsps)
    local scrsp::SubDataFrame = scrsps[i]
    if size(scrsp,1) > 1
      #first try to pick records by number missing
      nmissing::Vector{Int} = zeros(Int, size(scrsp,1))
      for f ∈ usedcrspfields
        nmissing .+= (ismissing).(scrsp[!,f])
      end
      minmissing::Int = minimum(nmissing)
      scrsp.keep .= nmissing .== minmissing

      #keep the top record iff it is identical to the others
      if sum(scrsp.keep) > 1
        sscrsp = view(scrsp, scrsp.keep, :)
        nchecks::Int = size(sscrsp,1)
        keep1::Bool = true
        for i ∈ 2:nchecks
          for f ∈ usedcrspfields
            #checks if the given used field is equal to the top row
            (coalesce(sscrsp[i,f],-99.0) ≠ coalesce(sscrsp[1,f], -99.0)) && (keep1 = false)
          end
        end

        sscrsp.keep .= false
        sscrsp.keep[1] = keep1
      end
    end
  end

  crsp = crsp[crsp.keep,:]

  #need to lag the marketcap
  crsp.Lmktcapdate = Vector{MDate}(undef, size(crsp,1))
  crsp.Lmktcap = Vector{MInt}(undef, size(crsp,1))

  sort!(crsp, [:permno, :date])
  scrsps = groupby(crsp, :permno)
  @mpar for i ∈ 1:length(scrsps)
    local scrsp::SubDataFrame = scrsps[i]

    Nscrsp::Int = size(scrsp,1)
    for i ∈ 2:Nscrsp
      local Lmktcap::MInt
      local Lmktcapdate::MDate

      #get the most recent market cap other than the current
      if !ismissing(scrsp[i-1, :mktcap])
        Lmktcap = scrsp[i-1, :mktcap]
        Lmktcapdate = scrsp[i-1, :date]
      else
        Lmktcap = scrsp[i-1, :Lmktcap]
        Lmktcapdate = scrsp[i-1, :Lmktcapdate]
      end

      #this si the check against stale data
      if (!ismissing(Lmktcapdate)) && (Lmktcapdate + time2stale ≥ scrsp[i,:date])
        scrsp[i, :Lmktcap] = Lmktcap
        scrsp[i, :Lmktcapdate] = Lmktcapdate
      end
    end
  end

  crsp.lLmktcap = (log).(crsp.Lmktcap)
  select!(crsp, crspcolumns)

  if testoutput
    d = Date(2015,12,31)
    testoutput && CSV.write("output\\crsp_$d.csv", crsp[(cd->cd==d).(crsp.date), :])
  end

  return crsp
end

#make tryparse handle missing values
#Base.tryparse(::Type{<:Any}, v::Missing)::Nothing = nothing

#parse a value where the type is known: Generic case
function parsecrsp(::Type{T}, s::String;
    parsedmissings = CRSP_PARSED_MISSINGS) where T

  #if it doesn't parse to the right type, set to missing
  v::Union{T,Missing} = something(tryparse(T, s), missing)

  #check against the list of missing codes
  (!ismissing(v)) && (v ∈ parsedmissings) && (v=missing)

  return v
end

#the date case
parsecrsp(::Type{Date}, s::String;
    crspdateformat::DateFormat = CRSP_DATE_FORMAT) = Dates.Date(s, crspdateformat)

#helper methods and special cases
parsecrsp(::Type{<:Any}, v::Missing) = missing
parsecrsp(::Type{Date}, i::Int) = parsecrsp(Date, "$i")

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA