text:
using FinancialModels

#----- Calibration -----#
include("NMTA.jl")
import HypothesisTests: ADFTest, ApproximateTwoSampleKSTest
import Statistics: mean, quantile, std
import StatsBase.kurtosis
import GLM: lm, coef
using ARCHModels, Polynomials, JLD
struct Moments # Moments of log-returns
    μ::Float64 # Mean
    σ::Float64 # Standard deviation
    κ::Float64 # Kurtosis
    ks::Float64 # Kolmogorov-Smirnov test statistic for the difference between distributions
    hurst::Float64 # Hurst exponent: hurst < 0.5 => mean reverting; hurst == 0.5 => random walk; hurst > 0.5 => momentum
    gph::Float64 # GPH estimator representing long-range dependence
    adf::Float64 # ADF statistic representing random walk property of returns
    garch::Float64 # GARCH paramaters representing short-range dependence
    hill::Float64 # Hill estimator
    function Moments(logreturns1::Vector{Float64}, logreturns2::Vector{Float64})
        μ = mean(logreturns1); σ = std(logreturns1); κ = kurtosis(logreturns1)
        ks = ApproximateTwoSampleKSTest(logreturns1, logreturns2).δ
        hurst = HurstExponent(logreturns1)
        gph = GPH(abs.(logreturns1))
        adf = ADFTest(logreturns1, :none, 0).stat
        garch = sum(coef(ARCHModels.fit(GARCH{1, 1}, logreturns1))[2:3])
        hill = HillEstimator(logreturns1[findall(x -> (x >= quantile(logreturns1, 0.95)) && (x > 0), logreturns1)], 50)
        new(μ, σ, κ, ks, hurst, gph, adf, garch, hill)
    end
end
function HurstExponent(x, d = 100)
    N = length(x)
    if mod(N, 2) != 0 x = push!(x, (x[N - 1] + x[N]) / 2); N += 1 end
    N₁ = N₀ = min(floor(0.99 * N), N-1); dv = Divisors(N₁, d)
    for i in (N₀ + 1):N
        dw = Divisors(i, d)
        if length(dw) > length(dv) N₁ = i; dv = copy(dw) end
    end
    OptN = Int(N₁); d = dv
    x = x[1:OptN]
    RSempirical = map(i -> RS(x, i), d)
    return coeffs(Polynomials.fit(Polynomial, log10.(d), log10.(RSempirical), 1))[2] # Hurst is slope of log-log linear fit
end
function Divisors(n, n₀)
    return filter(x -> mod(n, x) == 0, n₀:floor(n/2))
end
function RS(z, n)
    y = reshape(z, (Int(n), Int(length(z) / n)))
    μ = mean(y, dims = 1)
    σ = std(y, dims = 1)
    temp = cumsum(y .- μ, dims = 1)
    return mean((maximum(temp, dims = 1) - minimum(temp, dims = 1)) / σ)
end
function GPH(x, bandwidthExponent = 0.5)
    n = length(x); g = Int(trunc(n^bandwidthExponent))
    j = 1:g; kk = 1:(n - 1)
    w = 2 .* π .* j ./ n # x .-= mean(x)
    σ = sum(x .^ 2) / n
    Σ = map(k -> sum(x[1:(n - k)] .* x[(1 + k):n]) / n, kk)
    periodogram = map(i -> σ + 2 * sum(Σ .* cos.(w[i] .* kk)), j)
    indeces = j[findall(x -> x > 0, periodogram)]
    x_reg = 2 .* log.(2 .* sin.(w[indeces] ./ 2)); y_reg = log.(periodogram[indeces] ./ (2 * π))
    regression = lm(hcat(ones(length(x_reg)), x_reg), y_reg)
    return abs(coef(regression)[2])
end
function HillEstimator(x, iterations)
    N = length(x)
    logx = log.(x)
    L = minimum(x); R = maximum(x)
    α = 1 / ((sum(logx) / N) - log(L))
    for i in 1:iterations
        C = (log(L) * (L^(-α)) - log(R) * (R^(-α))) / ((L^(-α)) - (R^(-α)))
        D = (R^α * L^α * (log(L) - log(R))^2) / (L^α - R^α)^2
        α = α * (1 + (α * (sum(logx) / N) - α * C - 1) / (α^2 * D - 1))
    end
    return α
end
function MovingBlockBootstrap(midprice::Vector{Float64}, iterations::Int64 = 1000, windowsize::Int64 = 2000)
    logreturns = diff(log.(midprice))
    bootstrapmoments = Vector{Vector{Float64}}()
    Threads.@threads for i in 1:iterations
        indeces = rand(1:(length(logreturns) - windowsize + 1), Int(ceil(length(logreturns)/windowsize)))
        bootstrapreturns = Vector{Float64}()
        for index in indeces
            append!(bootstrapreturns, logreturns[index:(index  + windowsize - 1)])
        end
        moments = Moments(bootstrapreturns[1:length(logreturns)], logreturns)
        [moments.μ moments.σ moments.κ moments.ks moments.hurst moments.gph moments.adf moments.garch moments.hill]
    end
    W = inv(cov(bootstrapmoments))
    save("Data/W.jld", "W", W)
end
function WeightedSumofSquaredErrors(parameters::Parameters, replications::Int64, W::Array{Float64, 2}, empiricalmoments::Moments, empiricallogreturns::Vector{Float64}, gateway::TradingGateway)
    errormatrix = fill(0.0, (replications, 9))
    for i in 1:replications
        microprice = InjectSimulation(gateway, parameters, seed = i)
        if !isempty(microprice)
            filter!(x -> !isnan(x), microprice)
            logreturns = diff(log.(microprice))
            try
                simulatedmoments = Moments(logreturns, empiricallogreturns)
                errormatrix[i, :] = [simulatedmoments.μ-empiricalmoments.μ simulatedmoments.σ-empiricalmoments.σ simulatedmoments.κ-empiricalmoments.κ simulatedmoments.ks-empiricalmoments.ks simulatedmoments.hurst-empiricalmoments.hurst simulatedmoments.gph-empiricalmoments.gph simulatedmoments.adf-empiricalmoments.adf simulatedmoments.garch-empiricalmoments.garch simulatedmoments.hill-empiricalmoments.hill]
            catch e
                println(e)
                errormatrix[i, :] = errormatrix[i - 1, :]
            end
        else
            return Inf
        end
    end
    GC.gc() # Garbage collection
    errors = mean(errormatrix, dims = 1)
    return (errors * W * transpose(errors))[1]
end
#---------------------------------------------------------------------------------------------------

#----- Simulation -----#
# options = Options(time = 10000, show_trace = true, show_every = 100)
# results = Simulate(JacobLealParameters(), options = options)
#---------------------------------------------------------------------------------------------------

#----- Monte-Carlo Simulation -----#
# monte_carlo_relpications = 10
# simulations = Vector{Vector{Float64}}()#[Vector{Float64}() for r in 1:monte_carlo_relpications]
# Threads.@threads for r in 1:monte_carlo_relpications # Parallel for loop\
#     print(r)
#     options = Options(time = 1000, seed = r)
# 	push!(simulations, prices(Simulate(JacobLealParameters(), options = options)))
# end
#---------------------------------------------------------------------------------------------------

#----- Visualisation -----#
# using Plots
# using Distributions: fit, Normal
# using StatsBase: autocor, quantile
# using StatsPlots: histogram, qqplot!, density, density!
# function SavePlot(p, fn::String, save::Union{Bool, String} = false, format = "png")
#     if isa(save, String)
#         savefig(p, string(save, "/", fn, format))
#     elseif save
#         savefig(p, string(fn, ".", format))
#     end
# end
# function PlotSimulation(results::SimulationResults, save::Union{Bool, String} = false, format = "png")
#     p = prices(results)
#     log_returns = diff(log.(p))
#     l = @layout([a; b{0.2h}])
#     p1 = plot(p, linecolor = :black, ylabel = "pₜ", legend = false, xticks = false)
#     p2 = plot(log_returns, seriestype = :line, linecolor = :black, legend = false, ylabel = "Log-returns", xlabel = "t")
#     simulation_plot = plot(p1, p2, layout = l, link = :x, tickfontsize = 5, guidefontsize = 7)
#     SavePlot(simulation_plot, "JacobLealSimulation", save, format)
#     return simulation_plot
# end
# function PlotLogReturnDistribution(results::SimulationResults, save::Union{Bool, String} = false, format = "png")
#     log_returns = diff(log.(prices(results)))
#     fitted_normal = fit(Normal, log_returns)
#     distribution = histogram(log_returns, normalize = :pdf, fillcolor = :orange, linecolor = :orange, xlabel = "Log returns",
#     ylabel = "Probability Density", label = "Empirical", legend = :topright, legendfontsize = 5, fg_legend = :transparent,
#     guidefontsize = 7, tickfontsize = 5, ylim = (0, 150), xlim = (-0.04, 0.04))
#     plot!(distribution, fitted_normal, line = (:black, 2), label = "Fitted Normal")
#     qqplot!(distribution, Normal, log_returns, xlabel = "Normal theoretical quantiles", ylabel = "Sample quantiles", linecolor = :black,
#     guidefontsize = 7, tickfontsize = 5, xrotation = 30, yrotation = 30, marker = (:orange, stroke(:orange), 3), legend = false,
#     inset = (1, bbox(0.1, 0.03, 0.4, 0.4)), subplot = 2, title = "Normal QQ-plot", titlefontsize = 7)
#     SavePlot(distribution, "JacobLealDistribution", save, format)
#     return distribution
# end
# function PlotLogReturnAutoCorrelations(results::SimulationResults, save::Union{Bool, String} = false, format = "png")
#     log_returns = diff(log.(prices(results)))
#     auto_corr = autocor(log_returns, 1:500; demean = false)
#     abs_auto_corr = autocor(abs.(log_returns), 1:500; demean = false)
#     auto_corr_plot = plot(auto_corr, seriestype = [:sticks, :scatter], marker = (:orange, stroke(:orange), 3), linecolor = :black, xlabel = "Lag",
#     ylabel = "Autocorrelation", legend = false, guidefontsize = 7, tickfontsize = 5, ylim = (-0.4, 0.2))
#     plot!(auto_corr_plot, [1.96 / sqrt(length(log_returns)), -1.96 / sqrt(length(log_returns))], seriestype = :hline, line = (:dash, :black, 2))
#     plot!(auto_corr_plot, abs_auto_corr, seriestype = :scatter, marker = (:orange, stroke(:orange), 3), legend = false, xlabel = "Lag",
#     ylabel = "Autocorrelation", inset = (1, bbox(0.62, 0.5, 0.4, 0.4, :top)), subplot = 2, guidefontsize = 7, tickfontsize = 5, xrotation = 30,
#     yrotation = 30, title = "Absolute log-return autocorrelation", titlefontsize = 7)
#     SavePlot(auto_corr_plot, "JacobLealAutoCorrelation", save, format)
#     return auto_corr_plot
# end
# function PlotExtremeLogReturnDistribution(results::SimulationResults, save::Union{Bool, String} = false, format = "png")
#     log_returns = diff(log.(prices(results)))
#     upper_observations = log_returns[findall(x -> x >= quantile(log_returns, 0.95), log_returns)]; sort!(upper_observations)
#     upperxₘᵢₙ = minimum(upper_observations)
#     upperα = 1 + length(upper_observations) / sum(log.(upper_observations ./ upperxₘᵢₙ))
#     upper_theoretical_quantiles = map(i -> (1 - (i / length(upper_observations))) ^ (-1 / (upperα - 1)) * upperxₘᵢₙ, 1:length(upper_observations))
#     lower_observations = -log_returns[findall(x -> x <= quantile(log_returns, 0.05), log_returns)]; sort!(lower_observations)
#     lowerxₘᵢₙ = minimum(lower_observations)
#     lowerα = 1 + length(lower_observations) / sum(log.(lower_observations ./ lowerxₘᵢₙ))
#     lower_theoretical_quantiles = map(i -> (1 - (i / length(lower_observations))) ^ (-1 / (lowerα - 1)) * lowerxₘᵢₙ, 1:length(lower_observations))
#     distribution = density(upper_observations, seriestype = [:scatter, :line], marker = (:orange, stroke(:orange), :utriangle), linecolor = :orange,
#     xlabel = string("Log return extreme percentiles"), ylabel = "Density", label = string("Upper percentiles - α = ", round(upperα, digits = 3)),
#     legend = :topright, fg_legend = :transparent, guidefontsize = 7, tickfontsize = 5)
#     density!(distribution, lower_observations, seriestype = [:scatter, :line], marker = (:orange, stroke(:orange), :dtriangle), linecolor = :orange,
#     label = string("Lower percentiles - α = ", round(lowerα, digits = 3)))
#     plot!(distribution, hcat(upper_theoretical_quantiles, upper_theoretical_quantiles), hcat(upper_observations, upper_theoretical_quantiles),
#     scale = :log10, seriestype = [:scatter :line], inset = (1, bbox(0.6, 0.3, 0.34, 0.34, :top)), subplot = 2, guidefontsize = 7,
#     tickfontsize = 5, xrotation = 30, yrotation = 30, legend = :none, xlabel = "Power-Law Theoretical Quantiles", ylabel = "Sample Quantiles",
#     linecolor = :black, marker = (:orange, stroke(:orange), 3, [:utriangle :none]), fg_legend = :transparent, title = "Power-Law QQ-plot", titlefontsize = 7)
#     plot!(distribution, [lower_theoretical_quantiles lower_theoretical_quantiles], [lower_observations lower_theoretical_quantiles],
#     seriestype = [:scatter :line], subplot = 2, linecolor = :black, marker = (:orange, stroke(:orange), 3, [:dtriangle :none]))
#     SavePlot(distribution, "JacobLealExtremeDistribution", save, format)
#     return distribution
# end
# PlotSimulation(results, false)
# PlotLogReturnDistribution(results, false)
# PlotLogReturnAutoCorrelations(results, false)
# PlotExtremeLogReturnDistribution(results, false)
#---------------------------------------------------------------------------------------------------

max_stars_repo_name: IvanJericevich/FinancialModels

max_issues_repo_name: IvanJericevich/FinancialModels

max_forks_repo_name: IvanJericevich/FinancialModels


text:
module Romeo

include("_utils.jl")

include("lti/LTI.jl")
include("optimization/Optimization.jl")

end # module

ext: jl


text:
# Only works from within IJulia, since uses widgets

module RepoHistoryBrowser

export browse_history

using Interact

# Concatenate two HTML representations
# HTML is a type defined in Interact.jl
+(a::HTML, b::HTML) = html(string(a.value, b.value))

# Syntax for using git on another directory:
# git -C ~/Dropbox/projects/RepoBrowser/ log --oneline

# Get the copy of a filename from a particular commit with a given hash, in the given directory
function use_commit(directory::String, hash::String, filename::String)
    try
        run(`git -C $directory show $hash:$filename` |> "_tmp.jl")  # pipe output to a temporary file
    catch
        run(`echo "No file $filename in commit $hash"` |> "_tmp.jl" )
    end
end

# Use pygments to make the highlighted HTML version:
function generate_html(directory::String, hash::String, filename::String)
    use_commit(directory, hash, filename)
    code = readall(`pygmentize -O full -f html -l julia _tmp.jl`);
    html(code)
end

# Main function: browse the history of a file in a git repo in the given directory

function browse_history(directory::String, filename::String)
    commits = []

    try
        commits = reverse(readlines(`git -C $directory log --oneline`))
    catch
        error("Not possible to get the history for $directory. Is it a git repository?")
    end

    hashes = ASCIIString[split(commit)[1] for commit in commits];

    @manipulate for n in 1:length(hashes)
        html(commits[n]) + html("<p>") + generate_html(directory, hashes[n], filename)
    end

end

# browse in the current directory if no directory given
browse_history(filename::String) = browse_history(pwd(), filename)

end

max_stars_repo_name: dpsanders/RepoBrowser

max_issues_repo_name: dpsanders/RepoBrowser

max_forks_repo_name: dpsanders/RepoBrowser


max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
############################################################
# Power Layer
############################################################
@defstruct PowerLayer Layer (
  name :: AbstractString = "power",
  (power :: Number = 1, isreal(power)),
  (scale :: Number = 1, isreal(scale)),
  (shift :: Number = 0, isreal(shift)),
  (tops :: Vector{Symbol} = Symbol[], length(tops) > 0),
  (bottoms :: Vector{Symbol} = Symbol[], length(bottoms) == length(tops)),
)
@characterize_layer(PowerLayer,
  can_do_bp => true
)

type PowerLayerState <: LayerState
  layer      :: PowerLayer
  blobs      :: Vector{Blob}
  blobs_diff :: Vector{Blob}
end

function setup(backend::Backend, layer::PowerLayer, inputs::Vector{Blob}, diffs::Vector{Blob})
  blobs = Blob[make_blob(backend, eltype(x), size(x)) for x in inputs]
  blobs_diff = Array{Blob}(length(inputs))
  for i = 1:length(inputs)
    # if the bottom layer does not need back propagate, I don't need, either
    if isa(diffs[i], NullBlob)
      blobs_diff[i] = NullBlob()
    else
      blobs_diff[i] = make_blob(backend, eltype(inputs[i]), size(inputs[i]))
    end
  end
  state = PowerLayerState(layer, blobs, blobs_diff)
end
function shutdown(backend::Backend, state::PowerLayerState)
  map(destroy, state.blobs)
  map(destroy, state.blobs_diff)
end

function forward(backend::CPUBackend, state::PowerLayerState, inputs::Vector{Blob})
  for i = 1:length(inputs)
    input = inputs[i]
    output = state.blobs[i]

    # output = input
    copy!(output, input)

    # output *= scale
    if state.layer.scale != 1
      BLAS.scal!(length(output), convert(eltype(output),state.layer.scale), output.data, 1)
    end

    if state.layer.shift != 0
      # output += shift
      Vec.add_scal!(output.data, state.layer.shift)
    end

    # output = output ^ power
    if state.layer.power != 1
      if state.layer.power == 2
        Vec.mul!(output.data, output.data)
      else
        Vec.pow!(output.data, state.layer.power)
      end
    end
  end
end

function backward(backend::CPUBackend, state::PowerLayerState,
    inputs::Vector{Blob}, diffs::Vector{Blob})

  data_type = eltype(inputs[1])
  pow_scale = convert(data_type,state.layer.power * state.layer.scale)
  for i = 1:length(inputs)
    diff = diffs[i]
    if isa(diff, NullBlob)
      continue
    end
    if state.layer.power == 1 || state.layer.scale == 0
      # trivial case, derivative is constant
      fill!(diff, pow_scale)
    else
      input = inputs[i]
      output = state.blobs[i]

      erase!(diff)

      if state.layer.power == 2
        # dO/dI = 2 * scale * (scale * I + shift)
        #       = pow_scale * scale * I + pow_scale * shift
        BLAS.axpy!(length(input), convert(data_type, pow_scale*state.layer.scale),
            input.data, 1, diff.data, 1)
        if state.layer.shift != 0
          Vec.add_scal!(diff.data, pow_scale * state.layer.shift)
        end
      elseif state.layer.shift == 0
        # dO/dI = power * scale * (scale * I) ^ (power - 1)
        #       = power * O / I
        BLAS.axpy!(length(input), convert(data_type,state.layer.power),
            output.data, 1, diff.data, 1)
        Vec.div!(diff.data, input.data)
      else
        # general case
        # dO/dI = power * scale * (scale * I + shift) ^ (power - 1)
        #       = power * scale * O / (scale * I + shift)
        copy!(diff, input)
        if state.layer.scale != 1
          BLAS.scal!(length(diff), convert(data_type,state.layer.scale), diff.data, 1)
        end
        Vec.add_scal!(diff.data, state.layer.shift)
        Vec.div2!(output.data, diff.data)
        BLAS.scal!(length(diff), pow_scale, diff.data, 1)
      end
    end
    Vec.mul!(diff.data, state.blobs_diff[i].data)
  end
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using AsterReader

"""
Nodes are different order in Code Aster compared to ABAQUS. This is yet
incomplete mapping between the permutations. Most of this is still a
great mystery.

# References
- http://onelab.info/pipermail/gmsh/2008/003850.html
- http://caelinux.org/wiki/index.php/Proj:UNVConvert

"""
const med_connectivity = Dict{Symbol, Vector{Int}}(
    :Tet4   => [4,3,1,2],
    :Tet10  => [4,3,1,2,10,7,8,9,6,5],
    :Pyr5   => [1,4,3,2,5],
    :Wedge6 => [4,5,6,1,2,3],
    :Hex8   => [4,8,7,3,1,5,6,2],
    :Hex20  => [4,8,7,3,1,5,6,2,20,15,19,11,12,16,14,10,17,13,18,9],
    :Hex27  => [4,8,7,3,1,5,6,2,20,15,19,11,12,16,14,10,17,13,18,9,24,25,26,23,21,22,27])

"""
Map element names used in in Code Aster to element names used in in JuliaFEM
"""
const med_element_names = Dict{Symbol, Symbol}(
    :PO1 => :Poi1,
    :SE2 => :Seg2,
    :SE3 => :Seg3,
    :SE4 => :Seg4,
    :TR3 => :Tri3,
    :TR6 => :Tri6,
    :TR7 => :Tri7,
    :QU4 => :Quad4,
    :QU8 => :Quad8,
    :QU9 => :Quad9,
    :TE4 => :Tet4,
    :T10 => :Tet10,
    :PE6 => :Wedge6,
    :P15 => :Wedge15,
    :P18 => :Wedge18,
    :HE8 => :Hex8,
    :H20 => :Hex20,
    :H27 => :Hex27,
    :PY5 => :Pyr5,
    :P13 => :Pyr13)

"""
    aster_read_mesh(filename, mesh_name=nothing; reorder_element_connectivity=true)

Read code aster mesh from file and return `Mesh` structure.

If mesh file contains several meshes, a name of mesh must be given. By default,
elements are reordered so that they match to the conventions used in JuliaFEM.
"""
function aster_read_mesh(filename::String, mesh_name=nothing; reorder_element_connectivity=true)
    m = AsterReader.aster_read_mesh(filename, mesh_name)
    mesh = Mesh(m)
    for (elid, eltype) in mesh.element_types
        mesh.element_types[elid] = med_element_names[eltype]
    end
    if reorder_element_connectivity
        reorder_element_connectivity!(mesh, med_connectivity)
    end
    nnodes = length(mesh.nodes)
    nelements = length(mesh.elements)
    @info("Mesh parsed from Code Aster file $filename.")
    @info("Mesh contains $nnodes nodes and $nelements elements.")
    for (elset_name, elset_elids) in mesh.element_sets
        content = Dict{Symbol, Int}()
        for elid in elset_elids
            eltype = mesh.element_types[elid]
            content[eltype] = get(content, eltype, 0) + 1
        end
        s = join(("$v x $k" for (k, v) in content), ", ")
        nels = length(elset_elids)
        @info("Element set $elset_name contains $nels elements ($s).")
    end
    return mesh
end


text:
using DrWatson
@quickactivate :Continuation

using SparseArrays, Pardiso, LinearAlgebra
using PyCall, Printf

ddir = "/media/hdd/chimera/pathfollow_fftg001/"
#ddir = "../../chimera/pathfollow_fftg001/"

function readJac(filename)
    py"""
    import numpy as np
    import scipy.sparse as sp
    def read_jac(filename):
        mat = sp.coo_matrix(sp.load_npz(filename))
        return mat.row + 1, mat.col + 1, mat.data
    """
    sparse(py"read_jac"(filename)...)
end

function readJac(N, σ, α, γ)
    wdir = params_to_fname(N, σ, α, γ)
    if isfile(wdir * "/jac2.npz")
        readJac(wdir * "/jac2.npz")
    else
        readJac(wdir * "/jac1.npz")
    end
end

function float2str(fl::Real)
    replace(string(fl), "." => "")
end

function a2str(α)
    replace((@sprintf "a%.4f" α), "." => "")
end

function params_to_fname(N, σ, α, γ)
    ddir * "N$N" * "s040" * a2str(α) * "g" * float2str(γ)
end

N, σ, α, γ = 128, 0.4, 0.8073, 0.01
jac = readJac(N, σ, α, γ)

using PyPlot
pygui(true)
spy(jac)

diff = jac - transpose(jac)

maximum(diff)

ps = PardisoSolver()

val, vec = eigen(jac)

using BenchmarkTools

# 2.86 GB for N=64 ~ 7s
@btime lu(jac);
F = lu(jac);

detlog = sum(log10.(Complex.(diag(F.L)))) + sum(log10.(Complex.(diag(F.U))))

det(F)


text:
import Base.getindex

export NearestInterpolation,
	getindex,
	getindex_symbolic

abstract NearestInterpolation <: SimpleInterpolation

@inline function getindex{T<:Real}(interp::Type{NearestInterpolation}, x::T)
	return ((one(T), round(Int, x)),)
end

function getindex_symbolic(interp::Type{NearestInterpolation}, x::Union{Symbol, Expr})
	@gensym a b
	return :(), ((:(1), :(round(Int, $x))),)
end

max_stars_repo_name: JKrehl/Interpol

max_issues_repo_name: JKrehl/Interpol

max_forks_repo_name: JKrehl/Interpol


text:
@testset "PenaltyWrapper Inprocessing" begin
	X, y = @load_compas;
	model = @pipeline ContinuousEncoder @load(RandomForestClassifier, pkg=DecisionTree)
	wrappedModel = PenaltyWrapper(classifier=model, grp=:sex, measure=tpr, n_iters=10, lr=0.1, alpha=1)
	mach = machine(wrappedModel, X, y)
	fit!(mach)
	@test length(predict(mach, X)) == length(y)
	# TODO:test the fitresult of machines
end


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

mutable struct VisualizationCallback{SolutionVariables, VariableNames, PlotDataCreator, PlotCreator}
  interval::Int
  solution_variables::SolutionVariables
  variable_names::VariableNames
  show_mesh::Bool
  plot_data_creator::PlotDataCreator
  plot_creator::PlotCreator
  plot_arguments::Dict{Symbol,Any}
end

function Base.show(io::IO, cb::DiscreteCallback{Condition,Affect!}) where {Condition, Affect!<:VisualizationCallback}
  visualization_callback = cb.affect!
  @unpack interval, plot_arguments, solution_variables, variable_names, show_mesh, plot_creator, plot_data_creator = visualization_callback
  print(io, "VisualizationCallback(",
            "interval=", interval, ", ",
            "solution_variables=", solution_variables, ", ",
            "variable_names=", variable_names, ", ",
            "show_mesh=", show_mesh, ", ",
            "plot_data_creator=", plot_data_creator, ", ",
            "plot_creator=", plot_creator, ", ",
            "plot_arguments=", plot_arguments, ")")
end

function Base.show(io::IO, ::MIME"text/plain", cb::DiscreteCallback{Condition,Affect!}) where {Condition, Affect!<:VisualizationCallback}
  if get(io, :compact, false)
    show(io, cb)
  else
    visualization_callback = cb.affect!

    setup = [
             "interval" => visualization_callback.interval,
             "plot arguments" => visualization_callback.plot_arguments,
             "solution variables" => visualization_callback.solution_variables,
             "variable names" => visualization_callback.variable_names,
             "show mesh" => visualization_callback.show_mesh,
             "plot creator" => visualization_callback.plot_creator,
             "plot data creator" => visualization_callback.plot_data_creator,
            ]
    summary_box(io, "VisualizationCallback", setup)
  end
end

"""
    VisualizationCallback(; interval=0,
                            solution_variables=cons2prim,
                            variable_names=[],
                            show_mesh=false,
                            plot_data_creator=PlotData2D,
                            plot_creator=show_plot,
                            plot_arguments...)

Create a callback that visualizes results during a simulation, also known as *in-situ
visualization*.

!!! warning "Experimental implementation"
    This is an experimental feature and may change in any future releases.

The `interval` specifies the number of time step iterations after which a new plot is generated. The
available variables to plot are configured with the `solution_variables` parameter, which acts the
same way as for the [`SaveSolutionCallback`](@ref). The variables to be actually plotted can be
selected by providing a single string or a list of strings to `variable_names`, and if `show_mesh`
is `true`, an additional plot with the mesh will be generated.

To customize the generated figure, `plot_data_creator` allows to use different plot data types. With
`plot_creator` you can further specify an own function to visualize results, which must support the
same interface as the default implementation [`show_plot`](@ref). All remaining
keyword arguments are collected and passed as additional arguments to the plotting command.
"""
function VisualizationCallback(; interval=0,
                                 solution_variables=cons2prim,
                                 variable_names=[],
                                 show_mesh=false,
                                 plot_data_creator=PlotData2D,
                                 plot_creator=show_plot,
                                 plot_arguments...)
  mpi_isparallel() && error("this callback does not work in parallel yet")

  if variable_names isa String
    variable_names = String[variable_names]
  end

  visualization_callback = VisualizationCallback(interval,
                                                 solution_variables, variable_names, show_mesh,
                                                 plot_data_creator, plot_creator,
                                                 Dict{Symbol,Any}(plot_arguments))

  # Warn users if they create a visualization callback without having loaded the Plots package
  #
  # Note: This warning is added for convenience, as Plots is the only "officially" supported
  #       visualization package right now. However, in general nothing prevents anyone from using
  #       other packages such as Makie, Gadfly etc., given that appropriate `plot_creator`s are
  #       passed. This is also the reason why the visualization callback is not included via
  #       Requires.jl only when Plots is present.
  #       In the future, we should update/remove this warning if other plotting packages are
  #       starting to be used.
  if !(:Plots in names(@__MODULE__, all=true))
    @warn "Package `Plots` not loaded but required by `VisualizationCallback` to visualize results"
  end

  DiscreteCallback(visualization_callback, visualization_callback, # the first one is the condition, the second the affect!
                   save_positions=(false,false),
                   initialize=initialize!)
end

function initialize!(cb::DiscreteCallback{Condition,Affect!}, u, t, integrator) where {Condition, Affect!<:VisualizationCallback}
  visualization_callback = cb.affect!

  visualization_callback(integrator)

  return nothing
end

# this method is called to determine whether the callback should be activated
function (visualization_callback::VisualizationCallback)(u, t, integrator)
  @unpack interval = visualization_callback

  # With error-based step size control, some steps can be rejected. Thus,
  #   `integrator.iter >= integrator.destats.naccept`
  #    (total #steps)       (#accepted steps)
  # We need to check the number of accepted steps since callbacks are not
  # activated after a rejected step.
  return interval > 0 && ( (integrator.destats.naccept % interval == 0 &&
                           !(integrator.destats.naccept == 0 && integrator.iter > 0)) ||
                          isfinished(integrator))
end

# this method is called when the callback is activated
function (visualization_callback::VisualizationCallback)(integrator)
  u_ode = integrator.u
  semi = integrator.p
  @unpack plot_arguments, solution_variables, variable_names, show_mesh, plot_data_creator, plot_creator = visualization_callback

  # Extract plot data
  plot_data = plot_data_creator(u_ode, semi, solution_variables=solution_variables)

  # If variable names were not specified, plot everything
  if isempty(variable_names)
    variable_names = String[keys(plot_data)...]
  end

  # Create plot
  plot_creator(plot_data, variable_names;
               show_mesh=show_mesh, plot_arguments=plot_arguments,
               time=integrator.t, timestep=integrator.destats.naccept)

  # avoid re-evaluating possible FSAL stages
  u_modified!(integrator, false)
  return nothing
end

"""
    show_plot(plot_data, variable_names;
              show_mesh=true, plot_arguments=Dict{Symbol,Any}(),
              time=nothing, timestep=nothing)

Visualize the plot data object provided in `plot_data` and display result, plotting only the
variables in `variable_names` and, optionally, the mesh (if `show_mesh` is `true`).  Additionally,
`plot_arguments` will be unpacked and passed as keyword arguments to the `Plots.plot` command.

This function is the default `plot_creator` argument for the [`VisualizationCallback`](@ref).
`time` and `timestep` are currently unused by this function.

!!! warning "Experimental implementation"
    This is an experimental feature and may change in future releases.

See also: [`VisualizationCallback`](@ref), [`save_plot`](@ref)
"""
function show_plot(plot_data, variable_names;
                   show_mesh=true, plot_arguments=Dict{Symbol,Any}(),
                   time=nothing, timestep=nothing)
  # Gather subplots
  plots = []
  for v in variable_names
    push!(plots, Plots.plot(plot_data[v]; plot_arguments...))
  end
  if show_mesh
    push!(plots, Plots.plot(getmesh(plot_data); plot_arguments...))
  end

  # Determine layout
  cols = ceil(Int, sqrt(length(plots)))
  rows = div(length(plots), cols, RoundUp)
  layout = (rows, cols)

  # Show plot
  display(Plots.plot(plots..., layout=layout))
end

"""
    save_plot(plot_data, variable_names;
              show_mesh=true, plot_arguments=Dict{Symbol,Any}(),
              time=nothing, timestep=nothing)

Visualize the plot data object provided in `plot_data` and save result as a PNG file in the `out`
directory, plotting only the variables in `variable_names` and, optionally, the mesh (if `show_mesh`
is `true`).  Additionally, `plot_arguments` will be unpacked and passed as keyword arguments to the
`Plots.plot` command.

The `timestep` is used in the filename. `time` is currently unused by this function.

!!! warning "Experimental implementation"
    This is an experimental feature and may change in future releases.

See also: [`VisualizationCallback`](@ref), [`show_plot`](@ref)
"""
function save_plot(plot_data, variable_names;
                   show_mesh=true, plot_arguments=Dict{Symbol,Any}(),
                   time=nothing, timestep=nothing)
  # Gather subplots
  plots = []
  for v in variable_names
    push!(plots, Plots.plot(plot_data[v]; plot_arguments...))
  end
  if show_mesh
    push!(plots, Plots.plot(getmesh(plot_data); plot_arguments...))
  end

  # Determine layout
  cols = ceil(Int, sqrt(length(plots)))
  rows = div(length(plots), cols, RoundUp)
  layout = (rows, cols)

  # Create plot
  Plots.plot(plots..., layout=layout)

  # Determine filename and save plot
  filename = joinpath("out", @sprintf("solution_%06d.png", timestep))
  Plots.savefig(filename)
end

end # @muladd


text:
# Thanks @tlienart

using Franklin, JSON

using Markdown, Dates

include("youtube_videos.jl")

const DATEFMT = dateformat"yyyy-mm-dd HH:MMp"
const TZ = "America/New_York"

function hfun_plutonotebookpage(params)
    path = params[1]
    path_to_html = if endswith(path, ".jl")
        path[1:end-3] * ".html"
    else
        path * ".html"
    end

    # https://developer.mozilla.org/en-US/docs/Web/HTTP/Headers/Feature-Policy#directives

    return """

    <iframe width="100%" height="100%"
    src="$(path_to_html)"
    class="plutopage"
    frameborder="0"
    allow="accelerometer; ambient-light-sensor; autoplay; battery; camera; display-capture; document-domain; encrypted-media; execution-while-not-rendered; execution-while-out-of-viewport; fullscreen; geolocation; gyroscope; layout-animations; legacy-image-formats; magnetometer; microphone; midi; navigation-override; oversized-images; payment; picture-in-picture; publickey-credentials-get; sync-xhr; usb; wake-lock; screen-wake-lock; vr; web-share; xr-spatial-tracking"
    allowfullscreen>

    """
end

function hfun_doc(params)
    fname = join(params[1:max(1, length(params)-2)], " ")
    head = params[end-1]
    type = params[end]
    doc = eval(Meta.parse("@doc $fname"))
    txt = Markdown.plain(doc)
    # possibly further processing here
    body = Franklin.fd2html(txt, internal=true)
    return """

            $head

$type

$body

    """
end

function hfun_youtube(params)
    id = params[1]
    return """

    """
end

function hfun_showtime(params)
    id = params[1]
    str = locvar(id)
    if isnothing(str)
        @warn "Unknown datetime variable $str"
        return ""
    end
    try
        DateTime(str, DATEFMT)
    catch err
        @warn "There was an error parsing date $str, the format is yyyy-mm-dd HH:MMp (see ?DateFormat)"
        rethrow(err)
    end
end

function parse_duration(str)
    str = replace(str, r"^PT"=>"")
    hrex, mrex, srex = Regex.(string.("^([0-9]+)", ["H","M","S"]))

    t = 0
    hmatch = match(hrex, str)
    if !isnothing(hmatch)
        h = parse(Int, hmatch[1])
        t += 60*60*h
        str = replace(str, hrex=>"")
    end

    mmatch = match(mrex, str)
    if !isnothing(mmatch)
        m = parse(Int, mmatch[1])
        t += 60*m
        str = replace(str, mrex=>"")
    end

    smatch = match(srex, str)
    if !isnothing(smatch)
        s = parse(Int, smatch[1])
        t += s
        str = replace(str, srex=>"")
    end

    t
end

function hfun_go_live()
    seq = locvar("sequence")
    airtime = locvar("airtime")

    if isnothing(seq)
        @warn "airtime set, but no `sequence` variable not defined." *
        "sequence is an array of video IDs to play in order on this page"
    end

    vid_ids = [get(videos, s, s) for s in seq]

    f = tempname()
    # Get the duration of each video
    download("https://www.googleapis.com/youtube/v3/videos?id=$(join(vid_ids, ","))&part=contentDetails&key=AIzaSyDZhbWHc2PTEFTx173MaTgddnWCGPqdbB8", f)
    dict = JSON.parse(String(read(f)))

    durations = [parse_duration(video["contentDetails"]["duration"])
                 for video in dict["items"]]

    jrepr(x) = sprint(io->JSON.print(io, x))
    """

    """
end


text:
################################################################################
#
#    NfOrd/NfOrd.jl : Orders in absolute number fields
#
# This file is part of Hecke.
#
# Copyright (c) 2015, 2016: Claus Fieker, Tommy Hofmann
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# * Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
#
#  Copyright (C) 2015, 2016, 2017 Tommy Hofmann
#
################################################################################

export ==, +, basis, basis_mat, basis_mat_inv, discriminant, degree, den,
       gen_index, EquationOrder, index, isequation_order, isindex_divisor, lll,
       lll_basis, maximal_order, MaximalOrder, minkowski_mat, nf,
       norm_change_const, Order, parent, poverorder, pmaximal_overorder,
       ring_of_integers, signature, trace_matrix, different, codifferent,
       reduced_discriminant

################################################################################
#
#  Make NfOrd fully working Nemo ring
#
################################################################################

Nemo.parent_type(::Type{NfAbsOrdElem{S, T}}) where {S, T} = NfAbsOrd{S, T}

Nemo.elem_type(::Type{NfAbsOrd{S, T}}) where {S, T} = NfAbsOrdElem{S, T}

Nemo.elem_type(::NfAbsOrd{S, T}) where {S, T} = NfAbsOrdElem{S, T}

ideal_type(::NfAbsOrd{S, T}) where {S, T} = NfAbsOrdIdl{S, T}

ideal_type(::Type{NfAbsOrd{S, T}}) where {S, T} = NfAbsOrdIdl{S, T}

Nemo.show_minus_one(::Type{NfAbsOrdElem{S, T}}) where {S, T} = true

Nemo.base_ring(::NfAbsOrd) = Union{}

Nemo.needs_parentheses(::NfAbsOrdElem) = true

Nemo.isnegative(::NfAbsOrdElem) = false

################################################################################
#
#  Basic field access
#
################################################################################

doc"""
    nf(O::NfOrd) -> AnticNumberField

Returns the ambient number field of $\mathcal O$.
"""
@inline nf(O::NfAbsOrd) = O.nf

doc"""
    parent(O::NfOrd) -> NfOrdSet

Returns the parent of $\mathcal O$, that is, the set of orders of the ambient
number field.
"""
parent(O::NfOrd) = NfAbsOrdSet(nf(O), false)

doc"""
    isequation_order(O::NfOrd) -> Bool

Returns whether $\mathcal O$ is the equation order of the ambient number
field.
"""
@inline isequation_order(O::NfAbsOrd) = O.isequation_order

@inline ismaximal_known(O::NfAbsOrd) = O.ismaximal != 0

# The following function should actually do more!
@inline ismaximal(O::NfAbsOrd) = O.ismaximal == 1

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

doc"""
    degree(O::NfOrd) -> Int

Returns the degree of $\mathcal O$.
"""
degree(O::NfAbsOrd) = degree(O.nf)

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

function assure_has_basis(O::NfAbsOrd)
  if isdefined(O, :basis_ord)
    return nothing
  end
  b = O.basis_nf
  d = degree(O)
  B = Vector{elem_type(O)}(d)
  for i in 1:length(b)
    v = [fmpz(0) for j in 1:d]
    one!(v[i])
    B[i] = NfAbsOrdElem(O, b[i], v)
  end
  O.basis_ord = B
  return nothing
end

function assure_has_basis_mat(O::NfAbsOrd)
  if isdefined(O, :basis_mat)
    return nothing
  end
  A = O.basis_nf
  O.basis_mat = FakeFmpqMat(basis_mat(A))
  return nothing
end

function assure_has_basis_mat_inv(O::NfAbsOrd)
  if isdefined(O, :basis_mat_inv)
    return nothing
  end
  O.basis_mat_inv = inv(basis_mat(O, Val{false}))
  return nothing
end

################################################################################
#
#  Basis
#
################################################################################

function basis_ord(O::NfAbsOrd, copy::Type{Val{T}} = Val{true}) where T
  assure_has_basis(O)
  if copy == Val{true}
    return deepcopy(O.basis_ord)::Vector{elem_type(O)}
  else
    return O.basis_ord::Vector{elem_type(O)}
  end
end

doc"""
    basis(O::NfOrd) -> Array{NfOrdElem, 1}

Returns the $\mathbf Z$-basis of $\mathcal O$.
"""
@inline basis{T}(O::NfAbsOrd, copy::Type{Val{T}} = Val{true}) = basis_ord(O, copy)

doc"""
    basis(O::NfOrd, K::AnticNumberField) -> Array{nf_elem, 1}

Returns the $\mathbf Z$-basis elements of $\mathcal O$ as elements of the
ambient number field.
"""
function basis(O::NfAbsOrd, K::AnticNumberField)
  nf(O) != K && error()
  return deepcopy(O.basis_nf)
end

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

doc"""
    basis_mat(O::NfOrd) -> FakeFmpqMat

Returns the basis matrix of $\mathcal O$ with respect to the power basis
of the ambient number field.
"""
function basis_mat(O::NfAbsOrd, copy::Type{Val{T}} = Val{true}) where T
  assure_has_basis_mat(O)
  if copy == Val{true}
    return deepcopy(O.basis_mat)
  else
    return O.basis_mat
  end
end

doc"""
    basis_mat_inv(O::NfOrd) -> FakeFmpqMat

Returns the inverse of the basis matrix of $\mathcal O$.
"""
function basis_mat_inv(O::NfAbsOrd, copy::Type{Val{T}} = Val{true}) where T
  assure_has_basis_mat_inv(O)
  if copy == Val{true}
    return deepcopy(O.basis_mat_inv)
  else
    return O.basis_mat_inv
  end
end

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

function show(io::IO, S::NfOrdSet)
  print(io, "Set of orders of the number field ")
  print(io, S.nf)
end

function show(io::IO, O::NfAbsOrd)
  if ismaximal_known(O) && ismaximal(O)
    show_maximal(io, O)
  else
    show_gen(io, O)
  end
end

function show_gen(io::IO, O::NfAbsOrd)
  print(io, "Order of ")
  println(io, nf(O))
  print(io, "with Z-basis ")
  print(io, basis(O))
end

function show_maximal(io::IO, O::NfAbsOrd)
  print(io, "Maximal order of $(nf(O)) \nwith basis $(O.basis_nf)")
end

################################################################################
#
#  Discriminant
#
################################################################################

function assure_has_discriminant(O::NfAbsOrd)
  if isdefined(O, :disc)
    return nothing
  else
    if isequation_order(O) && issimple(nf(O))
      O.disc = numerator(discriminant(nf(O).pol))
    else
      O.disc = discriminant(basis(O))
    end
  end
  return nothing
end

doc"""
    discriminant(O::NfOrd) -> fmpz

Returns the discriminant of $\mathcal O$.
"""
function discriminant(O::NfAbsOrd)
  assure_has_discriminant(O)
  return O.disc
end

#TODO: compute differently in equation orders, this is the rres...
doc"""
   reduced_discriminant(O::NfOrd) -> fmpz
> Returns the reduced discriminant, ie. the largest elementary divisor of the
> trace matrix of $\mathcal O$.
"""
function reduced_discriminant(O::NfOrd)
  if isequation_order(O)
    Zx = PolynomialRing(FlintZZ, cached = false)[1]
    f = Zx(nf(O).pol)
    return rres(f, derivative(f))
  end
  return maximal_elementary_divisor(trace_matrix(O))
end

################################################################################
#
#  (Generalized) index
#
################################################################################

doc"""
    gen_index(O::NfOrd) -> fmpq

Returns the generalized index of $\mathcal O$ with respect to the equation
order of the ambient number field.
"""
function gen_index(O::NfAbsOrd)
  if isdefined(O, :gen_index)
    return deepcopy(O.gen_index)
  else
    O.gen_index = inv(det(basis_mat(O, Val{false})))#FlintQQ(O.basis_mat.den^degree(O), det(O.basis_mat.num))
    return deepcopy(O.gen_index)
  end
end

doc"""
    index(O::NfOrd) -> fmpz

Assuming that the order $\mathcal O$ contains the equation order
$\mathbf Z[\alpha]$ of the ambient number field, this function returns the
index $[ \mathcal O : \mathbf Z]$.
"""
function index(O::NfAbsOrd)
  if isdefined(O, :index)
    return deepcopy(O.index)
  else
    i = gen_index(O)
    !isone(denominator(i)) && error("Order does not contain the equation order")
    O.index = numerator(i)
    return deepcopy(O.index)
  end
end

################################################################################
#
#  Index divisor
#
################################################################################

doc"""
    isindex_divisor(O::NfOrd, d::fmpz) -> Bool
    isindex_divisor(O::NfOrd, d::Int) -> Bool

Returns whether $d$ is a divisor of the index of $\mathcal O$. It is assumed
that $\mathcal O$ contains the equation order of the ambient number field.
"""
function isindex_divisor(O::NfOrd, d::Union{fmpz, Integer})
  i = index(O)
  return i % d == 0
end

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

doc"""
    deepcopy(O::NfOrd) -> NfOrd

Makes a copy of $\mathcal O$.
"""
function Base.deepcopy_internal(O::NfAbsOrd{S, T}, dict::ObjectIdDict) where {S, T}
  z = NfAbsOrd{S, T}(O.nf)
  for x in fieldnames(O)
    if x != :nf && isdefined(O, x)
      setfield!(z, x, Base.deepcopy_internal(getfield(O, x), dict))
    end
  end
  return z
end

################################################################################
#
#  Signature
#
################################################################################

doc"""
    signature(O::NfOrd) -> Tuple{Int, Int}

Returns the signature of the ambient number field of $\mathcal O$.
"""
function signature(x::NfOrd)
  if x.signature[1] != -1
    return x.signature
  else
    x.signature = signature(nf(x))
    return x.signature
  end
end

################################################################################
#
#  Minkowski matrix
#
################################################################################

doc"""
    minkowski_mat(O::NfOrd, abs_tol::Int = 64) -> arb_mat

Returns the Minkowski matrix of $\mathcal O$.  Thus if $\mathcal O$ has degree
$d$, then the result is a matrix in $\operatorname{Mat}_{d\times d}(\mathbf
R)$. The entries of the matrix are real balls of type `arb` with radius less
then `2^-abs_tol`.
"""
function minkowski_mat(O::NfOrd, abs_tol::Int = 64)
  if isdefined(O, :minkowski_mat) && O.minkowski_mat[2] > abs_tol
    A = deepcopy(O.minkowski_mat[1])
  else
    T = Vector{Vector{arb}}(degree(O))
    B = O.basis_nf
    for i in 1:degree(O)
      T[i] = minkowski_map(B[i], abs_tol)
    end
    p = maximum([ prec(parent(T[i][j])) for i in 1:degree(O), j in 1:degree(O) ])
    M = ArbMatSpace(ArbField(p, false), degree(O), degree(O), false)()
    for i in 1:degree(O)
      for j in 1:degree(O)
        M[i, j] = T[i][j]
      end
    end
    O.minkowski_mat = (M, abs_tol)
    A = deepcopy(M)
  end
  return A
end

################################################################################
#
#  Inclusion of number field elements
#
################################################################################

# Check if a number field element is contained in O
# In this case, the second return value is the coefficient vector with respect
# to the basis of O
function _check_elem_in_order(a::T, O::NfAbsOrd{S, T},
                              short::Type{Val{U}} = Val{false}) where {S, T, U}
  assure_has_basis_mat_inv(O)
  t = O.tcontain
  elem_to_mat_row!(t.num, 1, t.den, a)
  t = mul!(t, t, O.basis_mat_inv)
  if short == Val{true}
    return isone(t.den)
  elseif short == Val{false}
    if !isone(t.den)
      return false, Vector{fmpz}()
    else
      v = Vector{fmpz}(degree(O))
      for i in 1:degree(O)
        v[i] = deepcopy(t.num[1, i])
      end
      return true, v
    end
  end
end

doc"""
    in(a::nf_elem, O::NfOrd) -> Bool

Checks whether $a$ lies in $\mathcal O$.
"""
function in(a::nf_elem, O::NfOrd)
  return _check_elem_in_order(a, O, Val{true})
end

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

doc"""
    denominator(a::nf_elem, O::NfOrd) -> fmpz

Returns the smallest positive integer $k$ such that $k \cdot a$ is contained in
$\mathcal O$.
"""
function denominator(a::nf_elem, O::NfOrd)
  assure_has_basis_mat_inv(O)
  M = O.tcontain
  elem_to_mat_row!(M.num, 1, M.den, a)
  M = mul!(M, M, O.basis_mat_inv)
  return deepcopy(M.den)
end

function integral_split(a::nf_elem, O::NfOrd)
  assure_has_basis_mat_inv(O)
  M = O.tcontain
  elem_to_mat_row!(M.num, 1, M.den, a)
  M = mul!(M, M, O.basis_mat_inv)
  return M.den, M.num
end

##################################3#############################################
#
#  Norm change constant
#
################################################################################

# For x = \sum_i x_i omega_i let |x|_1 = \sqrt(x_1^2 + ... + x_d^2).
# And let |x|_2 = sqrt(T_2(x))
# Then there exist c1, c2 such that
# |x|_2^2 <= c1 |x|_2^2, |x|_1^2 <= c2 |x|_1^2
# A suitable pair (c1, c2) can be determined using the Minkowski map/matrix
#
# Reference
# Fieker, Friedrichs
# On Reconstruction of Algebraic Numbers
# (in particular p. 288)
#
# TODO: Make this rigorous using interval arithmetic. The only problem is that
#       the characteristic polynomial might not be squarefree.
doc"""
    norm_change_const(O::NfOrd) -> (Float64, Float64)

Returns $(c_1, c_2) \in \mathbf R_{>0}^2$ such that for all
$x = \sum_{i=1}^d x_i \omega_i \in \mathcal O$ we have
$T_2(x) \leq c_1 \cdot \sum_i^d x_i^2$
and
$\sum_i^d x_i^2 \leq c_2 \cdot T_2(x)$,
where $(\omega_i)_i$ is the $\mathbf Z$-basis of $\mathcal O$.
"""
function norm_change_const(O::NfOrd)
  if O.norm_change_const[1] > 0
    return O.norm_change_const
  else
    d = degree(O)
    M = minkowski_mat(O, 64)
    # I don't think we have to swap rows,
    # since permutation matrices are orthogonal
    #r1, r2 = signature(O)
    #for i in 2:2:r2
    #  swap_rows!(M, r1 + i, r1 + 2*r2 - i + 1)
    #end

    M = M*M'

    N = Symmetric([ Float64(M[i, j]) for i in 1:rows(M), j in 1:cols(M) ])
    #forcing N to really be Symmetric helps julia - aparently
    r = sort(eigvals(N))
    if !(r[1] > 0)
      # more complicated methods are called for...
      m = ceil(Int, log(d)/log(2))
      m += m%2
      @assert iseven(m)
      l_max = root(trace(M^m), m) #an upper bound within a factor of 2
                                    #according to a paper by Victor Pan
                                    #https://doi.org/10.1016/0898-1221(90)90236-D
                                    #formula (1) and discussion
      pr = 128
      l_min = l_max
      if isodd(d) d+=1; end
      while true
        try
          M = inv(M)
          l_min = root(trace(M^d), d) #as above...
          if isfinite(l_min)
            z = (Float64(l_max), Float64(l_min))
            O.norm_change_const = z
            return z
          end
          M = minkowski_mat(O, pr)
          pr *= 2
        catch e  # should verify the correct error
          M = minkowski_mat(O, pr)
          pr *= 2
        end
      end
    end

    @assert r[1]>0
#    N = transpose(M)*M
#    N = MatrixSpace(AcbField(prec(base_ring(N))), rows(N), cols(N))(N)
#    chi = charpoly(PolynomialRing(base_ring(N), "x")[1], N)
#    return chi
#    r = roots(chi)
#    # I want upper bound for the largest and lower bound for the smallest root
#
#    tm = arf_struct(0, 0, 0, 0)
#    ccall((:arf_init, :libarb), Void, (Ptr{arf_struct}, ), &tm)
#    ccall((:arb_get_abs_ubound_arf, :libarb), Void, (Ptr{arf_struct}, Ptr{arb}), &tm, &real(r[end]))
#    # 3 is round to infinity
#    c1 = ccall((:arf_get_d, :libarb), Cdouble, (Ptr{arf_struct}, Cint), &tm, 3)
#
#    ccall((:arb_get_abs_ubound_arf, :libarb), Void, (Ptr{arf_struct}, Ptr{arb}), &tm, &(inv(real(r[1]))))
#    c2 = ccall((:arf_get_d, :libarb), Cdouble, (Ptr{arf_struct}, Cint), &tm, 3)
#
#    ccall((:arf_clear, :libarb), Void, (Ptr{arf_struct}, ), &tm)
#
#    z = (c1, c2)
    z = (r[end], inv(r[1]))
    O.norm_change_const = z
    return z
  end
end

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

doc"""
    Order(B::Array{nf_elem, 1}, check::Bool = true) -> NfOrd

Returns the order with $\mathbf Z$-basis $B$. If `check` is set, it is checked
whether $B$ defines an order.
"""
function Order(::S, a::Array{T, 1}, check::Bool = true,
               cache::Bool = true) where {S <: Union{AnticNumberField, NfAbsNS}, T}
  K = parent(a[1])
  if check
    b, bmat, bmat_inv, _ = defines_order(K, a)
    if !b
      error("The elements do not define an order")
    else
      return NfAbsOrd(K, bmat, bmat_inv, deepcopy(a), cache)
    end
  else
    return NfAbsOrd(deepcopy(a), cache)
  end
end

function Order(K, a::Vector, check::Bool = true,
               cache::Bool = true)
  local b
  try
    b = map(K, a)
  catch
    error("Cannot coerce elements from array into the number field")
  end
  return Order(K, b, check, cache)
end

doc"""
    Order(K::AnticNumberField, A::FakeFmpqMat, check::Bool = true) -> NfOrd

Returns the order which has basis matrix $A$ with respect to the power basis
of $K$. If `check` is set, it is checked whether $A$ defines an order.
"""
function Order(K::S, a::FakeFmpqMat, check::Bool = true,
               cache::Bool = false) where {S}
  if check
    b, ainv, d = defines_order(K, a)
    if !b
      error("The basis matrix does not define an order")
    else
      return NfAbsOrd(K, deepcopy(a), ainv, d, cache)
    end
  else
    return NfAbsOrd(K, deepcopy(a), cache)
  end
end

doc"""
    Order(K::AnticNumberField, A::fmpq_mat, check::Bool = true) -> NfOrd

Returns the order which has basis matrix $A$ with respect to the power basis
of $K$. If `check` is set, it is checked whether $A$ defines an order.
"""
function Order(K::S, a::fmpq_mat, check::Bool = true,
               cache::Bool = true) where {S}
  return Order(K, FakeFmpqMat(a), cache)
end

doc"""
    Order(K::AnticNumberField, A::fmpz_mat, check::Bool = true) -> NfOrd

Returns the order which has basis matrix $A$ with respect to the power basis
of $K$. If `check` is set, it is checked whether $A$ defines an order.
"""
function Order(K::S, a::fmpz_mat, check::Bool = true,
               cache::Bool = true) where {S}
  return Order(K, FakeFmpqMat(a), check, cache)
end

doc"""
    Order(A::NfOrdFracIdl) -> NfOrd

Returns the fractional ideal $A$ as an order of the ambient number field.
"""
function Order(a::NfOrdFracIdl, check::Bool = true, cache::Bool = true)
  return Order(nf(order(a)), basis_mat(a)*basis_mat(order(a)), check, cache)
end

doc"""
    EquationOrder(K::AnticNumberField) -> NfOrd

Returns the equation order of the number field $K$.
"""
function EquationOrder(K::T, cache::Bool = false) where {T}
  M = FakeFmpqMat(identity_matrix(FlintZZ, degree(K)))
  Minv = FakeFmpqMat(identity_matrix(FlintZZ, degree(K)))
  z = NfAbsOrd{T, elem_type(T)}(K, M, Minv, basis(K), cache)
  z.isequation_order = true
  return z
end

doc"""
   equation_order(M::NfOrd) -> NfOrd
> The equation order of then number field.
"""
equation_order(M::NfAbsOrd) = equation_order(nf(M))

function _order(K::S, elt::Array{T, 1}) where {S, T}
  o = one(K)

  if !(o in elt)
    push!(elt, o)
  end

  n = degree(K)

  closed = false

  dold = fmpq(0)

  # Since 1 is in elt, prods will contain all elements

  while !closed
    prods = T[]
    for u in elt
      for v in elt
        push!(prods, u * v)
      end
    end

    B = basis_mat(prods)
    H = sub(hnf(B), (rows(B) - degree(K) + 1):rows(B), 1:degree(K))

    # TODO: Just take the product of the diagonal
    dd = H.num[1, 1]
    for i in 2:degree(K)
      dd *= H.num[i, i]
    end
    d = dd//H.den^degree(K)

    #d = prod(H.num[i,i] for i=1:degree(K))//(H.den^degree(K))
    #d = det(H)

    if iszero(d)
      error("Elements do not define a module of full rank")
    end

    if dold == d
      closed = true
    else
      dold = d
      elt = T[]
      for i in 1:n
        push!(elt, elem_from_mat_row(K, H.num, i, H.den))
      end
    end
  end

  # Make an explicit check
  @hassert :NfOrd 1 defines_order(K, elt)[1]
  return Order(K, elt, false)
end

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

# This is the function which is used in dictionaries
function Base.isequal(R::NfOrd, S::NfOrd)
  return R === S
end

function ==(R::NfAbsOrd, S::NfAbsOrd)
  nf(R) != nf(S) && return false
  assure_has_basis_mat(R)
  assure_has_basis_mat(S)
  return R.basis_mat == S.basis_mat
end

function ==(R::NfAbsOrdSet, S::NfAbsOrdSet)
  return R.nf === S.nf
end

################################################################################
#
#  Trace matrix
#
################################################################################

doc"""
    trace_matrix(O::NfOrd) -> fmpz_mat

Returns the trace matrix of `\mathcal O`, that is, the matrix
$(\operatorname{tr}_{K/\mathbf Q}(b_i \cdot b_j))_{1 \leq i, j \leq d}$.
"""
function trace_matrix(O::NfAbsOrd)
  if isdefined(O, :trace_mat)
    return deepcopy(O.trace_mat)
  end
  K = nf(O)
  b = O.basis_nf
  n = degree(K)
  g = zero_matrix(FlintZZ, n, n)
  for i=1:n
    t = trace(b[i]^2)
    @assert isinteger(t)
    g[i, i] = numerator(t)
    for j in (i + 1):n
      t = trace(b[i]*b[j])
      @assert isinteger(t)
      g[i, j] = numerator(t)
      g[j, i] = numerator(t)
    end
  end
  O.trace_mat = g
  return deepcopy(g)
end

################################################################################
#
#  Addition of orders
#
################################################################################

doc"""
    +(R::NfOrd, S::NfOrd) -> NfOrd

Given two orders $R$, $S$ of $K$, this function returns the smallest order
containing both $R$ and $S$. It is assumed that $R$, $S$ contain the ambient
equation order and have coprime index.
"""
function +(a::NfAbsOrd, b::NfAbsOrd)
  nf(a) != nf(b) && error("Orders must have same ambient number field")
  if isone(gcd(index(a), index(b)))
    aB = basis_mat(a, Val{false})
    bB = basis_mat(b, Val{false})
    d = degree(a)
    m = zero_matrix(FlintZZ, 2*d, d)
    for i=1:d
      for j=1:d
        m[i,j]=bB.den*aB.num[i,j]
        m[i+d,j]= aB.den*bB.num[i,j]
      end
    end
    mat=_hnf_modular_eldiv(m, bB.den*aB.den, :lowerleft)
    c = sub(mat, d + 1:2*d, 1:d)
    O = Order(nf(a), FakeFmpqMat(c, aB.den*bB.den), false)
    O.primesofmaximality = unique(vcat(a.primesofmaximality, b.primesofmaximality))
    O.disc=gcd(discriminant(a), discriminant(b))
    if a.disc<0
      O.disc=-O.disc
    end
    return O
  else
    return _order(nf(a), vcat(a.basis_nf, b.basis_nf))
  end
end

################################################################################
#
#  p-overorder
#
################################################################################

function _poverorder(O::NfAbsOrd, p::fmpz)
  I= pradical(O, p)
  if isdefined(I, :princ_gen) && I.princ_gen==p
    return deepcopy(O)
  end
  R = ring_of_multipliers(I)
  return R
end

function _poverorder(O::NfAbsOrd, p::Integer)
  return _poverorder(O, fmpz(p))
end

doc"""
    poverorder(O::NfOrd, p::fmpz) -> NfOrd
    poverorder(O::NfOrd, p::Integer) -> NfOrd

This function tries to find an order that is locally larger than $\mathcal O$
at the prime $p$: If $p$ divides the index $[ \mathcal O_K : \mathcal O]$,
this function will return an order $R$ such that
$v_p([ \mathcal O_K : R]) < v_p([ \mathcal O_K : \mathcal O])$. Otherwise
$\mathcal O$ is returned.
"""
function poverorder(O::NfAbsOrd, p::fmpz)
  if isequation_order(O) && issimple(nf(O))
    return dedekind_poverorder(O, p)
  else
    return _poverorder(O, p)
  end
end

function poverorder(O::NfAbsOrd, p::Integer)
  return poverorder(O, fmpz(p))
end

################################################################################
#
#  p-maximal overorder
#
################################################################################

doc"""
    pmaximal_overorder(O::NfOrd, p::fmpz) -> NfOrd
    pmaximal_overorder(O::NfOrd, p::Integer) -> NfOrd

This function finds a $p$-maximal order $R$ containing $\mathcal O$. That is,
the index $[ \mathcal O_K : R]$ is not dividible by $p$.
"""
function pmaximal_overorder(O::NfAbsOrd, p::fmpz)
  @vprint :NfOrd 1 "computing p-maximal overorder for $p ... \n"
  if rem(discriminant(O), p^2) != 0
    push!(O.primesofmaximality, p)
    return O
  end

  d = discriminant(O)
  @vprint :NfOrd 1 "extending the order at $p for the first time ... \n"
  OO = poverorder(O, p)
  dd = discriminant(OO)
  v=valuation(dd,p)
  if v==0 || v==1
    push!(OO.primesofmaximality, p)
    return OO
  end
  i = 1
  while d != dd
    i += 1
    @vprint :NfOrd 1 "extending the order at $p for the $(i)th time ... \n"
    d = dd
    OO = poverorder(OO, p)
    dd = discriminant(OO)
    v = valuation(dd,p)
    if v==0 || v==1
      break
    end
  end
  push!(OO.primesofmaximality, p)
  return OO
end

function pmaximal_overorder(O::NfAbsOrd, p::Integer)
  return pmaximal_overorder(O, fmpz(p))
end

function MaximalOrder(O::NfOrd, primes::Array{fmpz, 1})
  OO = deepcopy(O)
  disc = abs(discriminant(O))
  for i in 1:length(primes)
    p = primes[i]
    (j, disc) = remove(disc, p)
    if j == 1
      continue
    end
    @vprint :NfOrd 1 "Computing p-maximal overorder for $p ..."
    O1 = pmaximal_overorder(O, p)
    if valuation(discriminant(O1), p)< valuation(discriminant(OO),p)
      OO += O1
    end
    if !(p in OO.primesofmaximality)
      push!(OO.primesofmaximality,p)
    end
    @vprint :NfOrd 1 "done\n"
  end
  return OO
end

doc"""
***
    MaximalOrder(K::NfOrd) -> NfOrd

Returns the maximal order of $K$.
"""
function MaximalOrder(O::NfAbsOrd)
  OO = deepcopy(O)
  @vtime :NfOrd fac = factor(Nemo.abs(discriminant(O)))
  for (p,j) in fac
    if j == 1
      continue
    end
    @vprint :NfOrd 1 "Computing p-maximal overorder for $p ..."
    O1 = pmaximal_overorder(O, p)
    if valuation(discriminant(O1), p)< valuation(discriminant(OO),p)
      OO += O1
    end
    @vprint :NfOrd 1 "done\n"
  end
  OO.ismaximal = 1
  return OO
end

###############################################################################
#
#  Some LLL-related functions
#
###############################################################################

# don't know what this is doing
#for totally real field, the T_2-Gram matrix is the trace matrix, hence exact.

function _lll_gram(M::NfOrd)
  K = nf(M)
  @assert istotally_real(K)
  g = trace_matrix(M)

  q,w = lll_gram_with_transform(g)
  On = NfOrd(K, FakeFmpqMat(w*basis_mat(M, Val{false}).num, denominator(basis_mat(M, Val{false}))))
  On.ismaximal = M.ismaximal
  return On
end

doc"""
    lll_basis(M::NfOrd) -> Array{nf_elem, 1}
> A basis for $m$ that is reduced using the LLL algorithm for the Minkowski metric.
"""
function lll_basis(M::NfOrd)
  I = ideal(M, parent(basis_mat(M, Val{false}).num)(1))
  return lll_basis(I)
end

doc"""
    lll(M::NfOrd) -> NfOrd
> The same order, but with the basis now being LLL reduced wrt. the Minkowski metric.
"""
function lll(M::NfOrd)
  K = nf(M)

  if istotally_real(K)
    return _lll_gram(M)
  end

  I = ideal(M, 1)

  prec = 100
  while true
    try
      q,w = lll(I, parent(basis_mat(M).num)(0), prec = prec)
      On = NfOrd(K, FakeFmpqMat(w*basis_mat(M, Val{false}).num, denominator(basis_mat(M, Val{false}))))
      On.ismaximal = M.ismaximal
      return On
    catch e
      if isa(e, LowPrecisionLLL)
        prec = Int(round(prec*1.2))
        if prec>1000
          error("precision too large in LLL");
        end
        continue;
      else
        rethrow(e)
      end
    end
  end
end

################################################################################
#
#  Constructors for users
#
################################################################################

doc"""
***
    MaximalOrder(K::AnticNumberField) -> NfOrd

Returns the maximal order of $K$.

# Example

```julia-repl
julia> Qx, xx = FlintQQ["x"];
julia> K, a = NumberField(x^3 + 2, "a");
julia> O = MaximalOrder(K);
```
"""
function MaximalOrder(K::AnticNumberField, cache::Bool = false)
  O = EquationOrder(K)
  @vprint :NfOrd 1 "Computing the maximal order ...\n"
  O = MaximalOrder(O)
  @vprint :NfOrd 1 "... done\n"
  M = NfOrd(K, basis_mat(O), cache)
  M.ismaximal = 1
  return M
end

doc"""
***
    MaximalOrder(K::AnticNumberField, primes::Array{fmpz, 1}) -> NfOrd

Assuming that ``primes`` contains all the prime numbers at which the equation
order $\mathbf{Z}[\alpha]$ of $K = \mathbf{Q}(\alpha)$ is not maximal,
this function returns the maximal order of $K$.
"""
function MaximalOrder(K::AnticNumberField, primes::Array{fmpz, 1})
  O = EquationOrder(K)
  @vprint :NfOrd 1 "Computing the maximal order ...\n"
  O = MaximalOrder(O, primes)
  @vprint :NfOrd 1 "... done\n"
  return NfOrd(K, basis_mat(O))
end

doc"""
***
    maximal_order(K::AnticNumberField) -> NfOrd
    ring_of_integers(K::AnticNumberField) -> NfOrd

Returns the maximal order of $K$.
"""
function maximal_order(K::AnticNumberField)
  try
    c = _get_maximal_order_of_nf(K)::NfOrd
    return c
  catch e
    if !isa(e, AccessorNotSetError)
      rethrow(e)
    end
    #O = MaximalOrder(K)::NfOrd
    O = new_maximal_order(K)::NfOrd
    O.ismaximal=1
    _set_maximal_order_of_nf(K, O)
    return O
  end
end

doc"""
***
    maximal_order(K::AnticNumberField, primes::Array{fmpz, 1}) -> NfOrd
    maximal_order(K::AnticNumberField, primes::Array{Integer, 1}) -> NfOrd

Assuming that ``primes`` contains all the prime numbers at which the equation
order $\mathbf{Z}[\alpha]$ of $K = \mathbf{Q}(\alpha)$ is not maximal
(e.g. ``primes`` may contain all prime divisors of the discriminant of
$\mathbf Z[\alpha]$), this function returns the maximal order of $K$.
"""
function maximal_order(K::AnticNumberField, primes::Array{fmpz, 1})
  O = MaximalOrder(K, primes)
  return O
end

maximal_order(K::AnticNumberField, primes::Array{T, 1}) where {T} =
  maximal_order(K, map(FlintZZ, primes))

doc"""
***
    ring_of_integers(K::AnticNumberField, primes::Array{fmpz, 1}) -> NfOrd
    ring_of_integers(K::AnticNumberField, primes::Array{Integer, 1}) -> NfOrd

Assuming that ``primes`` contains all the prime numbers at which the equation
order $\mathbf{Z}[\alpha]$ of $K = \mathbf{Q}(\alpha)$ is not maximal,
this function returns the maximal order of $K$.
"""
ring_of_integers(x...) = maximal_order(x...)

################################################################################
#
#  Check if something defines an order
#
################################################################################

# This is not optimizied for performance.
# If false, then this returns (false, garbage, garbage).
# If true, then this return (true, basis_mat, basis_mat_inv).
function defines_order(K::S, x::FakeFmpqMat) where {S}
  if rows(x) != degree(K) || cols(x) != degree(K)
    return false, x, Vector{elem_type(K)}()
  end
  local xinv
  try
    xinv = inv(x)
  catch
    return false, x, Vector{elem_type(K)}()
  end
  n = degree(K)
  B_K = basis(K)
  d = Vector{elem_type(K)}(n)
  # Construct the basis elements from the basis matrix
  for i in 1:n
    d[i] = elem_from_mat_row(K, x.num, i, x.den)
  end

  # Check if Z-module spanned by x is closed under multiplcation
  l = Vector{elem_type(K)}(n)
  for i in 1:degree(K)
    for j in 1:degree(K)
      l[j] = d[i]*d[j]
    end
    Ml = basis_mat(l)
    if !isone((Ml * xinv).den)
      return false, x, Vector{elem_type(K)}()
    end
  end
  # Check if 1 is contained in the Z-module
  Ml = basis_mat([one(K)])
  if !isone((Ml * xinv).den)
    return false, x, Vector{elem_type(K)}()
  end
  return true, xinv, d
end

function defines_order(K::S, A::Vector{T}) where {S, T}
  if length(A) != degree(K)
    return false, FakeFmpqMat(), FakeFmpqMat(), A
  else
    B = basis_mat(A)
    b, Binv, _ = defines_order(K, B)
    return b, B, Binv, A
  end
end

################################################################################
#
#  Different
#
################################################################################

doc"""
    different(x::NfOrdElem) -> NfOrdElem

> The different of $x$, ie. $0$ if x is not a primitive element, or
> $f'(x)$ for $f$ the minimal polynomial of $x$ otherwise.
"""
function different(x::NfOrdElem)
  if iszero(x)
    return x
  end
  f = minpoly(x)
  if degree(f) < degree(parent(x))
    return parent(x)(0)
  end
  return derivative(f)(x)
end

doc"""
    different(R::NfOrd) -> NfOrdIdl

> The differnt ideal of $R$, that is, the ideal generated by all differents
> of elements in $R$.
> For the maximal order, this is also the inverse ideal of the co-different.
"""
function different(R::NfOrd)
  D = ideal(R, different(R(gen(nf(R)))))
  d = abs(discriminant(R))
  while norm(D) != d
    #@show D, norm(D), d
    x = rand(R, -10:10)
    y = different(x)
    if !iszero(y)
      D += ideal(R, y)
    end
  end
  return D
end

doc"""
    codifferent(R::NfOrd) -> NfOrdIdeal
> The codiffernt ideal of $R$, ie. the trace-dual of $R$.
"""
function codifferent(R::NfOrd)
  t = trace_matrix(R)
  ti, d = pseudo_inv(t)
  return ideal(R, ti, true)//d
end

###############################################################################
#
#  Buchmann-Lenstra for the computation of maximal orders
#
###############################################################################

function new_maximal_order(K::AnticNumberField)
  O=EquationOrder(K)
  if degree(K)==1
    O.ismaximal=1
    return O
  end
  ds=discriminant(O)
  #First, factorization of the discriminant given by the snf of the trace matrix
  M = trace_matrix(O)
  l = coprime_base(_el_divs(M,ds))
  @vprint :NfOrd 1 "Factors of the discriminant: $l\n "
  l1=fmpz[]
  OO=O
  @vprint :NfOrd 1 "Trial division of the discriminant\n "
  for d in l
    fac = factor_trial_range(d)[1]
    rem= d
    for (p,v) in fac
      rem=divexact(rem, p^v)
    end
    @vprint :NfOrd 1 "Computing the maximal order at $(collect(keys(fac)))\n "
    O1=MaximalOrder(O, collect(keys(fac)))
    OO+=O1
    if abs(rem)!=1
      push!(l1,abs(rem))
    end
  end
  if isempty(l1)
    OO.ismaximal=1
    return OO
  end
  for i=1:length(l1)
    a,b = ispower(l1[i])
    if a>1
      l1[i]=b
    end
  end
  #=
  O1, l=DedekindComposite(O, deepcopy(l1))
  if discriminant(O1)!=discriminant(O)
    OO+=O1
    push!(l, discriminant(O1))
  end
  append!(l1,l)
  l1=coprime_base(l1)
  =#
  O1, Q=_TameOverorderBL(OO, l1)
  if !isempty(Q)
    @vprint :NfOrd 1 "I have to factor $Q\n "
    for el in Q
      d=factor(el).fac
      O1+=MaximalOrder(O1, collect(keys(d)))
    end
  end
  O1.ismaximal=1
  return O1

end

function DedekindComposite(O::NfOrd, l::Array{fmpz,1})

  @vprint :NfOrd 1 "Dedekind criterion with $rem \n"
  O1=O
  l1=fmpz[]
  while true
    @vprint :NfOrd 1 "$l"
    if isempty(l)
      break
    end
    div, f, OO = dedekind_test_composite(O,l[1])
    if div!=1
      @show "Divisor found"
      push!(l1, div)
      push!(l,div)
      l=coprime_base(l)
      continue
    end
    if f
      @show "Can't say"
      filter!(x-> x!=l[1], l)
      continue
    else
      @show "Enlarge!"
      @show index(O1)
      O1+=OO
      filter!(x-> x!=l[1], l)
      continue
    end
  end
  return O, l1
end

function _TameOverorderBL(O::NfOrd, lp::Array{fmpz,1})

  OO=O
  M=coprime_base(lp)
  Q=fmpz[]
  while !isempty(M)
    @vprint :NfOrd 1 M
    q=M[1]
    if isprime(q)
      OO1=pmaximal_overorder(O, q)
      if valuation(discriminant(OO1), q)<valuation(discriminant(OO), q)
        OO+=OO1
      end
      filter!(x-> x!=q, M)
      continue
    end
    OO, q1= _cycleBL(OO,q)
    if q1==1
      filter!(x->x!=q, M)
    elseif q1==q
      push!(Q,q)
      filter!(x-> x!=q, M)
    else
      push!(M, q1)
      push!(M, divexact(q,q1))
      M=coprime_base(M)
      continue
    end
  end
  if isempty(Q)
    OO.ismaximal=1
  end
  return OO, Q

end

function _qradical(O::NfOrd, q::fmpz)

  R=ResidueRing(FlintZZ, q, cached=false)
  #First, we compute the q-radical as the kernel of the trace matrix mod q.
  #By theory, this is free if q is prime; if I get a non free module, I have found a factor of q.
  @vprint :NfOrd 1 "radical computation\n "
  Mat=MatrixSpace(R, degree(O), degree(O), false)(trace_matrix(O))
  M=nullspace(Mat)[1]
  if iszero(M)
    @vprint "The radical is equal to the ideal generated by q"
    return fmpz(1), ideal(O,q)
  end
  if cols(M)!=0
    M=howell_form(transpose(M))
    for i=1:rows(M)
      if iszero_row(M,i)
        break
      end
      j=i
      while iszero(M[i,j])
        j+=1
      end
      if !isone(M[i,j])
        @vprint :NfOrd 1 "Split: $(M[i,j])"
        return lift(M[i,j]), NfOrdIdl()
      end
    end
  end
  # Now, we have the radical.
  # We want to compute the ring of multipliers.
  # So we need to lift the matrix.
  @vprint :NfOrd 1 "Computing hnf of basis matrix \n "
  MatIdeal=zero_matrix(FlintZZ, rows(M)+degree(O), cols(M))
  for i=1:rows(M)
    for j=1:degree(O)
      MatIdeal[i,j]=FlintZZ(M[i,j].data)
    end
  end
  for i=1:degree(O)
    MatIdeal[i+rows(M), i]=q
  end
  gens=[O(q)]
  for i=1:rows(M)
    if !iszero_row(MatIdeal, i)
      push!(gens, elem_from_mat_row(O, MatIdeal, i))
    end
  end
  M2=sub(Hecke._hnf_modular_eldiv(MatIdeal, fmpz(q),  :lowerleft), rows(M)+1:degree(O)+rows(M), 1:degree(O))
  I=NfOrdIdl(O, M2)
  I.gens=gens
  return fmpz(1), I
end

function _cycleBL(O::NfOrd, q::fmpz)

  q1,I=_qradical(O,q)
  if q1!=1
    return O,q1
  elseif isdefined(I, :princ_gens) && q==I.princ_gens
    return O, fmpz(1)
  end
  @vprint :NfOrd 1 "ring of multipliers\n"
  O1=ring_of_multipliers(I)
  @vprint :NfOrd 1 "ring of multipliers computed\n"
  if discriminant(O1)!=discriminant(O)
    if gcd(discriminant(O1), q)== 1
      return O1, fmpz(1)
    else
      return _cycleBL(O1,q)
    end
  end
  @vprint :NfOrd 1 "couldn't enlarge\n"
  # (I:I)=OO. Now, we want to prove tameness (or get a factor of q)
  # We have to test that (OO:a)/B is a free Z/qZ module.
  inva=colon(ideal(O,1),I, true)
  M1=basis_mat_inv(inva)
  @assert M1.den==1
  G1=_el_divs(M1.num,q)
  for i=1:length(G1)
    q1=gcd(q,G1[i])
    if q1!=q && q1!=1
      @vprint :NfOrd 1 "Found the factor $q1"
      return O,fmpz(q1)
    end
  end
  @vprint :NfOrd 1 "...is free\n"
  h=2
  ideals=Array{NfOrdIdl,1}(3)
  ideals[1]=I
  ideals[2]=I^2
  ideals[3]=I^3
  while true
    if h>degree(O)
      error("Not found!")
    end
    I1=(ideals[1]+ideal(O,q))*(ideals[3]+ideal(O,q))
    I2=(ideals[2]+ideal(O,q))^2
    M2=basis_mat(I2, Val{false})*basis_mat_inv(I1, Val{false})
    @assert M2.den==1
    G2=_el_divs(M2.num,q)
    if isempty(G2)
      h+=1
      ideals[1]=ideals[2]
      ideals[2]=ideals[3]
      ideals[3]=ideals[2]*I
      continue
    end
    for i=1:length(G2)
      q1=gcd(q,G2[i])
      if q1!=q && q1!=1
        return O, q1
      end
    end
    break
  end
  f,r = ispower(q,h)
  if f
    return O, r
  else
    return O, q
  end

end

function _el_divs(M::fmpz_mat, d::fmpz)
  M1=_hnf_modular_eldiv(M, d)
  while !isdiag(M1)
    M1= M1'
    hnf_modular_eldiv!(M1, d)
  end
  l=fmpz[]
  for j=1:rows(M1)
    if M1[j,j]!=1
      push!(l, M1[j,j])
    end
  end
  return l
end

function TameOverorderBL(O::NfOrd, lp::Array{fmpz,1}=fmpz[])

  # First, we hope that we can get a factorization of the discriminant by computing
  # the structure of the group OO^*/OO
  OO=O
  list=append!(_el_divs(trace_matrix(OO)),primes_up_to(degree(O)))
  l=coprime_base(list)
  #Some trivial things, maybe useless
  for i=1:length(l)
    a,b=ispower(l[i])
    if a>1
      l[i]=b
    end
    if isprime(l[i])
      @vprint :NfOrd 1 "pmaximal order at $(l[i])\n"
      OO1=pmaximal_overorder(O, l[i])
      if valuation(discriminant(OO1), l[i])<valuation(discriminant(OO), l[i])
        OO+=OO1
      end
      l[i]=0
    end
  end
  push!(l, discriminant(OO))
  append!(l,lp)
  filter!(x-> x!=0, l)
  for i=1:length(l)
    l[i]=abs(l[i])
  end
  M=coprime_base(l)
  Q=fmpz[]
  while !isempty(M)
    @vprint :NfOrd 1 M
    q=M[1]
    if isprime(q)
      OO1=pmaximal_overorder(O, q)
      if valuation(discriminant(OO1), q)<valuation(discriminant(OO), q)
        OO+=OO1
      end
      filter!(x-> x!=q, M)
      continue
    end
    OO, q1= _cycleBL(OO,q)
    if q1==1
      filter!(x->x!=q, M)
    elseif q1==q
      push!(Q,q)
      filter!(x-> x!=q, M)
    else
      push!(M, q1)
      push!(M, divexact(q,q1))
      M=coprime_base(M)
      continue
    end
  end
  if isempty(Q)
    OO.ismaximal=1
  end
  return OO, Q

end


text:
function get_data()
  d = []
  open("graph.txt", "r") do f
    #conditions = Array(Any, 100)
    for l in readlines(f)
      code = [strip(x) for x in split(l, ",")]
      code2 = []
      for y in code[2:end]
        zs = split(y)
        if length(zs) > 1
          push!(code2, [zs[1], parse(Float64, zs[2])])
        end
        #code2 = [split(y) for y in code[2:end]]
      end
      push!(d, [code[1], code2])
    end
  end
  return d
end

get_nums_from_str(str) = parse(Int, match(r"([0-9]+)", str).captures[1])

function data_rev!(data)
  for row in data
    row[1] = get_nums_from_str(row[1]) + 1
    #println(row[2])
    if length(row) != 1
      for el in row[2:end]
        el[1] = get_nums_from_str(el[1]) + 1
      end
    end
  end
end

#println(J0)

function set_C(data)
  #println(data)
  C = zeros(100, 100)
  for row in data
    a = row[1]
    if length(row) != 1
      for list in row[2:end]
        b = list[1]
        c = list[2]
        #println(a, b)
        C[a, b] = c
        C[b, a] = c
      end
    end
  end
  return C
end

function infplus(c, j)
  if j == 5000
    return 5000
  else
    return c + j
  end
end

function get_next_J(J, C)
  new_J = deepcopy(J)
  for i in 1:length(J)
    available_js = filter(j -> C[i, j] != 0, 1:length(J))
    #println(available_js)
    #println([(C[i, j] + J[j]) for j in available_js]...)
    #println([C[i, j] for j in available_js])
    #println([J[j] for j in available_js])
    #println([C[i, j] + J[j] for j in available_js])
    new_J[i] = minimum([infplus(C[i, j], J[j]) for j in available_js])
    #println(minimum([C[i, j] + J[j] for j in available_js]))
  end
  println(new_J)
  #println(new_J-J)
  return new_J
end
#println("hi")

data = get_data()
data_rev!(data)
#println(data)
J0 = ones(Float64, 100) .* 5000
J0[1] = 0#%%%%%%%%%%%%%%%%%%%%%%%%%%%
C = set_C(data)
#println(C)

J2, J = get_next_J(J0, C), J0
i = 1
while J != J2
  #println(J2-J)
  J2, J = get_next_J(J2, C), J2
  i += 1
  #println(i)
  if i == 50
    break
  end
end

println("finished by step $i")
#function data2matrix(data)
println(J2)


text:
# This file was generated, do not modify it. # hide
ŷ = predict(mach_pipe, rows=validation);
print(
    "Measurements:\n",
    "  brier loss: ", brier_loss(ŷ, y[validation]) |> mean, "\n",
    "  auc:        ", auc(ŷ, y[validation]),                "\n",
    "  accuracy:   ", accuracy(mode.(ŷ), y[validation])
)

max_stars_repo_name: alan-turing-institute/MLJTutorials

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
include("PrintMatrix.jl")
include("PsPot_HGH.jl")

include("../common/PWGrid_v02.jl")

import PyPlot
const plt = PyPlot

function test_main()
  #psp1 = PsPot_HGH(1, "H" , "LDA_HGH/1h.1.hgh")
  #psp2 = PsPot_HGH(2, "Ti", "LDA_HGH/22ti.12.hgh")
  #psp = PsPot_HGH(3, "Sm", "LDA_HGH/62sm.16.hgh")
  #psp4 = PsPot_HGH(4, "Pt", "LDA_HGH/78pt.18.hgh")
  psp = PsPot_HGH(1, "Ge", "LDA_HGH/32ge.4.hgh")

  info_PsPot_HGH(psp)

  #Ns = [32, 32, 32]
  #LatVecs = 18.0*diagm(ones(3))
  #pw = PWGrid( Ns, LatVecs )

  # use G-vectors for wavefunction
  #Gm = sqrt(pw.gvec.G2[pw.gvecw.idx_gw2r])

  #Vloc = eval_HGH_Vloc_G( psp, pw.gvec.G2, pw.Ω )
  #plt.clf()
  #plt.plot( sqrt(pw.gvec.G2[2:end]), Vloc[2:end], marker="o", label="Vloc" )
  #plt.grid()
  #plt.legend()
  #plt.xlim(0.0,2.0)
  #plt.savefig("Ge_Vloc.png", dpi=300)

  #Vprj_s = eval_HGH_proj_G( psp, 0, 1, Gm, pw.Ω )
  #Vprj_p = eval_HGH_proj_G( psp, 1, 1, Gm, pw.Ω )
  #Vprj_d = eval_HGH_proj_G( psp, 2, 1, Gm, pw.Ω )

  #plt.clf()
  #plt.plot( Gm, Vprj_s, marker="o", label="s-prj" )
  #plt.plot( Gm, Vprj_p, marker="o", label="p-prj" )
  #plt.plot( Gm, Vprj_d, marker="o", label="d-prj" )
  #plt.grid()
  #plt.legend()
  #plt.savefig("Ge_prj.png", dpi=300)

end

test_main()


text:
ccall((:openblas_get_num_threads64_, Base.libblas_name), Cint, ())

using LinearAlgebra
LinearAlgebra.BLAS.set_num_threads(4)

using DifferentialEquations
function rober(du,u,p,t)
  y₁,y₂,y₃ = u
  k₁,k₂,k₃ = p
  du[1] = -k₁*y₁+k₃*y₂*y₃
  du[2] =  k₁*y₁-k₂*y₂^2-k₃*y₂*y₃
  du[3] =  k₂*y₂^2
  nothing
end
prob = ODEProblem(rober,[1.0,0.0,0.0],(0.0,1e5),(0.04,3e7,1e4))
sol = solve(prob,Rosenbrock23())

using Plots
plot(sol, xscale=:log10, tspan=(1e-6, 1e5), layout=(3,1))

using BenchmarkTools
@btime solve(prob)

function rober_jac(J,u,p,t)
  y₁,y₂,y₃ = u
  k₁,k₂,k₃ = p
  J[1,1] = k₁ * -1
  J[2,1] = k₁
  J[3,1] = 0
  J[1,2] = y₃ * k₃
  J[2,2] = y₂ * k₂ * -2 + y₃ * k₃ * -1
  J[3,2] = y₂ * 2 * k₂
  J[1,3] = k₃ * y₂
  J[2,3] = k₃ * y₂ * -1
  J[3,3] = 0
  nothing
end
f = ODEFunction(rober, jac=rober_jac)
prob_jac = ODEProblem(f,[1.0,0.0,0.0],(0.0,1e5),(0.04,3e7,1e4))

@btime solve(prob_jac)

using ModelingToolkit
de = modelingtoolkitize(prob)
ModelingToolkit.generate_jacobian(de)[2] # Second is in-place

:((##MTIIPVar#376, u, p, t)->begin
          #= C:\Users\accou\.julia\packages\ModelingToolkit\czHtj\src\utils.jl:65 =#
          #= C:\Users\accou\.julia\packages\ModelingToolkit\czHtj\src\utils.jl:66 =#
          let (x₁, x₂, x₃, α₁, α₂, α₃) = (u[1], u[2], u[3], p[1], p[2], p[3])
              ##MTIIPVar#376[1] = α₁ * -1
              ##MTIIPVar#376[2] = α₁
              ##MTIIPVar#376[3] = 0
              ##MTIIPVar#376[4] = x₃ * α₃
              ##MTIIPVar#376[5] = x₂ * α₂ * -2 + x₃ * α₃ * -1
              ##MTIIPVar#376[6] = x₂ * 2 * α₂
              ##MTIIPVar#376[7] = α₃ * x₂
              ##MTIIPVar#376[8] = α₃ * x₂ * -1
              ##MTIIPVar#376[9] = 0
          end
          #= C:\Users\accou\.julia\packages\ModelingToolkit\czHtj\src\utils.jl:67 =#
          nothing
      end)

jac = eval(ModelingToolkit.generate_jacobian(de)[2])
f = ODEFunction(rober, jac=jac)
prob_jac = ODEProblem(f,[1.0,0.0,0.0],(0.0,1e5),(0.04,3e7,1e4))

I = [1,2,1,2,3,1,2]
J = [1,1,2,2,2,3,3]
using SparseArrays
jac_prototype = sparse(I,J,1.0)

f = ODEFunction(rober, jac=jac, jac_prototype=jac_prototype)
prob_jac = ODEProblem(f,[1.0,0.0,0.0],(0.0,1e5),(0.04,3e7,1e4))

const N = 32
const xyd_brusselator = range(0,stop=1,length=N)
brusselator_f(x, y, t) = (((x-0.3)^2 + (y-0.6)^2) <= 0.1^2) * (t >= 1.1) * 5.
limit(a, N) = a == N+1 ? 1 : a == 0 ? N : a
function brusselator_2d_loop(du, u, p, t)
  A, B, alpha, dx = p
  alpha = alpha/dx^2
  @inbounds for I in CartesianIndices((N, N))
    i, j = Tuple(I)
    x, y = xyd_brusselator[I[1]], xyd_brusselator[I[2]]
    ip1, im1, jp1, jm1 = limit(i+1, N), limit(i-1, N), limit(j+1, N), limit(j-1, N)
    du[i,j,1] = alpha*(u[im1,j,1] + u[ip1,j,1] + u[i,jp1,1] + u[i,jm1,1] - 4u[i,j,1]) +
                B + u[i,j,1]^2*u[i,j,2] - (A + 1)*u[i,j,1] + brusselator_f(x, y, t)
    du[i,j,2] = alpha*(u[im1,j,2] + u[ip1,j,2] + u[i,jp1,2] + u[i,jm1,2] - 4u[i,j,2]) +
                A*u[i,j,1] - u[i,j,1]^2*u[i,j,2]
    end
end
p = (3.4, 1., 10., step(xyd_brusselator))

using SparsityDetection, SparseArrays
input = rand(32,32,2)
output = similar(input)
sparsity_pattern = jacobian_sparsity(brusselator_2d_loop,output,input,p,0.0)
jac_sparsity = Float64.(sparse(sparsity_pattern))

using Plots
spy(jac_sparsity,markersize=1,colorbar=false,color=:deep)

f = ODEFunction(brusselator_2d_loop;jac_prototype=jac_sparsity)

function init_brusselator_2d(xyd)
  N = length(xyd)
  u = zeros(N, N, 2)
  for I in CartesianIndices((N, N))
    x = xyd[I[1]]
    y = xyd[I[2]]
    u[I,1] = 22*(y*(1-y))^(3/2)
    u[I,2] = 27*(x*(1-x))^(3/2)
  end
  u
end
u0 = init_brusselator_2d(xyd_brusselator)
prob_ode_brusselator_2d = ODEProblem(brusselator_2d_loop,
                                     u0,(0.,11.5),p)

prob_ode_brusselator_2d_sparse = ODEProblem(f,
                                     u0,(0.,11.5),p)

@btime solve(prob_ode_brusselator_2d,save_everystep=false)
@btime solve(prob_ode_brusselator_2d_sparse,save_everystep=false)

using SparseDiffTools
colorvec = matrix_colors(jac_sparsity)
@show maximum(colorvec)

f = ODEFunction(brusselator_2d_loop;jac_prototype=jac_sparsity,
                                    colorvec=colorvec)
prob_ode_brusselator_2d_sparse = ODEProblem(f,
                                     init_brusselator_2d(xyd_brusselator),
                                     (0.,11.5),p)
@btime solve(prob_ode_brusselator_2d_sparse,save_everystep=false)

@btime solve(prob_ode_brusselator_2d,TRBDF2(linsolve=LinSolveGMRES()),save_everystep=false)
@btime solve(prob_ode_brusselator_2d_sparse,TRBDF2(linsolve=LinSolveGMRES()),save_everystep=false)

using DiffEqOperators
Jv = JacVecOperator(brusselator_2d_loop,u0,p,0.0)

f = ODEFunction(brusselator_2d_loop;jac_prototype=Jv)
prob_ode_brusselator_2d_jacfree = ODEProblem(f,u0,(0.,11.5),p)
@btime solve(prob_ode_brusselator_2d_jacfree,TRBDF2(linsolve=LinSolveGMRES()),save_everystep=false)

using AlgebraicMultigrid
pc = aspreconditioner(ruge_stuben(jac_sparsity))
@btime solve(prob_ode_brusselator_2d_jacfree,TRBDF2(linsolve=LinSolveGMRES(Pl=pc)),save_everystep=false)

using Sundials
# Sparse Version
@btime solve(prob_ode_brusselator_2d_sparse,CVODE_BDF(),save_everystep=false)
# GMRES Version: Doesn't require any extra stuff!
@btime solve(prob_ode_brusselator_2d,CVODE_BDF(linear_solver=:GMRES),save_everystep=false)

using DifferentialEquations
function rober(du,u,p,t)
  y₁,y₂,y₃ = u
  k₁,k₂,k₃ = p
  du[1] = -k₁*y₁+k₃*y₂*y₃
  du[2] =  k₁*y₁-k₂*y₂^2-k₃*y₂*y₃
  du[3] =  y₁ + y₂ + y₃ - 1
  nothing
end
M = [1. 0  0
     0  1. 0
     0  0  0]
f = ODEFunction(rober,mass_matrix=M)
prob_mm = ODEProblem(f,[1.0,0.0,0.0],(0.0,1e5),(0.04,3e7,1e4))
sol = solve(prob_mm,Rodas5())

plot(sol, xscale=:log10, tspan=(1e-6, 1e5), layout=(3,1))

using SciMLTutorials
SciMLTutorials.tutorial_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])

max_stars_repo_name: SciML/SciMLTutorialsOutput

max_issues_repo_name: SciML/SciMLTutorialsOutput

max_forks_repo_name: SciML/SciMLTutorialsOutput


text:
# Dynamic Global Memory Allocation and Operations (B.21)

export malloc

@generated function malloc(sz::Csize_t)
    Context() do ctx
        T_pint8 = LLVM.PointerType(LLVM.Int8Type(ctx))
        T_size = convert(LLVMType, Csize_t; ctx)
        T_ptr = convert(LLVMType, Ptr{Cvoid}; ctx)

        # create function
        llvm_f, _ = create_function(T_ptr, [T_size])
        mod = LLVM.parent(llvm_f)

        # get the intrinsic
        # NOTE: LLVM doesn't have void*, Clang uses i8* for malloc too
        intr = LLVM.Function(mod, "malloc", LLVM.FunctionType(T_pint8, [T_size]))
        # should we attach some metadata here? julia.gc_alloc_obj has the following:
        #let attrs = function_attributes(intr)
        #    AllocSizeNumElemsNotPresent = reinterpret(Cuint, Cint(-1))
        #    packed_allocsize = Int64(1) << 32 | AllocSizeNumElemsNotPresent
        #    push!(attrs, EnumAttribute("allocsize", packed_allocsize; ctx))
        #end
        #let attrs = return_attributes(intr)
        #    push!(attrs, EnumAttribute("noalias", 0; ctx))
        #    push!(attrs, EnumAttribute("nonnull", 0; ctx))
        #end

        # generate IR
        Builder(ctx) do builder
            entry = BasicBlock(llvm_f, "entry"; ctx)
            position!(builder, entry)

            ptr = call!(builder, intr, [parameters(llvm_f)[1]])

            jlptr = ptrtoint!(builder, ptr, T_ptr)

            ret!(builder, jlptr)
        end

        call_function(llvm_f, Ptr{Cvoid}, Tuple{Csize_t}, :sz)
    end
end


text:
# ============================================================
# Extrapolation Functions
# ============================================================

function bc_extra!(ctr::AV, ng = 1::Integer; dirc)
    if Symbol(dirc) in (:xl, :xL)
        for i = 1:ng
            copy_ctr!(ctr[1-i], ctr[1])
        end
    elseif Symbol(dirc) in (:xr, :xR)
        nx = size(ctr, 1) - 2 * ng
        for i = 1:ng
            copy_ctr!(ctr[nx+i], ctr[nx])
        end
    end

    return nothing
end

function bc_extra!(ctr::AM, ng = 1::Integer; dirc)
    if Symbol(dirc) in (:xl, :xL)
        for j in axes(ctr, 2), i = 1:ng
            copy_ctr!(ctr[1-i, j], ctr[1, j])
        end
    elseif Symbol(dirc) in (:xr, :xR)
        nx = size(ctr, 1) - 2 * ng
        for j in axes(ctr, 2), i = 1:ng
            copy_ctr!(ctr[nx+i, j], ctr[nx, j])
        end
    elseif Symbol(dirc) in (:yl, :yL)
        for j = 1:ng, i in axes(ctr, 1)
            copy_ctr!(ctr[i, 1-j], ctr[i, 1])
        end
    elseif Symbol(dirc) in (:yr, :yR)
        ny = size(ctr, 2) - 2 * ng
        for j = 1:ng, i in axes(ctr, 1)
            copy_ctr!(ctr[i, ny+j], ctr[i, ny])
        end
    end

    return nothing
end


text:
using Pkg
Pkg.activate(".")
using DeutschlandDaten
using Dates
using DataFrames
using DataFramesMeta
using StatsBase
using Distributions
using StatsFuns

using StatsPlots
using Plots.PlotMeasures

deathdat = DeutschlandDaten.deaths_data()
popdata  = DeutschlandDaten.population()
popdatbl = DeutschlandDaten.population_bundesländer()

population_deaths(p, d=deathdat; deaths=tuple(), population=tuple(), kw...) =
    ( kw...,
      deaths...,
      population...,
      N=DeutschlandDaten.population(p,d; population..., kw...),
      D=DeutschlandDaten.deaths(d; deaths..., kw...)
      )

function dataprep!(mordat)
    @transform! mordat :geburt_start = :jahr .- rangestart.(:alter)
    @transform! mordat :kohorte = div.(:geburt_start, 10) * 10
    @transform! mordat :PD = :D ./ :N
    @transform! mordat :logitPD = logit.(:PD)
end

max_stars_repo_name: gkappler/DeutschlandDaten

max_issues_repo_name: gkappler/DeutschlandDaten

max_forks_repo_name: gkappler/DeutschlandDaten


text:
mutable struct Tree{T}
  head::T
  children::Vector{Tree{T}}
  Tree{T}(x::T, children = Tree{T}[]) where T =
    new{T}(x, children)
end

Tree(x, children = []) = Tree{typeof(x)}(x, children)

branch(xs) = length(xs) == 1 ? Tree(xs[1]) : Tree(xs[1], [branch(view(xs,2:lastindex(xs)))])

function mergechild!(children, child; (==) = Base.:(==), kws...)
  for x in children
    if x.head == child.head
      merge!(x, child; (==) = (==), kws...)
      return
    end
  end
  push!(children, child)
  return
end

function Base.merge!(a::Tree, b::Tree; merge = (a, b) -> a, kws...)
  a.head = merge(a.head, b.head)
  for child in b.children
    mergechild!(a.children, child; merge = merge, kws...)
  end
  return a
end

walk(x, inner, outer) = outer(Tree(x.head, map(inner, x.children)))
postwalk(f, x) = walk(x, x -> postwalk(f, x), f)
prewalk(f, x)  = walk(f(x), x -> prewalk(f, x), identity)


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

# The version of this JLL is decoupled from the upstream version.
# Whenever we package a new upstream release, we initially map its
# version X.Y.Z to X00.Y00.Z00 (i.e., multiply each component by 100).
# So for example version 2.6.3 would become 200.600.300.
#
# Moreover, all our packages using this JLL use `~` in their compat ranges.
#
# Together, this allows us to increment the patch level of the JLL for minor tweaks.
# If a rebuild of the JLL is needed which keeps the upstream version identical
# but breaks ABI compatibility for any reason, we can increment the minor version
# e.g. go from 200.600.300 to 200.601.300.
# To package prerelease versions, we can also adjust the minor version; e.g. we may
# map a prerelease of 2.7.0 to 200.690.000.
#
# There is currently no plan to change the major version, except when upstream itself
# changes its major version. It simply seemed sensible to apply the same transformation
# to all components.

name = "normaliz"
version = v"300.900.100"
upstream_version = v"3.9.1"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/Normaliz/Normaliz/releases/download/v$(upstream_version)/normaliz-$(upstream_version).tar.gz",
                  "ad5dbecc3ca3991bcd7b18774ebe2b68dae12ccca33c813ab29891beb85daa20")
]

# Bash recipe for building across all platforms
script = raw"""
cd normaliz-*
# avoid libtool problems
rm "${prefix}/lib/libgmpxx.la"
./configure --prefix=$prefix --host=$target --build=${MACHTYPE} --with-flint=$prefix --with-nauty=$prefix --with-gmp=$prefix CPPFLAGS=-I$prefix/include LDFLAGS=-L${libdir}
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
# windows build would require MPIR instead of GMP for 'long long'
platforms = supported_platforms(;experimental=true)
filter!(!Sys.iswindows, platforms)
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libnormaliz", :libnormaliz),
    ExecutableProduct("normaliz", :normaliz)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("GMP_jll", v"6.2.0"),
    Dependency("MPFR_jll", v"4.1.1"),
    Dependency("FLINT_jll"; compat = "~200.800.101"),
    Dependency("nauty_jll"; compat = "~2.6.13"),
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6", preferred_gcc_version=v"6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil