text:
export Apply

"""
    Apply{J, O} <: SFunc{Tuple{SFunc{J, O}, J}, O}

Apply represents an sfunc that takes two groups of arguments.  The first group is a single
argument, which is an sfunc to apply to the second group of arguments.

# Additional supported operators
- `support`
- `support_quality`
- `sample`
- `logcpdf`
- `compute_pi`
- `send_lambda`

# Type parameters
- `J`: the input type of the *sfunc* that may be applied; that *sfunc* is the input type of the `Apply`
- `O`: the output type of the *sfunc* that may be applied, which is also the output type of the `Apply`
"""
struct Apply{J <: Tuple, O} <: SFunc{Tuple{SFunc{J, O}, J}, O}
end

@impl begin
    struct ApplySupport end
    function support(::Apply{J,O},
                    parranges::NTuple{N,Vector},
                    size::Integer,
                    curr::Vector{<:O}) where {J<:Tuple,O,N}

        result = Vector{O}()
        for sf in parranges[1]
            append!(result, support(sf, (parranges[2],), size, curr))
        end
        return unique(result)
    end
end

@impl begin
    struct ApplySupportQuality end

    function support_quality(::Apply{J,O}, parranges) where {J,O}
        q = support_quality_rank(:CompleteSupport)
        for sf in parranges[1]
            imp = get_imp(MultiInterface.get_policy(), Support, sf, parranges[2], 0, O[])
            q = min(q, support_quality_rank(support_quality(imp, sf, parranges[2])))
        end
        return support_quality_from_rank(q)
    end
end

@impl begin
    struct ApplySample end
    function sample(::Apply{J,O}, input::Tuple{SFunc{J,O}, J})::O where {J<:Tuple,O}
        return sample(input[1], input[2])
    end
end

@impl begin
    struct ApplyLogcpdf end
    function logcpdf(::Apply{J,O}, i::Tuple{SFunc{J,O}, J}, o::O)::AbstractFloat where {J<:Tuple,O}
        return logcpdf(i[1], i[2], o)
    end
end

# WARNING: THIS LOGIC DOES NOT WORK WITH MORE THAN ONE PARENT
@impl begin
    struct ApplyComputePi end
    function compute_pi(::Apply{J,O},
                     range::__OptVec{<:O},
                     parranges::NTuple{N,Vector},
                     incoming_pis::Tuple)::Dist{<:O} where {N,J<:Tuple,O}

        sfrange = parranges[1]
        argsrange = parranges[2]
        sfpi = incoming_pis[1]
        argspi = incoming_pis[2]
        result = zeros(Float64, length(range))
        for sf in sfrange
            p1 = cpdf(sfpi, (), sf)
            p2 = compute_pi(sf, range, (argsrange,), (argspi,))
            p3 = [p1 * cpdf(p2, (), x) for x in range]
            result .+= p3
        end
        return Cat(range, result)
    end
end

# WARNING: THIS LOGIC DOES NOT WORK WITH MORE THAN ONE PARENT
@impl begin
    struct ApplySendLambda end

    function send_lambda(::Apply{J,O},
                       lambda::Score{<:O},
                       range::__OptVec{<:O},
                       parranges::NTuple{N,Vector},
                       incoming_pis::Tuple,
                       parent_idx::Integer)::Score where {N,J<:Tuple,O}
        @assert parent_idx == 1 || parent_idx == 2
        sfrange = parranges[1]
        argsrange = parranges[2]
        sfpi = incoming_pis[1]
        argspi = incoming_pis[2]

        if parent_idx == 2
            # For each x, we must sum over the sfunc argument compute P(y|x) for each possible sfunc
            result = Vector{Float64}()
            for args in argsrange
                resultpieces = Vector{Float64}()
                for sf in sfrange
                    sp = logcpdf(sfpi, (), sf)
                    for y in range
                        a = isa(args, Tuple) ? args : tuple(args)
                        push!(resultpieces, sp + logcpdf(sf, a, y) + get_log_score(lambda, y))
                    end
                end
                push!(result, logsumexp(resultpieces))
            end
            return LogScore(argsrange, result)

        else # parent_idx == 1
            # This is simpler; we must sum over the arguments, which is achieved by the embedded compute_pi
            result = Vector{Float64}()
            for sf in sfrange
                resultpieces = Vector{Float64}()
                ypi = compute_pi(sf, range, (argsrange,), (argspi,))
                for y in range
                    push!(resultpieces, logcpdf(ypi, (), y) + get_log_score(lambda, y))
                end
                push!(result, logsumexp(resultpieces))
            end
            return LogScore(sfrange, result)
        end
    end
end


text:
module MatrixEquations

const BlasFloat = Union{Float64,Float32,ComplexF64,ComplexF32}
const BlasReal = Union{Float64,Float32}
const BlasComplex = Union{ComplexF64,ComplexF32}

using LinearAlgebra
using LinearAlgebra: require_one_based_indexing
import LinearAlgebra: mul!
using LinearMaps

include("lapackutil.jl")
using .LapackUtil: tgsyl!, lanv2, ladiv, lag2, lacn2!, safemin, smlnum

export MatrixEquationsMaps
export utqu, utqu!, qrupdate!, rqupdate!, isschur, triu2vec, vec2triu, utnormalize!
export lanv2, ladiv, lag2, lacn2!
export lyapc, lyapd, lyapcs!, lyapds!
export plyapc, plyaps, plyapcs!, plyapd, plyapds!
export arec, ared, garec, gared
export sylvc, sylvd, gsylv, sylvcs!, sylvds!, gsylvs!
export sylvsys, dsylvsys, sylvsyss!, dsylvsyss!, tgsyl!
export sylvckr, sylvdkr, gsylvkr, sylvsyskr, dsylvsyskr
export opnorm1, opnorm1est, oprcondest, opsepest
export lyapop, invlyapop, sylvop, invsylvop, sylvsysop, invsylvsysop, trmatop

include("meutil.jl")
include("sylvester.jl")
include("lyapunov.jl")
include("riccati.jl")
include("sylvkr.jl")
include("plyapunov.jl")
include("meoperators.jl")
include("condest.jl")

end

ext: jl

max_stars_repo_name: andreasvarga/MatrixEquations

max_issues_repo_name: andreasvarga/MatrixEquations

max_forks_repo_name: andreasvarga/MatrixEquations


text:
# DTU02435 Decision Making under Uncertainty
# Assignment 3, Task 2
# Out-of-Sample Test of the Plan
# Edward J. Xu
# May 5th, 2019
push!(LOAD_PATH, "$(homedir())/Desktop/StochasticAirplaneRent")
cd("$(homedir())/Desktop/StochasticAirplaneRent")
using JuMP
# using GLPKMathProgInterface
# using Gurobi
using CPLEX
using CSV
using DataFrames

function processData()
    datf_cost = CSV.read("Data/OperationCosts.csv", header = true, delim = ',')
    datf_cap = CSV.read("Data/PlaneCapacity.csv", header = ["p", "cap"], delim = ',')
    datf_route = CSV.read("Data/DemandsRoute.csv", header = true, delim = ',')
    datf_costAdmin = CSV.read("Data/CostAdquisition.csv", header = ["p", "cost"], delim = ',')
    datf_demand = CSV.read("Data/OutofSample.csv", header = true, delim = ',')
    # Calculate the demand matrix
    mat4_demand = zeros(4, 10, 12, 500)
    for i = 1: 4
        for j = 1: 10
            for w = 1: 500
                for t = 1: 12
                    mat4_demand[i, j, t, w] = datf_demand[t, w+1] * datf_route[i, j+1]
                end
            end
        end
    end
    return (datf_cost, datf_cap.cap, datf_costAdmin.cost, mat4_demand)
end

function optim(vec_y, datf_cost, vec_cap, vec_costAdmin, mat_demand)
    model = Model(solver = CplexSolver(CPX_PARAM_SCRIND = 0))
    # model = Model(solver = GurobiSolver(Presolve=0))
    # model = Model(solver = GLPKSolverMIP())
    @variable(model, mat3_z[1:4, 1:10, 1:10] >= 0, Int)
    @variable(model, mat_m[1:4, 1:10] >= 0)
    @objective(model, Min, sum(datf_cost.cost[.&((datf_cost[:i] .== i), (datf_cost[:j] .== j),
        (datf_cost[:p] .== p))][1] * mat3_z[i, j, p] for p = 1:10, i = 1:4, j = 1:10) + sum(mat_m) * 1000)
    @constraint(model, [i = 1:4, j = 1:10], mat_demand[i, j] - sum(mat3_z[i, j, p] * vec_cap[p] for p = 1:10)
        <= mat_m[i, j])
    @constraint(model, [p = 1:10], sum(mat3_z[i, j, p] for i = 1:4, j = 1:10) <= vec_y[p])
    solve(model)
    result_obj = getobjectivevalue(model)
    return result_obj
end

function testOutOfSample(vec_y, datf_cost, vec_cap, vec_costAdmin, mat4_demand)
    mat_result_obj = zeros(500, 12)
    vec_result_obj = zeros(500)
    ## 1,  Begin Optimization
    for w = 1:500
        for t = 1:12
            mat_result_obj[w, t] = optim(vec_y, datf_cost, vec_cap, vec_costAdmin, mat4_demand[:, :, t, w])
        end
        vec_result_obj[w] = sum(mat_result_obj[w, tt] for tt = 1:12) + sum(vec_y[p] * vec_costAdmin[p] for p = 1: 10)
        println("Result: obj = $(vec_result_obj[w]), when w = $(w).")
    end
    return vec_result_obj
end

function main()
    (datf_cost, vec_cap, vec_costAdmin, mat4_demand) = processData()
    println("################################################################################\n",
            "############### 1/4, Out-of-Sample Test of Deterministic Program ###############\n",
            "################################################################################")
    vec_y_1 = [0.0, 1.0, 2.0, 2.0, 0.0, 12.0, 5.0, 4.0, 8.0, 8.0]
    timeStart = time()
    vec_result_obj_1 = testOutOfSample(vec_y_1, datf_cost, vec_cap, vec_costAdmin, mat4_demand)
    println("Elapsed time: $(time()-timeStart) seconds.")
    CSV.write("result/out-of-sample_1.csv", DataFrame(vec_result_obj_1'), writeheader = false, delim = ',')
    println("################################################################################\n",
            "################# 2/4, Out-of-Sample Test of Stochastic Program ################\n",
            "################################################################################")
    vec_y_2 = [0.0, 7.0, 2.0, 3.0, 0.0, 12.0, 4.0, 4.0, 8.0, 8.0]
    timeStart = time()
    vec_result_obj_2 = testOutOfSample(vec_y_2, datf_cost, vec_cap, vec_costAdmin, mat4_demand)
    println("Elapsed time: $(time()-timeStart) seconds.")
    CSV.write("result/out-of-sample_2.csv", DataFrame(vec_result_obj_2'), writeheader = false, delim = ',')
    println("################################################################################\n",
            "################## 3/4, Out-of-Sample Test of Robust Program ###################\n",
            "################################################################################")
    vec_y_3 = [0.0, 3.0, 1.0, 2.0, 0.0, 12.0, 2.0, 8.0, 8.0, 8.0]
    timeStart = time()
    vec_result_obj_3 = testOutOfSample(vec_y_3, datf_cost, vec_cap, vec_costAdmin, mat4_demand)
    println("Elapsed time: $(time()-timeStart) seconds.")
    CSV.write("result/out-of-sample_3.csv", DataFrame(vec_result_obj_3'), writeheader = false, delim = ',')
    println("################################################################################\n",
            "################################### 4/4, End ###################################\n",
            "################################################################################")
end

main()


text:
#Run all the tests together

require("Ito/test/calendars.jl")
require("Ito/test/statistics.jl")
require("Ito/test/integration.jl")
require("Ito/test/term_structure.jl")
require("Ito/test/currencies.jl")


text:
module FunctionDependencies

const dependency_table = Dict{Symbol,Symbol}(
    :√ => :sqrt,
    :adjoint => :conj,
    :< => :isless,
    :> => :<,
    :isgreater => :isless,
    :ismore => :>,
    :≥ => :(>=),
    :≤ => :(<=),
    :min => :isless,
    :max => :isless,
    :cmp => :isless,
    :isequal => :(==),
    :(!=) => :(==),
    :≠ => :(!=),
    :(!==) => :(===),
    :≢ => :(!==),
    :≡ => :(===),
    :⊻ => :xor,
    :⊼ => :nand,
    :⊽ => :nor,
    :% => :rem,
    :÷ => :div,
    :mod1 => :mod,
    :∈ => :in,
)

maybe_add_dependent_funccall!(funccalls::Set{Symbol}, ::Nothing) = funccalls
function maybe_add_dependent_funccall!(funccalls::Set{Symbol}, call)
    push!(funccalls, call)
    maybe_add_dependent_funccall!(funccalls, get(dependency_table, call, nothing))
    funccalls
end

function maybe_add_dependent_funccalls!(funccalls::Set{Symbol})
    calls_to_add = intersect(keys(dependency_table), funccalls)
    for call in calls_to_add
        maybe_add_dependent_funccall!(funccalls, dependency_table[call])
    end
    funccalls
end

end


text:
# A script for loading unregistered dependencies before CI action
using Pkg

pkg"dev https://github.com/kylejbrown17/GraphUtils.jl.git"


text:
module transitions

using conv_code
using encoder
using utils
using code_structs

export transition_list

#=
Esta função deve retornar algo do tipo Array{Transition},
em que o Transition.in_bit é um bit de entrada, o Transition.output
é a saída gerada pelo tal bit de entrada e Transition.next_state é
o próximo estado (notação binária)

O 'origin_state' correspondendo ao estado com:
	- todas as memórias nulas é 1
	- todas as memórias=1 é 2^(code.quant_mem)
=#
function transition_list(code::ConvolutionalCode, origin_state::Int64)
	if haskey(transitions_cache, code)
		return transitions_cache[code][origin_state]
	end

	transition_list = Array{Array{Transition, 1}, 1}(undef, 2^code.quant_mem)

	for i=0:2^(code.quant_mem) - 1
		transition_list[i+1] = []
		binary = dec2bin(i, code.quant_mem+1)
		binary[1] = 0
		push!(transition_list[i+1], Transition(0, current_out(code, binary), i >> 1 + 1))
		binary[1] = 1
		push!(transition_list[i+1], Transition(1, current_out(code, binary),
			2^(code.quant_mem-1) + (i >> 1) + 1
		))
	end

	transitions_cache[code] = transition_list

	return transitions_cache[code][origin_state]
end

transitions_cache = Dict()

end


text:
include("../gcc_common.jl")
build_and_upload_gcc(v"4.8.5")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
"""
    refine_triangulation_images(triang_vertices,
								triang_simplex_indices,
								splitting_rules,
                                image_vertices,
								split_indices = Int[])

Refine a triangulation (composed of simplices), also updating the images of the
simplices in the triangulation.

- `triang_vertices::Array{Float64, 2} is a
    (n_triangulation_vertices x embeddingdim) sized array where each row is a
    vertex of the triangulation.
`triang_simplex_indices::Array{Int, 2}` is an array of size
    (n_trinagulation_simplices x (embeddingdim+1)). Each row of simplex_indices contains
    the indices of the vertices (rows of the vertices array) furnishing the corresponding
    simplex.
`image_vertices::Array{Float64, 2} is a (n_triangulation_vertices x embeddingdim) sized
        array where each row is a vertex of the triangulation.
`split_indices::Array{Int, 1}` are the row numbers of simplex_indices indicating which
    simplices should be split.
"""
function refine_triangulation_images(triang_vertices, triang_simplex_indices,
                                    image_vertices, split_indices, k)

    if length(split_indices) == 0
        return triang_vertices, triang_simplex_indices, image_vertices
    end
    # The number of simplices to split
    n_split_simplices = length(split_indices)

    # The dimension of the space
    E = size(triang_vertices, 2)

    # Rules for forming the strictly new vertices of the subtriangulation
    splitting_rules = simplicial_subdivision(k, E)
    rules = splitting_rules[1]

    # Array where each row represents one of the new simplices in the splitted simplex
    subtriangulation = splitting_rules[2]

    # How many new vertices are created each split?
    n_newvertices_eachsplit = size(rules, 1)

    # We need an array that can accomodate all of them. Each row in this
    # array will be a new vertex. Stacks of n_newvertices_eachsplit * E arrays.
    # We have as many as we have simplices to split.
    new_vertices = zeros(n_newvertices_eachsplit * n_split_simplices, E)
    new_imagevertices = zeros(n_newvertices_eachsplit * n_split_simplices, E)

    # Fill the array by looping over all simplices that we want to split
    for i = 1:n_split_simplices
        # Figure out what the row indices corresponding to the ith simplex
        # must be. Marks the beginning of each of the simplex stacks in new_vertices
        ind = n_newvertices_eachsplit * (i -1)

        # Index of the simplex we need to split
        simplex_idx = split_indices[i]

        # Get the vertices of the simplex currently being splitted. Each of the
        # n_newvertices_eachsplit new vertices will be a linear combination
        # of these vertices. Each row is a vertex.
        vertices = triang_vertices[triang_simplex_indices[simplex_idx, :], :]
        imagevertices = image_vertices[triang_simplex_indices[simplex_idx, :], :]

        # Generate the strictly new vertices for each sub
        for j = 1:n_newvertices_eachsplit
            # The index for a particular new vertex of the ith new simplex
            ind_newvertex = ind + j

            # Compute the jth vertex of the subtriangulation of the ith simplex

            # Go the jth new subsimplex of the ith splitted simplex. The entries of this
            # vector
            jth_subsimplex = rules[j, :]

            # Pick the corresponding original vertices with indices contained in rules[j, :]
            original_vertices = vertices[rules[j, :], :]
            original_vertices_image = imagevertices[rules[j, :], :]

            new_vertices[ind_newvertex, :] = sum(original_vertices, dims=1) ./ k
            new_imagevertices[ind_newvertex, :] = sum(original_vertices_image, dims=1) ./ k

        end
    end

    # Find the unique new vertices
    new_vertices_noreps = unique(new_vertices, 1)
    new_imagevertices_noreps = Array{Float64}(size(new_vertices_noreps, 1), E)

    Ind = Vector{Int}(size(new_vertices, 1))
    for i = 1:size(new_vertices_noreps, 1)

        count = 0

        for j = 1:size(new_vertices, 1) # run over the possibly repeated vertices
            if new_vertices_noreps[i, :] == new_vertices[j, :]
                Ind[j] = i

                count = count + 1
                if count == 1
                    new_imagevertices_noreps[i, :] = new_imagevertices[j, :]
                end
            end
        end
    end

    # TEST:  assert all zeros.
    #@show unique(new_imagevertices, 1) - new_imagevertices_noreps

    # Combine old and newly introduced vertices
    num_vertices_beforesplit = size(triang_vertices, 1)
    triang_vertices = vcat(triang_vertices, new_vertices_noreps)
    image_vertices = vcat(image_vertices, new_imagevertices_noreps)

    # Update the Ind array, so that we start at the new vertices
    Ind = Ind + num_vertices_beforesplit
    num_simplices_each_split = size(subtriangulation, 1)

    # The subsimplices formed by the splitting. Each row contains E + 1 indices referencing
    # the vertices furnishing that particular subsimplex (now found in the updated
    # triang_vertices array).
    newtriangulation = Array{Float64}(num_simplices_each_split * n_split_simplices, E + 1)

    # For each simplex that we need to split,
    for i = 1:n_split_simplices
        # The beginning of the stack we need to fill.
        index = num_simplices_each_split * (i - 1)

        # Figure out what the row indices corresponding to the ith simplex
        # must be. Marks the beginning of each of the simplex stacks in new_vertices
        ind = n_newvertices_eachsplit * (i - 1)

        # Index of the simplex we need to split
        simplex_idx = split_indices[i]

        # Pick the indices of the original vertices. Should be a column vector.
        inds_original_vertices = triang_simplex_indices[simplex_idx, :]

        # Indices of the new vertices. Should be a column vector
        inds_new_vertices = Ind[(ind + 1):(ind + n_newvertices_eachsplit)]
        inds_all_vertices = vcat(inds_original_vertices, inds_new_vertices)

        # Go through each of the subsimplices formed by splitting that simplex
        for j = 1:num_simplices_each_split
            newtriangulation[index + j, :] = inds_all_vertices[subtriangulation[j, :]]
        end
    end

    # Indices of the simplices that are not split
    untouched_indices = complementary(split_indices, size(triang_simplex_indices, 1))

    triang_simplex_indices = round.(Int, vcat(triang_simplex_indices[untouched_indices, :],
                                  newtriangulation))

    return triang_vertices, triang_simplex_indices, image_vertices, untouched_indices
end


text:
# Simulate a classic beamformer with a circular array

# before running:
#  include("../src/DirectionFinding.jl")
#  using <favorite plotting program>

using DirectionFinding

begin
    fc = 1e9
    λ = 299.79e6/fc
    N₀ = 0.01
    a1 = pi/8
    ca = circulararray(11, 2λ, σ²=0.0001)
    s1 = SignalSource(a1, halfsine_train(1e6))
    sim = Simulation(ca, fc, N₀, s1)
    doa = cbf(sim, 2.00013e9, 100)
    println("Source Angle: $a1")
    println("Found: $(doa > 0 ? doa : doa+2π)")
end


text:
module GroundDataFVG
"""
Module for the download and processing of atmospheric data gathered by measuring stations located in Friuli Venezia Giulia, Italy.
"""
#=
Accesso risorse Friuli:
    https://www.dati.friuliveneziagiulia.it/resource/qp5k-6pvm.csv

    Qualità dell'aria:
        PM10:                       qp5k-6pvm
        PM2.5:                      d63p-pqpr
        Ozono:                      7vnx-28uy
        Monossido di carbonio:      t274-vki6
        Biossido di zolfo:          2zdv-x7g2
        Biossido di azoto:          ke9b-p6z2
        Pollini:                    svph-8w2g

    Meteo:
        Elenco sensori:             498i-2j88
        Previsioni e dati:          j654-ykm6 (Probabilmente non funzionante, il link è esterno e non ha il download)
        Dati stazioni:              4wxn-35av (Probabilmente non funzionante, il link è esterno e non ha il download)

Link:
    Dati stazioni:
        https://dev.meteo.fvg.it/xml/stazioni/GRA.xml
=#

using CombinedParsers
using CombinedParsers.Regexp
using CSV
using DataFrames
using Dates
using HTTP
using JSONTables

export getData,
       getRegionAttributes, getRegionIds, getRegionStationInfo

@syntax meteo_data = Repeat(
   "<",
   re"[^ >]+",
   Optional( re" [^ =]+=\"[^ \"]+\"" ),
   Optional( re" [^ =]+=\"[^\"]+\"" ),
   ">",
   Either( Numeric(Int64), Numeric(Float64), re"[^>]+" ),
   "</", re"[^>]+", ">"
)

"""
    getRegionAttributes( [ type::Symbol=:METEO ] )

Obtain the names of the columns of the region's dataframe required by `GroundData.createMap`'s `attributes` parameter to create `GroundData.standardize`'s `map` parameter.
"""
function getRegionAttributes( type::Symbol=:METEO )
    return type == :METEO ? [ :param, :unit, :value, nothing, :observation_time, :longitude, :latitude, :station_altitude, nothing, nothing, :rmh ] :
        type == :AIRQUALITY ? [ :parametro, :unita_misura, :value, nothing, :data_misura, :longitudine, :latitudine, nothing, :dati_insuff, nothing ] :
            throw(DomainError(type, "`type` must be either `:METEO` OR `:AIRQUALITY`."))
end

"""
    getRegionAttributes( [ type::Symbol=:METEO ] )

Obtain the names of the columns of the dataframe required for `GroundData.standardize`'s `bridge` parameter.
"""
function getRegionIds( type::Symbol=:METEO )
    return type == :METEO ? :nome :
        type == :AIRQUALITY ? nothing :
            throw(DomainError(type, "`type` must be either `:METEO` OR `:AIRQUALITY`."))
end

"""
    getRegionStationInfo( [ type::Symbol=:METEO  ] )

Obtain the names of the columns of the region's stations dataframe required by `GroundData.createMap`'s `attributes` parameter to be used in `GroundData.generateUuidsTable`.
"""
function getRegionStationsInfo( type::Symbol=:METEO )
    return type == :METEO ? [ nothing, :nome, :longitude, :latitude ] :
        type == :AIRQUALITY ? [ nothing, :ubicazione, :longitudine, :latitudine ] :
            throw(DomainError(type, "`type` must be either `:METEO` OR `:AIRQUALITY`."))
end

"""
    getMeteoStationsData()

Obtain a `DataFrame` describing the meteo stations in Friuli Venezia Giulia, Italy.
"""
function getMeteoStationsData()
    return CSV.read(split(@__DIR__, "src")[1]*"resources\\Ground stations data\\stazioni_meteoclimatiche-FVG.csv", DataFrame)
end

#   prec_type
#       0:nulla; 1:pioggia; 2:pioggia e neve; 3:neve
#   cloudiness
#       0:n.d.; 1:sereno; 2:poco nuvoloso; 3:variabile; 4:nuvoloso; 5:coperto
"""
    getMeteoData()

Obtain the data of the meteorological stations in Friuli Venezia Giulia as a `DataFrame`.
"""
function getMeteoData()
    resources = [ "ARI", "BAR", "BGG", "BIC", "BOA", "BOR", "BRU", "CAP", "CDP", "CER", "CHI", "CIV", "CMT", "COD", "COR",
                  "ENE", "FAG", "FOS", "FSP", "GEM", "GRA", "GRG", "GRM", "LAU", "LIG", "LSR", "MAT", "MGG", "MNF", "MUS",
                  "PAL", "PDA", "PIA", "213200", "POR", "PRD", "RIV", "SAN", "SGO", "SPN", "TAL", "TAR", "TOL", "TRI",
                  "UDI", "VIV", "ZON" ]
    #   prec_type = [ "nulla", "pioggia", "pioggia e neve", "neve" ]
    #   cloudiness = [ "n.d.", "sereno", "poco nuvoloso", "variabile", "nuvoloso", "coperto" ]
    data_str = String[]
    for res in resources
        try
            page = HTTP.get("https://dev.meteo.fvg.it/xml/stazioni/$res.xml")
            push!( data_str, String(page.body) )
        catch e
            if !isa( e, HTTP.ExceptionRequest.StatusError )
                throw(e)
            end
        end
    end
    data_split = @. replace( replace( getindex( split( data_str, r"</?meteo_data>" ), 2 ), r"\n *" => "" ), r"<!--[^-]+-->" => "" )
    data_parse = meteo_data.(data_split)
    vect = Dict[]
    for data in data_parse
        # Attributes of a single station, they are shared between all the parameters measured by the station
        others = [
                    Symbol( String( attribute[2] ) ) => attribute[6] isa Number ? attribute[6] : String( attribute[6] )
                    for attribute in data[1:6]
                 ]
        for attribute in data[7:end]
            # Parameters measured by each station
            dict = push!(
                       Dict(
                          :param => !ismissing( attribute[4] ) ? String( attribute[4][5] ) : String( attribute[2] ),
                          :value => attribute[6] isa Number ? attribute[6] : String( attribute[6] ),
                          :unit => !ismissing(attribute[3]) ? String( attribute[3][5] ) : missing
                       ),
                       others...
                   )
            push!( vect, dict )
        end
    end
    df = select!( DataFrame(vect), Not(2, 7) )
    rel_heights = [ length( split( param, " a " ) ) == 2 ? split( param, " a " )[2] : "0m" for param in df[:, :param] ]
    transform!(
        df,
        [:station_name] => ByRow( x -> x = uppercase(x) ) => :nome,
        [:observation_time] => ByRow( x -> DateTime( x[1:14], "dd/mm/yyyy H.M" ) ) => :observation_time
    )
    insertcols!( df, :rmh => rel_heights )
    return df
end

"""
    getAQData()

Obtain data on airquality gathered form measuring stations in Friuli Venezia Giulia.
"""
function getAQData()
    codes = [ "qp5k-6pvm" , "d63p-pqpr", "7vnx-28uy", "t274-vki6", "2zdv-x7g2", "ke9b-p6z2" ]
    params = [ "PM10", "PM2.5", "Ozono", "Monossido di carbonio", "Biossido di zolfo", "Biossido di Azoto" ]
    val_types = [ "media_giornaliera", "media_giornaliera", "media_oraria_max", "media_mobile_8h_max", "media_giornaliera", "media_oraria_max" ]
    data = [ DataFrame( CSV.File( HTTP.get( "https://www.dati.friuliveneziagiulia.it/resource/$code.csv" ).body ) ) for code in codes ]
    for (i, (df, param)) in enumerate(zip(data, params))
        rename!( df, val_types[i] => "value" )
        insertcols!(df, :type => val_types[i] )
        !in( "parametro", names(df) ) && insertcols!( df, :parametro => param )
    end
    dataframe = reduce( (x, y) -> vcat( x, y, cols=:intersect ), data )
    return dataframe
end

"""
    getData(; <keyword arguments> )

Obtain data of category `type` and source of category `kind`.

# Arguments
- `type::Symbol=:METEO`: defines the type of data to be downloaded, it must either be `:METEO` or `:AIRQUALITY`.
- `kind::Symbol=:STATIONS`: defines if the data to be downloaded has to regard information on the stations or their actual measurements, it must either be `:STATIONS` or `:SENSORS`.
"""
function getData(; type::Symbol=:METEO, kind::Symbol=:STATIONS )
    if type == :METEO
        if kind == :STATIONS
            return getMeteoStationsData()
        elseif kind == :SENSORS
            return getMeteoData()
        else
            throw(DomainError(kind, "`kind` must either be `:STATIONS` or `:SENSORS`."))
        end
    elseif type == :AIRQUALITY
        return getAQData()
    else
        throw(DomainError(type, "`type` must either be `:METEO` or `:AIRQUALITY`."))
    end
end

#   ressta = getData( type=:METEO, source=:STATIONS )
#   ressen = getData( type=:METEO, source=:SENSORS )
#   ressen = getData( type=:AIRQUALITY )

end # module

max_stars_repo_name: DavideFavaro/Tirocinio

max_issues_repo_name: DavideFavaro/Tirocinio

max_forks_repo_name: DavideFavaro/Tirocinio


text:
"""
`densityplot(x, y; nargs...)` → `Plot`

Description
============

Draws a density plot for the given points.
It uses the first parameter `x`
(which should be a vector or range) to denote
the horizontal position of each point,
and the second parameter `y`
(which should also be a vector or range)
as their vertical position.
This means that the two vectors have to have the same length.

Usage
======

    densityplot(x, y; title = "", name = "", width = 40, height = 15, border = :solid, xlim = [0, 0], ylim = [0, 0], margin = 3, padding = 1, color = :white, labels = true)

Arguments
==========

- **`x`** : The horizontal dimension for each point.

- **`y`** : The vertical dimension for each point.

- **`title`** : Text to display on the top of the plot.

- **`name`** : Annotation of the current drawing to displayed on the right

- **`width`** : Number of characters per row that should be used for plotting.

- **`height`** : Number of rows that should be used for plotting. Not applicable to `barplot`.

- **`border`** : The style of the bounding box of the plot. Supports `:solid`, `:bold`, `:dashed`, `:dotted`, `:ascii`, and `:none`.

- **`xlim`** : Plotting range for the x coordinate. `[0, 0]` stands for automatic.

- **`ylim`** : Plotting range for the y coordinate. `[0, 0]` stands for automatic.

- **`margin`** : Number of empty characters to the left of the whole plot.

- **`padding`** : Space of the left and right of the plot between the labels and the canvas.

- **`color`** : Color of the drawing. Can be any of `:blue`, `:red`, `:yellow`

- **`labels`** : Can be used to hide the labels by setting `labels=false`.

Returns
========

A plot object of type `Plot{DensityCanvas}`

Author(s)
==========

- Christof Stocker (Github: https://github.com/Evizero)

Examples
=========

    julia> densityplot(randn(1000), randn(1000), title = "Density Plot")

                         Density Plot
           ┌────────────────────────────────────────┐
       2.9 │                    ░                   │
           │                ░    ░  ░               │
           │                ░░░  ░ ░ ░      ░       │
           │             ░░░ ░░▒░▒░░░ ▒ ░           │
           │          ░░ ░░ ░▒░░▒▒▒▒░░░░░░░         │
           │         ░░ ░▒░░░░░▓▒▒░▒░▒░▓░▒░░░       │
           │        ░ ░░░░▓▓▒▓▒▒█▒▓▒▒▓▒▒░▒ ░ ░      │
           │           ░ ▒▒▒▒▓▓▒▓▓▓▓▒█▒▒░░░░░░ ░    │
           │     ░   ░░░░░░░▒▒▓░░▓▒█▒▒▓▓▒▒ ░    ░   │
           │           ░ ░░▒▒░░▒░░▒░░░░ ░░  ░ ░     │
           │         ░   ░ ░ ░░░░░░ ░░ ░░ ░         │
           │              ░  ░ ░ ░  ░ ░   ░░░       │
           │                ░░  ░                   │
           │                                        │
      -3.3 │                                        │
           └────────────────────────────────────────┘
           -3.4                                   2.9

See also
=========

`Plot`, `scatterplot`, `DensityCanvas`
"""
function densityplot(
        x::AbstractVector,
        y::AbstractVector;
        color::Symbol = :white,
        kw...)
    X = convert(Vector{Float64}, x)
    Y = convert(Vector{Float64}, y)
    new_plot = Plot(X, Y, DensityCanvas; grid = false, kw...)
    points!(new_plot, X, Y, color)
end

function densityplot!(
        plot::Plot{<:Canvas},
        x::AbstractVector,
        y::AbstractVector;
        color::Symbol = :white,
        kw...)
    X = convert(Vector{Float64}, x)
    Y = convert(Vector{Float64}, y)
    points!(plot, X, Y, color)
end


text:
export iauCpv
"""
Copy a position/velocity vector.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  vector/matrix support function.

Given:
   pv     double[2][3]    position/velocity vector to be copied

Returned:
   c      double[2][3]    copy

Called:
   iauCp        copy p-vector

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauCpv(double pv[2][3], double c[2][3])

function iauCpv(pv::AbstractMatrix{<:Real})

   # Allocate return value
   c = zeros(Float64, 3, 2)

   ccall((:iauCpv, libsofa_c), Cvoid,
         (Ptr{Cdouble}, Ptr{Cdouble}),
         convert(Array{Float64, 2}, pv'), c)

   return SMatrix{2,3}(c')
end


text:
using MPI, SeisNoise, CUDA, Serialization, Dates

function corr2serial(F::FFTData,arr::Array{FFTData,1},maxlag::Real,CORRDIR::String)
    N = length(arr)
    for ii = 1:N
        C = correlate(F,arr[ii],maxlag)
        stack!(C,allstack=true)
        outpath = joinpath(CORRDIR,F.name,C.name)
        serialize(outpath,cpu(C))
    end
    return nothing
end

function all2all(arr1::Array{FFTData,1},arr2::Array{FFTData,1},maxlag::Real,CORRDIR::String)
    N = length(arr1)
    for ii = 1:N
        corr2serial(arr1[ii],arr2,maxlag,CORRDIR)
    end
    return nothing
end

function one2all(arr::Array{FFTData,1},maxlag::Real,CORRDIR::String)
    N = length(arr)
    for ii = 1:N-1
        corr2serial(arr[ii],arr[ii+1:end],maxlag,CORRDIR)
    end
    return nothing
end

function freemem!(arr::Array{FFTData,1})
	N = length(arr)
	for ii = 1:N
		CUDA.unsafe_free!(arr[ii].fft)
	end
	return nothing
end

function loadFFT(files::Array{String,1})
    N = length(files)
    arr = Array{FFTData}(undef,N)
    for ii = 1:N
        arr[ii] = deserialize(files[ii])
    end
    return arr
end

function splitFFT(files::AbstractArray,Nper::Int)
    Nfiles = length(files)
    splits = ceil(Int,Nfiles/Nper)
    out = []
    for ii = 1:splits
        startind = (ii-1)*Nper + 1
        endind = min(startind + Nper - 1,Nfiles)
        push!(out,files[startind:endind])
    end
    return out
end

function splitCC(splits::AbstractArray,maxlag::Real,CORRDIR::String)

	# load first split on the GPU
	arr1 = loadFFT(splits[1]) .|> gpu
	N = length(splits)
    for ii = 1:N
        # xcorr all in first array
		println("Correlating split $ii $(now())")
		data = Channel{AbstractArray}(1)

		@sync begin
			@async one2all(arr1,maxlag,CORRDIR)
			@async put!(data,loadFFT(splits[N]))
		end

        # load other arrays
		arr2 = take!(data) .|> gpu
        for jj = N:-1:ii+1
			println("Correlating split $jj $(now())")

			data = Channel{AbstractArray}(1)
			# do compute & I/O at same time
			@sync begin
				@async all2all(arr1,arr2,maxlag,CORRDIR)
				@async begin
					loadind = jj > ii + 1 ? jj - 1 : jj
					put!(data,loadFFT(splits[loadind]))
				end
			end

			# transfer CPU arr to gpu
			freemem!(arr2)
			arr2 = take!(data) .|> gpu
        end

		# transfer last arr2 to arr1 on GPU
		freemem!(arr1)
		arr1 = deepcopy(arr2)
		freemem!(arr2)
    end
    return nothing
end


text:
"""
Created on Fri Jun 19 2020

@author: Yoann Pradat

Test clustering functions.
"""

using Distances
using Test
using VariantsNMF

@testset "build_conn" begin
    labels = [1, 0, 1, 2, 1]
    conn_true = [
        [1 0 1 0 1]
        [0 1 0 0 0]
        [1 0 1 0 1]
        [0 0 0 1 0]
        [1 0 1 0 1]
       ]

    conn = VariantsNMF.build_conn(labels)
    conn = convert(Matrix{Int64}, conn)
    @test conn == conn_true
end

@testset "kmeans_init" begin
    p_centroids = 96
    n_centroids = 16
    n_repeats = 400

    # 1. test on synthetic data where centroids in each repetition are identical
    # each cluster is then reduced to one point and the clustering is perfect
    centroids_true = rand(p_centroids, n_centroids)
    labels_true = repeat(collect(1:n_centroids), n_repeats)
    X = repeat(centroids_true, 1, n_repeats)

    # 1.1 stop when relative change in cost remains below 5e-2 over 10 iter
    centroids_ref, labels_ref, cost_ref = VariantsNMF.kmeans_init(
        X             = X,
        init          = collect(1:n_centroids),
        max_iter      = 100,
        dist_func     = CosineDist(),
        stopping_crit = :cost,
        stopping_tol  = 5e-2,
        stopping_iter = 10,
        verbose       = false
    )

    @test centroids_ref ≈ centroids_true atol = 1e-9
    @test cost_ref ≈ 0 atol = 1e-9
    @test labels_ref == labels_true

    # 1.2 stop when cosine distance between new centroids and ref centroids remains below 5e-2 for 10 iter
    centroids_ref, labels_ref, cost_ref = VariantsNMF.kmeans_init(
        X             = X,
        init          = collect(1:n_centroids),
        max_iter      = 100,
        dist_func     = CosineDist(),
        stopping_crit = :dist,
        stopping_tol  = 5e-2,
        stopping_iter = 10,
        verbose       = false
    )

    @test centroids_ref ≈ centroids_true atol = 1e-9
    @test cost_ref ≈ 0 atol = 1e-9
    @test labels_ref == labels_true

end


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

# ToDo: This is not a good test...

md = fmi2CreateModelDescription()
fmi2ModelDescriptionAddRealStateAndDerivative(md, "mass.s")
fmi2ModelDescriptionAddRealStateAndDerivative(md, "mass.v")
fmi2ModelDescriptionAddRealOutput(md, "mass.f")


text:
# region is an iterable subset of dimensions
# spec. an integer, range, tuple, or array

# inplace complex
function plan_fft!(X::CuArray{T,N}, region) where {T<:cufftComplexes,N}
    K = CUFFT_FORWARD
    inplace = true
    xtype = (T == cufftComplex) ? CUFFT_C2C : CUFFT_Z2Z

    pp = _mkplan(xtype, size(X), region)

    cCuFFTPlan{T,K,inplace,N}(pp, X, size(X), region, xtype)
end

function plan_bfft!(X::CuArray{T,N}, region) where {T<:cufftComplexes,N}
    K = CUFFT_INVERSE
    inplace = true
    xtype =  (T == cufftComplex) ? CUFFT_C2C : CUFFT_Z2Z

    pp = _mkplan(xtype, size(X), region)

    cCuFFTPlan{T,K,inplace,N}(pp, X, size(X), region, xtype)
end

# out-of-place complex
function plan_fft(X::CuArray{T,N}, region) where {T<:cufftComplexes,N}
    K = CUFFT_FORWARD
    xtype =  (T == cufftComplex) ? CUFFT_C2C : CUFFT_Z2Z
    inplace = false

    pp = _mkplan(xtype, size(X), region)

    cCuFFTPlan{T,K,inplace,N}(pp, X, size(X), region, xtype)
end

function plan_bfft(X::CuArray{T,N}, region) where {T<:cufftComplexes,N}
    K = CUFFT_INVERSE
    inplace = false
    xtype =  (T == cufftComplex) ? CUFFT_C2C : CUFFT_Z2Z

    pp = _mkplan(xtype, size(X), region)

    cCuFFTPlan{T,K,inplace,N}(pp, X, size(X), region, xtype)
end

# out-of-place real-to-complex
function plan_rfft(X::CuArray{T,N}, region) where {T<:cufftReals,N}
    K = CUFFT_FORWARD
    inplace = false
    xtype =  (T == cufftReal) ? CUFFT_R2C : CUFFT_D2Z

    pp = _mkplan(xtype, size(X), region)

    ydims = collect(size(X))
    ydims[region[1]] = div(ydims[region[1]],2)+1

    rCuFFTPlan{T,K,inplace,N}(pp, X, (ydims...,), region, xtype)
end

function plan_brfft(X::CuArray{T,N}, d::Integer, region::Any) where {T<:cufftComplexes,N}
    K = CUFFT_INVERSE
    inplace = false
    xtype =  (T == cufftComplex) ? CUFFT_C2R : CUFFT_Z2D
    ydims = collect(size(X))
    ydims[region[1]] = d

    pp = _mkplan(xtype, (ydims...,), region)

    rCuFFTPlan{T,K,inplace,N}(pp, X, (ydims...,), region, xtype)
end

# FIXME: plan_inv methods allocate needlessly (to provide type parameters)
# Perhaps use FakeArray types to avoid this.

function plan_inv(p::cCuFFTPlan{T,CUFFT_FORWARD,inplace,N}) where {T,N,inplace}
    X = CuArray{T}(undef, p.sz)
    pp = _mkplan(p.xtype, p.sz, p.region)
    ScaledPlan(cCuFFTPlan{T,CUFFT_INVERSE,inplace,N}(pp, X, p.sz, p.region,
                                                     p.xtype),
               normalization(X, p.region))
end

function plan_inv(p::cCuFFTPlan{T,CUFFT_INVERSE,inplace,N}) where {T,N,inplace}
    X = CuArray{T}(undef, p.sz)
    pp = _mkplan(p.xtype, p.sz, p.region)
    ScaledPlan(cCuFFTPlan{T,CUFFT_FORWARD,inplace,N}(pp, X, p.sz, p.region,
                                                     p.xtype),
               normalization(X, p.region))
end

function plan_inv(p::rCuFFTPlan{T,CUFFT_INVERSE,inplace,N}
                  ) where {T<:cufftComplexes,N,inplace}
    X = CuArray{real(T)}(undef, p.osz)
    Y = CuArray{T}(undef, p.sz)
    xtype = p.xtype == CUFFT_C2R ? CUFFT_R2C : CUFFT_D2Z
    pp = _mkplan(xtype, p.osz, p.region)
    ScaledPlan(rCuFFTPlan{real(T),CUFFT_FORWARD,inplace,N}(pp, X, p.sz, p.region,
                                                     xtype),
               normalization(X, p.region))
end

function plan_inv(p::rCuFFTPlan{T,CUFFT_FORWARD,inplace,N}
                  ) where {T<:cufftReals,N,inplace}
    X = CuArray{complex(T)}(undef, p.osz)
    Y = CuArray{T}(undef, p.sz)
    xtype = p.xtype == CUFFT_R2C ? CUFFT_C2R : CUFFT_Z2D
    pp = _mkplan(xtype, p.sz, p.region)
    ScaledPlan(rCuFFTPlan{complex(T),CUFFT_INVERSE,inplace,N}(pp, X, p.sz,
                                                              p.region, xtype),
               normalization(Y, p.region))
end

# The rest of the standard API

size(p::CuFFTPlan) = p.sz

function mul!(y::CuArray{Ty}, p::CuFFTPlan{T,K,false}, x::CuArray{T}
                  ) where {Ty,T,K}
    assert_applicable(p,x,y)
    unsafe_execute!(p,x,y)
    return y
end

function *(p::cCuFFTPlan{T,K,true,N}, x::CuArray{T,N}) where {T,K,N}
    assert_applicable(p,x)
    unsafe_execute!(p,x)
    x
end

function *(p::rCuFFTPlan{T,CUFFT_FORWARD,false,N}, x::CuArray{T,N}
           ) where {T<:cufftReals,N}
    @assert p.xtype ∈ [CUFFT_R2C,CUFFT_D2Z]
    y = CuArray{complex(T),N}(undef, p.osz)
    mul!(y,p,x)
    y
end

function *(p::rCuFFTPlan{T,CUFFT_INVERSE,false,N}, x::CuArray{T,N}
           ) where {T<:cufftComplexes,N}
    @assert p.xtype ∈ [CUFFT_C2R,CUFFT_Z2D]
    y = CuArray{real(T),N}(undef, p.osz)
    mul!(y,p,x)
    y
end

function *(p::cCuFFTPlan{T,K,false,N}, x::CuArray{T,N}) where {T,K,N}
    y = CuArray{T,N}(undef, p.osz)
    mul!(y,p,x)
    y
end


text:
module BaselineFunctions

	import Statistics

	export calculateBaseline

	function calculateBaseline(massAxis, avgSpectrum; baselinePointWidth = 0.3, threshold=0.5)

	  baselinePoints = ceil(massAxis[1])+0.6:baselinePointWidth:floor(massAxis[end]-1)-0.4
	  baselineValues = similar(baselinePoints)
	  baselineNoise = similar(baselinePoints)

	  for i=1:length(baselinePoints)
	    startIdx = searchsortedfirst(massAxis,baselinePoints[i]-baselinePointWidth)
	    endIdx = searchsortedfirst(massAxis,baselinePoints[i]+baselinePointWidth)
	    subSet = view(avgSpectrum, startIdx:endIdx)
	    threshold = Statistics.quantile(subSet,0.2)
	    baselineSamples = Array{Float64,1}()
	    #fill!(baselineSamples,0.0)
	    for j=1:length(subSet)
	      if subSet[j] <= threshold
		push!(baselineSamples,subSet[j])
	      end
	    end
	    baselineValues[i] = Statistics.mean(baselineSamples)
	    if length(baselineSamples) > 3
	      baselineNoise[i] = Statistics.std(baselineSamples)
	    else
	      baselineNoise[i] = baselineValues[i]
	    end
	    if (baselineValues[i] > 10000 || baselineValues == 0)
	      println("Strange Baseline at $(baselinePoints[i]):\n$baselineSamples")
	    end

	  end
	  return baselinePoints,baselineValues, baselineNoise
	end
end


text:
#creates a logarithmic range
function log_range(start, stop, points)
    r = LinRange(log(start), log(stop), points)
    r = exp.(r)
    return r
end

export log_range


text:
function dijkstra(net::Network{T}, snk::U) where {T<:Integer, U<:Integer}
    nnodes = numnodes(net)

    costs = fill(typemax(U), nnodes)
    childvec = zeros(T, nnodes)

    P = PriorityQueue{T,U}()
    P[snk] = costs[snk] = 0
    childvec[snk] = snk

    while !isempty(P)
        v = dequeue!(P)
        for u in inneighbors(net, v)
            lidx = idx(net, u, v)
            alt = costs[v] + net.links[lidx].l
            if (costs[u] > alt)
                P[u] = costs[u] = alt
                childvec[u] = v
            end
        end
    end
    childvec
end

function initchoices()
    # choice initialization
    DTC = zeros(T, size(trips)...);
    DTC[1,..] .= 1.
    #DTC[1:Tm,..] .= rand(size(DTC[1:Tm,..])...)
    #DTC ./= sum(DTC, dims=1)

    pathvecs = Dict(k => dijkstra(net, k) for k in snks);
    SR = Dict(d => Dict(i => zeros(T, nsinks, nclasses) for i in 1:2) for d in divs);
    for d in divs
        for (kno,k) in enumerate(snks)
            for c in 1:nclasses
                if (outneighbors(net, d)[1] == pathvecs[k][d])
                    SR[d][1][:,kno,c] .= 1.
                elseif (outneighbors(net, d)[2] == pathvecs[k][d])
                    SR[d][2][:,kno,c] .= 1.
                else
                    error()
                end
            end
        end
    end
    return (DTC, SR)
end


text:
using Documenter, VoronoiFVM, Literate

#
# Replace SOURCE_URL marker with github url of source
#
function replace_source_url(input,source_url)
    lines_in = collect(eachline(IOBuffer(input)))
    lines_out=IOBuffer()
    for line in lines_in
        println(lines_out,replace(line,"SOURCE_URL" => source_url))
    end
    return String(take!(lines_out))
end

function make_all()
    #
    # Generate Markdown pages from examples
    #
    example_jl_dir = joinpath(@__DIR__,"..","examples")
    example_md_dir  = joinpath(@__DIR__,"src","examples")

    for example_source in readdir(example_jl_dir)
        base,ext=splitext(example_source)
        if ext==".jl"
            source_url="https://github.com/j-fu/VoronoiFVM.jl/raw/master/examples/"*example_source
            preprocess(buffer)=replace_source_url(buffer,source_url)
            Literate.markdown(joinpath(@__DIR__,"..","examples",example_source),
                              example_md_dir,
                              documenter=false,
                              info=false,
                              preprocess=preprocess)
        end
    end

    generated_examples=vcat(["runexamples.md"],joinpath.("examples",readdir(example_md_dir)))

    makedocs(
        sitename="VoronoiFVM.jl",
        modules = [VoronoiFVM],
        clean = false,
        doctest = true,
        authors = "J. Fuhrmann",
        repo="https://github.com/j-fu/VoronoiFVM.jl",
        pages=[
            "Home"=>"index.md",
            "changes.md",
            "method.md",
            "API Documentation" => [
                "physics.md",
                "system.md",
                "solutions.md",
                "solver.md",
                "post.md",
                "misc.md",
                "allindex.md",
            ],
            "Examples" => generated_examples
        ]
    )

    rm(example_md_dir,recursive=true)

    if !isinteractive()
        deploydocs(repo = "github.com/j-fu/VoronoiFVM.jl.git")
    end
end

make_all()


text:
lattice = load_lattice(joinpath(@__DIR__, "./Mn3Si.in"))
atoms   = load_atoms(  joinpath(@__DIR__, "./Mn3Si.in"))
atoms   = attach_pseudos(atoms, Mn="hgh/pbe/mn-q15.hgh",
                                Si="hgh/pbe/si-q4.hgh")

# AFM ordering, configuration 5(c) in Jiang and Yang J. Alloys Compounds (2021)
@assert atoms[1][1].symbol == :Mn
@assert atoms[2][1].symbol == :Si
magnetic_moments = [atoms[1][1] => [5.0, -5.0, 5.0], atoms[2][1] => [0.0]]

cases = standard_cases()

config[@__DIR__] = (
    directory        = @__DIR__,
    lattice          = lattice,
    atoms            = atoms,
    magnetic_moments = magnetic_moments,
    #
    kwargs_model     = (smearing=Smearing.Gaussian(), temperature=0.01),
    Ecut             = 45,
    kgrid            = [13, 13, 13],
    supersampling    = 2.0,
    potential        = nothing,
    kwargs_scf       = (; ),
    mixing           = KerkerMixing(),
    cases            = cases,
)

max_stars_repo_name: mfherbst/supporting-adaptive-damping

max_issues_repo_name: mfherbst/supporting-adaptive-damping

max_forks_repo_name: mfherbst/supporting-adaptive-damping


text:
# comparison functions for HTML Nodes and Documents

# TODO right now hashing and equality completely ignore
# parents. I think this is *probably* appropriate but it deserves
# some more thought. There's an argument that two HTMLElements with
# the same contents and children but different parent pointers are not
# really equal. Perhaps an auxilliary equality function could be provided
# for this purpose?

# equality

import Base: ==, isequal, hash

isequal(x::HTMLDocument, y::HTMLDocument) =
    isequal(x.doctype,y.doctype) && isequal(x.root,y.root)

isequal(x::HTMLText,y::HTMLText) = isequal(x.text, y.text)

isequal(x::HTMLElement, y::HTMLElement) =
    isequal(x.attributes,y.attributes) && isequal(x.children,y.children)

==(x::HTMLDocument, y::HTMLDocument) =
    ==(x.doctype,y.doctype) && ==(x.root,y.root)

==(x::HTMLText,y::HTMLText) = ==(x.text, y.text)

==(x::HTMLElement, y::HTMLElement) =
    ==(x.attributes,y.attributes) && ==(x.children,y.children)

# hashing

function hash(doc::HTMLDocument)
    hash(hash(HTMLDocument),hash(hash(doc.doctype), hash(doc.root)))
end

function hash(elem::HTMLElement{T}) where {T}
    h = hash(HTMLElement)
    h = hash(h,hash(T))
    h = hash(h,hash(attrs(elem)))
    for child in children(elem)
        h = hash(h,hash(child))
    end
    return h
end

hash(t::HTMLText) = hash(hash(HTMLText),hash(t.text))


text:
using CellListMap
using StaticArrays
using NearestNeighbors
using Test

function nl_NN(x,y,r)
    balltree = BallTree(x)
    return inrange(balltree,y,r,true)
end

function compare_nb_lists(list_CL,list_NN;self=false)
    if self
        for (i,list_original) in pairs(list_NN)
            list = filter(!isequal(i),list_original)
            cl = filter(tup -> (tup[1] == i || tup[2] == i), list_CL)
            if length(cl) != length(list)
                @show i
                @show length(list), list
                @show length(cl), cl
                return false
            end
            for j in list
                length(findall(tup -> (tup[1] == j || tup[2] == j), cl)) == 1 || return false
            end
        end
    else
        for (i,list) in pairs(list_NN)
            cl = filter(tup -> tup[2] == i, list_CL)
            if length(cl) != length(list)
                @show i
                @show length(list), list
                @show length(cl), cl
                return false
            end
            for j in list
                length(findall(tup -> tup[1] == j, cl)) == 1 || return false
            end
        end
    end
    return true
end

@testset "neighborlist" begin

    r = 0.1

    for N in [2,3]

        #
        # Using vectors as input
        #

        # With y smaller than x
        x = [ rand(SVector{N,Float64}) for _ in 1:1000 ]
        y = [ rand(SVector{N,Float64}) for _ in 1:500 ]

        nb = nl_NN(x,x,r)
        cl = CellListMap.neighborlist(x,r)
        @test compare_nb_lists(cl,nb,self=true)

        nb = nl_NN(x,y,r)
        cl = CellListMap.neighborlist(x,y,r,autoswap=false)
        @test compare_nb_lists(cl,nb,self=false)
        cl = CellListMap.neighborlist(x,y,r,autoswap=true)
        @test compare_nb_lists(cl,nb,self=false)

        # with x smaller than y
        x = [ rand(SVector{N,Float64}) for _ in 1:500 ]
        y = [ rand(SVector{N,Float64}) for _ in 1:1000 ]
        nb = nl_NN(x,y,r)
        cl = CellListMap.neighborlist(x,y,r,autoswap=false)
        @test compare_nb_lists(cl,nb,self=false)
        cl = CellListMap.neighborlist(x,y,r,autoswap=true)
        @test compare_nb_lists(cl,nb,self=false)

        # Using matrices as input
        x = rand(N,1000)
        y = rand(N,500)

        nb = nl_NN(x,x,r)
        cl = CellListMap.neighborlist(x,r)
        @test compare_nb_lists(cl,nb,self=true)

        nb = nl_NN(x,y,r)
        cl = CellListMap.neighborlist(x,y,r,autoswap=false)
        @test compare_nb_lists(cl,nb,self=false)
        cl = CellListMap.neighborlist(x,y,r,autoswap=true)
        @test compare_nb_lists(cl,nb,self=false)

        # with x smaller than y
        x = rand(N,500)
        y = rand(N,1000)
        nb = nl_NN(x,y,r)
        cl = CellListMap.neighborlist(x,y,r,autoswap=false)
        @test compare_nb_lists(cl,nb,self=false)
        cl = CellListMap.neighborlist(x,y,r,autoswap=true)
        @test compare_nb_lists(cl,nb,self=false)

        # Check random coordinates to test the limits more thoroughly
        check_random_NN = true
        for i in 1:500
            x = rand(SVector{N,Float64},100); y = rand(SVector{N,Float64},50);
            nb = nl_NN(x,y,r); cl = CellListMap.neighborlist(x,y,r,autoswap=false);
            check_random_NN = compare_nb_lists(cl,nb,self=false)
        end
        @test check_random_NN

        # with different types
        x = rand(Float32,N,500)
        y = rand(Float32,N,1000)
        nb = nl_NN(x,y,r)
        cl = CellListMap.neighborlist(x,y,r,autoswap=false)
        @test compare_nb_lists(cl,nb,self=false)
        cl = CellListMap.neighborlist(x,y,r,autoswap=true)
        @test compare_nb_lists(cl,nb,self=false)

    end

end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule ASL_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("ASL")
JLLWrappers.@generate_main_file("ASL", UUID("ae81ac8f-d209-56e5-92de-9978fef736f9"))
end  # module ASL_jll


text:
# ---
# title: 581. Shortest Unsorted Continuous Subarray
# id: problem581
# author: Indigo
# date: 2021-06-26
# difficulty: Medium
# categories: Array
# link: <https://leetcode.com/problems/shortest-unsorted-continuous-subarray/description/>
# hidden: true
# ---
#
# Given an integer array `nums`, you need to find one **continuous subarray**
# that if you only sort this subarray in ascending order, then the whole array
# will be sorted in ascending order.
#
# Return  _the shortest such subarray and output its length_.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [2,6,4,8,10,9,15]
#     Output: 5
#     Explanation: You need to sort [6, 4, 8, 10, 9] in ascending order to make the whole array sorted in ascending order.
#
#
# **Example 2:**
#
#
#
#     Input: nums = [1,2,3,4]
#     Output: 0
#
#
# **Example 3:**
#
#
#
#     Input: nums = [1]
#     Output: 0
#
#
#
#
# **Constraints:**
#
#   * `1 <= nums.length <= 104`
#   * `-105 <= nums[i] <= 105`
#
#
## @lc code=start
using LeetCode

function find_unsorted_subarray(nums::Vector{Int})
    stk = Int[]
    l, r = length(nums) + 1, 1
    for i in eachindex(nums)
        while !isempty(stk) && nums[stk[end]] > nums[i]
            l = min(l, pop!(stk))
        end
        push!(stk, i)
    end
    empty!(stk)
    for i in length(nums):-1:1
        while !isempty(stk) && nums[stk[end]] < nums[i]
            r = max(r, pop!(stk))
        end
        push!(stk, i)
    end
    return r - l > 0 ? r - l + 1 : 0
end
## @lc code=end


text:
struct Bundle{S<:System,P,F}
    produce::Produce{P}
    recursive::Bool
    filter::F
end

Bundle(p::Produce{P,V}, ops::AbstractString) where {P,V} = begin
    recursive = false
    filter = nothing
    index = 0
    for op in split(ops, "/")
        if op == "*"
            # collecting children only at the current level
            recursive = false
        elseif op == "**"
            # collecting all children recursively
            recursive = true
        else
            #TODO: support generic indexing function?
            filter = op
        end
    end
    S = eltype(p)
    F = typeof(filter)
    Bundle{S,P,F}(p, recursive, filter)
end

Base.collect(b::Bundle{S}) where {S<:System} = begin
    p = getfield(b, :produce)
    v = collect(p)
    if getfield(b, :recursive)
        l = S[]
        #TODO: possibly reduce overhead by reusing calculated values in child nodes
        g(V::Vector{<:System}) = for s in V; g(s) end
        g(s::System) = (push!(l, s); g(value(getfield(s, p.name))))
        g(::Nothing) = nothing
        g(v)
    else
        l = copy(v)
    end
    f = getfield(b, :filter)
    if !isnothing(f)
        filter!(s -> value(s[f]), l)
    end
    l
end

Base.getindex(s::Produce, ops::AbstractString) = Bundle(s, ops)

struct Bunch{V}
    it::Base.Generator
end

Base.iterate(b::Bunch, i...) = iterate(getfield(b, :it), i...)
Base.length(b::Bunch) = length(getfield(b, :it))
Base.eltype(::Type{<:Bunch{<:State{V}}}) where V = V

Base.getproperty(b::Bundle{S}, p::Symbol) where {S<:System} = (value(getfield(x, p)) for x in value(b)) |> Bunch{vartype(S, p)}
Base.getproperty(b::Bunch{S}, p::Symbol) where {S<:System} = (value(getfield(x, p)) for x in b) |> Bunch{vartype(S, p)}
Base.getindex(b::Bundle, i::AbstractString) = getproperty(b, Symbol(i))
Base.getindex(b::Bunch, i::AbstractString) = getproperty(b, Symbol(i))

value(b::Bundle) = collect(b)
#TODO: also make final value() based on generator, but then need sum(x; init=0) in Julia 1.6 for empty generator
#value(b::Bunch) = (value(v) for v in b)
value(b::Bunch) = collect(b)

Base.getindex(b::Bundle) = value(b)
Base.getindex(b::Bunch) = value(b)

Base.adjoint(b::Bundle) = value(b)
Base.adjoint(b::Bunch) = value(b)


text:
"""
    nonlinear_gradient_correction(x)
    nonlinear_gradient_correction!(x, val)

Returns `Bool` stating if the image saved has been corrected for gradient
nonlinearities by the scanner sequence. Default is `false`.
"""
@defprop NonlinearGradientCorrection{:nonlinear_gradient_correction}::Bool

"""
    pulse_sequence(x)
    pulse_sequence!(x, val)

General description of the pulse sequence used for the scan (i.e. MPRAGE,
Gradient Echo EPI, Spin Echo EPI, Multiband gradient echo EPI).
"""
@defprop PulseSequence{:pulse_sequence}::String

"""
    pulse_sequence_details(x)
    pulse_sequence_details!(x, val)

Information beyond pulse sequence type that identifies the specific pulse
sequence used (i.e. "Standard Siemens Sequence distributed with the VB17
software," "Siemens WIP ### version #.##," or "Sequence written by X using a
version compiled on MM/DD/YYYY").
"""
@defprop PulseSequenceDetails{:pulse_sequence_details}::String

"""
    pulse_seqtype(x)
    pulse_seqtype!(x, val)

A general description of the pulse sequence used for the scan (i.e. MPRAGE,
Gradient Echo EPI, Spin Echo EPI, Multiband gradient echo EPI).
"""
@defprop PulseSequenceType{:pulse_seqtype}::String

"""
    scanning_sequence(x)
    scanning_sequence!(x, val)

Description of the type of sequence data acquired.
"""
@defprop ScanningSequence{:scanning_sequence}::String

"""
    sequence_name(x)
    sequence_name!(x, val)

Manufacturer’s designation of the sequence name.
"""
@defprop SequenceName{:sequence_name}::String

"""
    sequence_variant(x)
    sequence_variant!(x, val)

Variant of the `scanning_sequence` property.
"""
@defprop SequenceVariant{:sequence_variant}::String
"""
    SequenceMetadata

Metadata structure for general MRI sequence information.

## Supported Properties
$(description_list(nonlinear_gradient_correction, pulse_sequence, pulse_sequence_details,
                   pulse_sequence_type, scanning_sequence, sequence_name, sequence_variant))
"""
struct SequenceMetadata
    nonlinear_gradient_correction::Bool
    pulse_sequence::String
    pulse_sequence_details::String
    pulse_sequence_type::String
    scanning_sequence::String
    sequence_name::String
    sequence_variant::String
end

@properties SequenceMetadata begin
    nonlinear_gradient_correction(self) => :nonlinear_gradient_correction
    pulse_sequence(self) => :pulse_sequence
    pulse_sequence_details(self) => :pulse_sequence_details
    pulse_seqtype(self) => :pulse_sequence_type
    scanning_sequence(self) => :scanning_sequence
    sequence_name(self) => :sequence_name
    sequence_variant(self) => :sequence_variant
end


text:
struct CartesianFDMContext{N,P,R,T,D,S,M}
    top::P
    n::NTuple{N,Int}
    ρ::R
    τ::NTuple{2,NTuple{N,T}}
    δ::NTuple{2,NTuple{N,D}}
    σ::NTuple{2,NTuple{N,S}}
    μ::Tuple{M,NTuple{N,M}}
end

function cartesianfdmcontext(top, n)
    ρ = spdiagm(0 => ones(Bool, prod(n)))
    τ = backwardshiftmatrices(top, n), forwardshiftmatrices(top, n)
    δ = Ref(ρ) .- τ[1], τ[2] .- Ref(ρ)
    σ = Ref(ρ) .+ τ[1], τ[2] .+ Ref(ρ)
    μ = maskmatrices(top, n)
    CartesianFDMContext(top, n, ρ, τ, δ, σ, μ)
end

"""

"""
_kron(opn) = opn

_kron(opn...) = kron(reverse(opn)...)

_kron(opn::NTuple{1}, eye::NTuple{1}) = opn

_kron(opn::NTuple{2}, eye::NTuple{2}) =
    kron(eye[2], opn[1]),
    kron(opn[2], eye[1])

_kron(opn::NTuple{3}, eye::NTuple{3}) =
    kron(eye[3], eye[2], opn[1]),
    kron(eye[3], opn[2], eye[1]),
    kron(opn[3], eye[2], eye[1])

"""

"""
function forwardshiftmatrices(top, n)
    opn = _forwardshift.(top, n)
    eye = I.(n)
    _kron(opn, eye)
end

_forwardshift(::Topology, n::Int) =
    spdiagm(1 => ones(Bool, n-1))

_forwardshift(::Periodic, n::Int) =
    spdiagm(1-n => ones(Bool, 1),
            1 => ones(Bool, n-1))

"""

"""
function backwardshiftmatrices(top, n)
    opn = _backwardshift.(top, n)
    eye = I.(n)
    _kron(opn, eye)
end

_backwardshift(::Topology, n::Int) =
    spdiagm(-1 => ones(Bool, n-1))

_backwardshift(::Periodic, n::Int) =
    spdiagm(-1 => ones(Bool, n-1),
            n-1 => ones(Bool, 1))

"""

"""
function maskmatrices(top, n)
    opn = _mask.(top, n)
    μ⁰ = _kron(opn...)

    opn = _normal.(top, n)
    eye = _tangent.(top, n)
    μ¹ = _kron(opn, eye)

    μ⁰, μ¹
end

_mask(::Topology, n::Int) =
    spdiagm(0 => [i ≠ n for i in 1:n])

_mask(::Periodic, n::Int) =
    spdiagm(0 => ones(Bool, n))

_normal(::Topology, n::Int) =
    spdiagm(0 => [i ≠ 1 && i ≠ n for i in 1:n])

_normal(::Periodic, n::Int) =
    spdiagm(0 => ones(Bool, n))

_tangent(::Topology, n::Int) =
    spdiagm(0 => [i ≠ n for i in 1:n])

_tangent(::Periodic, n::Int) =
    spdiagm(0 => ones(Bool, n))


text:
## utility methods

XXX(msg) = println(msg); XXX() = XXX("define me")

function members(object::PyObject)
    convert(Vector{(String,PyObject)},
            pycall(PyCall.inspect["getmembers"], PyObject, object))
end
has_member(object::PyObject, key) = contains(pluck(members(object),1), key)

pluck(xs, key) = [ x[key] for x = xs ]


text:
f = :(f(x, y = 3; a, b = 5) = 50)

doc_f = :(Core.@doc "I document function f" f)

f2 = :(function f(x)
    println(2)
    return 5
end)

s = :(mutable struct HandleType
    handle::Ptr{Nothing}
end)

doc_s = :(Core.@doc "I document structure s" s)

c = :(const sym = val)

e = :(@enum myenum a b c d)
e2 = :(@cenum myotherenum::Int begin
    a = 1
    b = 2
    c = 3
    d = 4
end)

@testset "Expressions" begin
    include("exprs/utils.jl")
    include("exprs/deconstruct.jl")
    include("exprs/reconstruct.jl")

    @test is_broadcast(:(f.(x, y)))
    @test !is_broadcast(:(f(x, y)))
    @test !is_broadcast(:(x.y))

    @test broadcast_ex(:(f(x, y)), true) == :(f.(x, y))
    @test broadcast_ex(:(f(x, y)), false) == :(f(x, y))
    @test prettify(concat_exs(:(x = 3), :(y = 4))) == prettify(:(x = 3; y = 4))
    @test prettify(concat_exs(:(
        begin
            x = 3
        end
    ), :(y = 4))) == prettify(:(x = 3; y = 4))
end


text:
@testset "PopulationGrouping" begin
  @testset "ExampleGrouping" begin
    num_groups = 100
    num_individuals = 100
    rng = MersenneTwister(13)
    group_ids = rand(rng, 1:num_groups, num_individuals)
    @test minimum(group_ids) > 0
    @test maximum(group_ids) <= num_groups

    grouping = MocosSim.PopulationGrouping(group_ids, num_groups)

    @test MocosSim.numgroups(grouping) == num_groups

    @testset "Every group has the right people" begin
      for group_id in 1:num_groups
        group = MocosSim.getgroup(grouping, group_id)

        for person_id in group
          @test group_ids[person_id] == group_id
        end
      end
    end

    @testset "Every person is in the right group" begin
      for person_id in 1:num_individuals
        group_id = group_ids[person_id]
        group = MocosSim.getgroup(grouping, group_id)
        @test person_id in group
      end
    end

    @testset "Group sizes" begin
      sizes = MocosSim.groupsizes(grouping)
      @test length(sizes) == num_groups
      for i in 1:MocosSim.numgroups(grouping)
        @test length(MocosSim.getgroup(grouping, i)) == sizes[i]
      end
    end
  end
end


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

#ToDo: Move functionality of reshape_variate and reshape_variates to package ValueShapes.

reshape_variate(shape::Missing, v::Any) = v

function reshape_variate(shape::AbstractValueShape, v::Any)
    v_shape = valshape(v)
    if !(v_shape <= shape)
        throw(ArgumentError("Shape of variate doesn't match target variate shape, with variate of type $(typeof(v)) and expected shape $(shape)"))
    end
    v
end

function reshape_variate(shape::ArrayShape{<:Real,1}, v::Any)
    unshaped_v = unshaped(v)::AbstractVector{<:Real}
    reshape_variate(shape, unshaped_v)
end

function reshape_variate(shape::AbstractValueShape, v::AbstractVector{<:Real})
    _reshape_realvec(shape, v)
end

function reshape_variate(shape::ArrayShape{<:Real,1}, v::AbstractVector{<:Real})
    _reshape_realvec(shape, v)
end

function _reshape_realvec(shape::AbstractValueShape, v::AbstractVector{<:Real})
    ndof = length(eachindex(v))
    ndof_expected = totalndof(shape)
    if ndof != ndof_expected
        throw(ArgumentError("Invalid length ($ndof) of variate, target shape  $(shape) has $ndof_expected degrees of freedom"))
    end
    shape(v)
end

reshape_variates(shape::Missing, vs::AbstractVector{<:Any}) = v

function reshape_variates(shape::AbstractValueShape, vs::AbstractVector{<:Any})
    v_elshape = elshape(vs)
    if !(v_elshape <= shape)
        throw(ArgumentError("Shape of variates doesn't match target variate shape, with variates of type $(eltype(vs)) and expected element shape $(shape)"))
    end
    vs
end

function reshape_variates(shape::ArrayShape{<:Real,1}, vs::AbstractVector{<:Any})
    unshaped_vs = unshaped.(vs)::AbstractVector{<:AbstractVector{<:Real}}
    reshape_variates(shape, unshaped_vs)
end

function reshape_variates(shape::AbstractValueShape, vs::AbstractVector{<:AbstractVector{<:Real}})
    _reshape_realvecs(shape, vs)
end

function reshape_variates(shape::ArrayShape{<:Real,1}, vs::AbstractVector{<:AbstractVector{<:Real}})
    _reshape_realvecs(shape, vs)
end

function _reshape_realvecs(shape::AbstractValueShape, vs::AbstractVector{<:AbstractVector{<:Real}})
    ndof = first(innersize(vs))
    ndof_expected = totalndof(shape)
    if ndof != ndof_expected
        throw(ArgumentError("Invalid length ($ndof) of variates, target shape $(shape) has $ndof_expected degrees of freedom"))
    end
    shape.(vs)
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

rst(x) = sprint(rst, x)

function rst(io::IO, content::Vector)
    isempty(content) && return
    for md in content[1:end-1]
        rst(io, md)
        println(io)
    end
    rst(io, content[end])
end

rst(io::IO, md::MD) = rst(io, md.content)

function rst(io::IO, header::Header{l}) where l
    s = rstinline(header.text)
    println(io, s)
    println(io, string("*=-~:.^"[l])^length(s))
    println(io)
end

function rst(io::IO, code::Code)
    if code.language == "jldoctest"
        println(io, ".. doctest::\n")
    elseif code.language != "rst"
        println(io, ".. code-block:: julia\n")
    end
    for l in lines(code.code)
        println(io, "    ", l)
    end
end

function rst(io::IO, p::Paragraph)
    rstinline(io, p.content)
    println(io)
end

function rst(io::IO, list::List)
    for (i, item) in enumerate(list.items)
        bullet = isordered(list) ? "$(i + list.ordered - 1). " : "* "
        print(io, bullet)
        lines = split(rstrip(sprint(rst, item)), '\n')
        for (n, line) in enumerate(lines)
            print(io, (n == 1 || isempty(line)) ? "" : " "^length(bullet), line)
            n < length(lines) && println(io)
        end
        println(io)
    end
end

function rst(io::IO, q::BlockQuote)
    s = sprint(rst, q.content)
    for line in split(rstrip(s), "\n")
        println(io, "    ", line)
    end
    println(io)
end

function rst(io::IO, f::Footnote)
    print(io, ".. [", f.id, "]")
    s = sprint(rst, f.text)
    lines = split(rstrip(s), "\n")
    # Single line footnotes are printed on the same line as their label
    # rather than taking up an additional line.
    if length(lines) == 1
        println(io, " ", lines[1])
    else
        println(io)
        for line in lines
            println(io, isempty(line) ? "" : "   ", rstrip(line))
        end
        println(io)
    end
end

function rst(io::IO, md::Admonition)
    s = sprint(rst, md.content)
    title = md.title == ucfirst(md.category) ? "" : md.title
    println(io, ".. ", md.category, "::", isempty(title) ? "" : " $title")
    for line in split(rstrip(s), "\n")
        println(io, isempty(line) ? "" : "   ", line)
    end
    println(io)
end

function rst(io::IO, md::HorizontalRule)
    println(io, "–" ^ 5)
end

function rst(io::IO, l::LaTeX)
    println(io, ".. math::\n")
    for line in lines(l.formula)
        println(io, "    ", line)
    end
end

rst(io::IO, md) = show(io, "text/rst", md)

# Inline elements

rstinline(x) = sprint(rstinline, x)

function rstinline(io::IO, md...)
    wasCode = false
    for el in md
        wasCode && isa(el, AbstractString) && !Base.startswith(el, " ") && print(io, "\\ ")
        wasCode = (isa(el, Code) || isa(el, LaTeX) || isa(el, Link)) && (wasCode = true)
        rstinline(io, el)
    end
end

rstinline(io::IO, md::Vector) = !isempty(md) && rstinline(io, md...)

# rstinline(io::IO, md::Image) = rstinline(io, ".. image:: ", md.url)

function rstinline(io::IO, md::Link)
    if ismatch(r":(func|obj|ref|exc|class|const|data):`\.*", md.url)
        rstinline(io, md.url)
    else
        rstinline(io, "`", md.text, " <", md.url, ">`_")
    end
end

rstinline(io::IO, f::Footnote) = print(io, "[", f.id, "]_")

rstescape(s) = replace(s, "\\", "\\\\")

rstinline(io::IO, s::AbstractString) = print(io, rstescape(s))

rstinline(io::IO, md::Bold) = rstinline(io, "**", md.text, "**")

rstinline(io::IO, md::Italic) = rstinline(io, "*", md.text, "*")

rstinline(io::IO, md::Code) = print(io, "``", md.code, "``")

rstinline(io::IO, br::LineBreak) = println(io)

rstinline(io::IO, l::LaTeX) = print(io, ":math:`", l.formula, "`")

rstinline(io::IO, x) = show(io, MIME"text/rst"(), x)

# show

Base.show(io::IO, ::MIME"text/rst", md::MD) = rst(io, md)

max_stars_repo_name: xorJane/julia

max_issues_repo_name: JuliaLangEs/julieta

max_forks_repo_name: JuliaLangEs/julieta


text:
using Meganet, BenchmarkTools

history = Pkg.dir("Meganet")*"//benchmarks//micro//bm_doubleSymLayer.jld"

TYPE = Float64

npixel = 500
nex = 1000

nchannel = 3
K      = getDenseKernel(TYPE,[4000,2000])
Bin    = randn(TYPE,nFeatOut(K),4)
Bout   = randn(TYPE,nFeatIn(K),3)
nLayer = getTVNormLayer(TYPE,[nex,4])
L      = getDoubleSymLayer(TYPE,K,nLayer,Bin,Bout)

theta = initTheta(L)

Y     = randn(TYPE,nFeatIn(L),nex)

# Warmup
Yout2,Yout2,tmp2 = apply(L,theta,Y,true)

@code_warntype apply(L,theta,Y,true)

trial = @benchmark apply(L,theta,Y,true);

Meganet.updatehistory!(history, trial)
hist = JLD.load(history, "hist")
judge(hist)


text:
export iscommutative, trred_matrix, any_order, pmaximal_overorder, phereditary_overorder, ismaximal

elem_type(::AlgAssRelOrd{S, T, U}) where {S, T, U} = AlgAssRelOrdElem{S, T, U}

elem_type(::Type{AlgAssRelOrd{S, T, U}}) where {S, T, U} = AlgAssRelOrdElem{S, T, U}

ideal_type(::AlgAssRelOrd{S, T, U}) where {S, T, U} = AlgAssRelOrdIdl{S, T, U}

ideal_type(::Type{AlgAssRelOrd{S, T, U}}) where {S, T, U} = AlgAssRelOrdIdl{S, T, U}

@doc Markdown.doc"""
    algebra(O::AlgAssRelOrd) -> AbsAlgAss

Returns the algebra which contains $O$.
"""
algebra(O::AlgAssRelOrd) = O.algebra

_algebra(O::AlgAssRelOrd) = algebra(O)

@doc Markdown.doc"""
    base_ring(O::AlgAssRelOrd) -> Union{NfAbsOrd, NfRelOrd}

Returns an order $R$ in the base ring of the algebra of $O$, such that $O$ is an $R$-order.
"""
base_ring(O::AlgAssRelOrd) = order(basis_pmatrix(O, copy = false).coeffs[1])

@doc Markdown.doc"""
    iscommutative(O::AlgAssRelOrd) -> Bool

Returns `true` if $O$ is a commutative ring and `false` otherwise.
"""
iscommutative(O::AlgAssRelOrd) = iscommutative(algebra(O))

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

@doc Markdown.doc"""
    Order(A::AbsAlgAss{<: NumFieldElem}, M::Generic.Mat{<: NumFieldElem})
      -> AlgAssRelOrd

Returns the order of $A$ with basis matrix $M$.
"""
function Order(A::AbsAlgAss{S}, M::Generic.Mat{S}) where S <: NumFieldElem
  return AlgAssRelOrd{S, fractional_ideal_type(order_type(base_ring(A))), typeof(A)}(A, deepcopy(M))
end

@doc Markdown.doc"""
    Order(A::AbsAlgAss{<: NumFieldElem}, M::PMat{<: NumFieldElem, T})
      -> AlgAssRelOrd

Returns the order of $A$ with basis pseudo-matrix $M$.
"""
function Order(A::AbsAlgAss{S}, M::PMat{S, T}) where { S <: NumFieldElem, T }
  return AlgAssRelOrd{S, T, typeof(A)}(A, deepcopy(M))
end

@doc Markdown.doc"""
    Order(A::AbsAlgAss{<: NumFieldElem}, B::Vector{<: AbsAlgAssElem{ <: NumFieldElem}})
      -> AlgAssRelOrd

Returns the order of $A$ with basis $B$.
"""
function Order(A::AbsAlgAss{S}, B::Vector{ <: AbsAlgAssElem{S} }) where { S <: NumFieldElem }
  @assert length(B) == dim(A)
  M = zero_matrix(base_ring(A), dim(A), dim(A))
  for i = 1:dim(A)
    elem_to_mat_row!(M, i, B[i])
  end
  return Order(A, M)
end
################################################################################
#
#  "Assure" functions for fields
#
################################################################################

function assure_has_basis_pmatrix(O::AlgAssRelOrd{S, T, U}) where {S, T, U}
  if isdefined(O, :basis_pmatrix)
    return nothing
  end
  if !isdefined(O, :pseudo_basis)
    error("No pseudo_basis and no basis_pmatrix defined.")
  end
  pb = pseudo_basis(O, copy = false)
  A = algebra(O)
  M = zero_matrix(base_ring(A), degree(O), degree(O))
  C = Vector{T}()
  for i = 1:degree(O)
    elem_to_mat_row!(M, i, pb[i][1])
    push!(C, deepcopy(pb[i][2]))
  end
  O.basis_pmatrix = PseudoMatrix(M, C)
  return nothing
end

function assure_has_pseudo_basis(O::AlgAssRelOrd{S, T, U}) where {S, T, U}
  if isdefined(O, :pseudo_basis)
    return nothing
  end
  if !isdefined(O, :basis_pmatrix)
    error("No pseudo_basis and no basis_pmatrix defined.")
  end
  P = basis_pmatrix(O, copy = false)
  A = algebra(O)
  pseudo_basis = Vector{Tuple{elem_type(A), T}}()
  for i = 1:degree(O)
    a = elem_from_mat_row(A, P.matrix, i)
    push!(pseudo_basis, (a, deepcopy(P.coeffs[i])))
  end
  O.pseudo_basis = pseudo_basis
  return nothing
end

function assure_has_basis_matrix(O::AlgAssRelOrd)
  if isdefined(O, :basis_matrix)
    return nothing
  end
  O.basis_matrix = basis_pmatrix(O).matrix
  return nothing
end

function assure_has_basis_mat_inv(O::AlgAssRelOrd)
  if isdefined(O, :basis_mat_inv)
    return nothing
  end
  O.basis_mat_inv = inv(basis_matrix(O, copy = false))
  return nothing
end

function assure_has_inv_coeff_ideals(O::AlgAssRelOrd)
  if isdefined(O, :inv_coeff_ideals)
    return nothing
  end
  pb = pseudo_basis(O, copy = false)
  O.inv_coeff_ideals = [ inv(pb[i][2]) for i in 1:degree(O) ]
  return nothing
end

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

@doc Markdown.doc"""
    pseudo_basis(O::AlgAssRelOrd; copy::Bool = true)

Returns the pseudo basis of $O$, i. e. a vector $v$ of pairs $(e_i, a_i)$ such
that $O = \bigoplus_i a_i e_i$, where $e_i$ is an element of `algebra(O)`
and $a_i$ a fractional ideal of `base_ring(O)`.
"""
function pseudo_basis(O::AlgAssRelOrd{S, T, U}; copy::Bool = true) where {S, T, U}
  assure_has_pseudo_basis(O)
  if copy
    return deepcopy(O.pseudo_basis)::Vector{Tuple{elem_type(U), T}}
  else
    return O.pseudo_basis::Vector{Tuple{elem_type(U), T}}
  end
end

@doc Markdown.doc"""
    basis_pmatrix(O::AlgAssRelOrd; copy::Bool = true) -> PMat

Returns the basis pseudo-matrix of $O$.
"""
function basis_pmatrix(O::AlgAssRelOrd; copy::Bool = true)
  assure_has_basis_pmatrix(O)
  if copy
    return deepcopy(O.basis_pmatrix)
  else
    return O.basis_pmatrix
  end
end

function inv_coeff_ideals(O::AlgAssRelOrd; copy::Bool = true)
  assure_has_inv_coeff_ideals(O)
  if copy
    return deepcopy(O.inv_coeff_ideals)
  else
    return O.inv_coeff_ideals
  end
end

# Returns a basis of O as Z-module
function absolute_basis(O::AlgAssRelOrd)
  pb = pseudo_basis(O, copy = false)
  res = Vector{elem_type(O)}()
  for i = 1:degree(O)
    for b in absolute_basis(pb[i][2])
      push!(res, O(b*pb[i][1]))
    end
  end
  return res
end

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

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

Returns the basis matrix of $O$, that is the basis pseudo-matrix of $O$ without
the coefficient ideals.
"""
function basis_matrix(O::AlgAssRelOrd; copy::Bool = true)
  assure_has_basis_matrix(O)
  if copy
    return deepcopy(O.basis_matrix)
  else
    return O.basis_matrix
  end
end

@doc Markdown.doc"""
    basis_mat_inv(O::AlgAssRelOrd; copy::Bool = true) -> MatElem

Returns the inverse of the basis matrix of $O$.
"""
function basis_mat_inv(O::AlgAssRelOrd; copy::Bool = true)
  assure_has_basis_mat_inv(O)
  if copy
    return deepcopy(O.basis_mat_inv)
  else
    return O.basis_mat_inv
  end
end

################################################################################
#
#  Degree
#
################################################################################

@doc Markdown.doc"""
    degree(O::AlgAssRelOrd) -> Int

Returns the dimension of the algebra containing $O$.
"""
function degree(O::AlgAssRelOrd)
  return dim(algebra(O))
end

################################################################################
#
#  Inclusion of algebra elements
#
################################################################################

function _check_elem_in_order(a::AbsAlgAssElem{S}, O::AlgAssRelOrd{S, T, V}, short::Type{Val{U}} = Val{false}) where {S, T, U, V}
  t = zero_matrix(base_ring(algebra(O)), 1, degree(O))
  elem_to_mat_row!(t, 1, a)
  t = t*basis_mat_inv(O, copy = false)
  b_pmat = basis_pmatrix(O, copy = false)
  if short == Val{true}
    for i = 1:degree(O)
      if !(t[1, i] in b_pmat.coeffs[i])
        return false
      end
    end
    return true
  else
    for i = 1:degree(O)
      if !(t[1, i] in b_pmat.coeffs[i])
        return false, Vector{S}()
      end
    end
    v = Vector{S}(undef, degree(O))
    for i = 1:degree(O)
      v[i] = deepcopy(t[1, i])
    end
    return true, v
  end
end

@doc Markdown.doc"""
    in(a::AbsAlgAssElem, O::AlgAssRelOrd) -> Bool

Returns `true` if the algebra element $a$ is in $O$ and `false` otherwise.
"""
function in(a::AbsAlgAssElem{S}, O::AlgAssRelOrd{S, T, U}) where {S, T, U}
  return _check_elem_in_order(a, O, Val{true})
end

################################################################################
#
#  Denominator in an order
#
################################################################################

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

Returns $d\in \mathbb Z$ such that $d \cdot a \in O$.
"""
function denominator(a::AbsAlgAssElem, O::AlgAssRelOrd)
  t = zero_matrix(base_ring(algebra(O)), 1, degree(O))
  elem_to_mat_row!(t, 1, a)
  t = t*basis_mat_inv(O, copy = false)
  d = fmpz(1)
  inv_coeff = inv_coeff_ideals(O, copy = false)
  for i = 1:degree(O)
    tt = inv_coeff[i]*t[1, i]
    tt = simplify(tt)
    d = lcm(d, denominator(tt))
  end
  return d
end

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

@doc Markdown.doc"""
    rand(O::AlgAssRelOrd, B::Int) -> AlgAssRelOrdElem

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

################################################################################
#
#  Print
#
################################################################################

function show(io::IO, O::AlgAssRelOrd)
  compact = get(io, :compact, false)
  if compact
    print(io, "Order of ")
    show(IOContext(io, :compact => true), algebra(O))
  else
    print(io, "Order of ")
    println(io, algebra(O))
    print(io, "with pseudo-basis ")
    pb = pseudo_basis(O, copy = false)
    for i = 1:degree(O)
      print(io, "\n(")
      show(IOContext(io, :compact => true), pb[i][1])
      print(io, ", ")
      show(IOContext(io, :compact => true), pb[i][2])
      print(io, ")")
    end
  end
end

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

@doc Markdown.doc"""
    ==(R::AlgAssRelOrd, S::AlgAssRelOrd) -> Bool

Returns `true` if $R$ and $S$ are equal and `false` otherwise.
"""
function ==(R::AlgAssRelOrd, S::AlgAssRelOrd)
  algebra(R) != algebra(S) && return false
  return basis_pmatrix(R, copy = false) == basis_pmatrix(S, copy = false)
end

################################################################################
#
#  Discriminant and Reduced Trace Matrix
#
################################################################################

@doc Markdown.doc"""
    trred_matrix(O::AlgssRelOrd) -> MatElem

Returns the reduced trace matrix $M$ of $O$, i. e. `M[i, j] = trred(b[i]*b[j])`,
where $b$ is a basis of $O$.
"""
function trred_matrix(O::AlgAssRelOrd)
  if isdefined(O, :trred_matrix)
    return deepcopy(O.trred_matrix)
  end
  A = algebra(O)
  b = pseudo_basis(O, copy = false)
  d = dim(A)
  M = zero_matrix(base_ring(A), d, d)
  for i = 1:d
    t = trred(b[i][1]*b[i][1])
    M[i, i] = t
    for j = i + 1:d
      t = trred(b[i][1]*b[j][1])
      M[i, j] = t
      M[j, i] = t
    end
  end
  O.trred_matrix = M
  return deepcopy(M)
end

@doc Markdown.doc"""
    discriminant(O::AlgssRelOrd)

Returns the discriminant of $O$.
"""
function discriminant(O::AlgAssRelOrd{S, T, U}) where {S, T, U}
  if isdefined(O, :disc)
    return O.disc::ideal_type(order_type(parent_type(S)))
  end
  d = det(trred_matrix(O))
  pb = pseudo_basis(O, copy = false)
  a = pb[1][2]^2
  for i = 2:degree(O)
    a *= pb[i][2]^2
  end
  disc = d*a
  simplify(disc)
  O.disc = numerator(disc)
  return deepcopy(O.disc)::ideal_type(order_type(parent_type(S)))
end

################################################################################
#
#  Maximal Order
#
################################################################################

function ismaximal_order_known(A::AbsAlgAss{T}) where { T <: NumFieldElem }
  return isdefined(A, :maximal_order)
end

@doc Markdown.doc"""
    maximal_order(A::AbsAlgAss{ <: NumFieldElem }) -> AlgAssRelOrd

Returns a maximal $R$-order of $A$ where $R$ is the maximal order of the base ring of $A$.
"""
function maximal_order(A::AbsAlgAss{T}) where { T <: NumFieldElem }
  if isdefined(A, :maximal_order)
    return A.maximal_order::order_type(A)
  end

  # So far ..._absolute is usually faster for linear, quadratic and cubic base fields,
  # but of course there are exceptions.
  # Feel free to adjust this if-condition.
  if base_field(base_ring(A)) == FlintQQ && degree(base_ring(A)) <= 3
    O = maximal_order_via_absolute(A)
  else
    O = maximal_order_via_relative(A)
  end
  A.maximal_order = O
  return O::order_type(A)
end

function maximal_order_via_absolute(A::AbsAlgAss{T}) where { T <: NumFieldElem }
  B, BtoA = AlgAss(A)
  C, CtoB = restrict_scalars(B, FlintQQ)
  OC = maximal_order(C)
  M = zero_matrix(base_ring(A), degree(OC), dim(A))
  for i = 1:degree(OC)
    elem_to_mat_row!(M, i, BtoA(CtoB(elem_in_algebra(basis(OC, copy = false)[i], copy = false))))
  end
  PM = sub(pseudo_hnf(PseudoMatrix(M), :lowerleft, true), (degree(OC) - dim(A) + 1):degree(OC), 1:dim(A))
  O = Order(A, PM)
  O.ismaximal = 1
  return O
end

function maximal_order_via_relative(A::AbsAlgAss{T}) where { T <: NumFieldElem }
  O = any_order(A)
  return maximal_order(O)
end

@doc Markdown.doc"""
    maximal_order(O::AlgAssRelOrd) -> AlgAssRelOrd

Returns a maximal order of the algebra of $O$ containing itself.
"""
function maximal_order(O::AlgAssRelOrd{S, T, U}) where {S, T, U}
  A = algebra(O)

  if isdefined(A, :maximal_order)
    # Check whether O \subseteq OO
    OO = A.maximal_order::AlgAssRelOrd{S, T, U}
    if _spans_subset_of_pseudohnf(basis_pmatrix(O, copy = false), basis_pmatrix(OO, copy = false), :lowerleft)
      return OO
    end
  end

  d = discriminant(O)
  fac = factor(d)

  OO = O
  for (p, e) in fac
    if e == 1
      continue
    end
    OO += pmaximal_overorder(O, p)
  end
  OO.ismaximal = 1

  if !isdefined(A, :maximal_order)
    A.maximal_order = OO
  end
  return OO::AlgAssRelOrd{S, T, U}
end

@doc Markdown.doc"""
    any_order(A::AbsAlgAss{ <: NumFieldElem }) -> AlgAssRelOrd

Returns any $R$-order of $A$ where $R$ is the maximal order of the base ring of $A$.
"""
function any_order(A::AbsAlgAss{T}) where { T <: NumFieldElem }
  K = base_ring(A)
  return any_order(A, maximal_order(K))
end

@doc Markdown.doc"""
    any_order(A::AbsAlgAss{ <: NumFieldElem}, R::Union{ NfAbsOrd, NfRelOrd })
      -> AlgAssRelOrd

Returns any $R$-order of $A$.
"""
function any_order(A::AbsAlgAss{T}, R::Union{ NfAbsOrd, NfRelOrd }) where { T <: NumFieldElem }
  K = base_ring(A)
  d = _denominator_of_mult_table(A, R)

  M = vcat(zero_matrix(K, 1, dim(A)), d*identity_matrix(K, dim(A)))
  oneA = one(A)
  for i = 1:dim(A)
    M[1, i] = deepcopy(coefficients(oneA, copy = false)[i])
  end
  PM = PseudoMatrix(M)
  PM = pseudo_hnf(PM, :lowerleft, true)
  O = Order(A, sub(PM, 2:dim(A) + 1, 1:dim(A)))
  return O
end

function _denominator_of_mult_table(A::AbsAlgAss{T}, R::Union{ NfAbsOrd, NfRelOrd }) where { T <: NumFieldElem }
  l = denominator(multiplication_table(A, copy = false)[1, 1, 1], R)
  for i = 1:dim(A)
    for j = 1:dim(A)
      for k = 1:dim(A)
        l = lcm(l, denominator(multiplication_table(A, copy = false)[i, j, k], R))
      end
    end
  end
  return l
end

_denominator_of_mult_table(A::AlgGrp{T}, R::Union{ NfAbsOrd, NfRelOrd }) where { T <: NumFieldElem } = fmpz(1)

# TODO: This is type unstable
# Requires that O is maximal and A = K^(n\times n) for a number field K.
# Computes a maximal order of type
#  (  O    ...   O  a^-1 )
#  (  :          :  :    )
#  (  O    ...   O  a^-1 )
#  (  a    ...   a  O    )
# for an ideal a of O.
# See Bley, Johnston "Computing generators of free modules over orders in group
# algebras", Prop. 5.1.
function _simple_maximal_order(O::AlgAssRelOrd, make_free::Bool = true, with_trafo::Type{Val{T}} = Val{false}) where T
  A = algebra(O)
  @assert A isa AlgMat
  n = degree(A)
  K = coefficient_ring(A)

  # Build a matrix with the first columns of basis elements of O
  M = zero_matrix(K, dim(A), n)
  for i = 1:dim(A)
    b = matrix(pseudo_basis(O, copy = false)[i][1], copy = false)
    for j = 1:n
      M[i, j] = deepcopy(b[j, 1])
    end
  end
  PM = PseudoMatrix(M, [pseudo_basis(O, copy = false)[i][2] for i in 1:dim(A)])
  PM = pseudo_hnf(PM, :upperright)

  M = sub(PM.matrix, 1:n, 1:n)
  PM = PseudoMatrix(M, PM.coeffs[1:n])
  U = similar(PM.matrix, 0, 0)
  steinitz_form!(PM, U, false)

  a = PM.coeffs[end]
  a = simplify(a)
  fl = false

  if make_free
    fl, beta = isprincipal(a)
  end

  if fl
    mul_row!(PM.matrix, nrows(PM.matrix), beta)
    a = K(1) * base_ring(PM)
  else
    d = denominator(a)
    if !isone(d)
      mul_row!(PM.matrix, nrows(PM.matrix), K(1//d))
    end
    a = a * d
  end

  M = transpose(PM.matrix)
  iM = inv(M)
  N = zero_matrix(K, dim(A), dim(A))
  for i = 1:dim(A)
    elem_to_mat_row!(N, i, iM*pseudo_basis(O, copy = false)[i][1]*M)
  end

  PN = PseudoMatrix(N, deepcopy(basis_pmatrix(O, copy = false).coeffs))
  PN = pseudo_hnf(PN, :lowerleft)

  niceorder = Order(A, PN)
  niceorder.isnice = true
  niceorder.nice_order_ideal = a

  if with_trafo == Val{true}
    return niceorder, A(iM)
  else
    return niceorder
  end
end

@doc Markdown.doc"""
    nice_order(O::AlgAssRelOrd) -> AlgAssRelOrd, AlgElem

Given a maximal order $O$ in a full matrix algebra over a number field, return a
nice maximal order $R$ and element $a$ such that $a O a^-1 = R$.
"""
function nice_order(O::AlgAssRelOrd{S, T, U}; cached::Bool = true) where {S, T, U}
  if cached && isdefined(O, :nice_order)
    return O.nice_order::Tuple{typeof(O), elem_type(U)}
  else
    sO, A = _simple_maximal_order(O, true, Val{true})
    if cached
      O.nice_order = sO, A
    end
    return sO::typeof(O), A::elem_type(U)
  end
end

function nice_order_ideal(O::AlgAssRelOrd)
  !O.isnice && error(throw("Order must be nice"))
  return O.nice_order_ideal
end

@doc Markdown.doc"""
    ismaximal(O::AlgAssRelOrd) -> Bool

Returns `true` if $O$ is a maximal order and `false` otherwise.
"""
function ismaximal(O::AlgAssRelOrd)
  if O.ismaximal == 1
    return true
  end
  if O.ismaximal == 2
    return false
  end

  A = algebra(O)
  d = discriminant(O)
  if isdefined(A, :maximal_order)
    if d == discriminant(maximal_order(A))
      O.ismaximal = 1
      return true
    else
      O.ismaximal = 2
      return false
    end
  end

  fac = factor(d)

  for (p, e) in fac
    if e == 1
      continue
    end

    d2 = discriminant(pmaximal_overorder(O, p))
    if d != d2
      O.ismaximal = 2
      return false
    end
  end
  O.ismaximal = 1
  return true
end

ismaximal_known(O::AlgAssRelOrd) = O.ismaximal != 0

################################################################################
#
#  p-hereditary / p-maximal overorders
#
################################################################################

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

Returns an order $O'$ containing $O$ such that the localization $O'_p$ is
hereditary where $p$ is a prime ideal of the base ring of $O$.
"""
function phereditary_overorder(O::AlgAssRelOrd, p::Union{ NfAbsOrdIdl, NfRelOrdIdl }; return_pradical::Type{Val{T}} = Val{false}) where T
  d = discriminant(O)
  prad = pradical(O, p)
  OO = left_order(prad)
  dd = discriminant(OO)
  while d != dd
    d = dd
    prad = pradical(OO, p)
    OO = left_order(prad)
    dd = discriminant(OO)
    if valuation(dd, p) < 2
      break
    end
  end
  if return_pradical == Val{true}
    return OO, prad
  else
    return OO
  end
end

# See Friedrichs: "Berechnung von Maximalordnungen über Dedekindringen", Algorithmus 3.16
function _pmaximal_overorder(O::AlgAssRelOrd, p::Union{ NfAbsOrdIdl, NfRelOrdIdl })
  return _pmaximal_overorder(O, pradical(O, p), p)
end

function _pmaximal_overorder(O::AlgAssRelOrd, prad::AlgAssRelOrdIdl, p::Union{ NfAbsOrdIdl, NfRelOrdIdl }; strict_containment::Bool = false)
  d = discriminant(O)
  primes = _prime_ideals_over(O, prad, p, strict_containment = strict_containment)
  for P in primes
    OO = left_order(P)
    dd = discriminant(OO)
    if valuation(dd, p) < 2
      return OO
    end
    if d != dd
      return _pmaximal_overorder(OO, p)
    end
  end
  return O
end

@doc Markdown.doc"""
    pmaximal_overorder(O::AlgAssRelOrd, p::Union{ NfAbsOrdIdl, NfRelOrdIdl })
      -> AlgAssRelOrd

Returns an order $O'$ containing $O$ such that the index $(O'':O')$ of any maximal
order $O''$ containing $O$ is not divisible by $p$ where $p$ is a prime ideal
of the base ring of $O$.
"""
function pmaximal_overorder(O::AlgAssRelOrd, p::Union{ NfAbsOrdIdl, NfRelOrdIdl })
  O, prad = phereditary_overorder(O, p, return_pradical = Val{true})
  return _pmaximal_overorder(O, prad, p, strict_containment = true)
end

################################################################################
#
#  Addition
#
################################################################################

function +(a::AlgAssRelOrd{S, T, U}, b::AlgAssRelOrd{S, T, U}) where { S, T, U}
  @assert algebra(a) === algebra(b)
  aB = basis_pmatrix(a, copy = false)
  bB = basis_pmatrix(b, copy = false)
  d = degree(a)
  PM = sub(pseudo_hnf(vcat(aB, bB), :lowerleft, true), d + 1:2*d, 1:d)
  return Order(algebra(a), PM)
end

################################################################################
#
#  Units of quotients
#
################################################################################

# Computes a generating system of U in O, where U is a set of representatives of
# the image of the projection map \pi:O^\times -> (O/g*O)^\times.
# Assumes that O is a maximal order in Mat_{n\times n}(K).
# See Bley, Johnson: "Computing generators of free modules over orders in
# group algebras", section 6.
function enum_units(O::AlgAssRelOrd, g::NfAbsOrdIdl)
  A = algebra(O)
  @assert A isa AlgMat
  @assert degree(A)^2 == dim(A)

  n = degree(A)

  K = base_ring(A)
  OK = base_ring(O)
  L = _simple_maximal_order(O)
  a = deepcopy(basis_pmatrix(L, copy = false).coeffs[end - 1])
  ai = deepcopy(basis_pmatrix(L, copy = false).coeffs[n])

  gensOKg = Vector{elem_type(K)}()
  for b in basis(OK)
    bmod = mod(b, g)
    if iszero(bmod)
      continue
    end
    push!(gensOKg, elem_in_nf(bmod))
  end

  if isone(a)
    gensinvag = gensOKg
  else
    gensinvag = Vector{elem_type(K)}()
    aig = ai*g
    for b in basis(ai)
      bmod = mod(b, aig)
      if iszero(bmod)
        continue
      end
      push!(gensinvag, bmod)
    end
  end

  if isone(a)
    gensag = gensOKg
  else
    gensag = Vector{elem_type(K)}()
    ag = a*g
    for b in basis(a)
      bmod = mod(b, ag)
      if iszero(bmod)
        continue
      end
      push!(gensag, bmod)
    end
  end

  result = Vector{elem_type(L)}()
  n1 = n - 1
  # n \nmid i, j or n \mid i, j
  for i = 1:n1
    for j = 1:n1
      if j == i
        continue
      end
      for x in gensOKg
        E = identity_matrix(K, n)
        E[i, j] = deepcopy(x)
        push!(result, L(A(E)))
      end
    end
  end

  # n \nmid i and n \mid j
  for i = 1:n1
    for x in gensag
      E = identity_matrix(K, n)
      E[i, n] = deepcopy(x)
      push!(result, L(A(E)))
    end
  end

  # n \mid i and n \nmid j
  for j = 1:n1
    for x in gensinvag
      E = identity_matrix(K, n)
      E[n, j] = deepcopy(x)
      push!(result, L(A(E)))
    end
  end

  U, mU = unit_group(OK)
  for i = 1:ngens(U)
    x = elem_in_nf(mU(U[i]))
    E = identity_matrix(K, n)
    E[1, 1] = x
    push!(result, L(A(E)))
  end
  return result
end

################################################################################
#
#  Ramifying primes
#
################################################################################

# Returns a vector of tuples (p, m, k) where p is prime ideal at which the
# algebra ramifies, m its local index (Schur index) and k its local capacity.
# See Reiner: "Maximal order" Theorem 32.1
function ramified_prime_ideals(O::AlgAssRelOrd)
  A = algebra(O)
  @assert issimple(A)
  n2 = dim(A)
  n = isqrt(n2)
  @assert n^2 == n2

  d = discriminant(O)
  facd = factor(d)
  result = Vector{Tuple{ideal_type(base_ring(O)), Int, Int}}()
  for (p, e) in facd
    k = divexact(n2 - e, n)
    m = divexact(n, k)
    push!(result, (p, m, k))
    @assert m*k == n
  end

  return result
end

################################################################################
#
#  "All" maximal orders
#
################################################################################

# Returns a vector containing a system of representatives of the maximal orders
# of A with respect to conjugation, that is, any maximal order of A is conjugated
# to one of them and no two returned orders are conjugated.
# Only works for algebras fulfilling the Eichler condition.
representatives_of_maximal_orders(A::AlgAss{nf_elem}) = representatives_of_maximal_orders(maximal_order(A))

function representatives_of_maximal_orders(O::AlgAssRelOrd)
  A = algebra(O)
  @assert issimple(A)
  @assert iseichler(A)
  @assert ismaximal(O)
  K = base_ring(A)
  OK = base_ring(O)
  n2 = dim(A)
  n = isqrt(n2)
  @assert n^2 == n2

  inf_plc = ramified_infinite_places(A)

  R, mR = ray_class_group(OK(1)*OK, inf_plc)
  S, mS = snf(R)
  if order(S) == 1
    return [ O ]
  end

  ram_primes = ramified_prime_ideals(O)

  U = Vector{elem_type(S)}()
  for i = 1:ngens(S)
    push!(U, n*S[i])
  end
  for (p, m, k) in ram_primes
    push!(U, k*(mS\(mR\p)))
  end

  SU, mSU = quo(S, U)
  if order(SU) == 1
    return [ O ]
  end

  # Each element of SU corresponds now to a maximal order.
  # We have to find a prime ideal in each of these classes.
  reps_found = Set{elem_type(SU)}()
  primes = Vector{ideal_type(OK)}()
  push!(reps_found, id(SU))
  P = PrimeIdealsSet(OK, 1, -1, indexdivisors = false, ramified = false)
  for p in P
    g = mSU(mS\(mR\p))
    if g in reps_found
      continue
    end

    push!(reps_found, g)
    push!(primes, p)

    if length(reps_found) == order(SU)
      break
    end
  end

  # For each of the prime ideals compute a maximal ideal with left order O.
  # Then the right orders of these form a system of representatives.
  max_ideals = Vector{ideal_type(O)}()
  for i = 1:length(primes)
    push!(max_ideals, maximal_integral_ideal(O, primes[i], :left))
  end

  result = Vector{typeof(O)}()
  push!(result, O)
  for i = 1:length(max_ideals)
    push!(result, right_order(max_ideals[i]))
  end

  return result
end