text:
using StaticArrays
using CSV

"""
Represents an atom or molecule, irespective of its charge.
"""
struct Formula
    #supports up to triatomic molecules, can be trivially extended.
    #Unlike tuples SVectors support sorting
    atoms::SVector{3, UInt8}

    function Formula(Z::Integer)
        @assert 1 <= Z <= Natoms
        new([0x00, 0x00, Z])
    end
    function Formula(Zs::AbstractVector{<:Integer})
        l = length(Zs)
        if l == 0
            throw(ArgumentError("Can't construct an empty Formula"))
        elseif l == 1
            Zs = [0 ; 0 ; Zs]
        elseif l == 2
            Zs = [0 ; Zs]
        elseif l > 3
            throw(ArgumentError("Can't construct Formula with more than three atoms"))
        end
        @assert(issorted(Zs))
        new(Zs)
    end

    """
        Formula(code::String)

        Construct a Formula from an encoded string form.  This can be a MOOG-style numeric code, i.e.
        "0801" for OH, or an atomic or molecular symbol, i.e. "FeH", "Li", or "C2".
    """
    function Formula(code::AbstractString)
        if code in atomic_symbols #quick-parse single elements
            return new([0, 0, atomic_numbers[code]])
        end

        #handle numeric codes, e.g. 0801 -> OH
        if all(isdigit(c) for c in code)
            if length(code) <= 2
                return Formula(parse(Int,code))
            elseif length(code) <= 4
                if length(code) == 3
                    code = "0"*code
                end
                el1 = parse(Int, code[1:2])
                el2 = parse(Int, code[3:4])
                return new([0x00, min(el1, el2), max(el1, el2)])
            else
                throw(ArgumentError("numeric codes for molecules with more than 4 chars like " *
                                    "$(code) are not supported"))
            end
        end
        #otherwise, code should be "OH", "FeH", "Li", "C2", etc.
        inds = filter(1:length(code)) do i
            isdigit(code[i]) || isuppercase(code[i])
        end
        push!(inds, length(code)+1)
        subcode = map(1:(length(inds)-1)) do j
            code[inds[j]:inds[j+1]-1]
        end

        atoms = UInt8[]
        for s in subcode
            num = tryparse(Int, s)
            if num isa Int
                previous = atoms[end]
                for _ in 1:(num-1)
                    push!(atoms, previous)
                end
            else
                push!(atoms, atomic_numbers[s])
            end
        end

        sort!(atoms)
        Formula(atoms)
    end
end

function get_atoms(f::Formula)
    if f.atoms[2] == 0
        view(f.atoms,3:3)
    elseif f.atoms[1] == 0
        view(f.atoms,2:3)
    else
        view(f.atoms, 1:3)
    end
end

#pretty-print lines in REPL and jupyter notebooks
function Base.show(io::IO, m::MIME"text/plain", f::Formula)
    print(io, *([atomic_symbols[i] for i in f.atoms if i != 0]...))
end

"""
    ismolecule(f::Formula)

`true` when `f` is composed of more than one atom
"""
ismolecule(f::Formula) = f.atoms[2] != 0

"""
    get_mass(f::Formula)

Returns the mass [g] of `f`.
"""
function get_mass(f::Formula)
    sum(atomic_masses[a] for a in get_atoms(f))
end

const roman_numerals = ["I", "II", "III", "IV", "V", "VI", "VII", "VIII", "IX", "X"]
"""
Represents an atom or molecule (a `Formula`) with a particular number of electrons (regardless of
their configuration).
"""
struct Species
    formula::Formula
    charge::Int
end

"""
    Species(code::AbstractString)

Parse the "species code" as it is often specified in linelists and return a the "astronomy"
notation. 01.00 → "H_I", 02.01 → "He_II", 02.1000 → "He_II", 0608 → "CO_I", etc.

To parse at compile time, use the `species` string macro, i.e. `species"H I"`.
"""
function Species(code::AbstractString)
    code = strip(code, ['0', ' '])
    toks = split(code, [' ', '.', '_'])
    if length(toks) > 2
        throw(ArgumentError(code * " isn't a valid species code"))
    end
    formula = Formula(toks[1])
    charge = if length(toks) == 1 || length(toks[2]) == 0
        0 #no charge specified -> assume neutral
    else
        charge = findfirst(toks[2] .== roman_numerals)
        charge = (charge isa Int ? charge : parse(Int, toks[2]))
        #if this is a MOOG-style numeric code, the charge is correct, otherwise subtract 1
        if tryparse(Float64, code) == nothing
            charge -= 1
        end
        charge
    end
    Species(formula, charge)
end

#used to contruct Species at compile time and avoid parsing in hot loops
macro species_str(s)
    Species(s)
end

#pretty-print lines in REPL and jupyter notebooks
function Base.show(io::IO, m::MIME"text/plain", s::Species)
    show(io, m, s.formula)
    print(io, " ", get_roman_numeral(s))
end

ismolecule(s::Species) = ismolecule(s.formula)
get_mass(s::Species) = get_mass(s.formula)
get_atoms(s::Species) = get_atoms(s.formula)
get_roman_numeral(s::Species) = roman_numerals[s.charge+1]

#This type represents an individual line.
struct Line{F}
    wl::F                     #cm
    log_gf::F                 #unitless
    species::Species
    E_lower::F                #eV (also called the excitation potential)
    gamma_rad::F              #s^-1
    gamma_stark::F            #s^-1
    vdW::Union{F, Tuple{F,F}} #either Γ_vdW [s^-1] per electron or (σ, α) from ABO theory

    @doc """
        Line(wl::F, log_gf::F, species::Species, E_lower::F,
             gamma_rad::Union{F, Missing}=missing, gamma_stark::Union{F, Missing}=missing,
             vdw::Union{F, Tuple{F, F}, Missing}, missing) where F <: Real

    Construct a `Line`.  If any of `gamma_rad`, `gamma_stark`, or `vdW` are `missing`, guess them.
    `vdW` may be log(Γ_vdW) (assumed if negative), Γ_vdW (assumed if 0 < `vdW` < 1), or packed ABO
    parameters (assumed if `vdW` > 1).  It may also be passed as a Tuple, `(σ, α)`.
    """
    function Line(wl::F, log_gf::F, species::Species, E_lower::F,
                  gamma_rad::Union{F, Missing}=missing, gamma_stark::Union{F, Missing}=missing,
                  vdW::Union{F, Tuple{F, F}, Missing}=missing) where F <: Real
        if ismissing(gamma_stark) || ismissing(vdW)
            gamma_stark_approx, vdW_approx = approximate_gammas(wl, species, E_lower)
            if ismissing(gamma_stark)
                gamma_stark = gamma_stark_approx
            end
            if ismissing(vdW)
                vdW = vdW_approx
            end
        end
        if ismissing(gamma_rad)
            gamma_rad = approximate_radiative_gamma(wl, log_gf)
        end

        if vdW isa F
            if vdW < 0 #if vdW is negative, assume it's log(Γ_vdW)
                vdW = 10^vdW
            elseif vdW > 1 #if it's > 1 assume it's packed ABO params
                vdW = (floor(vdW) * bohr_radius_cgs * bohr_radius_cgs, vdW - floor(vdW))
            end
        end

        new{F}(wl, log_gf, species, E_lower, gamma_rad, gamma_stark, vdW)
    end
end

#pretty-print lines in REPL and jupyter notebooks
function Base.show(io::IO, m::MIME"text/plain", line::Line)
    show(io, m, line.species)
    print(io, " ", round(line.wl*1e8, digits=6), " Å")
end

"""
    approximate_radiative_gamma(wl, log_gf)

Approximate radiate broadening parameter.
"""
function approximate_radiative_gamma(wl, log_gf)
    e = electron_charge_cgs
    m = electron_mass_cgs
    c = c_cgs
    8π^2 * e^2 / (m * c * wl^2) * 10^log_gf
end

"""
    approximate_gammas(wl, species, E_lower; ionization_energies=Korg.ionization_energies)

A simplified form of the Unsoeld (1955) approximation for van der Waals broadening and the
[Cowley 1971](https://ui.adsabs.harvard.edu/abs/1971Obs....91..139C/abstract) approximation for
Stark broadening, evaluated at 10,000 K.
Used for atomic lines with no vdW and stark broadening info in the linelist.

Returns `(γ_stark`, `log10(γ_vdW))` in Hz, where these are the per-perturber quantities.
For autoionizing lines (those for which E_upper > χ), Returns 0.0 for γ_vdW.

In the calculation of `n*²`, uses the approximation that
``\\overbar{r^2} = 5/2 {n^*}^4 / Z^2``
which neglects the dependence on the angular momentum quantum number, l, in the the form given by
[Warner 1967](https://ui.adsabs.harvard.edu/abs/1967MNRAS.136..381W/abstract) (the earliest english
work reporting the Unsoeld result).
"""
function approximate_gammas(wl, species, E_lower; ionization_energies=ionization_energies)
    Z = species.charge + 1 #Z is ionization stage, not atomic number
    if ismolecule(species) || Z > 3
        return 0.0,0.0
    end
    χ = ionization_energies[get_atoms(species.formula)[1]][Z]
    c = c_cgs
    h = hplanck_eV
    k = kboltz_cgs
    E_upper = E_lower + (h * c / wl)

    #It's not obvious to me which Rydberg constant to use here, and below in Δrbar2.  The sources
    #are not entirely clear. It doesn't make a big difference.
    nstar4_upper = (Z^2 * RydbergH_eV / (χ - E_upper))^2
    #I'm not actually able to reproduce Crowley 1971 equation 7 (his simplified form) from equation
    #5, but these match the values in the Turbospectrum source, so they are probably correct.
    #The constants here were calculated assuming that "v" is the mean (not modal) electron speed
    if Z == 1
        γstark = 2.25910152e-7 * nstar4_upper #Cowley (1971) equation 5 evaluated at T=10,000 K
    else
        #Cowley (1971) equation 6 @ T=10,000 K (n.b. the constant is 12/5 * that above)
        γstark = 5.42184365e-7 * nstar4_upper / (Z + 1)^2
    end

    Δrbar2 = (5/2) * Rydberg_eV^2 * Z^2 * (1/(χ - E_upper)^2 - 1/(χ - E_lower)^2)
    if χ < E_upper
        γvdW = 0.0
    else
        # (log) γ_vdW From R J Rutten's course notes.
        # Equations 11.29 and 11.30 from Gray 2005 are equivalent
        γvdW = 6.33 + 0.4log10(Δrbar2) + 0.3log10(10_000) + log10(k)
    end

    γstark, γvdW
end

"""
    read_linelist(filename; format="vald", isotopic_abundances=Korg.isotopic_abundances)

Parse a linelist file, returning a vector of [`Line`](@ref)s.

Pass `format="kurucz"` for a [Kurucz linelist](http://kurucz.harvard.edu/linelists.html)
(`format=kurucz_vac` if it uses vacuum wavelengths; Be warned that Korg will not assume that
wavelengths are vacuum below 2000 Å),`format="vald"` for a
[VALD](http://vald.astro.uu.se/~vald/php/vald.php) linelist, and `format="moog"` for a MOOG linelist
(doesn't support broadening parameters or dissociation energies).

VALD linelists (the default and preferred format) can be either "short" or "long" format,
"extract all" or "extract stellar".  Air wavelengths will automatically be converted into vacuum
wavelengths, and energy levels will be automatically converted from cm``^{-1}`` to eV.

When they are not pre-scaled by isotopic abundace (which VALD does by default), Korg will
automatically adjust the `log_gf` of each line according to `isotopic_abundances`, which defaults
to the values from [NIST](https://www.nist.gov/pml/atomic-weights-and-isotopic-compositions-relative-atomic-masses).
To use custom isotopic abundances, just pass `isotopic_abundances` as a dictionary mapping
`(atomic number, atomic weight)` pairs to abundances between 0 and 1.
"""
function read_linelist(fname::String; format="vald", isotopic_abundances::Dict=isotopic_abundances
                      ) :: Vector{Line}
    format = lowercase(format)
    linelist = open(fname) do f
        if format == "kurucz"
            parse_kurucz_linelist(f; vacuum=false)
        elseif format == "kurucz_vac"
            parse_kurucz_linelist(f; vacuum=true)
        elseif format == "vald"
            parse_vald_linelist(f, isotopic_abundances)
        elseif format == "moog"
            parse_moog_linelist(f)
        else
            throw(ArgumentError("$(format) is not a supported linelist format"))
        end
    end

    filter!(linelist) do line #filter triply+ ionized and hydrogen lines
        (0 <= line.species.charge <= 2) && (line.species != species"H_I")
    end

    #ensure linelist is sorted
    if !issorted(linelist, by=l->l.wl)
        sort!(linelist, by=l->l.wl)
    end

    linelist
end

#used to handle missing gammas in vald and kurucz lineslist parsers
expOrMissing(x) = x == 0.0 ? missing : 10.0^x
idOrMissing(x) = x == 0.0 ? missing : x

function parse_kurucz_linelist(f; vacuum=false)
    lines = Line[]
    for row in eachline(f)
        row == "" && continue #skip empty lines

        #some linelists have a missing column in the wavelenth region
        if length(row) == 159
            row = " " * row
        end

        #kurucz provides wavenumbers for "level 1" and "level 2", which is which is
        #determined by parity
        E_levels = map((row[25:36], row[53:64])) do s
            #abs because Kurucz multiplies predicted values by -1
            abs(parse(Float64,s)) * c_cgs * hplanck_eV
        end

        wl_transform = vacuum ? identity : air_to_vacuum

        push!(lines, Line(wl_transform(parse(Float64, row[1:11])*1e-7), #convert from nm to cm
                     parse(Float64, row[12:18]),
                     Species(row[19:24]),
                     min(E_levels...),
                     expOrMissing(parse(Float64, row[81:86])),
                     expOrMissing(parse(Float64, row[87:92])),
                     idOrMissing(parse(Float64, row[93:98]))))
    end
    lines
end

function parse_vald_linelist(f, isotopic_abundances)
    lines = filter!(collect(eachline(f))) do line
        length(line) > 0 && line[1] != '#' #remove comments and empty lines
    end

    lines = replace.(lines, "'"=>"\"") #replace single quotes with double

    # is this an "extract all" or an "extract stellar" linelist?
    extractall = !occursin(r"^\s+\d", lines[1])
    firstline = extractall ? 3 : 4
    header = lines[firstline - 1]

    scale_isotopes = any(startswith.(lines, "* oscillator strengths were NOT scaled "))
    if !scale_isotopes && !any(startswith.(lines,"* oscillator strengths were scaled "))
        throw(ArgumentError("Can't parse linelist.  Can't detect whether log(gf)s are scaled by "*
                            "isotopic abundance."))
    end

    #we take the linelist to be long-format when the second line after the header starts with a
    #space or a quote followed a space.  In some linelists the quotes are there, but in others
    #they are not.
    shortformat = !(occursin(r"^\"? ", lines[firstline + 1]))
    body = lines[firstline : (shortformat ? 1 : 4) : end]
    body = body[1 : findfirst(l->l[1]!='\"' || !isuppercase(l[2]), body)-1]

    CSVheader = if shortformat && extractall
        ["species", "wl", "E_low", "loggf", "gamma_rad", "gamma_stark", "gamma_vdW", "lande",
         "reference"]
    elseif shortformat #extract stellar
        ["species", "wl", "E_low", "Vmic", "loggf", "gamma_rad", "gamma_stark", "gamma_vdW",
         "lande", "depth", "reference"]
    else #long format (extract all or extract stellar)
        ["species", "wl", "loggf", "E_low", "J_lo", "E_up", "J_up", "lower_lande", "upper_lande",
         "mean_lande", "gamma_rad", "gamma_stark", "gamma_vdW"]
    end
    body = CSV.File(reduce(vcat, codeunits.(body.*"\n")), header=CSVheader, delim=',',
                    silencewarnings=true)

    E_low = if contains(header, "cm") #convert E_low to eV if necessary
        body.E_low * c_cgs * hplanck_eV
    elseif contains(header, "eV")
        body.E_low
    else
        throw(ArgumentError( "Can't parse linelist.  Can't determine energy units: " * E_col))
    end

    wl = if contains(header, "air") #convert wls to vacuum if necessary
        air_to_vacuum.(body.wl)
    elseif contains(header, "vac")
        body.wl
    else
        throw(ArgumentError( "Can't parse linelist.  Can't determine vac/air wls: " * header))
    end

    loggf = body.loggf .+ if scale_isotopes
        refs = if !shortformat #the references are on different lines
            lines[firstline+3 .+ ((0:length(body)-1) .* 4)]
        else #references are in the last column
            body.reference
        end

        map(refs) do ref
            #find things that look like (16)O or (64)Ni in reference string
            regexp = r"\((?<isotope>\d\d?\d?)\)(?<elem>\p{Lu}\p{Ll}?)"
            #add up the adjustments to log(gf) from isotopic abundances (zero if no info is present)
            log_probs = map(findall(regexp, ref)) do r
                m = match(regexp, ref[r])
                log10(isotopic_abundances[atomic_numbers[m["elem"]], parse(Int64, m["isotope"])])
            end
            sum([0 ; log_probs])
        end
    else
        0
    end

    Line.(wl * 1e-8, loggf, Species.(body.species), E_low, expOrMissing.(body.gamma_rad),
          expOrMissing.(body.gamma_stark), idOrMissing.(body.gamma_vdW))
end

#todo support moog linelists with broadening parameters?
function parse_moog_linelist(f)
    lines = collect(eachline(f))
    #moog format requires blank first line
    linelist = map(lines[2:end]) do line
        toks = split(line)
        Line(parse(Float64, toks[1]) * 1e-8, #convert Å to cm
             parse(Float64, toks[4]),
             Species(toks[2]),
             parse(Float64, toks[3]))
    end
    linelist
end


text:
using IntervalOptimization
using Test

@testset "IntervalOptimization.jl" begin
    # Write your tests here.
end

ext: jl


text:
# Test Benders Algorithm for Mixed Integer Linear Programming
# Edward J. Xu, edxu96@outlook.com
# April 5th, 2019
# ---------------------------------------------------------------------------
push!(LOAD_PATH, "$(homedir())/Desktop/Benders Decomposition, MILP")
cd("$(homedir())/Desktop/Benders Decomposition, MILP")
using Benders
# ---------------------------------------------------------------------------
# Test Example 1: mat_a and mat_b being column vector
n_x = 1
n_y = 1
vec_min_y = hcat([0])
vec_max_y = hcat([10])
vec_c = hcat([5])
vec_f = hcat([-3])
vec_b = hcat([4; 0; -13])
mat_a = hcat([1; 2; 1])
mat_b = hcat([2; -1; -3])
Benders.milp(n_x = n_x,
             n_y = n_y,
             vec_min_y = vec_min_y,
             vec_max_y = vec_max_y,
             vec_c = vec_c,
             vec_f = vec_f,
             vec_b = vec_b,
             mat_a = mat_a,
             mat_b = mat_b,
             epsilon = 0,
             timesIterationMax = 5)
# ---------------------------------------------------------------------------
# Test Example 2: mat_a and mat_b being matrix
n_x_2 = 2
n_y_2 = 2
vec_min_y_2 = hcat([0; 0])
vec_max_y_2 = hcat([10; 10])
vec_c_2 = hcat([5; 3])
vec_f_2 = hcat([-3; 1])
vec_b_2 = hcat([4; 0; -13])
mat_a_2 = [1 3; 2 1; 1 -5]
mat_b_2 = [2 -4; -1 2; -3 1]
Benders.milp(n_x = n_x_2,
             n_y = n_y_2,
             vec_min_y = vec_min_y_2,
             vec_max_y = vec_max_y_2,
             vec_c = vec_c_2,
             vec_f = vec_f_2,
             vec_b = vec_b_2,
             mat_a = mat_a_2,
             mat_b = mat_b_2,
             epsilon = 0.000001,
             timesIterationMax = 10)
# ---------------------------------------------------------------------------
# Test Example 3: sub and ray problem
n_x_3 = 2
n_y_3 = 2
vec_min_y_3 = hcat([0; 0])
vec_max_y_3 = hcat([2; 2])
vec_c_3 = hcat([2; 6])
vec_f_3 = hcat([2; 3])
vec_b_3 = hcat([5; 4])
mat_a_3 = [-1 2; 1 -3]
mat_b_3 = [3 -1; 2 2]
Benders.milp(n_x = n_x_3,
             n_y = n_y_3,
             vec_min_y = vec_min_y_3,
             vec_max_y = vec_max_y_3,
             vec_c = vec_c_3,
             vec_f = vec_f_3,
             vec_b = vec_b_3,
             mat_a = mat_a_3,
             mat_b = mat_b_3,
             epsilon = 0.000001,
             timesIterationMax = 10)

max_stars_repo_name: eitt/BendersDecomposition


text:
using DiscreteElasticRods
DER = DiscreteElasticRods
using Zygote
using ForwardDiff
using FiniteDiff
using ReverseDiff
using BenchmarkTools
using IterativeSolvers
using SparseDiffTools
using DelimitedFiles

struct AHesVec{F,uType}
    f::F
    u::uType
end

Base.size(L::AHesVec) = (length(L.u),length(L.u))
Base.size(L::AHesVec,i::Int) = length(L.u)
Base.:*(L::AHesVec,x::AbstractVector) =  SparseDiffTools.autoback_hesvec(L.f, L.u, x)

import LinearAlgebra
function LinearAlgebra.mul!(du::AbstractVector,L::AHesVec,v::AbstractVector)
    du .= L*v
end

# Rod sizes
n_array = [10, 20, 30, 40, 50]
grad_array = zeros(3, length(n_array))
hess_array = zeros(5, length(n_array))
minres_array = zeros(2, length(n_array))

#for k = 1:length(n_array)
k = 1
    # Setup an energy functional
n = n_array[k]
ref_rod = DER.random_rod(n)
ref_strains = DER.full_kinematics(ref_rod)
rod_props = DER.elliptical_stiffness(10, 0.3, 0.01, 0.04)

function total_energy(q)
    Δx = reshape(q[1:3*n+6], 3, n+2)
    Δθ = reshape(q[3*n+7:end], 1, n+1)
    rod_coords = DER.rod_delta(Δx, Δθ)

    new_rod = DER.rod_update(ref_rod, rod_coords)
    strains = DER.full_kinematics(new_rod)
    DER.elastic_energy(ref_strains, strains, rod_props)
end

function total_energy_array(q)
    Δx = reshape(q[1:3*n+6], 3, n+2)
    Δθ = reshape(q[3*n+7:end], 1, n+1)
    x, d = DER.rod_update(ref_rod.x, ref_rod.d, Δx, Δθ)
    l, κ, τ = DER.full_kinematics(x, d)
    strains = DER.rod_strains(l, κ, τ)
    DER.elastic_energy(ref_strains, strains, rod_props)
end

nq = 4*n+7
q = rand(nq)
# tape = ReverseDiff.GradientTape(total_energy_array, q)
# compiled_tape = ReverseDiff.compile(tape)
#
# grad_reverse!(out, q) = ReverseDiff.gradient!(out, compiled_tape, q)
# grad_forward(q) = ForwardDiff.gradient(total_energy, q)
# grad_zygote(q) = first(Zygote.gradient(total_energy, q))
#
# # Timing the Gradients
# out = similar(q)
# grad_reverse!(out,q)
# grad_forward(q)
# grad_zygote(q)
# rbench = @benchmark grad_reverse!($out, $q)
# fbench = @benchmark grad_forward($q)
# zbench = @benchmark grad_zygote($q)
#
# grad_array[1,k] = minimum(fbench).time/1e9
# grad_array[2,k] = minimum(rbench).time/1e9
# grad_array[3,k] = minimum(zbench).time/1e9

# q = zeros(nq)
# v = rand(nq)
# nh = @benchmark num_hesvec(total_energy, $q, $v)
# nah = @benchmark numauto_hesvec(total_energy, $q, $v)
# anh = @benchmark autonum_hesvec(total_energy, $q, $v)
# nbh = @benchmark numback_hesvec(total_energy, $q, $v)
# abh = @benchmark autoback_hesvec(total_energy, $q, $v)
#
# hess_array[1,k] = minimum(nh).time/1e9
# hess_array[2,k] = minimum(nah).time/1e9
# hess_array[3,k] = minimum(anh).time/1e9
# hess_array[4,k] = minimum(nbh).time/1e9
# hess_array[5,k] = minimum(abh).time/1e9

#writedlm("hessian_timing.txt", hess_array, ',')

grd = rand(nq)

Hvec1 = HesVec(total_energy_array, q; autodiff=true)
out = zeros(length(q))
mna = @benchmark minres!($out, $Hvec1, $grd)

AHVec = AHesVec(total_energy, q)
out = zeros(length(q))
mab = @benchmark minres!($out, $AHVec, $grd)

minres_array[1,k] = minimum(mna).time/1e9
minres_array[2,k] = minimum(mab).time/1e9
writedlm("minres_timing.txt", minres_array, ',')

print("Benchmarked n = $n\n")
#end

writedlm("gradient_timing.txt", grad_array, ',')
writedlm("hessian_timing.txt", hess_array, ',')
writedlm("minres_timing.txt", minres_array, ',')

# ------------------------------------------------------------------------------
# Minres Timing
# ------------------------------------------------------------------------------


text:
# ------------------------------------------------------------------------------------------
# # Visual Q&A Demo
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Setup
# 1-Adds required packages to Julia.
# 2-Downloads sample data and a pretrained model.
# ------------------------------------------------------------------------------------------

include("demosetup.jl")

# ------------------------------------------------------------------------------------------
# ## Initialization
# 1-Loads the sample demo data (image features,questions,vocabulary).
# 2-Loads the pretrained model.
# ------------------------------------------------------------------------------------------

include("src/newmacnetwork.jl")
feats,qstsns,(w2i,a2i,i2w,i2a) = loadDemoData("data/demo/");
_,wrun,r,_,p = loadmodel("models/macnet.jld";onlywrun=true);
if !(typeof(first(wrun)) <: atype)
    wrun = map(atype,wrun);
end;

# ------------------------------------------------------------------------------------------
# ## Sample Data
# 1-Randomly selects (question,image) pair from the sample data
# 2-Make predictions for the question and checks whether the prediction is correct
# ------------------------------------------------------------------------------------------

rnd        = rand(1:length(qstsns))
inst       = qstsns[rnd]
feat       = feats[:,:,:,rnd:rnd]
question   = Array{Int}(inst[2])
answer     = inst[3];
family     = inst[4];
results,prediction = singlerun(wrun,r,feat,question;p=p);
answer==prediction

img = load("data/demo/CLEVR_v1.0/images/val/$(inst[1])")

textq  = i2w[question];
println("Question: ",join(textq," "))
texta  = i2a[answer];
println("Answer: $(texta)\nPrediction: $(i2a[prediction]) ")

# ------------------------------------------------------------------------------------------
# ## User Data
# You can enter your own question about the image and test whether the prediction is correct
# ------------------------------------------------------------------------------------------

userinput = readline(STDIN)
words = split(userinput) # tokenize(userinput)
question = [get!(w2i,wr,1) for wr in words]
results,prediction = singlerun(wrun,r,feat,question;p=p);
println("Question: $(join(i2w[question]," "))")
println("Prediction: $(i2a[prediction])")

# ------------------------------------------------------------------------------------------
# ## Visualize
# `visualize` function visualizes attention maps for each time step of the mac network
# ------------------------------------------------------------------------------------------

visualize(img,results;p=p)

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
# Package initiation hook
function __init__()
    # Unfortunately, the default GR backend does not work well with
    # offline plotting (see this issue https://github.com/JuliaPlots/Plots.jl/issues/2127).
    # So, we will use PyPlot backend instead.
    @static if Sys.islinux()
        pyplot()
    end
end

function help_message()
    act_commands = collect(keys(filter(c -> c.second, ACTIVE_COMMANDS)))
    names = collect(values(filter(c -> ACTIVE_COMMANDS[c.first], COMMANDS_NAMES)))

    commands = mapreduce(c -> string(c) * "\n", *, act_commands)
    opt = mapreduce(c -> string(c) * "\n", *, filter(c -> c ∈ names, OPT_SERVICES_LIST))
    str = """
    ```
    HoJBot accepts the following commands:
    ```
    $commands
    ```
    The following services are opt-in. Please check the related help command (`service help` for any `service` below).
    ```
    $opt
    """
    return str
end

function start_bot(;
    commands=ACTIVE_COMMANDS,
    handlers=HANDLERS_LIST,
    run_duration = Minute(365 * 24 * 60),  # run for a very long time by default
)
    @info "Starting bot... command prefix = $COMMAND_PREFIX"
    global client = Client(
        ENV["HOJBOT_DISCORD_TOKEN"];
        presence=(game=(name="HoJ", type=AT_GAME),),
        prefix=COMMAND_PREFIX,
    )
    init_handlers!(client, handlers)
    init_commands!(client, commands)
    # add_help!(client; help = help_message())
    warm_up_enabled() && warm_up()
    open(client)
    auto_shutdown(client, run_duration, "SHUTDOWN")
    wait(client)
end

function init_handlers!(client::Client, handlers)
    for (symbol, event, active) in handlers
        active && add_handler!(client, event, (c, e) -> handler(c, e, Val(symbol)))
    end
end

function init_commands!(client::Client, commands)
    for (com, active) in commands
        active && add_command!(client, com, (c, m) -> commander(c, m, COMMANDS_NAMES[com]))
    end
end

"""
    auto_shutdown(run_duration::TimePeriod)

Run a background process to track the program's run time and exit
the program when it has exceeded the specified `run_duration` or
when a file exists at `trigger_path`.
"""
function auto_shutdown(c::Client, run_duration::TimePeriod, trigger_path::AbstractString="")
    start_time = now()
    @async while true
        if now() > start_time + run_duration
            @info "Times up! The bot is shutting down automatically."
            shutdown_gracefully(c)
        end
        if length(trigger_path) > 0 && isfile(trigger_path)
            @info "The bot is shutting down via trigger path `$trigger_path`."
            rm(trigger_path)
            shutdown_gracefully(c)
        end
        sleep(5)
    end
end

function shutdown_gracefully(c::Client)
    try
        close(c)
    catch ex
        @warn "Unable to close client connection" ex
    end
    sleep(1)
    exit(0)
end

commander(c::Client, m::Message, service) = begin
    @info "commander" c m service
    commander(c, m, Val(service))
end
help_commander(c::Client, m::Message, service) = help_commander(c, m, Val(service))

function get_opt(username, discriminator)
    user = username * "_" * discriminator
    path = joinpath(pwd(), "data", "opt", user)
    !isfile(path) && write(path, "{}")
    return JSON.parsefile(path; dicttype=LittleDict)
end

function get_opt!(username, discriminator, service)
    get!(get_opt(username, discriminator), string(service), false)
end

function set_opt!(username, discriminator, service, value)
    opt = get_opt(username, discriminator)
    opt[string(service)] = value
    user = username * "_" * discriminator
    path = joinpath(pwd(), "data", "opt", user)
    write(path, json(opt))
end

function opt_in(c::Client, m::Message, service)
    username = m.author.username
    discriminator = m.author.discriminator
    if service ∉ OPT_SERVICES_LIST
        reply(
            c,
            m,
            """
            Oh, @$username#$discriminator, the `$(string(service))` service either does not exist or accept opt-in/opt-out requests!
            """,
        )
    elseif get_opt!(username, discriminator, service)
        reply(
            c,
            m,
            """
            Oh, @$username#$discriminator, you already suscribed to the `$(string(service))` service!
            """,
        )
    else
        set_opt!(username, discriminator, service, true)
        reply(
            c,
            m,
            """
            Thanks @$username#$discriminator for joining our `$(string(service))` service!
            """,
        )
    end
end

function opt_out(c::Client, m::Message, service)
    username = m.author.username
    discriminator = m.author.discriminator
    if service ∉ OPT_SERVICES_LIST
        reply(
            c,
            m,
            """
            Oh, @$username#$discriminator, the `$(string(service))` service either does not exist or accept opt-in/opt-out requests!
            """,
        )
    elseif !get_opt!(username, discriminator, service)
        reply(
            c,
            m,
            """
            Oh, @$username#$discriminator, you haven't suscribed to the `$(string(service))` service!
            """,
        )
    else
        set_opt!(username, discriminator, service, false)
        reply(
            c,
            m,
            """
            It is sad that you're leaving our `$(string(service))` service, @$username#$discriminator. We hope you will come back soon and enjoy other HoJBot stuff!
            """,
        )
    end
end

function warm_up_enabled()
    return get(ENV, "ENABLE_WARM_UP", "1") in ("1", "Y", "YES")
end

# TODO use SnoopCompile to find precompile methods
function warm_up()
    @info "Warming up..."
    Threads.@spawn begin
        dummy_user_id = UInt64(0)
        elapsed = @elapsed try
            symbol = "AAPL"
            ig_get_quote(symbol)
            ig_save_portfolio(
                dummy_user_id, IgPortfolio(100, [IgHolding(symbol, 100, today(), 130)])
            )
            ig_load_portfolio(dummy_user_id)
            ig_ranking_table(Client("hey"))

            from_date, to_date = Date(2020, 1, 1), Date(2020, 12, 31)
            df = ig_historical_prices(symbol, from_date, to_date)
            ig_chart(symbol, df.Date, df."Adj Close")
        catch ex
            @error "Warm up error: " ex
            Base.showerror(stdout, ex, catch_backtrace())
        finally
            ig_remove_game(dummy_user_id)
        end
        @info "Completed warm up in $elapsed seconds"
    end
end


text:
module TensorTrain

using LinearAlgebra, Random
using ..Auxiliary

include("TensorTrainFactor.jl")
include("TensorTrainFactorization.jl")

include("TT.jl")

end


text:
module Exceptions

export @ArgumentError, @APIUsageError, @ArgumentEvaluationError, @FileError, @IncoherentArgumentError, @InsecureRecognitionException, @KeywordArgumentError
export                  APIUsageError,  ArgumentEvaluationError,  FileError,  IncoherentArgumentError,  InsecureRecognitionException,  KeywordArgumentError

macro ArgumentError(msg) esc(:(throw(ArgumentError($msg)))) end
macro APIUsageError(msg) esc(:(throw(APIUsageError($msg)))) end
macro ArgumentEvaluationError(msg) esc(:(throw(ArgumentEvaluationError($msg)))) end
macro FileError(msg) esc(:(throw(APIUsageError($msg)))) end
macro IncoherentArgumentError(msg) esc(:(throw(IncoherentArgumentError($msg)))) end
macro InsecureRecognitionException(msg) esc(:(throw(InsecureRecognitionException($msg)))) end
macro KeywordArgumentError(msg) esc(:(throw(KeywordArgumentError($msg)))) end

struct APIUsageError <: Exception
    msg::String
end
Base.showerror(io::IO, e::APIUsageError) = print(io, "APIUsageError: ", e.msg)

struct ArgumentEvaluationError <: Exception
    msg::String
end
Base.showerror(io::IO, e::ArgumentEvaluationError) = print(io, "ArgumentEvaluationError: ", e.msg)

struct FileError <: Exception
    msg::String
end
Base.showerror(io::IO, e::FileError) = print(io, "FileError: ", e.msg)

struct InsecureRecognitionException <: Exception
    msg::String
end
Base.showerror(io::IO, e::InsecureRecognitionException) = print(io, "InsecureRecognitionException: ", e.msg)

struct IncoherentArgumentError <: Exception
    msg::String
end
Base.showerror(io::IO, e::IncoherentArgumentError) = print(io, "IncoherentArgumentError: ", e.msg)

struct KeywordArgumentError <: Exception
    msg::String
end
Base.showerror(io::IO, e::KeywordArgumentError) = print(io, "KeywordArgumentError: ", e.msg)

end # Module Exceptions


text:
# ## http://prosite.expasy.org/scanprosite/scanprosite-doc.html#pattern_syntax
# ## Pattern syntax:

# 1 The standard IUPAC one-letter codes for the amino acids are used in PROSITE.
# 2 The symbol `x' is used for a position where any amino acid is accepted.
# 3 Ambiguities are indicated by listing the acceptable amino acids for a given position, between square brackets `[ ]'. For example: [ALT] stands for Ala or Leu or Thr.
# 4 Ambiguities are also indicated by listing between a pair of curly brackets `{ }' the amino acids that are not accepted at a given position. For example: {AM} stands for any amino acid except Ala and Met.
# 5 Each element in a pattern is separated from its neighbor by a `-'.
# 6 Repetition of an element of the pattern can be indicated by following that element with a numerical value or, if it is a gap ('x'), by a numerical range between parentheses.
#     Examples:
#         x(3) corresponds to x-x-x
#         x(2,4) corresponds to x-x or x-x-x or x-x-x-x
#         A(3) corresponds to A-A-A
#     Note: You can only use a range with 'x', i.e. A(2,4) is not a valid pattern element.
# 7 When a pattern is restricted to either the N- or C-terminal of a sequence, that pattern either starts with a `<' symbol or respectively ends with a `>' symbol. In some rare cases (e.g. PS00267 or PS00539), '>' can also occur inside square brackets for the C-terminal element. 'F-[GSTV]-P-R-L-[G>]' means that either 'F-[GSTV]-P-R-L-G' or 'F-[GSTV]-P-R-L>' are considered.

# ## Extended syntax allowed on BioSeq:

# If your pattern consists of one-letter amino acid codes only, without any ambiguous residues, you need not specify the '-', i.e. you can directly copy/paste peptide sequences into the text field.
# Example: M-A-S-K-E can be written as MASKE.
# To search all sequences which do not contain a certain amino acid, e.g Cys, you can use <{C}*>.

const _prosite_alpha = IntSet(Int['A','C','D','E','F','G','H','I','K','L','M','N','P','Q','R','S','T','V','W','Y',
		              '[',']','{','}','(',')','<','>',',','x','-','*',
		              '1','2','3','4','5','6','7','8','9','0'])

const _prosite_dict = UInt8[ i for i in 1:255 ]
_prosite_dict[Int('(')] = '{'
_prosite_dict[Int(')')] = '}'
_prosite_dict[Int('<')] = '^'
_prosite_dict[Int('>')] = '$'

# Function for the { } syntax for excluding amino acids [ See pattern syntax 4 ]
function _exclude(s::ASCIIString)
  if search(s,'x')!=0
    throw("All possible amino acids are being excluded")
  end
  possibles = IntSet(Int[ i for i in AMINO_20_UPPERCASE ])
  setdiff!(possibles,s.data)
  ascii(UInt8[ i for i in possibles])
end

function _element_parser(elem::ASCIIString)
  firstChar = Char(elem.data[1])
  if length(elem) == 1
    if 'x' == firstChar
      return("[ACDEFGHIKLMNPQRSTVWY]") # [ See Pattern syntax 2 ]
    else
      if in(firstChar, AMINO_20_UPPERCASE)
	return( elem )
      else
	throw("$firstChar isn't in xACDEFGHIKLMNPQRSTVWY")
      end
    end
  else
    if firstChar == '[' && search(elem,'>')!=0
      cap = match( r"\[(.+)>\]$" , elem )
      return("($( cap.captures[1] ).*)?\$") # [ See rare cases on pattern syntax 7 ]
    elseif firstChar!='x' && search(elem,',')!=0
      throw("You can only use a range with 'x'") # [ See Note on Pattern Syntax 6 ]
    elseif search(elem,'{')!=0 # [ See pattern syntax 4 ]
      cap = match( r"(<*)\{(.+)\}([\(>]*.*)" , elem )
      elem = string(cap.captures[1],"[",_exclude(ascii(cap.captures[2])),"]",cap.captures[3])
    else
      elem = replace(elem,'x',"[ACDEFGHIKLMNPQRSTVWY]")
      return(ascii(_prosite_dict[elem.data])) # Replace ( ) to { } and < > to ^ $
    end
  end
end

# Prosite pattern string to regex string
function _prosite2regex(str::ASCIIString)
  # Split elements  [ See Pattern syntax 5 ]
  elements = split(str,'-')
  if length(elements)==1
    # Without '-' can be something like <{C}*> or a sequence [ See extended syntax ]
    if str[1] == '<'
      values = _exclude( ascii( match( r"<\{(.+)\}\*>", str ).captures[1] ) )
      return( "^[$(values)]+\$" )
    else
      return( _element_parser( ascii(str) ) )
    end
  else
    return( string( [ _element_parser( ascii(elem) ) for elem in elements ]... ) )
  end
end

macro prosite_str(pattern, flags...)
  check(pattern,_prosite_alpha)
  out = _prosite2regex(pattern)
  Regex(out,flags...)
end


text:
# overview over all variables used

meteo_UT	# unixtime meteodata
meteo_DT	# datetime meteodata
SWd		# shortwave radiation downwards
WS		# wind strength
WD		# wind direction
WDstdv		# standard deviation of wind direction
Temp 		# temperature at the peak
RH  		# relative humidity at the peak

TEMPstation
RHstation
PRESSUREstation

BC_UT 		# unixtime black carbon
BC_DT 		# datetime black carbon
BC_sum 		# total black carbon

ACSM_UT 	# unixtime ACSM
ACSM_DT		# datetime ACSM
ACSM_org	# total organics in particle
ACSM_sulf	# total sulfates in particle
ACSM_nitr	# total nitrate in particle
ACSM_ammon	# total ammonia in particle
ACSM_chlor	# total chlorine species in particle

SMPS_DT 	# datetime SMPS
SMPS_UT		# unixtime SMPS
SMPS_binsize 	# binsizes SMPS
SMPS_totconc	# total concentration of particles from SMPS
SMPS_binconc	# dNdlog(Dp) SMPS

no3cims_DT	# datetime Nitrate CIMS
no3cims_UT	# unixtime Nitrate CIMS
Nitrophenol 	# Nitrophenol from Nitrate CIMS
SO5trace	# SO5- trace from Nitrate CIMS
MSAtrace 	# MSA from Nitrate CIMS

SRRs		# Source Receptor Relationship Timetraces
SRRnames	# Names of SRRs

ptr3Traces	# all Traces from PTR3
ptr3Times	# Times from PTR3
ptr3Peaks	# all Peaks from PTR3


text:
using InformationDistances
using Documenter

DocMeta.setdocmeta!(InformationDistances, :DocTestSetup, :(using InformationDistances); recursive=true)

makedocs(;
    modules=[InformationDistances],
    authors="Simon Schoelly <sischoel@gmail.com> and contributors",
    repo="https://github.com/simonschoelly/InformationDistances.jl/blob/{commit}{path}#L{line}",
    sitename="InformationDistances.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://simonschoelly.github.io/InformationDistances.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
    doctest=false
)

deploydocs(;
    repo="github.com/simonschoelly/InformationDistances.jl",
)


text:
immutable SwapIndex{S <: AbstractString, TP <: TenorPeriod, B <: BusinessCalendar, C <: BusinessDayConvention, DC <: DayCount, T <: TermStructure} <: InterestRateIndex
  familyName::S
  tenor::TP
  fixingDays::Int
  currency::Currency
  fixingCalendar::B
  fixedLegTenor::Dates.Period
  fixedLegConvention::C
  fixedLegDayCount::DC
  discount::T
  iborIndex::IborIndex
  exogenousDiscount::Bool
  lastSwap::VanillaSwap
  lastFixingDate::Date

  SwapIndex(familyName::S, tenor::TP, fixingDays::Int, currency::Currency, fixingCalendar::B, fixedLegTenor::Dates.Period, fixedLegConvention::C, fixedLegDayCount::DC,
            discount::T, iborIndex::IborIndex, exogenousDiscount::Bool = true) =
            new{S, TP, B, C, DC, T}(familyName, tenor, fixingDays, currency, fixingCalendar, fixedLegTenor, fixedLegConvention, fixedLegDayCount, discount, iborIndex, exogenousDiscount)
end

SwapIndex{S <: AbstractString, TP <: TenorPeriod, B <: BusinessCalendar, C <: BusinessDayConvention, DC <: DayCount, T <: TermStructure}(familyName::S, tenor::TP, settlementDays::Int, currency::Currency, fixingCalendar::B,
          fixedLegTenor::Dates.Period, fixedLegConvention::C, fixedLegDayCount::DC, discount::T, iborIndex::IborIndex) = SwapIndex{S, TP, B, C, DC, T}(familyName, tenor, settlementDays, currency, fixingCalendar, fixedLegTenor,
          fixedLegConvention, fixedLegDayCount, discount, iborIndex)

function EuriborSwapIsdaFixA(tenor::TenorPeriod, forwardingTS::YieldTermStructure, discTS::YieldTermStructure)
  ibor = Int(tenor.period) > 1 ? euribor_index(TenorPeriod(Dates.Month(6)), forwardingTS) : euribor_index(TenorPeriod(Dates.Month(3)), forwardingTS)

  return SwapIndex("EuriborSwapIsdaFixA", tenor, 2, EURCurrency(), TargetCalendar(), Dates.Year(1), ModifiedFollowing(), BondThirty360(), discTS, ibor)
end


text:
using JuMP

##### Simple LP
# Example taken directly from https://jump.readthedocs.org/en/latest/quickstart.html

println("***** Problem 1 *****")

m = Model()
@variable(m, 0 <= x <= 2 )
@variable(m, 0 <= y <= 30 )

@objective(m, Max, 5x + 3*y )
@constraint(m, 1x + 5y <= 3.0 )

print(m)

status = solve(m)

println("Objective value: ", getobjectivevalue(m))
println("x = ", getvalue(x))
println("y = ", getvalue(y))

##### Another simple LP with logical conditions in constraint definition
# Example from https://groups.google.com/forum/#!topic/julia-opt/pMHz-9YHN2o
# Note: This example only works with Julia-v0.4 and later

println(" ")
println("***** Problem 2 *****")

m = Model()

I = 1:5

@variable(m, 0 <= x[I] <= 10)
@variable(m, 0 <= y[I] <= 10)

# This constraint is only considered when i + j are less than or equal to 3
# Note: Only supported by Julia-v0.4 and later
@constraint(m, constr[i=I,j=I; i+j <= 3], x[i] - y[j] == 1)

@objective(m, Min, sum{x[i] + y[i], i=I})

print(m)

status = solve(m)

println(" ")
println("Objective value: ", getobjectivevalue(m))
println("x = ", getvalue(x))
println("y = ", getvalue(y))


text:
"""
Takes a CSV file and makes it into a dataframe with reasonable header, etc.
"""
function importCSV(filename;header=3)
    #Read File as Dataframe
    df = CSV.read(filename,DataFrame,header=header)
    return df
end

"""
Cleans a MACCOR file based on everything wrong I see so far
"""
function clean_df(df)
    df = convertTestTime(df)
    df = convertStepTime(df)
    df = cleanGarbageCols(df)
end

"""
Converts a dataframe with raw time vector to usable date-time format
"""
function convertTestTime(df;start=DateTime(1,1,1,0,0))
    testTime = df.TestTime
    testTimes = split.(testTime)
    newTestTimes = Array{DateTime,1}(undef,length(testTime))
    for (i,time) in enumerate(testTimes)
        testday = time[1]
        #Since d is one letter, the day is 1:end-1
        numTestDays = testday[1:end-1]
        numTestDays = parse(Int,numTestDays)
        testtime = time[2]
        testtime = split(testtime,(':','.'))
        testtime = parse.(Int,testtime)
        testtime = start+Day(numTestDays)+Hour(testtime[1])+Minute(testtime[2])+Second(testtime[3])+Millisecond(testtime[4]*10)
        newTestTimes[i] = testtime
    end
    df.TestTime = newTestTimes
    return df
end

"""This converts the step times of the testing DF's to reasonable times"""
function convertStepTime(df;start=DateTime(1,1,1,0,0))
    stepTime = df.StepTime
    stepTimes = split.(stepTime)
    newStepTimes = Array{DateTime,1}(undef,length(stepTime))
    for (i,time) in enumerate(stepTimes)
        stepday = time[1]
        #Since d is one letter, the day is 1:end-1
        numStepDays = stepday[1:end-1]
        numStepDays = parse(Int,numStepDays)
        steptime = time[2]
        steptime = split(steptime,(':','.'))
        steptime = parse.(Int,steptime)
        steptime = start+Day(numStepDays)+Hour(steptime[1])+Minute(steptime[2])+Second(steptime[3])+Millisecond(steptime[4]*10)
        newStepTimes[i] = steptime
    end
    df.StepTime = newStepTimes
    return df
end

"""This gets rid of cols that are useless"""
function cleanGarbageCols(df)
    df = deepcopy(df) #Don't want to mutate here
    cols = names(df)
    for col in cols
        #My Current Best guess heuristic on what a garbage column looks like, LOL
        if any(ismissing.(df[:,col]))
            @warn "Deleted $col"
            select!(df,Not(col))
            continue
        end
        if all(df[:,col].==0)
            select!(df,Not(col))
        end
    end
    return df
end


text:
using PyPlot

n = 500
n = 500
x = randn(n)  # N(0, 1)
x = exp(x)  # Map x to lognormal
y = randn(n) + 2.0  # N(2, 1)
z = randn(n) + 4.0  # N(4, 1)

fig, ax = subplots()
ax[:boxplot]([x y z])
ax[:set_xticks]((1, 2, 3))
ax[:set_ylim](-2, 14)
ax[:set_xticklabels]((L"$X$", L"$Y$", L"$Z$"), fontsize=16)
plt.show()


text:
using GLFW, ModernGL

# whenever the window size changed (by OS or user resize) this callback function executes
function framebuffer_size_callback(window::GLFW.Window, width::Cint, height::Cint)
    # make sure the viewport matches the new window dimensions; note that width and
    # height will be significantly larger than specified on retina displays.
	glViewport(0, 0, width, height)
end

# process all input: query GLFW whether relevant keys are pressed/released this frame and react accordingly
function processInput(window::GLFW.Window)
    GLFW.GetKey(window, GLFW.KEY_ESCAPE) == GLFW.PRESS && GLFW.SetWindowShouldClose(window, true)
end

# misc. config
const SCR_WIDTH = 800
const SCR_HEIGHT = 600

const vertexShaderSource = """
	#version 330 core
	layout (location = 0) in vec3 aPos;
	void main(void)
	{
	    gl_Position = vec4(aPos.x, aPos.y, aPos.z, 1.0);
	}"""

const fragmentShaderSource = """
	#version 330 core
	out vec4 FragColor;
	void main(void)
	{
	    FragColor = vec4(1.0f, 0.5f, 0.2f, 1.0f);
	}"""

GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, 3)
GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, 3)
GLFW.WindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE)
GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, GL_TRUE)

# create window
window = GLFW.CreateWindow(SCR_WIDTH, SCR_HEIGHT, "LearnOpenGL")
window == C_NULL && error("Failed to create GLFW window.")
GLFW.MakeContextCurrent(window)
GLFW.SetFramebufferSizeCallback(window, framebuffer_size_callback)

# build and compile shader program
vertexShader = glCreateShader(GL_VERTEX_SHADER)
glShaderSource(vertexShader, 1, Ptr{GLchar}[pointer(vertexShaderSource)], C_NULL)
glCompileShader(vertexShader)
successRef = Ref{GLint}(-1)
glGetShaderiv(vertexShader, GL_COMPILE_STATUS, successRef)
if successRef[] != GL_TRUE
    infoLog = Vector{GLchar}(512)
    glGetShaderInfoLog(vertexShader, 512, C_NULL, infoLog)
    error("ERROR::SHADER::VERTEX::COMPILATION_FAILED\n", String(infoLog))
end

fragmentShader = glCreateShader(GL_FRAGMENT_SHADER)
glShaderSource(fragmentShader, 1, Ptr{GLchar}[pointer(fragmentShaderSource)], C_NULL)
glCompileShader(fragmentShader)
glGetShaderiv(vertexShader, GL_COMPILE_STATUS, successRef)
if successRef[] != GL_TRUE
    infoLog = Vector{GLchar}(512)
    glGetShaderInfoLog(fragmentShader, 512, C_NULL, infoLog)
    error("ERROR::SHADER::FRAGMENT::COMPILATION_FAILED\n", String(infoLog))
end

# link shaders
shaderProgram = glCreateProgram()
glAttachShader(shaderProgram, vertexShader)
glAttachShader(shaderProgram, fragmentShader)
glLinkProgram(shaderProgram)
glGetProgramiv(shaderProgram, GL_LINK_STATUS, successRef)
if successRef[] != GL_TRUE
    infoLog = Vector{GLchar}(512)
    glGetProgramInfoLog(shaderProgram, 512, C_NULL, infoLog)
    error("ERROR::SHADER::PROGRAM::LINKING_FAILED\n", String(infoLog))
end
glDeleteShader(vertexShader)
glDeleteShader(fragmentShader)

# set up vertex data (and buffer(s)) and configure vertex attributes
firstTriangle = GLfloat[-0.9,  0.5, 0.0,    # left
                         0.0, -0.5, 0.0,    # right
                         0.45, 0.5, 0.0]    # top
secondTriangle = GLfloat[0.0, -0.5, 0.0,    # left
						 0.9, -0.5, 0.0,    # right
					     0.45, 0.5, 0.0]    # top
vbos = zeros(GLuint, 2)
vaos = zeros(GLuint, 2)
glGenVertexArrays(2, vaos)
glGenBuffers(2, vbos)
# first triangle setup
glBindVertexArray(vaos[1])
glBindBuffer(GL_ARRAY_BUFFER, vbos[1])
glBufferData(GL_ARRAY_BUFFER, sizeof(firstTriangle), firstTriangle, GL_STATIC_DRAW)
glVertexAttribPointer(0, 3, GL_FLOAT, GL_FALSE, 3 * sizeof(GLfloat), Ptr{Cvoid}(0))
glEnableVertexAttribArray(0)
# second triangle setup
glBindVertexArray(vaos[2])
glBindBuffer(GL_ARRAY_BUFFER, vbos[2])
glBufferData(GL_ARRAY_BUFFER, sizeof(secondTriangle), secondTriangle, GL_STATIC_DRAW)
glVertexAttribPointer(0, 3, GL_FLOAT, GL_FALSE, 3 * sizeof(GLfloat), Ptr{Cvoid}(0))
glEnableVertexAttribArray(0)

# uncomment this call to draw in wireframe polygons.
# glPolygonMode(GL_FRONT_AND_BACK, GL_LINE)

# render loop
while !GLFW.WindowShouldClose(window)
    processInput(window)

    # render
    glClearColor(0.2, 0.3, 0.3, 1.0)
    glClear(GL_COLOR_BUFFER_BIT)

	glUseProgram(shaderProgram)
	# draw first triangle
    glBindVertexArray(vaos[1])
	glDrawArrays(GL_TRIANGLES, 0, 3)

	# draw  second triangle
    glBindVertexArray(vaos[2])
	glDrawArrays(GL_TRIANGLES, 0, 3)

    # swap buffers and poll IO events
    GLFW.SwapBuffers(window)
    GLFW.PollEvents()
end

# optional: de-allocate all resources once they've outlived their purpose
glDeleteVertexArrays(2, vaos)
glDeleteBuffers(2, vbos)

GLFW.DestroyWindow(window)

max_stars_repo_name: Gnimuc/Videre

max_issues_repo_name: Gnimuc/Videre

max_forks_repo_name: Gnimuc/Videre


text:
using PyPlot
using PyCall
using DelimitedFiles
using Statistics
mpl = pyimport("tikzplotlib")

res = readdlm("result.txt")
d = Dict()
for i = 1:size(res,1)
    if haskey(d, res[i,1])
        push!(d[res[i,1]], res[i,2])
    else
        d[res[i,1]] = [res[i,2]]
    end
end

arr = [0.0,0.002,0.004,0.006,0.008,0.01,0.03,0.05,0.06, 0.07, 0.075,0.08,0.09,0.1]
s = zeros(length(arr))
m = zeros(length(arr))

for (ix,k) in enumerate(arr)
    s[ix] = std(d[k])
    m[ix] = mean(d[k])
end

semilogy(arr, m)
fill_between(arr, m - s, m+s, alpha=0.5)
xlabel("\$\\sigma_{\\mathrm{noise}}\$")
ylabel("Error")
mpl.save("poroerror.tex")

max_stars_repo_name: ADCMEMarket/ADCMEImages

max_issues_repo_name: ADCMEMarket/ADCMEImages

max_forks_repo_name: ADCMEMarket/ADCMEImages


text:
using Images, Plots
include("util.jl")

function parse_line(line::String)::Dict{Symbol,UnitRange{Int64}}
    parts = split(line, ", ")
    split_parts = [split(part, "=") for part in parts]
    x = 1:1
    y = 1:1
    for part in split_parts
        if part[1] == "x"
            if occursin("..", part[2])
                sides = split(part[2], "..")
                x = parse(Int64, sides[1]):parse(Int64, sides[2])
            else
                x = parse(Int64, part[2]):parse(Int64,part[2])
            end
        elseif part[1] == "y"
            if occursin("..", part[2])
                sides = split(part[2], "..")
                y = parse(Int64, sides[1]):parse(Int64,sides[2])
            else
                y = parse(Int64, part[2]):parse(Int64,part[2])
            end
        end
    end
    return Dict(:x => x, :y => y)
end

function get_bounds(parsed_lines::Array{Dict{Symbol,UnitRange{Int64}}})::Array{Int64}
    xmin = minimum([minimum(d[:x]) for d in parsed_lines])
    xmax = maximum([maximum(d[:x]) for d in parsed_lines])
    ymin = minimum([minimum(d[:y]) for d in parsed_lines])
    ymax = maximum([maximum(d[:y]) for d in parsed_lines])
    return [xmin-1, xmax+1, ymin, ymax]
end

function make_grid(parsed_lines::Array{Dict{Symbol,UnitRange{Int64}}})::Array{Char,2}
    bounds = get_bounds(parsed_lines)
    grid = map(x->'.', zeros(bounds[4],bounds[2]))
    for line in parsed_lines
        for x in line[:x]
            for y in line[:y]
                grid[y,x] = '#'
            end
        end
    end
    return grid
end

function read_input()
    lines = getlines("inputs/17.txt")
    parsed_lines = parse_line.(lines)
    grid = make_grid(parsed_lines)
    return parsed_lines, grid
end

function hasWall(x::Int64,y::Int64,grid::Array{Char,2},direction::Int64)
    if grid[y,x+direction] == '#' && grid[y+1,x] in "#~"
        return true
    elseif grid[y,x+direction] == '.'
        return false
    elseif grid[y+1,x] in ".|"
        return false
    else
        return hasWall(x+direction,y,grid,direction)
    end
end

function hasBothWalls(x::Int64,y::Int64,grid::Array{Char,2})
    if hasWall(x,y,grid,1) && hasWall(x,y,grid,-1)
        return true
    else
        return false
    end
end

function fillside(x::Int64,y::Int64,grid::Array{Char,2},direction::Int64)
    if grid[y,x] == '#'
        return
    else
        grid[y,x] = '~'
        fillside(x+direction,y,grid,direction)
    end
end

function fillLevel(x::Int64,y::Int64,grid::Array{Char,2})
    fillside(x,y,grid,1)
    fillside(x,y,grid,-1)
end

function ffill(x::Int64,y::Int64,grid::Array{Char,2},bounds::Array{Int64,1})
    #println("plotting")
    #plotgrid(grid)
    #sleep(0.025)
    if y < bounds[4]
        if grid[y+1,x] == '.'
            grid[y+1,x] = '|'
            ffill(x,y+1,grid,bounds)
        end
        if grid[y+1,x] in "~#" && grid[y,x+1] == '.'
            grid[y,x+1] = '|'
            ffill(x+1,y,grid,bounds)
        end
        if grid[y+1,x] in "~#" && grid[y,x-1] == '.'
            grid[y,x-1] = '|'
            ffill(x-1,y,grid,bounds)
        end
        if hasBothWalls(x,y,grid)
            fillLevel(x,y,grid)
        end
    end
    return
end

function color_cell(x::Char)::RGB
    if x == '~'
        return RGB(0,0,1)
    elseif x == '|'
        return RGB(0,0,0.5)
    elseif x == '#'
        return RGB(1,0,0)
    else
        return RGB(1,1,1)
    end
end

function plotgrid(grid)
    plt = plot(color_cell.(grid[1:40,480:520]))
    gui(plt)
end

parsed_lines, grid = read_input()
bounds = get_bounds(parsed_lines)
ffill(500,1,grid,bounds)
i = color_cell.(grid[1:40,480:520])

function part_1()
    count = 0
    for row in 1:size(grid)[1]
        for col in 1:size(grid)[2]
            if grid[row,col] in "~|"
                count += 1
            end
        end
    end
    return count
end

c = part_1()
println(c)

function part_2()
    count = 0
    for row in 1:size(grid)[1]
        for col in 1:size(grid)[2]
            if grid[row,col] in "~"
                count += 1
            end
        end
    end
    return count
end

d = part_2()
println(d)


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

"""
    Domain{Dim,T}

A domain is an indexable collection of geometries (i.e. mesh) or points
where each element can be mapped to a feature vector. For example, a
collection of polygonal areas representing the states of a country can
be seen as a domain. In this case, the features can be the name of the
state, the total population of the state, or any other quantity stored
in each polygon.
"""
abstract type Domain{Dim,T} end

"""
    getindex(domain, ind)

Return the `ind`-th element in the `domain`.
"""
Base.getindex(domain::Domain, ind::Int)

"""
    nelements(domain)

Return the number of elements in the `domain`.
"""
function nelements end

# ----------
# FALLBACKS
# ----------

"""
    embeddim(domain)

Return the number of dimensions of the space where the `domain` is embedded.
"""
embeddim(::Type{<:Domain{Dim,T}}) where {Dim,T} = Dim
embeddim(domain::Domain) = embeddim(typeof(domain))

"""
    coordtype(domain)

Return the machine type of each coordinate used to describe the `domain`.
"""
coordtype(::Type{<:Domain{Dim,T}}) where {Dim,T} = T
coordtype(domain::Domain) = coordtype(typeof(domain))

"""
    coordinates!(buff, domain, ind)

Compute the coordinates `buff` of the centroid of the `ind`-th element in
the `domain` in place.
"""
function coordinates!(buff, domain::Domain{Dim}, ind::Int) where {Dim}
  x = coordinates(centroid(domain[ind]))
  @inbounds for i in 1:Dim
    buff[i] = x[i]
  end
  buff
end

function coordinates!(buff, domain::Domain, inds::AbstractVector{Int})
  for j in eachindex(inds)
    coordinates!(view(buff,:,j), domain, inds[j])
  end
  buff
end

function coordinates(domain::Domain, ind::Int)
  buff = MVector{embeddim(domain),coordtype(domain)}(undef)
  coordinates!(buff, domain, ind)
end

function coordinates(domain::Domain, inds::AbstractVector{Int})
  buff = Matrix{coordtype(domain)}(undef, embeddim(domain), length(inds))
  coordinates!(buff, domain, inds)
end

Base.eltype(domain::Domain) = eltype([domain[i] for i in 1:nelements(domain)])

# -----------
# IO METHODS
# -----------

function Base.show(io::IO, domain::Domain{Dim,T}) where {Dim,T}
  nelm = nelements(domain)
  name = nameof(typeof(domain))
  print(io, "$nelm $name{$Dim,$T}")
end

function Base.show(io::IO, ::MIME"text/plain", domain::Domain)
  println(io, domain)
  N = nelements(domain)
  I, J = N > 10 ? (5, N-4) : (N, N+1)
  lines = [["  └─$(domain[i])" for i in 1:I]
           (N > 10 ? ["  ⋮"] : [])
           ["  └─$(domain[i])" for i in J:N]]
  print(io, join(lines, "\n"))
end


text:
using MatrixNetworks

function _kcliques(A::SparseMatrixCSC{Int64,Int64}, k::Int64)
    A = min.(A, 1)
    A -= spdiagm(diag(A))
    # "global" variables for the recursive algorithm
    n = size(A, 2)
    adjacency_lists = Vector{Vector{Int64}}(n)
    for v = 1:n
        adjacency_lists[v] = filter(x -> x != v, find(A[:, v]))
    end
    C = Int64[]
    labels = k * ones(Int64, n)
    counts = zeros(Int64, n)

    # Get all neighbors of a node with a given label
    function neighbors(node::Int64, label::Int64)
        index = 1
        for nbr in adjacency_lists[node]
            nbr_label = labels[nbr]
            if nbr_label == label; index += 1; end
            if nbr_label  > label; break; end
        end
        return adjacency_lists[node][1:(index - 1)]
    end

    # Get the number of neighbors of a node with a given label
    function num_neighbors(node::Int64, label::Int64)
        index = 1
        for nbr in adjacency_lists[node]
            nbr_label = labels[nbr]
            if nbr_label == label; index += 1; end
            if nbr_label  > label; break; end
        end
        return index - 1
    end

    # Get the nodes in U in sorted order that need to be processed
    function nodes_to_process(U::Vector{Int64}, label::Int64)
        sorted_nbrs = sort([(num_neighbors(v, label), v) for v in U], rev=true)
        return [x[2] for x in sorted_nbrs if x[1] > 1]
    end

    # Adjust adjacency lists of neighbors of nodes in U to put neighbors with a
    # given label in the front.
    function move_front(U::Vector{Int64}, label::Int64)
        for v in U
            front_nbrs = Int64[]
            back_nbrs  = Int64[]
            for nbr in adjacency_lists[v]
                if labels[nbr] == label; push!(front_nbrs, nbr)
                else                     push!(back_nbrs, nbr)
                end
            end
            num_front = length(front_nbrs)
            total = length(adjacency_lists[v])
            adjacency_lists[v][1:num_front] = front_nbrs
            adjacency_lists[v][(num_front + 1):total] = back_nbrs
        end
    end

    # Adjust adjacency lists of neighbors of nodes in U to put node vi directly
    # after nodes with a given label.
    function move_after(U::Vector{Int64}, vi::Int64, label::Int64)
        for v in U
            front_nbrs = Int64[]
            for nbr in adjacency_lists[v]
                if nbr == vi; continue; end
                if labels[nbr] <= label; push!(front_nbrs, nbr)
                else                     break
                end
            end
            num_front = length(front_nbrs)
            adjacency_lists[v][1:num_front] = front_nbrs
            adjacency_lists[v][num_front + 1] = vi
        end
    end

    # Process all current cliques, which each consist of:
    #   1. the nodes in the global stack C
    #   2. a node in U
    #   3. a neighbor of u with label 2
    function process_cliques(U::Vector{Int64})
        C_cnt = 0
        for v in U
            v_cnt = 0
            for nbr in neighbors(v, 2)
                # Only count (nbr, v) edge once.
                if nbr > v
                    counts[nbr] += 1
                    v_cnt += 1
                end
            end
            counts[v] += v_cnt
            C_cnt += v_cnt
        end
        counts[C] += C_cnt
    end

    # Recursive algorithm of Chiba and Nishizeki.
    function rcliques(U::Vector{Int64}, r::Int64)
        # Check base case
        if r == 2
            process_cliques(U)
            return
        end

        for (ind, v) in enumerate(nodes_to_process(U, r))
            if r == k
                print(@sprintf("%d of %d (%d-cliques) \r", ind, length(U), r))
                flush(STDOUT)
                if ind == length(U); println(""); end
            end
            # Get neighborhood of v
            Up = neighbors(v, r)
            labels[Up] = r - 1
            move_front(Up, r - 1)

            # Recurse on neighborhood of node v
            push!(C, v)
            rcliques(Up, r - 1)
            pop!(C)

            # Restore neighborhood
            labels[Up] = r

            # Eliminate v
            labels[v] = r + 1
            move_after(Up, v, r)
        end
    end

    rcliques(find(vec(sum(A, 2)) .>= (k - 1)), k)
    return counts
end

"""
Count k-cliques in the graph A using the Chiba-Nishizeki algorithm.

Returns a vector c, where c[v] is the number of k-cliques containing node v
in the graph.
"""
function kcliques(A::SparseMatrixCSC{Int64,Int64}, k::Int64)
    A = min.(A, 1)
    A -= spdiagm(diag(A))
    clique_counts = zeros(Int64, size(A, 1))
    (d, rt) = corenums(A)
    inds = find(d .>= (k - 1))
    counts = _kcliques(A[inds, inds], k)
    clique_counts[inds] = counts
    return clique_counts
end


text:
add_verbose_scope(:LLL)
add_assert_scope(:LLL)

################################################################################
#
#  Special functions for real fields and quadratic fields
#
################################################################################

function _lll_gram(A::NfOrdIdl)
  K = nf(order(A))
  @assert istotally_real(K)
  g = trace_matrix(A)
  @hassert :LLL 1 !iszero(det(g))
  @hassert :LLL 1 isposdef(g)
  l, t = lll_gram_with_transform(g)
  return FakeFmpqMat(l, fmpz(1)), t::fmpz_mat
end

function _lll_quad(A::NfOrdIdl)
  K = nf(order(A))
  @assert degree(K) == 2 && discriminant(order(A)) < 0
  b = basis(A)
  a1 = 2*numerator(norm(b[1]))
  a2 = 2*numerator(norm(b[2]))
  a12 = numerator(trace(b[1] * conjugate_quad(K(b[2]))))
  g = matrix(FlintZZ, 2, 2, [a1, a12, a12, a2])
  @hassert :LLL 1 isposdef(g)
  l, t = lll_gram_with_transform(g)
  return FakeFmpqMat(l, fmpz(1)), t::fmpz_mat
end

function _lll_CM(A::NfOrdIdl)
  OK = order(A)
  @vprint :LLL 3 "Reduction\n"
  M = _minkowski_matrix_CM(OK)
  @vtime :LLL 3 BM, T = lll_with_transform(basis_matrix(A, copy = false), lll_ctx(0.3, 0.51))
  g = BM*M*transpose(BM)
  @hassert :LLL 1 isposdef(g)
  @vtime :LLL 3 l, t = lll_gram_with_transform(g)
  return FakeFmpqMat(l, fmpz(1)), t*T::fmpz_mat
end

################################################################################
#
#  lll for ideals
#
################################################################################

function lll(A::NfOrdIdl, v::fmpz_mat = zero_matrix(FlintZZ, 1, 1); prec::Int = 100)

  K = nf(order(A))

  if iszero(v)
    if istotally_real(K)
      #in this case the gram-matrix of the minkowski lattice is the trace-matrix
      #which is exact.
      return _lll_gram(A)
    elseif degree(K) == 2 && discriminant(order(A)) < 0
       #in this case the gram-matrix of the minkowski lattice is related to the
      #trace-matrix which is exact.
      return _lll_quad(A)
    end
    if iscm_field_known(K) || isautomorphisms_known(K)
      if iscm_field(K)[1]
        return _lll_CM(A)
      end
    end
  end
  #General case. _lll could fail, we need a try and catch in a loop
  local t::fmpz_mat
  local l::FakeFmpqMat
  while true
    if prec > 2^18
      error("Something wrong in short_elem")
    end
    try
      l, t = _lll(A, v, prec = prec)
      break
    catch e
      if !(e isa LowPrecisionLLL || e isa InexactError)
        rethrow(e)
      end
    end
    prec = 2 * prec
  end
  return l, t

end

function _lll(A::NfOrdIdl, v::fmpz_mat = zero_matrix(FlintZZ, 1, 1); prec::Int = 100)
  K = nf(order(A))
  n = degree(order(A))
  prec = max(prec, 4*n)

  ctx1 = lll_ctx(0.5, 0.51)
  l, t1 = lll_with_transform(basis_matrix(A, copy = false), ctx1)

  if iszero(v)
    d = minkowski_gram_mat_scaled(order(A), prec)
    ccall((:fmpz_mat_mul, libflint), Nothing, (Ref{fmpz_mat}, Ref{fmpz_mat}, Ref{fmpz_mat}), d, d, transpose(l))
    ccall((:fmpz_mat_mul, libflint), Nothing, (Ref{fmpz_mat}, Ref{fmpz_mat}, Ref{fmpz_mat}), d, l, d)
    g = zero_matrix(FlintZZ, n, n)
    den = fmpz(1)
    sv = fmpz(0)
  else
    c = minkowski_matrix(nf(order(A)), prec) ## careful: current iteration
                                          ## c is NOT a copy, so don't change.
    b = l*basis_matrix(order(A), copy = false)

    rt_c = roots_ctx(K)
    if !isdefined(rt_c, :cache_z1)
      rt_c.cache_z1 = zero_matrix(FlintZZ, n, n)
      rt_c.cache_z2 = zero_matrix(FlintZZ, n, n)
    end

    d::fmpz_mat = rt_c.cache_z1
    g::fmpz_mat = rt_c.cache_z2

    round_scale!(g, c, prec)

    @v_do :ClassGroup 2 println("using inf val", v)
    c = deepcopy(c)
    mult_by_2pow_diag!(c, v)
    sv = max(fmpz(0), sum(v[1,i] for i=1:ncols(l)))

    round_scale!(d, c, prec)
    ccall((:fmpz_mat_mul, libflint), Nothing, (Ref{fmpz_mat}, Ref{fmpz_mat},  Ref{fmpz_mat}), g, (b.num), d)
    den = b.den

    ccall((:fmpz_mat_gram, libflint), Nothing, (Ref{fmpz_mat}, Ref{fmpz_mat}), d, g)
    shift!(d, -prec)
  end

  prec = div(prec, 2)
  shift!(d, -prec)  #TODO: remove?

  for i=1:n
    fmpz_mat_entry_add_ui!(d, i, i, UInt(nrows(d)))
  end

  ctx = Nemo.lll_ctx(0.99, 0.51, :gram)

  ccall((:fmpz_mat_one, libflint), Nothing, (Ref{fmpz_mat}, ), g)
  ccall((:fmpz_lll, libflint), Nothing, (Ref{fmpz_mat}, Ref{fmpz_mat}, Ref{Nemo.lll_ctx}), d, g, ctx)

  l, t = d, g
  ## test if entries in l are small enough, if not: increase precision
  ## or signal that prec was too low
  @v_do :ClassGroup 2 printstyled("lll basis length profile\n", color=:green);
  @v_do :ClassGroup 2 for i = 1:nrows(l)
    print(Float64(div(l[i,i], fmpz(2)^prec*den*den)*1.0), " : ")
  end
  @v_do :ClassGroup 2 println("")
  if nbits(maximum(abs, t)) >  div(prec, 2)
    @v_do :ClassGroup 2 printstyled("lll trafo too large\n", color=:red);
    throw(LowPrecisionLLL())
  end
  ## lattice has lattice disc = order_disc * norm^2
  ## lll needs to yield a basis sth
  ## l[1,1] = |b_i|^2 <= 2^((n-1)/2) disc^(1/n)
  ## and prod(l[i,i]) <= 2^(n(n-1)/2) disc
  n = nrows(l)
  disc = abs(discriminant(order(A)))*norm(A)^2 * den^(2*n) * fmpz(2)^(2*sv)
  di = root(disc, n)+1
  di *= fmpz(2)^(div(n+1,2)) * fmpz(2)^prec

  if compare_index(l, 1, 1, di) > 0
    @v_do :ClassGroup 2 printstyled("LLL basis too large\n", color = :red);
    @v_do :ClassGroup 3 println("bound is ", di, " value at ", 1, " is ", l[1,1]);
    throw(LowPrecisionLLL())
  end
  pr = prod_diagonal(l)
  if pr > fmpz(2)^(div(n*(n-1), 2)) * disc * fmpz(2)^(n*prec)
    @v_do :ClassGroup 2 printstyled("LLL basis too large\n", color = :red);
    @v_do :ClassGroup 2 println("prod too large: ", pr, " > 2^(n(n-1)/2) disc = ", fmpz(2)^(div(n*(n-1), 2)) * disc * fmpz(2)^(n*prec));
    throw(LowPrecisionLLL())
  end

  return FakeFmpqMat(deepcopy(l), fmpz(2)^prec), t*t1
end

###############################################################################
#
#  LLL computation for orders
#
###############################################################################

@doc Markdown.doc"""
    lll(M::NfAbsOrd) -> NfAbsOrd

The same order, but with the basis now being LLL reduced wrt. the Minkowski metric.
"""
function lll(M::NfAbsOrd; prec::Int = 100)

  if isdefined(M, :lllO)
    return M.lllO::typeof(M)
  end
  K = nf(M)

  if istotally_real(K)
    On =  _lll_gram(M)
    M.lllO = On
    return On::typeof(M)
  end

  if degree(K) == 2
    On = _lll_quad(M)
    M.lllO = On
    return On::typeof(M)
  end

  if iscm_field_known(K) || isautomorphisms_known(K)
    fl, f_conj = iscm_field(K)
    if fl
      On = _lll_CM(M)
      M.lllO = On
      return On::typeof(M)
    end
  end

  return _lll(M, prec)
end

#for totally real field, the T_2-Gram matrix is the trace matrix, hence exact.
function _lll_gram(M::NfAbsOrd)
  K = nf(M)
  @assert istotally_real(K)
  g = trace_matrix(M)
  w = lll_gram_with_transform(g)[2]
  On = NfAbsOrd(K, w*basis_matrix(M, copy = false))
  On.ismaximal = M.ismaximal
  if isdefined(M, :index)
    On.index = M.index
  end
  if isdefined(M, :disc)
    On.disc = M.disc
  end
  if isdefined(M, :gen_index)
    On.gen_index = M.gen_index
  end
  return On
end

function _minkowski_matrix_CM(M::NfAbsOrd)
  if isdefined(M,  :minkowski_gram_CMfields)
    return M.minkowski_gram_CMfields
  end
  B = basis(M, nf(M))
  g = _exact_minkowski_matrix(B)
  M.minkowski_gram_CMfields = g
  return g
end

function _exact_minkowski_matrix(B::Vector{T}) where T <: NumFieldElem
  K = parent(B[1])
  if istotally_real(K)
    return trace_matrix(B)
  else
    return _minkowski_via_approximation(B)
  end
end

function _minkowski_via_approximation(B::Vector{T}) where T <: NumFieldElem
  K = parent(B[1])
  n = length(B)
  g = zero_matrix(FlintZZ, n, n)
  prec = 16
  imgs = Vector{Vector{arb}}(undef, n)
  for i = 1:n
    imgs[i] = minkowski_map(B[i], prec)
  end
  i = 1
  t = arb()
  while i <= n
    j = i
    while j <= n
      el = imgs[i][1]*imgs[j][1]
      for k = 2:n
        mul!(t, imgs[i][k], imgs[j][k])
        add!(el, el, t)
      end
      fl, r = unique_integer(el)
      if fl
        g[i, j] = r
        g[j, i] = r
        j += 1
      else
        prec *= 2
        for k = i:n
          imgs[k] = minkowski_map(B[k], prec)
        end
      end
    end
    i += 1
  end
  return g
end

function trace_matrix(b::Vector{T}) where T <: NumFieldElem
  K = parent(b[1])
  n = absolute_degree(K)
  g = zero_matrix(FlintZZ, n, n)
  aux = K()
  for i = 1:n
    mul!(aux, b[i], b[i])
    t = absolute_tr(aux)
    @assert isinteger(t)
    g[i, i] = numerator(t)
    for j in (i + 1):n
      mul!(aux, b[i], b[j])
      t = absolute_tr(aux)
      @assert isinteger(t)
      g[i, j] = numerator(t)
      g[j, i] = numerator(t)
    end
  end
  return g
end

function _lll_CM(M::NfAbsOrd)
  K = nf(M)
  g = _minkowski_matrix_CM(M)
  @vprint :LLL 1 "Now LLL\n"
  @hassert :LLL 1 isposdef(g)
  w = lll_gram_with_transform(g)[2]
  On = NfAbsOrd(K, w*basis_matrix(M, copy = false))
  On.ismaximal = M.ismaximal
  if isdefined(M, :index)
    On.index = M.index
  end
  if isdefined(M, :disc)
    On.disc = M.disc
  end
  if isdefined(M, :gen_index)
    On.gen_index = M.gen_index
  end
  return On
end

function _lll_quad(M::NfAbsOrd)
  K = nf(M)
  b = basis(M)
  a1 = 2*numerator(norm(b[1]))
  a2 = 2*numerator(norm(b[2]))
  a12 = numerator(trace(b[1] * conjugate_quad(K(b[2]))))
  g = matrix(FlintZZ, 2, 2, fmpz[a1, a12, a12, a2])
  @hassert :ClassGroup 1 isposdef(g)
  w = lll_gram_with_transform(g)[2]
  On = NfAbsOrd(K, w*basis_matrix(M, copy = false))
  On.ismaximal = M.ismaximal
  if isdefined(M, :index)
    On.index = M.index
  end
  if isdefined(M, :disc)
    On.disc = M.disc
  end
  if isdefined(M, :gen_index)
    On.gen_index = M.gen_index
  end
  return On
end

function _lll(M::NfAbsOrd, prec::Int)

  K = nf(M)
  n = degree(K)
  M1 = _ordering_by_T2(M)
  prec = max(prec, 10*n)
  prec = max(prec, 100 + 25*div(degree(M), 3) + Int(round(log(abs(discriminant(M))))))

  if n > 10
    if n > 100
      prec, M1 = lll_precomputation(M1, prec)
    end
    prec, M1 = lll_precomputation(M1, prec)
  end
  M1, prec = _lll_with_parameters(M1, (0.73, 0.51), prec)
  M1 = _lll_with_parameters(M1, (0.99, 0.51), prec)[1]
  M.lllO = M1
  return M1
end

function _ordering_by_T2(M::NfAbsOrd, prec::Int = 32)

  K = nf(M)
  B = basis(M, K)
  ints = fmpz[lower_bound(t2(x, prec), fmpz) for x in B]
  p = sortperm(ints)
  On = NfAbsOrd(B[p])
  On.ismaximal = M.ismaximal
  if isdefined(M, :index)
    On.index = M.index
  end
  if isdefined(M, :disc)
    On.disc = M.disc
  end
  if isdefined(M, :gen_index)
    On.gen_index = M.gen_index
  end
  return On
end

function subsets_it(n::Int, k::Int)
  if n == k
    return (Int[i for i = 1:n])
  end
  if k == 1
    return ([Int[i] for i = 1:n])
  end
  res = subsets_it(n-1, k-1)
  res1 = (push!(x, n) for x in res)
  res2 = subsets_it(n-1, k)
  res3 = [res1, res2]
  return (x for y in res3 for x in y)
end

#Inefficient, but at least it works.
function subsets(n::Int, k::Int)
  if n == k
    return Vector{Int}[Int[i for i = 1:n]]
  end
  if k == 1
    return Vector{Int}[Int[i] for i = 1:n]
  end
  res = subsets(n-1, k-1)
  for x in res
    push!(x, n)
  end
  append!(res, subsets(n-1, k))
  return res
end

function subsets(v::Vector{T}, k::Int) where T
  indices = subsets(length(v), k)
  res = Vector{T}[]
  for i in indices
    si = Vector{T}(undef, k)
    ind = 1
    for j in i
      si[ind] = v[j]
      ind += 1
    end
    push!(res, si)
  end
  return res
end

function _has_trivial_intersection(v::Vector{Vector{Int}}, V::Vector{Vector{Int}})
  for z in v
    for w in V
      if !isempty(intersect(z, w))
        return false
      end
    end
  end
  return true
end

function lll_precomputation(M::NfAbsOrd, prec::Int, nblocks::Int = 4)
  n = degree(M)
  K = nf(M)
  dimension_blocks = div(n, nblocks)
  blocks = Vector{Int}[]
  for i = 1:nblocks-1
    int = (dimension_blocks*(i-1)+1):dimension_blocks*i
    push!(blocks, collect(int))
  end
  int = (dimension_blocks*(nblocks-1)+1):n
  push!(blocks, collect(int))
  g = identity_matrix(FlintZZ, n)
  new_prec = prec
  to_do = subsets(blocks, 2)
  done = falses(length(to_do))
  blocks_selection = Vector{Int}[]
  while !all(done)
    block = 1
    while block < length(to_do)+1 && (done[block] || !_has_trivial_intersection(to_do[block], blocks_selection))
      block += 1
    end
    if block == length(to_do)+1
      blocks_selection = Vector{Int}[]
      On = NfAbsOrd(K, g*basis_matrix(M, copy = false))
      On.ismaximal = M.ismaximal
      if isdefined(M, :index)
      On.index = M.index
      end
      if isdefined(M, :disc)
        On.disc = M.disc
      end
      if isdefined(M, :gen_index)
        On.gen_index = M.gen_index
      end
      M = On
      g = identity_matrix(FlintZZ, n)
      continue
    end
    indices = vcat(to_do[block][1], to_do[block][2])
    @vprint :LLL 3 "Simplifying block $(block)\n"
    prec, g1 = _lll_sublattice(M, indices, prec = prec)
    _copy_matrix_into_matrix(g, indices, indices, g1)
    done[block] = true
    push!(blocks_selection, indices)
    @vprint :LLL 3 "Precision: $(prec)\n"
  end
  @vprint :LLL 3 "Precomputation finished with precision $(prec)\n"
  return prec, M
end

function _lll_sublattice(M::NfAbsOrd, u::Vector{Int}; prec = 100)
  K = nf(M)
  n = degree(M)
  l = length(u)
  @vprint :LLL 3 "Block of dimension $(l)\n"
  prec = max(prec, 10*n)
  local g::fmpz_mat

  bas = basis(M, K)[u]
  profile_sub = nbits(prod(Hecke.upper_bound(fmpz, t2(x)) for x in bas))
  @vprint :LLL 3 "Starting with profile $(profile_sub)\n"
  while true
    local d::fmpz_mat
    @vprint :LLL 3 "Computing Minkowski matrix\n"
    while true
      @vprint :LLL 3 "Precision: $(prec)\n"
      try
        d = minkowski_gram_mat_scaled(M, prec)
        break
      catch e
        prec = prec*2
      end
    end
    @vprint :LLL 3 "Minkowski matrix computed\n"
    g = identity_matrix(FlintZZ, l)
    d1 = sub(d, u, u)
    prec = div(prec, 2)
    shift!(d1, -prec)  #TODO: remove?
    for i=1:l
      fmpz_mat_entry_add_ui!(d1, i, i, UInt(l))
    end
    ctx = Nemo.lll_ctx(0.99, 0.51, :gram)
    @vtime :LLL 3 ccall((:fmpz_lll, libflint), Nothing, (Ref{fmpz_mat}, Ref{fmpz_mat}, Ref{Nemo.lll_ctx}), d1, g, ctx)

    if nbits(maximum(abs, g)) <= div(prec, 2)
      prec *= 2
      break
    end
    @vprint :LLL 3 "Still in the loop\n"
    prec *= 4
  end
  @vprint :LLL 3 "Computing the profile of the new basis \n"
  new_basis = g*basis_matrix(bas, FakeFmpqMat)
  els = elem_type(K)[elem_from_mat_row(K, new_basis.num, i, new_basis.den) for i = 1:nrows(new_basis)]
  new_profile = nbits(prod(Hecke.upper_bound(fmpz, t2(x)) for x in els))
  if new_profile <= profile_sub
    @vprint :LLL 3 "Output a better basis!\n"
    @vprint :LLL 3 "New profile: $(new_profile)\n"
    return prec, g
  else
    @vprint :LLL 3 "Output the same basis :(\n"
    return prec, identity_matrix(FlintZZ, l)
  end
end

function _lll_with_parameters(M::NfAbsOrd, parameters::Tuple{Float64, Float64}, prec; steps::Int = -1)

  K = nf(M)
  n = degree(M)
  prec = max(prec, 10*n)
  disc = abs(discriminant(M))
  local g::fmpz_mat
  local d::fmpz_mat
  ctx = Nemo.lll_ctx(parameters[1], parameters[2], :gram)
  dM = sum(nbits(Hecke.upper_bound(fmpz, t2(x))) for x in basis(M, K))
  @vprint :LLL 1 "Input profile: $(dM)\n"
  @vprint :LLL 1 "Target profile: $(nbits(disc^2)+divexact(n*(n-1), 2)) \n"
  att = 0
  while steps == -1 || att < steps
    att += 1
    if att > 3
      @vprint :LLL "Having a hard time computing a LLL basis"
    end
    @vprint :LLL 3 "Attempt number : $(att)\n"
    while true
      try
        d = minkowski_gram_mat_scaled(M, prec)
        break
      catch e
        prec = prec*2
      end
    end

    @vprint :LLL 3 "Minkowski matrix computed\n"
    diag_d = prod_diagonal(d)
    g = identity_matrix(FlintZZ, n)

    prec = div(prec, 2)
    shift!(d, -prec)  #TODO: remove?

    for i=1:n
      fmpz_mat_entry_add_ui!(d, i, i, UInt(nrows(d)))
    end
    @vtime :LLL 3 ccall((:fmpz_lll, libflint), Nothing, (Ref{fmpz_mat}, Ref{fmpz_mat}, Ref{Nemo.lll_ctx}), d, g, ctx)

    if nbits(maximum(abs, g)) <= div(prec, 2)
      fl = true
      disc = abs(discriminant(M))
      di = root(disc, n)+1
      di *= fmpz(2)^(div(n+1,2)+prec)
      if compare_index(d, 1, 1, di) > 0
        fl = false
      else
        pr = prod_diagonal(d)
        if pr > fmpz(2)^(div(n*(n-1), 2)+(n*prec)) * disc
          fl = false
        end
      end
      if fl
        break
      end
    end
    On = NfAbsOrd(K, g*basis_matrix(M, copy = false))
    On.ismaximal = M.ismaximal
    if isdefined(M, :index)
      On.index = M.index
    end
    if isdefined(M, :disc)
      On.disc = M.disc
    end
    if isdefined(M, :gen_index)
      On.gen_index = M.gen_index
    end
    prec *= 2
    dOn = nbits(prod(Hecke.upper_bound(fmpz, t2(x)) for x in basis(On, K)))
    if dOn < dM
      @vprint :LLL 3 "I use the transformation\n"
      @vprint :LLL 3 "New profile: $(dOn)\n"
      M = On
      dM = dOn
      prec = Int(floor(prec*1.5))
    else
      prec *= 2
    end
    if att == steps
      return M, prec
    end
    @vprint :LLL 3 "Still in the loop\n"
  end
  On = NfAbsOrd(K, g*basis_matrix(M, copy = false))
  On.ismaximal = M.ismaximal
  if isdefined(M, :index)
    On.index = M.index
  end
  if isdefined(M, :disc)
    On.disc = M.disc
  end
  if isdefined(M, :gen_index)
    On.gen_index = M.gen_index
  end
  @vprint :LLL 3 "LLL finished with parameters $(parameters)\n"
  return On, prec
end

################################################################################
#
#  Applications
#
################################################################################

@doc Markdown.doc"""
    lll_basis(M::NumFieldOrd) -> Vector{NumFieldElem}

A basis for $M$ that is reduced using the LLL algorithm for the Minkowski metric.
"""
function lll_basis(M::NfAbsOrd)
  M1 = lll(M)
  return basis(M1, nf(M1))
end

@doc Markdown.doc"""
    lll_basis(I::NumFieldOrdIdl) -> Vector{NumFieldElem}

A basis for $I$ that is reduced using the LLL algorithm for the Minkowski metric.
"""
function lll_basis(A::NfOrdIdl, v::fmpz_mat = zero_matrix(FlintZZ, 1, 1); prec::Int = 100)
  L, T = lll(A, v, prec=prec)
  S = FakeFmpqMat(T)*basis_matrix(A, copy = false)*basis_matrix(order(A), copy = false)
  K = nf(order(A))
  nS = numerator(S)
  dS = denominator(S)
  q = nf_elem[elem_from_mat_row(K, nS, i, dS) for i=1:degree(K)]
  return q
end

function lll_basis(A::NfOrdFracIdl, v::fmpz_mat = zero_matrix(FlintZZ, 1, 1); prec::Int = 100)
  assure_has_numerator_and_denominator(A)
  L, T = lll(A.num, v, prec=prec)
  S = FakeFmpqMat(T)*basis_matrix(A.num)*basis_matrix(order(A))
  K = nf(order(A))
  nS = numerator(S)
  dS = denominator(S)
  q = nf_elem[elem_from_mat_row(K, nS, i, dS*A.den) for i=1:degree(K)]
  return q
end

################################################################################
#
#  Short element in ideals
#
################################################################################

function short_elem(A::NfOrdFracIdl,
                v::fmpz_mat=zero_matrix(FlintZZ, 1,1); prec::Int = 100)
  assure_has_numerator_and_denominator(A)
  return divexact(short_elem(A.num, v, prec = prec), A.den)
end

function short_elem(A::NfOrdIdl,
                v::fmpz_mat = zero_matrix(FlintZZ, 1,1); prec::Int = 100)
  K = nf(order(A))
  t = lll(A, v, prec = prec)[2]
  w = view(t, 1:1, 1:ncols(t))
  mul!(w, w, basis_matrix(A, copy = false))
  c = w*basis_matrix(order(A), copy = false)
  q = elem_from_mat_row(K, c.num, 1, c.den)
  return q
end

################################################################################
#
#  Reduction of ideals
#
################################################################################

function reduce_ideal(A::NfOrdIdl)
  B = inv(A)
  b = short_elem(B)
  C = b*A
  simplify(C)
  @assert C.den == 1
  return C.num, b
end

function reduce_product(A::NfOrdIdl, B::NfOrdIdl)
  I = inv(A)
  J = inv(B)
  @vtime :LLL 3 bIJ = _lll_product_basis(I.num, J.num)
  pp = NfOrdIdl(order(A), bIJ)
  @vtime :LLL 3 b = divexact(short_elem(pp), I.den * J.den)
  AB = A*B
  C = b*AB
  simplify(C)
  @assert isone(C.den)
  return C.num, b
end

function reduce_ideal(A::NfOrdFracIdl)
  B = inv(A)
  b = short_elem(B.num)
  C = divexact(b, B.den)*A
  simplify(C)
  @assert C.den == 1
  return C.num, divexact(b, B.den)
end

@doc Markdown.doc"""
    reduce_ideal(A::FacElem{NfOrdIdl}) -> NfOrdIdl, FacElem{nf_elem}

Computes $B$ and $\alpha$ in factored form, such that $\alpha B = A$.
"""
function reduce_ideal(I::FacElem{NfOrdIdl, NfOrdIdlSet})
  @assert !isempty(I.fac)
  O = order(first(keys(I.fac)))
  K = nf(O)
  fst = true
  a = FacElem(Dict{nf_elem, fmpz}(one(K) => fmpz(1)))
  A = ideal(O, 1)
  for (k, v) = I.fac
    @assert order(k) === O
    if iszero(v)
      continue
    end
    if fst
      A, a = power_reduce(k, v)
      @hassert :PID_Test 1 (v>0 ? fractional_ideal(O, k)^Int(v) : inv(k)^Int(-v)) == A*evaluate(a)
      fst = false
    else
      B, b = power_reduce(k, v)
      mul!(a, a, b)
      @hassert :PID_Test (v>0 ? fractional_ideal(O, k)^Int(v) : inv(k)^Int(-v)) == B*evaluate(b)
      if norm(A, copy = false)*norm(B, copy = false) > abs(discriminant(O))
        A, c = reduce_product(A, B)
        add_to_key!(a.fac, c, -1)
      else
        A = A*B
      end
    end
  end
  @hassert :PID_Test 1 A*evaluate(a) == evaluate(I)
  return A, a
end

# The bound should be sqrt(disc) (something from LLL)
@doc Markdown.doc"""
    power_reduce(A::NfOrdIdl, e::fmpz) -> NfOrdIdl, FacElem{nf_elem}

Computes $B$ and $\alpha$ in factored form, such that $\alpha B = A^e$
$B$ has small norm.
"""
function power_reduce(A::NfOrdIdl, e::fmpz)
  O = order(A)
  K= nf(O)
  if norm(A, copy = false) > abs(discriminant(O))
    A1, a = reduce_ideal(A)
    @hassert :PID_Test 1 a*A == A1
    A = A1
    al = FacElem(Dict(a=>-e))
  else
    al = FacElem(Dict(K(1) => fmpz(1)))
  end

  #we have A_orig^e = (A*a)^e = A^e*a^e = A^e*al and A is now small

  if e < 0
    B = inv(A)
    A = numerator(B)
    add_to_key!(al.fac, K(denominator(B)), fmpz(e))
    e = -e
  end

  if isone(e)
    return A, al
  end
  # A^e = A^(e/2)^2 A or A^(e/2)^2
  # al * A^old^(e/2) = A_new
  C, cl = power_reduce(A, div(e, 2))
  @hassert :PID_Test 1 C*evaluate(cl) == A^Int(div(e, 2))
  mul!(al, al, cl^2)
  if norm(C, copy = false)^2 > abs(discriminant(O))
    @vtime :CompactPresentation :4 C2, a = reduce_product(C, C)
    add_to_key!(al.fac, a, -1)
  else
    C2 = C^2
  end

  if isodd(e)
    if norm(A, copy = false)*norm(C2, copy = false) > abs(discriminant(O))
      @vtime :CompactPresentation :4 A1, a = reduce_product(C2, A)
      A = A1
      add_to_key!(al.fac, a, -1)
    else
      A = C2*A
    end
  else
    A = C2
  end
  return A, al
end

function new_power_reduce(A::NfOrdIdl, e::fmpz)
  O = order(A)
  if iszero(e)
    return ideal(O, 1)
  end
  K = nf(O)
  al = FacElem(Dict(K(1) => fmpz(1)))
  if e < 0
    A1 = inv(A)
    A = A1.num
    add_to_key!(al.fac, K(A1.den), e)
    e = -e
  end

  if norm(A)^e <= abs(discriminant(O))
    return A^e, al
  end
  Ainvtot = inv(A)
  Ainv = Ainvtot.num
  d = Ainvtot.den
  res = _new_power_reduce(A, e, Ainv, d)
  mul!(al, al, res[2])
  return res[1], al
end

function _new_power_reduce(A::NfOrdIdl, e::fmpz, Ainv::NfOrdIdl, d::fmpz)
  #Ainv//d is the inverse of A
  #We want to reduce A^e
  O = order(A)
  K = nf(O)
  bdisc = abs(discriminant(O))
  if norm(A) > bdisc
    a = divexact(basis(Ainv)[1].elem_in_nf, d)
    A1 = a*A
    simplify(A1)
    A = A1.num
    al = FacElem(Dict(a=>-e))
    #I need to update Ainv!
    Ainvtot = inv(A)
    Ainv = Ainvtot.num
    d = Ainvtot.den
    newb = lll(Ainv)[2]
    mul!(newb, newb, basis_matrix(Ainv, copy = false))
    Ainv.basis_matrix = newb
  else
    al = FacElem(Dict(K(1) => fmpz(1)))
  end

  if isone(e)
    return A, al, Ainv, d
  end

  # A^e = A^(e/2)^2 A or A^(e/2)^2
  # al * A^old^(e/2) = A_new
  C, cl, Cinv, dCinv = _new_power_reduce(A, div(e, 2), Ainv, d)
  @hassert :PID_Test 1 C*evaluate(cl) == A^Int(div(e, 2))
  mul!(al, al, cl^2)
  if norm(C)^2 > bdisc
    a = divexact(short_elem_product(Cinv, Cinv), dCinv^2)
    C21 = a*(C^2)
    simplify(C21)
    C2 = C21.num
    C2invtot = inv(C2)
    C2inv = C2invtot.num
    C2d = C2invtot.den
    newb = lll(C2inv)[2]
    mul!(newb, newb, basis_matrix(C2inv, copy = false))
    C2inv.basis_matrix = newb
    add_to_key!(al.fac, a, -1)
  else
    C2 = C^2
    basis_IJ = _lll_product_basis(Cinv, Cinv)
    IJ = NfOrdIdl(O, basis_IJ)
    newb = lll(IJ)[2]
    mul!(newb, newb, basis_matrix(IJ, copy = false))
    IJ.basis_matrix = newb
    C2inv = IJ
    C2d = dCinv*dCinv
  end

  if isodd(e)
    if norm(A)*norm(C2) > bdisc
      a = divexact(short_elem_product(C2inv, Ainv), C2d*d)
      A1 = a*C2*A
      simplify(A1)
      A = A1.num
      add_to_key!(al.fac, a, -1)
      Ainvtot = inv(A)
      d = Ainvtot.den
      Ainv = Ainvtot.num
      bnew = lll(Ainv)[2]
      mul!(bnew, bnew, basis_matrix(Ainv, copy = false))
      Ainv.basis_matrix = bnew
    else
      A = C2*A
      basis_IJ = _lll_product_basis(C2inv, Ainv)
      IJ = NfOrdIdl(O, basis_IJ)
      newb = lll(IJ)[2]
      mul!(newb, newb, basis_matrix(IJ, copy = false))
      IJ.basis_matrix = newb
      Ainv = IJ
      d = d*C2d
    end
  else
    A = C2
    d = C2d
    Ainv = C2inv
  end
  return A, al, Ainv, d
end

function short_elem_product(A::NfOrdIdl, B::NfOrdIdl)
  return lll_basis_product(A, B)[1]
end

################################################################################
#
#  Short basis of product of ideals
#
################################################################################

# Claus and Tommy:
# We express the basis of IJ in terms of the basis of I (I contains IJ)
# Then we compute the lll of the matrix of the coordinates. This way we get a
# better basis to start the computation of LLL
#We compute the hnf to have a guaranteed bound on the entries
function _lll_product_basis(I::NfOrdIdl, J::NfOrdIdl)
  A = lll(I)[2]
  mul!(A, A, basis_matrix(I, copy = false))
  IJ = I*J
  C = basis_matrix(IJ, copy = false)
  @vtime :LLL 3 iA, de = pseudo_inv(A)
  mul!(iA, C, iA)
  divexact!(iA, iA, de)
  hnf_modular_eldiv!(iA, minimum(J))
  @vtime :LLL 3 lll!(iA, lll_ctx(0.3, 0.51))
  @vtime :LLL 3 lll!(iA)
  mul!(iA, iA, A)
  return iA
end

function lll_basis_product(I::NfOrdIdl, J::NfOrdIdl)

  basis_IJ = _lll_product_basis(I, J)
  IJ = NfOrdIdl(order(I), basis_IJ)
  res = lll_basis(IJ)
  return res
end


text:
#=ДАНО: Робот - в произвольной клетке ограниченного прямоугольной рамкой поля без внутренних перегородок и маркеров.

РЕЗУЛЬТАТ: Робот - в исходном положении в центре косого креста (в форме X) из маркеров.

(Подсказка: решение будет подобно решению задачи 1, если направление премещения Робота задавать кортежами пары значений типа HorizonSide)=#

function x_mark(r::Robot)
    for side in ((Nord,Ost),(Nord,West),(Sud,West),(Sud,Ost))
        putmarkers!(r,side)
        move_by_markers!(r,inverse(side))
    end
    putmarker!(r)
end

function putmarkers!(r::Robot,side::NTuple{2,HorizonSide})
    while isborder(r,side)==false
        move!(r,side)
        putmarker!(r)
    end
end
#то же самое что и первое задание;отличие в измененной входной переменной(side->NTuple{2,HorizonSide})
HorizonSideRobots.isborder(r::Robot,side::NTuple{2,HorizonSide}) = (isborder(r,side[1]) || isborder(r,side[2]))

HorizonSideRobots.move!(r::Robot, side::NTuple{2,HorizonSide}) = for s in side move!(r,s) end
move_by_markers!(r::Robot,side::NTuple{2,HorizonSide}) = while ismarker(r) move!(r,side) end

inverse(side::NTuple{2,HorizonSide}) = (inverse(side[1]),inverse(side[2]))
inverse(side::HorizonSide)=HorizonSide(mod(Int(side)+2, 4))
#=NTuple{2,HorizonSide}-кортеж,где первый элемент-кол-во элементов,а второй тип элемента-стороны=#

max_stars_repo_name: ilya-design/mirea-progs

max_issues_repo_name: ilya-design/mirea-progs

max_forks_repo_name: ilya-design/mirea-progs


text:
function get_a(obj::FieldDSSTThirdBodyContext)
    return jcall(obj, "getA", RealFieldElement, ())
end

function get_alpha(obj::FieldDSSTThirdBodyContext)
    return jcall(obj, "getAlpha", RealFieldElement, ())
end

function get_ao_r3_pow(obj::FieldDSSTThirdBodyContext)
    return jcall(obj, "getAoR3Pow", Vector{RealFieldElement}, ())
end

function get_bb(obj::FieldDSSTThirdBodyContext)
    return jcall(obj, "getBB", RealFieldElement, ())
end

function get_bbb(obj::FieldDSSTThirdBodyContext)
    return jcall(obj, "getBBB", RealFieldElement, ())
end

function get_beta(obj::FieldDSSTThirdBodyContext)
    return jcall(obj, "getBeta", RealFieldElement, ())
end

function get_bo_a(obj::FieldDSSTThirdBodyContext)
    return jcall(obj, "getBoA", RealFieldElement, ())
end

function get_bo_a_bpo(obj::FieldDSSTThirdBodyContext)
    return jcall(obj, "getBoABpo", RealFieldElement, ())
end

function get_gamma(obj::FieldDSSTThirdBodyContext)
    return jcall(obj, "getGamma", RealFieldElement, ())
end

function get_hxxx(obj::FieldDSSTThirdBodyContext)
    return jcall(obj, "getHXXX", RealFieldElement, ())
end

function get_kxxx(obj::FieldDSSTThirdBodyContext)
    return jcall(obj, "getKXXX", RealFieldElement, ())
end

function get_m2ao_a(obj::FieldDSSTThirdBodyContext)
    return jcall(obj, "getM2aoA", RealFieldElement, ())
end

function get_m_co2ab(obj::FieldDSSTThirdBodyContext)
    return jcall(obj, "getMCo2AB", RealFieldElement, ())
end

function get_max_ar3_pow(obj::FieldDSSTThirdBodyContext)
    return jcall(obj, "getMaxAR3Pow", jint, ())
end

function get_max_ecc_pow(obj::FieldDSSTThirdBodyContext)
    return jcall(obj, "getMaxEccPow", jint, ())
end

function get_max_freq_f(obj::FieldDSSTThirdBodyContext)
    return jcall(obj, "getMaxFreqF", jint, ())
end

function get_mean_motion(obj::FieldDSSTThirdBodyContext)
    return jcall(obj, "getMeanMotion", RealFieldElement, ())
end

function get_muo_r3(obj::FieldDSSTThirdBodyContext)
    return jcall(obj, "getMuoR3", RealFieldElement, ())
end

function get_oo_ab(obj::FieldDSSTThirdBodyContext)
    return jcall(obj, "getOoAB", RealFieldElement, ())
end

function get_qns(obj::FieldDSSTThirdBodyContext)
    return jcall(obj, "getQns", Vector{Vector{RealFieldElement}}, ())
end

function get_x(obj::FieldDSSTThirdBodyContext)
    return jcall(obj, "getX", RealFieldElement, ())
end

function getb(obj::FieldDSSTThirdBodyContext)
    return jcall(obj, "getb", RealFieldElement, ())
end


text:
# ---
# title: 472. Concatenated Words
# id: problem472
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Dynamic Programming, Depth-first Search, Trie
# link: <https://leetcode.com/problems/concatenated-words/description/>
# hidden: true
# ---
#
# Given a list of words ( **without duplicates** ), please write a program that
# returns all concatenated words in the given list of words.
#
# A concatenated word is defined as a string that is comprised entirely of at
# least two shorter words in the given array.
#
# **Example:**
#
#
#
#     Input: ["cat","cats","catsdogcats","dog","dogcatsdog","hippopotamuses","rat","ratcatdogcat"]
#
#     Output: ["catsdogcats","dogcatsdog","ratcatdogcat"]
#
#     Explanation: "catsdogcats" can be concatenated by "cats", "dog" and "cats";
#      "dogcatsdog" can be concatenated by "dog", "cats" and "dog";
#     "ratcatdogcat" can be concatenated by "rat", "cat", "dog" and "cat".
#
#
# **Note:**
#
#   1. The number of elements of the given array will not exceed `10,000 `
#   2. The length sum of elements in the given array will not exceed `600,000`.
#   3. All the input string will only include lower case letters.
#   4. The returned elements order does not matter.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
module CartesianGeneticProgramming

import YAML
using Cambrian
import JSON

include("functions.jl")
include("config.jl")
include("individual.jl")
include("process.jl")
include("mutation.jl")
include("crossover.jl")
include("evolution.jl")

end


text:
function second_moment(q::Normal)
    return abs2(mean(q)) + var(q)
end

function second_moment(q::Normal, y::Real)
    return abs2(mean(q) - y) + var(q)
end

# This is not exactly an approximation but corresponds anyway
# to the expectation of the function σ(f)
function approx_expected_logistic(μ, c)
    lower, upper = LogExpFunctions._logistic_bounds(μ)
    return μ < lower ? zero(μ) : (μ > upper ? one(μ) : exp(μ / 2) * sech(c / 2) / 2)
end


text:
module VTKDataIO

__precompile__(false)

using VTKDataTypes
using WriteVTK
using PyCall
using LightXML
using IterTools
using Printf

if get(ENV, "JULIA_REGISTRYCI_AUTOMERGE", "false") == "false"
    @eval begin
        const vtkns = pyimport_conda("vtk.util.numpy_support", "vtk")
        const vtk = pyimport_conda("vtk", "vtk")
        const np = pyimport_conda("numpy", "numpy")

        include("vtkreaders.jl")
        include("vtkwriters.jl")
        include("vtpwriter.jl")
        include("stlIO.jl")
        include("extra_readers.jl")
        include("PyVTK.jl")
        include("visualize.jl")

        export read_vtk, write_vtk, read_stl, write_stl, valid_to_write, PyVTK,
            _VTKDataTypes, visualize, write_x3d, write_ply, read_obj, read_ply, visualize_3ds
    end
end

end


text:
module Utils

using CUDA

function include_tests(m, dir = @__DIR__)
    for file in readdir(dir)
        if match(r"^test_.*\.jl$", file) !== nothing
            Base.include(m, joinpath(dir, file))
        end
    end
end

should_test_gpu() =
    lowercase(get(ENV, "JULIA_PKGEVAL", "false")) != "true" &&
    lowercase(get(ENV, "CUDAFOLDS_JL_TEST_GPU", "true")) == "true"

"""
    unsafe_free_all!(namespace::Module)

Call `CUDA.unsafe_free!` on all `CuArray`s in `namespace`.
"""
function unsafe_free_all!(namespace::Module)
    for n in names(namespace; all = true)
        x = try
            getfield(namespace, n)
        catch err
            err isa UndefVarError || rethrow()
            continue
        end
        if x isa CuArray
            @debug "Free: `$namespace.$n`"
            CUDA.unsafe_free!(x)
        end
    end
end

end  # module