text:
# using QuadGK
# using Plots
# using GSL
using Cuba

push!(LOAD_PATH, ENV["HOME"]*"/Documents/amplitude_analysis/modules")

using isobars
using masses
using DalitzPlotAnalysis

function phsp(s)
    function integrand(x, f)
        # complex path
        #  end points:
        s1_i = 4mπ2;
        s1_f = (sqrt(s)-mπ)^2 - 1e-8 # otherwise, rounding error
        s1_m = real(s1_f)*one(s);
        #  variables:
        s1 = (x[1] < 0.5) ? s1_i + 2x[1]*(s1_m-s1_i) : s1_m + 2(x[1]-0.5)*(s1_f-s1_m)
        #  jacobian
        res = one(1.0im)
        res *= (x[1] < 0.5) ? 2(s1_m-s1_i) : 2(s1_f-s1_m);
        # angular variables
        τ1 = (2*x[2]-1,0.0,2*x[3]-1,π*(2*x[4]-1))
        τ3 = fill(zero(s),4);
        s3,τ3[1],τ3[2],τ3[3],τ3[4] = change_basis(s1,τ1[1],τ1[2],τ1[3],τ1[4],mπ2,mπ2,mπ2,s);
        [ if !isfinite(τi)
           print("s = ",s,":\n",s1,", ",τ1,"\n ",s3,", ",τ3,"\n\n")
           end for τi in τ3 ]
        # the angular part of the rho-pi S-wave
        # For the rho pi S-wave, it reads:
        #    cosθ1*cosθ23-sinθ1*sinθ23*cosϕ23
        ang = [
            begin
                cI = τ[1]; sI = sqrt(1.0-cI^2);
                cV = τ[3]; sV = sqrt(1.0-cV^2);
                sqrt(3.)*(cI*cV-sI*sV*cos(τ[4]))
            end for τ in [τ1,τ3] ]
        # value
        amp  = f1_I(s1)*ang[1]-f1_I(s3)*ang[2]
        ampc = f1_II(s1)*ang[1]-f1_II(s3)*ang[2]
        res *= amp*ampc/2
        # intagra jacobian factors
        res *= sqrt(λ(s,s1,mπ2)*λ(s1,mπ2,mπ2))/(s*s1)
        f[1],f[2] = reim(res)
    end
    result = cuhre(integrand, 4, 2)
    complex(result[1]...) / (2*π) / (8*π)^2
end

function phsp_sim(s)
    function integrand(x, f)
        # complex path
        #  end points:
        s1_i = 4mπ2;
        s1_f = (sqrt(s)-mπ)^2 - 1e-8 # otherwise, rounding error
        s1_m = real(s1_f)*one(s);
        #  variables:
        s1 = (x[1] < 0.5) ? s1_i + 2x[1]*(s1_m-s1_i) : s1_m + 2(x[1]-0.5)*(s1_f-s1_m)
        #  jacobian
        res = one(1.0im)
        res *= (x[1] < 0.5) ? 2(s1_m-s1_i) : 2(s1_f-s1_m);
        # value
        amp  = f1_I(s1)
        ampc = f1_II(s1)
        res *= amp*ampc
        # intagra jacobian factors
        res *= sqrt(λ(s,s1,mπ2)*λ(s1,mπ2,mπ2))/(s*s1)
        f[1],f[2] = reim(res)
    end
    result = cuhre(integrand, 1, 2)
    complex(result[1]...) / (2*π) / (8*π)^2
end

# phsp_sim(1.2+0.1im)

let fout_name="/tmp/phsp_sym_50.txt", func=phsp, nVx = 50, nVy = 50
    sxv = linspace(0.8, 2.0, nVx)
    syv = linspace(-1.0, 0.1, nVy)[end:-1:1]
    close(open(fout_name,"w")) # to clean the file
    for (j, sy) in enumerate(syv), (i, sx) in enumerate(sxv)
        v = func(sx+1im*sy)
        file = open(fout_name,"a")
        writedlm(file, [sx sy real(v) imag(v)])
        close(file)
        perc = 100(nVx*(j-1)+i)/(nVx*nVy)
        println("\%$perc done: ($i/$nVx, $j/$nVy), y = $sy")
    end
end

# data = readdlm("/tmp/phsp.txt")
#
# matr = hcat([data[(50(i-1)+1):(50i),3]+1im*data[(50(i-1)+1):(50i),4] for i in 1:50]...)
#
# contour(real(matr))
# contour!(imag(matr))
#
# heatmap(real(matr))
# heatmap(imag(matr))
#


text:
struct BudgetedElementaryPathInstance{T} <: CombinatorialInstance
  graph::AbstractGraph{T}
  rewards::Dict{Edge{T}, Float64}
  weights::Dict{Edge{T}, Int}
  src::T
  dst::T

  budget::Int # weights * solution >= budget

  function BudgetedElementaryPathInstance(graph::AbstractGraph{T}, rewards::Dict{Edge{T}, Float64},
                                          weights::Dict{Edge{T}, Int}, src::T, dst::T;
                                          budget::Union{Nothing, Int}=nothing,
                                          max_weight::Union{Nothing, Int}=nothing) where T
    # Error checking.
    if src == dst
      error("Source node is the same as destination node.")
    end

    if ! isnothing(budget) && budget < 0
      error("Budget is present and negative.")
    end

    if any(collect(values(weights)) .< 0)
      error("At least a weight is negative.")
    end

    if length(rewards) != length(weights)
      error("Not the same number of values and weights; these two vectors must have the same size.")
    end

    # Complete the optional parameters.
    if isnothing(max_weight) # Not useful to store, it is just used to compute the max budget if it is not given.
      max_weight = maximum(collect(values(weights)))
    end

    if isnothing(budget)
      d = length(values) # Dimension of the problem.
      budget = d * max_weight
    end

    # Return a new instance.
    new{T}(graph, rewards, weights, src, dst, budget)
  end
end

graph(i::BudgetedElementaryPathInstance{T}) where T = i.graph
rewards(i::BudgetedElementaryPathInstance{T}) where T = i.rewards
weights(i::BudgetedElementaryPathInstance{T}) where T = i.weights
src(i::BudgetedElementaryPathInstance{T}) where T = i.src
dst(i::BudgetedElementaryPathInstance{T}) where T = i.dst
budget(i::BudgetedElementaryPathInstance{T}) where T = i.budget

# dimension(i::BudgetedElementaryPathInstance{T}) where T = ne(graph(i))
reward(i::BudgetedElementaryPathInstance{T}, u::T, v::T) where T = rewards(i)[Edge(u, v)]
weight(i::BudgetedElementaryPathInstance{T}, u::T, v::T) where T = weights(i)[Edge(u, v)]

struct BudgetedElementaryPathSolution{T} <: CombinatorialSolution
  instance::BudgetedElementaryPathInstance{T}
  path::Vector{Edge{T}}
  states::Dict{Tuple{T, Int}, Float64}
  solutions::Dict{Tuple{T, Int}, Vector{Edge{T}}}
end

function paths_all_budgets(s::BudgetedElementaryPathSolution{T}, max_budget::Int) where T
  if max_budget > budget(s.instance)
    @warn "The asked maximum budget $max_budget is higher than the instance budget $(budget(s.instance)). Therefore, some values have not been computed and are unavailable."
  end

  mb = min(max_budget, budget(s.instance))
  return Dict{Int, Vector{Edge{T}}}(
    budget => s.solutions[s.instance.dst, budget] for budget in 0:mb)
end

function paths_all_budgets_as_tuples(s::BudgetedElementaryPathSolution{T}, max_budget::Int) where T
  if max_budget > budget(s.instance)
    @warn "The asked maximum budget $max_budget is higher than the instance budget $(budget(s.instance)). Therefore, some values have not been computed and are unavailable."
  end

  mb = min(max_budget, budget(s.instance))
  return Dict{Int, Vector{Tuple{T, T}}}(
    budget => [(src(e), dst(e)) for e in s.solutions[s.instance.dst, budget]]
    for budget in 0:mb)
end

ext: jl


text:
sp = " "
num1 = split("one two three four five six seven eight nine", sp)
num2 = split("eleven twelve thirteen fourteen fifteen sixteen seventeen eighteen nineteen", sp)
num3 = split("ten twenty thirty forty fifty sixty seventy eighty ninety", sp)
num4 = split("hundred thousand", sp)
conj = "and"

function n2w(i)
  words = []
  if i == 1000
    push!(words, num1[1], num4[2])
  else
    if i >= 100
      push!(words, num1[Int(floor(i/100))], num4[1])
      if i % 100 != 0
        push!(words, conj)
      end
    end
    if i % 100 != 0
      if 11 <= i%100 <= 19
        push!(words, num2[Int(floor(i%100%10))])
      else
        if i % 100 >= 10
          push!(words, num3[Int(floor(i%100/10))])
        end
        if i % 10 != 0
          push!(words, num1[Int(floor(i%10))])
        end
      end
    end
  end
  return words
end

total = 0
for i = 1:1000
  global total += sum(map(x -> length(x), n2w(i)))
end
println(total)

max_stars_repo_name: sujimodern/project-euler-jl

max_issues_repo_name: sujimodern/project-euler-jl

max_forks_repo_name: sujimodern/project-euler-jl


text:
const SEARCH_CACHE = Dict{String,Any}()
"""
    synonyms(query)

Given a chemical search query, return the synonyms associated to that query.
This function doesn't have any cache.

# Examples:
```repl
synonyms("water")
```
"""
function synonyms(query)
    compound_id,key = search_chemical_id(detect_query(query))
    return __synonyms(compound_id,key)
end

function __synonyms(idx,key)
    db,sdb = DATA_DB[key]
    return sdb.list[findall(isequal(idx),sdb.index)]
end

function build_result(idx,key)
    if idx == -1
        return missing
    end
    db,sdb = DATA_DB[key]
    pubchemid= db.pubchemid[idx]
    CAS= db.CAS[idx]
    formula= db.formula[idx]
    MW= db.MW[idx]
    smiles= db.smiles[idx]
    InChI= db.InChI[idx]
    InChI_key= db.InChI_key[idx]
    iupac_name= db.iupac_name[idx]
    common_name= db.common_name[idx]
    return (;pubchemid, CAS, formula, MW, smiles, InChI, InChI_key, iupac_name, common_name)
end

function detect_query(id::String)
    if is_element(id)
        return ElementQuery(id)
    elseif is_cas(id)
        return CASQuery(id)
    elseif is_pubchemid(id)
        return PubChemIDQuery(id)
    elseif is_inchi(id)
        return InChIQuery(id)
    elseif is_inchikey(id)
        return InChIKeyQuery(id)
    elseif is_smiles(id)
        return SMILESQuery(id)
    else
        return AnyQuery(id)
    end
end

function detect_query(id::Int)
        return PubChemIDQuery(id)
end

function detect_query(id::Tuple{Integer,Integer,Integer})
    return CASQuery(id)
end

function db_iteration_order(DB)
    #generating keys iteration order: user, short database, long database
    dbnames = [:short,:long]
    for k in keys(DATA_DB)
        if !(k in (:short,:long))
            push!(dbnames,k)
        end
    end
    return dbnames
end

function search_chemical(query,cache=SEARCH_CACHE)
    if cache !== nothing
        #the base is that two chemicals with different casing should be the same.
        if query isa NTuple{3,Integer}
            ID = tuple_to_casstr(query)
        else
            ID = string(query)
        end

        if haskey(cache,ID)
            return cache[ID]
        end

        normalized_id = Unicode.normalize(ID,casefold = true,stripmark=true)
        if haskey(cache,normalized_id)
            return cache[normalized_id]
        end

        compound_id,key = search_chemical_id(detect_query(query))
        res = build_result(compound_id,key)      #return db.common_name[compound_id]
        cache[ID] = res
        cache[normalized_id] = res
    else
        compound_id,key = search_chemical_id(detect_query(query))
        return build_result(compound_id,key)      #return db.common_name[compound_id]
    end
end

function search_chemical_id(ID::AnyQuery;skip_common_name = false,try_strategies = true)::Tuple{Int,Symbol}
    #skip_common_name skips search on the common_name col
    compound_id = -1
    id = value(ID)
    search_done = false
    _keys = db_iteration_order(DATA_DB)
    for key in _keys
    db,sdb,sortdb = DATA_DB[key]
        if !skip_common_name
            searchvec = view(db.common_name,sortdb.common_name_sort)
            idx_sort = searchsorted(searchvec,id)
            #@show idx_sort
            if length(idx_sort) == 1 #found an element
                idx = only(sortdb.common_name_sort[idx_sort])
                compound_id =idx
                search_done = true
            end
        end
        if !search_done
            idx = searchsorted(sdb.list,id)
            if length(idx) == 1 #found an element
                compound_id = only(sdb.index[idx])
                search_done = true
            end
        end
        #found in db,returning
        if search_done
            return compound_id,key
        end
    end

    if !search_done
        if !try_strategies #bail out here if requested
            return -1,:not_found
        end
    end
    #==
    result not found, trying same strategies as present in CalebBell/Chemicals
    #strategy 1: trying without spaces and dashs.
    ==#
    _ids = Vector{String}(undef,7)
    _ids[1] = Unicode.normalize(id,casefold = true,stripmark=true)
    _ids[2] = replace(_ids[1]," "=>"")
    _ids[3] = replace(_ids[2],"-"=>"")
    _ids[4] = replace(id," "=>"")
    _ids[5] = replace(_ids[4],"-"=>"")

    #those matches find chemicals of the form n-name
    #or 1-name
    _ids[6] = begin
        if occursin(r"1-[A-Za-z]+",_ids[1]) |  occursin(r"n-[A-Za-z]+",_ids[1])
            chop(id,head=2,tail=0)
        else
            _ids[6] = id
        end
    end
    _ids[7] = Unicode.normalize(_ids[6],casefold = true,stripmark=true)
    _ids = unique!(_ids)
    _ids = setdiff!(_ids,[id])

    for _id in _ids
        compound_id,key = search_chemical_id(AnyQuery(_id))
        if compound_id !== -1
            search_done = true
            break
        end
    end

    #strategy 2: trying to match in the form 'water (H2O)'
    if !search_done
        re = r"\w+\s+\([\s\w]+\)"
        if occursin(re,id)
            _id = id |> z->replace(z,")"=>"") |> z->split(z,"(") .|> strip
            id1,id2 = first(id1),last(id2)
            compound_id1,key1 = search_chemical_id(AnyQuery(id1))
            compound_id2,key2 = search_chemical_id(AnyQuery(id2))
            if (compound_id1 == compound_id2) & (key1==key2)
                search_done = true
                compound_id = compound_id2
                key = key2
            end
        end
    end

    #if something worked, return here, else, return not found
    if search_done
       return compound_id,key
    else
        return -1,:not_found
    end
end

function search_chemical_id(ID::CASQuery)::Tuple{Int,Symbol}
    id = cas(ID)
    compound_id,key = search_id_impl(id,:CAS)
    if compound_id != -1
        return compound_id,key
    end
    return search_chemical_id(AnyQuery(value(ID)),skip_common_name=true,try_strategies=false)
end

function search_chemical_id(ID::InChIKeyQuery)::Tuple{Int,Symbol}
    id = value(ID)
    return search_id_impl(id,:InChI_key)
end

function search_chemical_id(ID::PubChemIDQuery)::Tuple{Int,Symbol}
    id = id_num(ID)
    return search_id_impl(id,:pubchemid)
end

function search_chemical_id(ID::InChIQuery)::Tuple{Int,Symbol}
    id = value(ID)
    search_id_impl(id,:InChI)
end

function search_chemical_id(ID::SMILESQuery)::Tuple{Int,Symbol}
    id = value(ID)
    search_id_impl(id,:smiles)
end

arrowtype(::Type{Int}) = Arrow.Primitive{Int,Vector{Int}}
arrowtype(::Type{String}) = Arrow.List{String, Int32, Vector{UInt8}}
arrowtype(::Type{Tuple{Int32,Int16,Int16}}) = Arrow.Struct{Tuple{Int32, Int16, Int16}, Tuple{Arrow.Primitive{Int32, Vector{Int32}}, Arrow.Primitive{Int16, Vector{Int16}}, Arrow.Primitive{Int16, Vector{Int16}}}}

function search_id_impl(id::T,k::Symbol)::Tuple{Int,Symbol} where {T}
    arrowT = arrowtype(T)
    return search_id_impl(id,k,arrowT)
end

function search_id_impl(id::T,sym::Symbol,::Type{A})::Tuple{Int,Symbol} where {T,A}
    compound_id::Int = -1
    search_done = false
    dbnames = db_iteration_order(DATA_DB)::Vector{Symbol}
    for dbname in dbnames
        db,_,sortdb = DATA_DB[dbname]
        dbcol = getproperty(db,sym)::A
        sort_sym = Symbol(sym,:_sort)
        dbidx = getproperty(sortdb,sort_sym)::Arrow.Primitive{Int,Vector{Int}}
        searchvec = view(dbcol,dbidx)
        idxs = searchsorted(searchvec,id)
        #@show idx_sort
        if length(idxs) == 1 #found and element
            compound_id = only(dbidx[idxs])::Int
            search_done = true
        else
            throw("Search is not unique, multiple matches found for $id in database $dbname, on the $sym column")
        end
        if search_done
            return compound_id,dbname
        end
    end
    return compound_id,sym #to identify where it fails
end


text:
""" Supertype for all generation technologies"""
abstract type Generator <: StaticInjection end
const Generators = Array{<:Generator, 1}

""" Supertype for all Hydropower generation technologies"""
abstract type HydroGen <: Generator end

"""
Supertype for all renewable generation technologies
Requires the implementation of `get_rating`and `get_power_factor` methods
"""
abstract type RenewableGen <: Generator end

""" Supertype for all Thermal generation technologies"""
abstract type ThermalGen <: Generator end

function IS.get_limits(
    valid_range::Union{NamedTuple{(:min, :max)}, NamedTuple{(:max, :min)}},
    unused::T,
) where {T <: Generator}
    # Gets min and max value defined for a field,
    # e.g. "valid_range": {"min":-1.571, "max":1.571}.
    return (min = valid_range.min, max = valid_range.max, zero = 0.0)
end

"""
Return the max active power for the Renewable Generation calculated as the rating * power_factor
"""
function get_max_active_power(d::T) where {T <: RenewableGen}
    return get_rating(d) * get_power_factor(d)
end

"""
Return the max reactive power for the Renewable Generation calculated as the rating * power_factor
"""
function get_max_reactive_power(d::T) where {T <: RenewableGen}
    return get_rating(d) * sin(acos(get_power_factor(d)))
end


text:
@testset "ellipse minor axis length" begin

    test_image_minor_axes_lengths = ([0], [5.773503], [1.154701],
                        [1.154701, 1.154701, 1.154701, 1.154701], [5.982882],
                        [5.537749], [3.651484], [3.944053], [1.154701, 1.154701], [5.341660],
                        [5.617433], [5.272805], [6.323095], [6.070635], [5.033223],
                        [4.570858, 3.651484, 5.163978], [10.105787],
                        [5.341433, 4.435154],
                        [11.796250 , 4.953113], [11.012302, 3.909324])

    for T in (Int, Bool, Gray{Bool}, Gray{N0f8}, Gray{N0f16}, Gray{N0f32}, Gray{Float64})
        for i = 1:20
            test_image = eval(Symbol("test_image_$(i)"))
            labels = Images.label_components(test_image,trues(3,3))
            t = measure_components((RegionEllipse(),), labels)
            if i > 1
                semi_axes = select(t,:semi_axes)
                a = first.(collect(semi_axes)) * 2
                @test all(isapprox(a, test_image_minor_axes_lengths[i]; atol = 1e-5))
            end
        end
    end
end

@testset "ellipse major axis length" begin

    test_image_major_axes_lengths = ([0], [5.773503], [1.154701],
                        [1.154701, 1.154701, 1.154701, 1.154701], [5.982882],
                        [5.537749], [3.651484], [4.760952], [4.618802, 4.618802], [5.773503],
                        [5.925463], [5.925463], [6.461763], [6.218253], [6.110101],
                        [9.988690, 3.651484, 5.163978], [10.479679],
                        [8.610212, 5.171661],
                        [12.170053 , 5.773503], [12.078644, 6.426807])

    for T in (Int, Bool, Gray{Bool}, Gray{N0f8}, Gray{N0f16}, Gray{N0f32}, Gray{Float64})
        for i = 1:20
            test_image = eval(Symbol("test_image_$(i)"))
            labels = Images.label_components(test_image,trues(3,3))
            t = measure_components((RegionEllipse(),), labels)
            if i > 1
                semi_axes = select(t,:semi_axes)
                b = last.(collect(semi_axes)) * 2
                @test all(isapprox(b, test_image_major_axes_lengths[i]; atol = 1e-5))
            end
        end
    end
end

@testset "ellipse centroid" begin

    test_image_centroids = ([(0,0)], [(3, 3)],
                        [(3, 3)],
                        [(2, 2), (4, 2), (2, 4), (4, 4)],
                        [(3, 3)],
                        [(3, 3)], [(3, 3)],
                        [(3.333333, 3)], [(2, 2.5), (4, 3.5)],
                        [(3, 3)],
                        [(3, 3)], [(2.888889, 3)],
                        [(3.052632, 3.052632)], [(2.904762, 3.095238)],
                        [(3.5, 3)],
                        [(3.25, 5.25), (8, 3), (7.5, 7.5)],
                        [(5.333333, 5.444444)],
                        [(3.833333, 4.333333), (8.5, 7.625)],
                        [(5.909091, 5.500000) , (6, 5.5)],
                        [(4.633333, 5.633333), (7.8, 5.3)])

    for T in (Int, Bool, Gray{Bool}, Gray{N0f8}, Gray{N0f16}, Gray{N0f32}, Gray{Float64})
        for i = 1:20
            test_image = eval(Symbol("test_image_$(i)"))
            labels = Images.label_components(test_image,trues(3,3))
            t = measure_components((RegionEllipse(),), labels)
            if i > 1
                semi_axes = select(t,:centroid)
                c = collect.(semi_axes)
                @test all(isapprox(c, collect.(test_image_centroids[i]); atol = 1e-5))
            end
        end
    end
end

@testset "eccentricity" begin

    test_ellipse_image_eccentricities = ([3.650024149988857e-8],
                                         [0.7323338372368249],
                                         [0.8694282672026615],
                                         [5.960464477539063e-8],
                                         [0.9428090415820629])

    for T in (Int, Bool, Gray{Bool}, Gray{N0f8}, Gray{N0f16}, Gray{N0f32}, Gray{Float64})
        j = 1
        for i = 10:14
            test_image = eval(Symbol("test_ellipse_$(i)"))
            labels = Images.label_components(test_image,trues(3,3))
            t = measure_components((RegionEllipse(),), labels)
            e = select(t,:eccentricity)
            @test isapprox(e, test_ellipse_image_eccentricities[j]; atol = 1e-5)
            j += 1
        end
    end
end

# # TODO Revisit tests once we change the convention with which we represent the
# # ellipse orientation.
# @testset "ellipse orientation" begin
#
#     # test_image_orientations = ([0], [0],
#     #                     [0],
#     #                     [0, 0, 0, 0],
#     #                     [0],
#     #                     [0], [0],
#     #                     [0], [90],
#     #                     [0],
#     #                     [0], [0],
#     #                     [45], [-45],
#     #                     [0],
#     #                     [6.595305, 0, 0],
#     #                     [85.571021],
#     #                     [38.523618 , -28.997308],
#     #                     [90.000000 , 90.000000],
#     #                     [-2.443658, -62.871225])
#
#     for T in (Int, Bool, Gray{Bool}, Gray{N0f8}, Gray{N0f16}, Gray{N0f32}, Gray{Float64})
#         for i = 1:9
#             test_image = eval(Symbol("test_ellipse_$(i)"))
#             labels = Images.label_components(test_image,trues(3,3))
#             t = measure_components((RegionEllipse(),), labels)
#             θ = select(t,:orientation)
#             println(" $i => ")
#             Base.display(θ)
#             #Base.display(test_image_orientations[i])
#             println(" \n ")
#             #@test all(isapprox(θ, test_image_orientations[i]); atol = 5)
#         end
#     end
# end


text:
# ---
# title: 1323. Maximum 69 Number
# id: problem1323
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Math
# link: <https://leetcode.com/problems/maximum-69-number/description/>
# hidden: true
# ---
#
# Given a positive integer `num` consisting only of digits 6 and 9.
#
# Return the maximum number you can get by changing **at most** one digit (6
# becomes 9, and 9 becomes 6).
#
#
#
# **Example 1:**
#
#
#
#     Input: num = 9669
#     Output: 9969
#     Explanation:
#     Changing the first digit results in 6669.
#     Changing the second digit results in 9969.
#     Changing the third digit results in 9699.
#     Changing the fourth digit results in 9666.
#     The maximum number is 9969.
#
#
# **Example 2:**
#
#
#
#     Input: num = 9996
#     Output: 9999
#     Explanation: Changing the last digit 6 to 9 results in the maximum number.
#
# **Example 3:**
#
#
#
#     Input: num = 9999
#     Output: 9999
#     Explanation: It is better not to apply any change.
#
#
#
# **Constraints:**
#
#   * `1 <= num <= 10^4`
#   * `num`'s digits are 6 or 9.
#
#
## @lc code=start
using LeetCode

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


text:
x =3965421
y = 569705405152896
#operadores
+x
-x
x + y
x - y
x * y
x / y
x ÷ y
x \ y
x ^ y
x % y
2 % 3
#not
!true
3*2/12
-x+2

#Bitwise
~x #bitwise not
bitstring(x)
bitstring(y)

bitstring(a)
bitstring(b)
bitstring(c)
a = x & y #bitwise and
b = x | y #bitwise or
c = x ⊻ y #bitwise xor
x >>> 5
bitstring(x >>> 5)  #logical shift right
(x >> 5) < x
bitstring(x >> 5)#arithmetic shift right
(x << 5) > x

bitstring(x << 5) #log/arithm shift left

123 ⊻ 234
xor(123, 234)
~UInt32(123)

#updating operator
x += 1
# += -= /= \= ÷= %= ^= &= |= ⊻= >>>= >>= <<=

#"dot" operators
#elemento por elemento
[1 2 3] ^ 3
[1 2 3] .^3
[1,2,3] .* [2,3,4]
⊗(A,B) = kron(A,B) #produto de Kronecker
[1 2 ; 3 4] ⊗ [0 5 ; 6 7]
[1 2 ; 3 4] .⊗ [0 5 ; 6 7]
#Comparações
# ==  !=,≠  <  <=,≤  >  >=,≥
1≤2
Inf>NaN
NaN == NaN
NaN < NaN
[1 NaN] == [1 NaN]
isequal(NaN, NaN)
isfinite(x)
isinf(x)
isnan(x)
-0.0 == 0.0
isequal(-0.0, 0.0)

v(x) = (println(x); x)
v(1) < v(2) ≤ v(3)

0 .< [0 .5 1.5 .8] .< 1
sin.([1 2 3])

#ordem de precedencia na pag 51(71)
Base.operator_precedence(:+), Base.operator_precedence(:*), Base.operator_precedence(:.)
Base.operator_precedence(:(=))
Base.operator_precedence(:sin)

Base.operator_associativity(:-), Base.operator_associativity(:+)
Base.operator_associativity(:^)
Base.operator_associativity(:sin)

#Note that symbols such as :sin return precedence 0.
#This value represents invalid operators and not operators
#of lowest precedence.
#Similarly, such operators are assigned associativity :none.

#Conversions
Int8(127.0)
Int8(128)
Int8(3.14)
127 % Int8
128 % Int8
round(Int8, 127.4)
round(Int8, 127.6)
 #ver Conversion and Promotion
typeof(round(2.5))
x=7.568

#Divisions
floor(x)
ceil(x)
trunc(x)
y=2.0
x/y
div(x,y) #truncated division
fld(x,y)
cld(x,y)
rem(x,y)
x == div(x,y)y + rem(x,y)
mod(x,y)
x == fld(x,y)y + mod(x,y)
mod1(x,y)
mod2pi(x)
0 <= mod2pi(x) < 2π
divrem(x,y)
fldmod(x,y)
x=49;y=56
gcd(x,y) #MDC
lcm(x,y) #MMC
x=-49
#Sinal abs
abs(x)
abs2(x)
sign(x)
signbit(x)
copysign(y,x)
flipsign(x,y)
#x com sinal de x*y

sin(x)
cos(x)
tan(x)
cot(x)
sec(x)
csc(x)

sqrt(abs(x))
√abs(x)
cbrt(abs(x))
∛abs(x)
hypot(3,4)
exp(1)
expm1(1) #e^x - 1
ldexp(1.0,8)
log(2)
log(ℯ,2)
log2(2)
log10(10)
log1p(0.0052)
#log(1+x) for x near zero
exponent(128.0)
significand(128.0)

#sinh cosh tanh coth sech csch
#asin acos atan acot asec acsc
#asinh acosh atanh acoth asech acsch
#sinc cosc
#These are all single-argument functions, with atan also accepting two
#arguments corresponding to a traditional atan2 function.

#Additionally, sinpi(x) and cospi(x) are provided for more accurate
#computations of sin(pi*x) and cos(pi*x) respectively.

#In order to compute trigonometric functions with degrees instead of radians,
#suffix the function with d. For example, sind(x) computes the sine of x
#where x is specified in degrees.

#Mais funções matemàticas especiais em SpecialFunctions.jl

#endcap

max_stars_repo_name: Mutanne/hiworld

max_issues_repo_name: Mutanne/hiworld

max_forks_repo_name: Mutanne/hiworld


text:
using Sockets
using JSON
import Base.run, Base.close

struct CompletionServer
    s::Base.LibuvServer
    debugmod::Bool
end

function CompletionServer(host::IPAddr, port::Integer, debugmod::Bool)
    CompletionServer(listen(host, port), debugmod)
end

function CompletionServer(port::Integer, debugmod::Bool)
    CompletionServer(listen(port), debugmod)
end

function close(server::CompletionServer)
    close(server.s)
end

function cmdparser(io, cmd)
    opt = cmd[1]
    args = cmd[2:end]
    if opt == "-c"
        writecompletion(io, args...)
    elseif opt == "-e"
        writeeval(io, args...)
    end
    nothing
end

function run(server::CompletionServer)
    stream = accept(server.s)
    while isopen(stream)
        try
            cmd = JSON.Parser.parse(stream)
            if server.debugmod
                println(cmd)
            end
            cmdparser(stream, cmd)
        catch error
            println(error)
            stream = accept(server.s)
        end
    end
    close(server)
end

function serverstart(host::AbstractString, port::Integer, debugmod::Bool=false)
    if host == "localhost"
        host = Sockets.localhost
    else
        host = IPv4(host)
    end
    server = CompletionServer(host, port, debugmod)
    run(server)
end


text:
abstract type MaxwellOperator3D <: IntegralOperator end
abstract type MaxwellOperator3DReg <: MaxwellOperator3D end

struct KernelValsMaxwell3D{T,U,P,Q}
    "gamma = im * wavenumber"
    gamma::U
    vect::P
    dist::T
    green::U
    gradgreen::Q
end

const inv_4pi = 1/(4pi)
function kernelvals(biop::MaxwellOperator3D, p, q)

    γ = biop.gamma
    r = cartesian(p) - cartesian(q)
    R = norm(r)
    γR = γ*R

    inv_R = 1/R

    expn = exp(-γR)
    green = expn * inv_R * inv_4pi
    gradgreen = -(γ + inv_R) * green * inv_R * r

    KernelValsMaxwell3D(γ, r, R, green, gradgreen)
end

function kernelvals(kernel::MaxwellOperator3DReg, p, q)

    γ = kernel.gamma
    r = p.cart - q.cart
    R = norm(r)
    γR = γ*R

    Exp = exp(-γ*R)
    green = (Exp - 1 + γR - 0.5*γR^2) / (4pi*R)
    gradgreen = ( - (γR + 1)*Exp + (1 - 0.5*γR^2) ) * (r/R^3) / (4π)

    KernelValsMaxwell3D(γ, r, R, green, gradgreen)
end

struct MWSingleLayer3D{T,U} <: MaxwellOperator3D
  gamma::T
  α::U
  β::U
end

MWSingleLayer3D(gamma)  = MWSingleLayer3D(gamma, -gamma, -1/(gamma))
MWWeaklySingular(gamma) = MWSingleLayer3D(gamma, 1, 0)
MWHyperSingular(gamma)  = MWSingleLayer3D(gamma, 0, 1)

export Maxwell3D

struct MWSingleLayer3DReg{T,U} <: MaxwellOperator3DReg
    gamma::T
    α::U
    β::U
end

struct MWSingleLayer3DSng{T,U} <: MaxwellOperator3D
    gamma::T
    α::U
    β::U
end

scalartype(op::MaxwellOperator3D) = typeof(op.gamma)

regularpart(op::MWSingleLayer3D) = MWSingleLayer3DReg(op.gamma, op.α, op.β)
singularpart(op::MWSingleLayer3D) = MWSingleLayer3DSng(op.gamma, op.α, op.β)

function _legendre(n,a,b)
    x, w = FastGaussQuadrature.gausslegendre(n)
    w .*= (b-a)/2
    x .= (x.+1)/2*(b-a).+a
    collect(zip(x,w))
end

defaultquadstrat(op::MaxwellOperator3D, tfs, bfs) = DoubleNumWiltonSauterQStrat(2,3,6,7,5,5,4,3)

function quaddata(op::MaxwellOperator3D,
    test_local_space::RefSpace, trial_local_space::RefSpace,
    test_charts, trial_charts, qs::DoubleNumWiltonSauterQStrat)

    tqd = quadpoints(test_local_space,  test_charts,  (qs.outer_rule_far,qs.outer_rule_near))
    bqd = quadpoints(trial_local_space, trial_charts, (qs.inner_rule_far,qs.inner_rule_near))
    leg = (
      _legendre(qs.sauter_schwab_common_vert,0,1),
      _legendre(qs.sauter_schwab_common_edge,0,1),
      _legendre(qs.sauter_schwab_common_face,0,1),)

    # High accuracy rules (use them e.g. in LF MFIE scenarios)
    # tqd = quadpoints(test_local_space, test_charts, (8,8))
    # bqd = quadpoints(trial_local_space, trial_charts, (8,9))
    # leg = (_legendre(8,a,b), _legendre(10,a,b), _legendre(5,a,b),)

    return (tpoints=tqd, bpoints=bqd, gausslegendre=leg)
end

# use Union type so this code can be shared between the operator
# and its regular part.
MWSL3DGen = Union{MWSingleLayer3D,MWSingleLayer3DReg}
function integrand(biop::MWSL3DGen, kerneldata, tvals, tgeo, bvals, bgeo)

  gx = tvals[1]
  fy = bvals[1]

  dgx = tvals[2]
  dfy = bvals[2]

  G = kerneldata.green
  γ = kerneldata.gamma

  α = biop.α
  β = biop.β

  t = (α * dot(gx, fy) + β * (dgx*dfy)) * G
end

struct MWDoubleLayer3D{T} <: MaxwellOperator3D
  gamma::T
end

struct MWDoubleLayer3DSng{T} <: MaxwellOperator3D
    gamma::T
end

struct MWDoubleLayer3DReg{T} <: MaxwellOperator3DReg
    gamma::T
end

regularpart(op::MWDoubleLayer3D) = MWDoubleLayer3DReg(op.gamma)
singularpart(op::MWDoubleLayer3D) = MWDoubleLayer3DSng(op.gamma)

const MWDL3DGen = Union{MWDoubleLayer3D,MWDoubleLayer3DReg}
function integrand(biop::MWDL3DGen, kerneldata, tvals, tgeo, bvals, bgeo)
    g = tvals[1]
    f = bvals[1]
    ∇G = kerneldata.gradgreen
    (f × g) ⋅ ∇G
end

# quadrule(op::MaxwellOperator3D, g::RTRefSpace, f::RTRefSpace, i, τ, j, σ, qd) = qrss(op, g, f, i, τ, j, σ, qd)

# function qrss(op, g, f, i, τ, j, σ, qd)
function quadrule(op::MaxwellOperator3D, g::RTRefSpace, f::RTRefSpace,  i, τ, j, σ, qd,
      qs::DoubleNumWiltonSauterQStrat)

    hits = 0
    dtol = 1.0e3 * eps(eltype(eltype(τ.vertices)))
    dmin2 = floatmax(eltype(eltype(τ.vertices)))
    for t in τ.vertices
        for s in σ.vertices
            d2 = LinearAlgebra.norm_sqr(t-s)
            dmin2 = min(dmin2, d2)
            hits += (d2 < dtol)
        end
    end

    hits == 3 && return SauterSchwabQuadrature.CommonFace(qd.gausslegendre[3])
    hits == 2 && return SauterSchwabQuadrature.CommonEdge(qd.gausslegendre[2])
    hits == 1 && return SauterSchwabQuadrature.CommonVertex(qd.gausslegendre[1])

    h2 = volume(σ)
    xtol2 = 0.2 * 0.2
    k2 = abs2(op.gamma)
    max(dmin2*k2, dmin2/16h2) < xtol2 && return WiltonSERule(
        qd.tpoints[2,i],
        DoubleQuadRule(
            qd.tpoints[2,i],
            qd.bpoints[2,j],),)
    return DoubleQuadRule(
        qd.tpoints[1,i],
        qd.bpoints[1,j],)
end

function quadrule(op::MaxwellOperator3D, g::BDMRefSpace, f::BDMRefSpace,  i, τ, j, σ, qd,
  qs::DoubleNumWiltonSauterQStrat)

  hits = 0
  dtol = 1.0e3 * eps(eltype(eltype(τ.vertices)))
  dmin2 = floatmax(eltype(eltype(τ.vertices)))
  for t in τ.vertices
      for s in σ.vertices
          d2 = LinearAlgebra.norm_sqr(t-s)
          dmin2 = min(dmin2, d2)
          hits += (d2 < dtol)
      end
  end

  hits == 3 && return SauterSchwabQuadrature.CommonFace(qd.gausslegendre[3])
  hits == 2 && return SauterSchwabQuadrature.CommonEdge(qd.gausslegendre[2])
  hits == 1 && return SauterSchwabQuadrature.CommonVertex(qd.gausslegendre[1])

  h2 = volume(σ)
  xtol2 = 0.2 * 0.2
  k2 = abs2(op.gamma)
  # max(dmin2*k2, dmin2/16h2) < xtol2 && return WiltonSERule(
  #     qd.tpoints[2,i],
  #     DoubleQuadRule(
  #         qd.tpoints[2,i],
  #         qd.bpoints[2,j],),)
  return DoubleQuadRule(
      qd.tpoints[1,i],
      qd.bpoints[1,j],)
end

function qrib(op::MaxwellOperator3D, g::RTRefSpace, f::RTRefSpace, i, τ, j, σ, qd)
  # defines coincidence of points
  dtol = 1.0e3 * eps(eltype(eltype(τ.vertices)))

  # decides on whether to use singularity extraction
  xtol = 0.2

  k = norm(op.gamma)

  hits = 0
  xmin = xtol
  for t in τ.vertices
    for s in σ.vertices
      d = norm(t-s)
      xmin = min(xmin, k*d)
      if d < dtol
        hits +=1
        break
      end
    end
  end

    hits == 3   && return BogaertSelfPatchStrategy(5)
    hits == 2   && return BogaertEdgePatchStrategy(8, 4)
    hits == 1   && return BogaertPointPatchStrategy(2, 3)
    rmin = xmin/k
    xmin < xtol && return WiltonSERule(
      qd.tpoints[1,i],
      DoubleQuadRule(
        qd.tpoints[2,i],
        qd.bpoints[2,j],
      ),
    )
    return DoubleQuadRule(
      qd.tpoints[1,i],
      qd.bpoints[1,j],
    )

end

function qrdf(op::MaxwellOperator3D, g::RTRefSpace, f::RTRefSpace, i, τ, j, σ, qd)
  # defines coincidence of points
  dtol = 1.0e3 * eps(eltype(eltype(τ.vertices)))

  # decides on whether to use singularity extraction
  xtol = 0.2

  k = norm(op.gamma)

  hits = 0
  xmin = xtol
  for t in τ.vertices
    for s in σ.vertices
      d = norm(t-s)
      xmin = min(xmin, k*d)
      if d < dtol
        hits +=1
        break
      end
    end
  end

  xmin < xtol && return WiltonSERule(
    qd.tpoints[1,i],
    DoubleQuadRule(
      qd.tpoints[2,i],
      qd.bpoints[2,j],
    ),
  )
  return DoubleQuadRule(
    qd.tpoints[1,i],
    qd.bpoints[1,j],
  )

end


text:
using DistributedFactorGraphs

struct NavAbilityBlobStore{T} <: AbstractBlobStore{T}

end


text:
"""
Loads dependencies and bootstraps a Genie app.
"""
module Genie

push!(LOAD_PATH, @__DIR__)

include(joinpath(@__DIR__, "Configuration.jl"))
include(joinpath(@__DIR__, "constants.jl"))

config = Configuration.Settings(app_env = Configuration.DEV)

using Revise

isfile(joinpath(CONFIG_PATH, "plugins.jl")) && include(joinpath(CONFIG_PATH, "plugins.jl"))

push!(LOAD_PATH,  joinpath(@__DIR__, "cache_adapters"),
                  joinpath(@__DIR__, "session_adapters"),
                  RESOURCES_PATH, HELPERS_PATH)

include(joinpath(@__DIR__, "genie_types.jl"))

include("Loggers.jl")
include("Inflector.jl")
include("Util.jl")
include("FileTemplates.jl")
include("Generator.jl")
include("Tester.jl")
include("Encryption.jl")
include("Cookies.jl")
include("Sessions.jl")
include("Input.jl")
include("Flax.jl")
include("Renderer.jl")
include("Router.jl")
include("Helpers.jl")
include("WebChannels.jl")
include("AppServer.jl")
include("Commands.jl")
include("Cache.jl")
include("Responses.jl")

using .Loggers
using .Inflector, .Util
using .FileTemplates, .Generator, .Tester, .Encryption, .Cookies, .Sessions
using .Input, .Renderer, .Router, .Helpers, .AppServer, .Commands
using .Flax, .AppServer

include(joinpath(@__DIR__, "REPL.jl"))

"""
    run() :: Nothing

Runs the Genie app by parsing the command line args and invoking the corresponding actions.
"""
function run() :: Nothing
  Commands.execute(Genie.config)

  nothing
end

"""
    newapp(path = "."; db_support = false, skip_dependencies = false, autostart = true) :: Nothing

Scaffolds a new Genie app.
"""
function newapp(path = "."; db_support = false, autostart = true) :: Nothing
  REPL.newapp(path, db_support = db_support, autostart = autostart)
end
const new_app = newapp

"""
"""
function loadapp(path = "."; autostart = false) :: Nothing
  REPL.loadapp(path, autostart = autostart)
end

"""
"""
function startup(port::Int = 8000, host::String = Genie.config.server_host;
                  ws_port::Int = port + 1, async::Bool = ! Genie.config.run_as_server,
                  verbose::Bool = false, ratelimit::Rational{Int} = 10_000//1)
  AppServer.startup(port, host, ws_port = ws_port, async = async, verbose = verbose, ratelimit = ratelimit)
end
const startapp = startup

end


text:
module TessellationWrapper

using JavaCall

const AbsoluteDate = @jimport org.orekit.time.AbsoluteDate
const Action = @jimport org.hipparchus.ode.events.Action
const AlongTrackAiming = @jimport org.orekit.models.earth.tessellation.AlongTrackAiming
const BSPTree = @jimport org.hipparchus.geometry.partitioning.BSPTree
const BSPTreeVisitor_Order = @jimport org.hipparchus.geometry.partitioning.BSPTreeVisitor$Order
const Class = @jimport java.lang.Class
const ConstantAzimuthAiming = @jimport org.orekit.models.earth.tessellation.ConstantAzimuthAiming
const Direction = @jimport org.orekit.models.earth.tessellation.Direction
const DivertedSingularityAiming = @jimport org.orekit.models.earth.tessellation.DivertedSingularityAiming
const EllipsoidTessellator = @jimport org.orekit.models.earth.tessellation.EllipsoidTessellator
const Enum = @jimport java.lang.Enum
const EventDetector = @jimport org.orekit.propagation.events.EventDetector
const EventHandler = @jimport org.orekit.propagation.events.handlers.EventHandler
const GeodeticPoint = @jimport org.orekit.bodies.GeodeticPoint
const HalfTrackSampler = @jimport org.orekit.models.earth.tessellation.HalfTrackSampler
const HalfTrackSpanHandler = @jimport org.orekit.models.earth.tessellation.HalfTrackSpanHandler
const InsideFinder = @jimport org.orekit.models.earth.tessellation.InsideFinder
const JString = @jimport java.lang.String
const LatitudeExtremumDetector = @jimport org.orekit.propagation.events.LatitudeExtremumDetector
const List = @jimport java.util.List
const Mesh = @jimport org.orekit.models.earth.tessellation.Mesh
const Mesh_Node = @jimport org.orekit.models.earth.tessellation.Mesh$Node
const Object = @jimport java.lang.Object
const OneAxisEllipsoid = @jimport org.orekit.bodies.OneAxisEllipsoid
const Optional = @jimport java.util.Optional
const Orbit = @jimport org.orekit.orbits.Orbit
const OrekitFixedStepHandler = @jimport org.orekit.propagation.sampling.OrekitFixedStepHandler
const S2Point = @jimport org.hipparchus.geometry.spherical.twod.S2Point
const SpacecraftState = @jimport org.orekit.propagation.SpacecraftState
const SphericalPolygonsSet = @jimport org.hipparchus.geometry.spherical.twod.SphericalPolygonsSet
const Tile = @jimport org.orekit.models.earth.tessellation.Tile
const TileAiming = @jimport org.orekit.models.earth.tessellation.TileAiming
const Vector3D = @jimport org.hipparchus.geometry.euclidean.threed.Vector3D

export AlongTrackAiming
export ConstantAzimuthAiming
export Direction
export DivertedSingularityAiming
export EllipsoidTessellator
export HalfTrackSampler
export HalfTrackSpanHandler
export InsideFinder
export Mesh
export Tile
export TileAiming
export add_node
export along_tile_direction
export build_simple_zone
export compare_to
export describe_constable
export equals
export event_occurred
export get_across_gap
export get_along_gap
export get_center
export get_class
export get_closest_existing_node
export get_coverage
export get_declaring_class
export get_end
export get_half_track
export get_inside_nodes
export get_inside_point
export get_interpolated_point
export get_max_across_index
export get_max_along_index
export get_min_across_index
export get_min_along_index
export get_node
export get_number_of_nodes
export get_singular_points
export get_start
export get_taxicab_boundary
export get_vertices
export handle_step
export hash_code
export init
export motion
export name
export neighbor_across_index
export neighbor_along_index
export next
export notify
export notify_all
export ordinal
export reset_state
export sample
export tessellate
export to_string
export value_of
export values
export visit_internal_node
export visit_leaf_node
export visit_order
export wait

include("along_track_aiming.jl")
include("constant_azimuth_aiming.jl")
include("direction.jl")
include("diverted_singularity_aiming.jl")
include("ellipsoid_tessellator.jl")
include("half_track_sampler.jl")
include("half_track_span_handler.jl")
include("inside_finder.jl")
include("mesh.jl")
include("tile.jl")
include("tile_aiming.jl")

end


text:
# File for generating data for Figure 2

using Pkg
using Random, Distributions
using ProgressMeter
using DataFrames
using Gadfly
using StatsBase

include("/home/christian/UniMA/EMMA/src/RFR.jl")
include("/home/christian/UniMA/EMMA/src/cross_val.jl")
include("/home/christian/UniMA/EMMA/src/aux_functions.jl")

# Friedman function
function friedman(x::Matrix, errors::Matrix)

    Y = 10 .* sin.(π .* x[:,1] .* x[:,2]) .+ 20 .* (x[:,3] .- 0.5).^2 .+ 10 .* x[:,4] + 5 .* x[:,5] .+ errors

    return Y
end

# Make data function with parameters
function make_data(n, d, func, σ)

    x_train = rand(Uniform(0, 1), n, d)
    x_test = rand(Uniform(0, 1), n, d)

    d = Normal(0, σ)
    td = truncated(d, -Inf, Inf)

    errors_train = rand(td, n, 1)
    errors_test = zeros(n, 1)

    if func=="friedman"
        y_train = friedman(x_train, errors_train)
        y_test = friedman(x_test, errors_test)
    else
        error("Provide function to compute Y")
    end

    return x_train, x_test, y_train, y_test
end

# Set seed for reproducability
Random.seed!(68151)

# Set number of runs and allocate space
n_runs = 100
result_arr = Array{Float64}(undef, 0, 2)
σ = 8 # Setting error term variance

@showprogress for run in 1:n_runs

    # Setting number of observations and dimensions
    n = 2000
    d = 5 # Only taking strong variables, we are in subsection: One dimension

    x_train, x_test, y_train, y_test = make_data(n, d, "friedman", σ)
    a_list = collect(LinRange(0, 30, 31))

    # Make parameter dictionary
    d1 = Dict{Symbol, Vector{Float64}}(
        :max_features => [d],
        :n_trees => [100],
        :α => [0.0])

    rf = RFR(param_dict = d1)
    fit!(rf, x_train, y_train)

    for i in 1:length(rf.trees)

        # Save tree depth and P_(t_L)
        res = zeros(length(rf.trees[i].depth_list), 2)
        res[:,1] = rf.trees[i].depth_list
        res[:,2] = rf.trees[i].pl

        result_arr = vcat(result_arr, res)
    end
end

# Only considering non-terminal nodes
plot_data = result_arr[result_arr[:,2].!=-2, :]
# Calculating λ = 4 * P(t_L) * (1 - P(t_L))
plot_data[:,2] = abs.(plot_data[:,2] .- 0.5)
plot_data[:,2] = 4 .* (plot_data[:,2]) .* (1 .- plot_data[:,2])

# Make dataframe and save results
plot_data = DataFrame(plot_data, :auto)
rename!(plot_data, ["depth", "deviation_from_median"])

jldsave("data/deep_sigma8.jld2"; plot_data)

max_stars_repo_name: christianhilscher/emma

max_issues_repo_name: christianhilscher/emma

max_forks_repo_name: christianhilscher/emma


text:
using DelimitedFiles

C = Vector{Matrix{Float64}}(undef, 6)
for j = 1:size(C,2)
    C[j] = readdlm("data/u$j.csv",',')
end

Fun(coefficients(C[1][:,1] + im*C[1][:,2], Chebyshev(), Legendre())


text:
using PyPlot
using Convex
using PyCall
using HDF5, JLD
@pyimport matplotlib2tikz

#using CenterOfMass
include("CenterOfMass.jl")
include("setup_parameters.jl")

interior_q(x) = norm(x) <= 1

masses = [1.0]
resolution = 0.02
prior = initialize_prior(circle_ps, resolution, interior_q, masses)

critical_forces = get_critical_values(circle_ws, attachment_ws[5], prior, interior_q)

circle_polygon = hcat(circle_ps..., circle_ps[1])
#plot(circle_polygon[1,:][:], circle_polygon[2,:][:], color = "k", alpha = 1, linestyle = "solid")

xs = get_range(prior, 1)
ys = get_range(prior, 2)

xg = repmat(xs', length(ys), 1)
yg = repmat(ys, 1, length(xs))

critical_forces = max(critical_forces, zeros(size(critical_forces)))
plot_surface(xg, yg, critical_forces[:,:], cmap="BuPu", rstride=2, cstride=2, alpha=0.8, linewidth=0.25)

#extent= [minimum(xs), maximum(xs), minimum(ys), maximum(ys)]

#imshow(critical_forces[:,:]', cmap= "BuPu", vmin=minimum(critical_forces), vmax=maximum(critical_forces),
  #extent = extent, interpolation="nearest", origin="lower")

max_stars_repo_name: mcorah/CenterOfMass

max_issues_repo_name: mcorah/CenterOfMass

max_forks_repo_name: mcorah/CenterOfMass


text:
export pairwise

import Base: show

"""
    pairwise([ initial ])

Creates a pairwise operator, which groups pairs of consecutive emissions together and emits them as a tuple of two values.
Accepts optional initial seed value to start pairing from.

```jldoctest
using Rocket

source = from(1:5) |> pairwise()

subscribe!(source, logger())
;

# output

[LogActor] Data: (1, 2)
[LogActor] Data: (2, 3)
[LogActor] Data: (3, 4)
[LogActor] Data: (4, 5)
[LogActor] Completed
```

```jldoctest
using Rocket

source = from(1:5) |> pairwise(0)

subscribe!(source, logger())
;

# output

[LogActor] Data: (0, 1)
[LogActor] Data: (1, 2)
[LogActor] Data: (2, 3)
[LogActor] Data: (3, 4)
[LogActor] Data: (4, 5)
[LogActor] Completed
```

See also: [`AbstractOperator`](@ref), [`InferableOperator`](@ref), [`logger`](@ref)
"""
pairwise()        = PairwiseOperator()
pairwise(initial) = start_with(initial) + pairwise()

struct PairwiseOperator <: InferableOperator end

operator_right(::PairwiseOperator, ::Type{L}) where L = Tuple{L, L}

function on_call!(::Type{L}, ::Type{ Tuple{L, L} }, ::PairwiseOperator, source) where L
    return proxy(Tuple{L, L}, source, PairwiseProxy{L}())
end

struct PairwiseProxy{L} <: ActorProxy end

actor_proxy!(::Type{ Tuple{L, L} }, ::PairwiseProxy{L}, actor::A) where { L, A } = PairwiseActor{L, A}(actor, nothing)

mutable struct PairwiseActor{L, A} <: Actor{L}
    actor    :: A
    previous :: Union{Nothing, L}
end

function on_next!(actor::PairwiseActor{L}, data::L) where L
    if actor.previous !== nothing
        next!(actor.actor, (actor.previous, data))
    end
    actor.previous = data
end

on_error!(actor::PairwiseActor, err) = error!(actor.actor, err)
on_complete!(actor::PairwiseActor)   = complete!(actor.actor)

Base.show(io::IO, ::PairwiseOperator)          = print(io, "PairwiseOperator()")
Base.show(io::IO, ::PairwiseProxy)             = print(io, "PairwiseProxy()")
Base.show(io::IO, ::PairwiseActor{L})  where L = print(io, "PairwiseActor($L)")


text:
struct Highway
    n_lanes::Int
    lanewidth::Float64
    len::Float64
end


text:
add_verbose_scope(:AbExt)
add_assert_scope(:AbExt)

add_verbose_scope(:MaxAbExt)

export abelian_fields, abelian_normal_extensions, abelian_extensions

###############################################################################
#
#  Abelian extensions
#
###############################################################################

function abelian_fields(O::Union{FlintIntegerRing, FlintRationalField},
                            gtype::Vector{Int}, discriminant_bound::fmpz;
                            only_real::Bool = false,
                            tame::Bool = false)

  Qx, x = PolynomialRing(FlintQQ, "x", cached = false)
  K, _ = NumberField(x - 1, "a", cached = false)
  OK = maximal_order(K)
  l = abelian_fields(OK, gtype, discriminant_bound,
                         only_real = only_real,
                         tame = tame)
  return l
end

function abelian_fields(gtype::Vector{Int}, conds::Vector{Int}, absolute_discriminant_bound::fmpz; only_real::Bool = false)
  K = rationals_as_number_field()[1]
  O = maximal_order(K)
  gtype = map(Int, snf(abelian_group(gtype))[1].snf)
  n = prod(gtype)

  #Getting conductors
  @vprint :AbExt 1 "Number of conductors: $(length(conds)) \n"
  fields = ClassField{MapRayClassGrp, GrpAbFinGenMap}[]

  #Now, the big loop
  fun = (x, y) -> quo(x, y, false)[2]
  for (i, k) in enumerate(conds)
    @vprint :AbExt 1 "Conductor: $k \n"
    @vprint :AbExt 1 "Left: $(length(conds) - i)\n"
    r, mr = Hecke.ray_class_groupQQ(O, k, !only_real, gtype[end])
    if !has_quotient(r, gtype)
      continue
    end
    ls = subgroups(r, quotype = gtype, fun = fun)
    for s in ls
      C = ray_class_field(mr, s)
      if Hecke._is_conductor_minQQ(C, n) && Hecke.discriminant_conductorQQ(O, C, k, absolute_discriminant_bound)
        @vprint :AbExt 1 "New Field \n"
        push!(fields, C)
      end
    end
  end
  return fields
end

function abelian_fields(O::NfOrd, gtype::Vector{Int}, absolute_discriminant_bound::fmpz; only_real::Bool = false, only_complex::Bool = false, tame::Bool = false)
  K = nf(O)
  @assert degree(K)==1
  gtype = map(Int, snf(abelian_group(gtype))[1].snf)
  n = prod(gtype)
  real = only_real

  expo = lcm(gtype)

  #Getting conductors
  l_conductors = conductorsQQ(O, gtype, absolute_discriminant_bound, tame)
  @vprint :AbExt 1 "Number of conductors: $(length(l_conductors)) \n"
  fields = ClassField{MapRayClassGrp, GrpAbFinGenMap}[]

  #Now, the big loop
  for (i, k) in enumerate(l_conductors)
    @vprint :AbExt 1 "Conductor: $k \n"
    @vprint :AbExt 1 "Left: $(length(l_conductors) - i)\n"
    r, mr = Hecke.ray_class_groupQQ(O, k, !real, expo)
    if !has_quotient(r, gtype)
      continue
    end
    ls = subgroups(r,quotype = gtype, fun = (x, y) -> quo(x, y, false)[2])
    for s in ls
      C = ray_class_field(mr, s)
      if Hecke._is_conductor_minQQ(C,n) && Hecke.discriminant_conductorQQ(O,C,k,absolute_discriminant_bound)
        if only_complex
          rC, sC = signature(C)
          if !iszero(rC)
            continue
          end
        end
        @vprint :AbExt 1 "New Field \n"
        push!(fields, C)
      end
    end
  end
  return fields
end

function abelian_normal_extensions(K::AnticNumberField, gtype::Vector{Int}, absolute_discriminant_bound::fmpz; only_real::Bool = false, only_complex::Bool = false, tame::Bool = false, absolute_galois_group::Tuple{Int, Int} = (0, 0))

  @assert !(only_real && only_complex)
  O = maximal_order(K)
  d = degree(K)
  if d == 1
    return abelian_fields(O, gtype, absolute_discriminant_bound, only_real = only_real, tame = tame)
  end
  gtype = map(Int, snf(abelian_group(gtype))[1].snf)
  n = prod(gtype)
  inf_plc = InfPlc[]

  fields = ClassField{MapRayClassGrp, GrpAbFinGenMap}[]
  bound = div(absolute_discriminant_bound, abs(discriminant(O))^n)

  if iszero(bound)
    return fields
  end

  if iseven(n) && !only_real
    inf_plc = real_places(K)
  end

  Aut = automorphisms(K)
  @assert length(Aut) == d #The field is normal over Q
  gs = Hecke.small_generating_set(Aut)

  expo = lcm(gtype)
  Cl, mCl = class_group(O)
  cgrp = !iscoprime(n, order(Cl))
  allow_cache!(mCl)

  #Getting conductors
  l_conductors = conductors(O, gtype, bound, tame)
  @vprint :AbExt 1 "Number of conductors: $(length(l_conductors)) \n"

  ctx = rayclassgrp_ctx(O, expo)
  #Now, the big loop
  for (i, k) in enumerate(l_conductors)
    @vprint :AbExt 1 "Conductor: $k \n"
    @vprint :AbExt 1 "Left: $(length(l_conductors) - i)\n"
    r, mr = ray_class_group_quo(O, k[1], k[2], inf_plc, ctx)
    if !has_quotient(r, gtype)
      continue
    end
    act = induce_action(mr, gs)
    ls = stable_subgroups_for_abexts(r, act, gtype)
    for s in ls
      s::GrpAbFinGenMap
      @hassert :AbExt 1 order(codomain(s)) == n
      C = ray_class_field(mr, s)::ClassField{MapRayClassGrp, GrpAbFinGenMap}
      if _is_conductor_min_normal(C) && discriminant_conductor(C, bound)
        if only_complex
          rC, sC = signature(C)
          if iszero(rC)
            continue
          end
        end
        if absolute_galois_group != (0, 0)
          autabs = absolute_automorphism_group(C, gs)
          G = generic_group(autabs, *)[1]
          id_G = find_small_group(G)
          if id_G == absolute_galois_group
            @vprint :AbExt 1 "New Field \n"
            push!(fields, C)
          end
        else
          @vprint :AbExt 1 "New Field \n"
          push!(fields, C)
        end
      end
    end
  end
  return fields
end

################################################################################
#
#  Abelian extensions of arbitrary fields
#
################################################################################

function abelian_extensions(K::AnticNumberField, gtype::Vector{Int}, absolute_discriminant_bound::fmpz; absolutely_distinct::Bool = false, ramified_at_inf_plc::Tuple{Bool, Vector{InfPlc}} = (false, InfPlc[]), only_tame::Bool = false)

  OK = maximal_order(K)
  gtype = map(Int, snf(abelian_group(gtype))[1].snf)
  n = prod(gtype)
  inf_plc = InfPlc[]

  fields = ClassField{MapRayClassGrp, GrpAbFinGenMap}[]
  bound = div(absolute_discriminant_bound, abs(discriminant(OK))^n)
  if iszero(bound)
    return fields
  end

  inf_plc = real_places(K)
  if ramified_at_inf_plc[1]
    inf_plc = ramified_at_inf_plc[2]
  end
  expo = gtype[end]
  auts = automorphisms(K)
  gens_auts = small_generating_set(auts)
  if isone(length(auts))
    absolutely_distinct = false
  end
  Cl, mCl = class_group(OK)
  cgrp = !iscoprime(n, order(Cl))
  allow_cache!(mCl)

  #Getting conductors
  l_conductors = conductors_generic(K, gtype, absolute_discriminant_bound, only_tame = only_tame)
  if absolutely_distinct
    l_conductors = _sieve_conjugates(auts, l_conductors)
  end
  @vprint :AbExt 1 "Number of conductors: $(length(l_conductors)) \n"

  ctx = rayclassgrp_ctx(OK, expo)
  fsub = (x, y) -> quo(x, y, false)[2]
  fsub_distinct = (x, y) -> (quo(x, y, false)[2], sub(x, y, false)[2])
  #Now, the big loop
  for (i, k) in enumerate(l_conductors)
    @vprint :AbExt 1 "Left: $(length(l_conductors) - i)\n"
    r, mr = ray_class_group_quo(OK, k, inf_plc, ctx)
    if !has_quotient(r, gtype)
      continue
    end
    if absolutely_distinct && _isstable(gens_auts, k)
      act = induce_action(mr, gens_auts)
      full_action = closure(act, *)
      ls_with_emb = subgroups(r, quotype = gtype, fun = fsub_distinct)
      _closure = Vector{GrpAbFinGenMap}()
      ls = Vector{GrpAbFinGenMap}()
      for (proj, emb) in ls_with_emb
        found = false
        for j = 1:length(_closure)
          if _issubset(emb, _closure[j])
            found = true
            break
          end
        end
        if !found
          push!(ls, proj)
          for mp in full_action
            push!(_closure, emb*mp)
          end
        end
      end
    else
      ls1 = subgroups(r, quotype = gtype, fun = fsub)
      ls = collect(ls1)::Vector{GrpAbFinGenMap}
    end
    new_fields = ClassField{MapRayClassGrp, GrpAbFinGenMap}[]
    for s in ls
      @hassert :AbExt 1 order(codomain(s)) == n
      C = ray_class_field(mr, s)::ClassField{MapRayClassGrp, GrpAbFinGenMap}
      cC = conductor(C)
      if ramified_at_inf_plc[1]
        if Set(cC[2]) != Set(ramified_at_inf_plc[2])
          continue
        end
      end
      if cC[1] == mr.defining_modulus[1] && norm(discriminant(C)) <= bound
        @vprint :AbExt 1 "New Field \n"
        push!(new_fields, C)
      end
    end
    append!(fields, new_fields)
  end
  return fields
end

function _isstable(auts::Vector{NfToNfMor}, d::Dict{NfOrdIdl, Int})
  if isempty(d)
    return true
  end
  OK = order(first(keys(d)))
  #CAREFUL: BASE FIELD NOT NORMAL.
  #NEED TO ACT WITH automorphisms
  primes = Set{fmpz}(minimum(x, copy = false) for x in keys(d))
  for p in primes
    lP = prime_decomposition(OK, p)
    prime_ideals = NfOrdIdl[x[1] for x in lP]
    perms = [induce_action(prime_ideals, f) for f in auts]
    orbs = orbits(perms)
    for i = 1:length(orbs)
      P = prime_ideals[orbs[i][1]]
      if !haskey(d, P)
        for j = 2:length(orbs[i])
          if haskey(d, prime_ideals[orbs[i][j]])
            return false
          end
        end
      else
        e = d[P]
        for j = 2:length(orbs[i])
          if !haskey(d, prime_ideals[orbs[i][j]]) || d[prime_ideals[orbs[i][j]]] != e
            return false
          end
        end
      end
    end
  end
  return true
end

function _image(cache, auts, I, i)
  pos = Base.ht_keyindex(cache[i], I)
  if pos != -1
    return cache[i].vals[pos]
  end
  img = auts[i](I)
  cache[i][I] = img
  return img
end

function _sieve_conjugates(auts::Vector{NfToNfMor}, conds::Vector{Dict{NfOrdIdl, Int}})
  if isone(length(auts))
    return conds
  end
  closure = Set{Dict{NfOrdIdl, Int}}()
  reps = Vector{Dict{NfOrdIdl, Int}}()
  cache = Vector{Dict{NfOrdIdl, NfOrdIdl}}(undef, length(auts))
  for i = 1:length(cache)
    cache[i] = Dict{NfOrdIdl, NfOrdIdl}()
  end
  for j = 1:length(conds)
    if conds[j] in closure
      continue
    end
    push!(reps, conds[j])
    for i = 1:length(auts)
      push!(closure, _induce_image(auts, i, conds[j], cache))
    end
  end
  return reps
end

function _induce_image(auts::Vector{NfToNfMor}, i::Int, cond::Dict{NfOrdIdl, Int}, cache::Vector{Dict{NfOrdIdl, NfOrdIdl}})
  res = Dict{NfOrdIdl, Int}()
  for (k, v) in cond
    res[_image(cache, auts, k, i)] = v
  end
  return res
end

################################################################################
#
#  Valuation bounds for discriminants
#
################################################################################

function valuation_bound_discriminant(n::Int, p::Union{Integer, fmpz})
  # First compute the p-adic expansion of n
  S = Vector{typeof(p)}()
	q = typeof(p)(n)
  q, r = divrem(q, p)
  push!(S, r)
  while q >= p
    q, r = divrem(q, p)
    push!(S, r)
  end

	if !iszero(q)
		push!(S, q)
	end

	@assert sum(S[i + 1] * p^i for i in 0:length(S)-1) == n

	b = zero(typeof(p))

	for i in 0:length(S) - 1
		b = b + S[i + 1] * (i + 1) * p^i
		if !iszero(S[i + 1])
			b = b - 1
		end
	end

  return b
end

###########################################################################
#
#  Some useful functions
#
###########################################################################

#This function gets a quotient of the ray class group and the action on
# the ray class group
# In output, we get the quotient group as a ZpnGModule

function _action_on_quo(mq::GrpAbFinGenMap, act::Array{GrpAbFinGenMap,1})

  q=mq.header.codomain
  S,mS=snf(q)
  n=Int(S.snf[end])
  R=ResidueField(FlintZZ, n, cached=false)
  quo_action=Array{nmod_mat,1}(undef, length(act))
  for s=1:length(act)
    quo_action[s]= change_base_ring(mS.map*act[i].map*mS.imap, R)
  end
  return ZpnGModule(S, quo_action)

end

###############################################################################
#
#  Quadratic Extensions of Q
#
###############################################################################

function quadratic_fields(bound::Int; tame::Bool=false, real::Bool=false, complex::Bool=false, with_autos::Type{Val{T}}=Val{false}) where T

  @assert !(real && complex)
  Qx,x=PolynomialRing(FlintQQ, "x")
  sqf=squarefree_up_to(bound)
  if real
    deleteat!(sqf,1)
  elseif complex
    sqf=Int[-i for i in sqf]
  else
    sqf= vcat(sqf[2:end], Int[-i for i in sqf])
  end
  if tame
    filter!( x -> mod(x,4)==1, sqf)
    return ( number_field(x^2-x+divexact(1-i,4), cached=false, check = false)[1] for i in sqf)
  end
  final_list=Int[]
  for i=1:length(sqf)
    if abs(sqf[i]*4)< bound
      @views push!(final_list,sqf[i])
      continue
    end
    if mod(sqf[i],4)==1
      @views push!(final_list,sqf[i])
    end
  end
  return ( mod(i,4)!=1 ? number_field(x^2-i, cached=false, check = false)[1] : number_field(x^2-x+divexact(1-i,4), cached = false, check = false)[1] for i in final_list)

end

function _quad_ext(bound::Int, only_real::Bool = false; unramified_outside::Vector{fmpz} = fmpz[])

  Qx, x = PolynomialRing(FlintQQ, "x", cached = false)
  K = NumberField(x-1, cached = false, check = false)[1]
  sqf = squarefree_up_to(bound, prime_base = unramified_outside)
  final_list = Int[]
  for i=2:length(sqf)
    if abs(sqf[i]*4)< bound
      @views push!(final_list, sqf[i])
      continue
    end
    if mod(sqf[i], 4) == 1
      @views push!(final_list, sqf[i])
    end
  end
  if !only_real
    for i = 1:length(sqf)
      if abs(sqf[i]*4)< bound
        @views push!(final_list, -sqf[i])
        continue
      end
      if mod(sqf[i], 4) == 3
        @views push!(final_list, -sqf[i])
      end
    end
  end
  fields_list = Array{Tuple{AnticNumberField, Vector{NfToNfMor}, Vector{NfToNfMor}}, 1}(undef, length(final_list))
  for i = 1:length(final_list)
    if mod(final_list[i],4) != 1
      cp = Vector{fmpz}(undef, 3)
      cp[1] = -final_list[i]
      cp[2] = 0
      cp[3] = 1
      L, gL = number_field(Qx(cp), cached=false, check = false)
      auts = NfToNfMor[hom(L, L, -gL, check = false)]
      emb = NfToNfMor[hom(K, L, one(L), check = false)]
      fields_list[i] = (L, auts, emb)
    else
      cp = Vector{fmpz}(undef, 3)
      cp[1] = divexact(1-final_list[i], 4)
      cp[2] = -1
      cp[3] = 1
      L, gL = number_field(Qx(cp), cached=false, check = false)
      auts = NfToNfMor[hom(L, L, 1-gL, check = false)]
      emb = NfToNfMor[hom(K, L, one(L), check = false)]
      fields_list[i] = (L, auts, emb)
    end
  end
  return fields_list

end

###############################################################################
#
#  C2 x C2 extensions of Q
#
###############################################################################

function C22_extensions(bound::Int)

  Qx, x=PolynomialRing(FlintZZ, "x")
  K, _=NumberField(x-1, cached = false)
  Kx,x=PolynomialRing(K,"x", cached=false)
  b1=ceil(Int,Base.sqrt(bound))
  n=2*b1+1
  pairs = _find_pairs(bound)
  return (_ext(Kx,x,i,j) for (i, j) in pairs)

end

function _ext(Ox,x,i,j)

  y=mod(i,4)
  pol1 = x^2
  if y != 1
    pol1 = x^2-i
  else
    pol1 = x^2-x+divexact(1-i,4)
  end
  y=mod(j,4)
  pol2=Ox(1)
  if y!=1
    pol2=x^2-j
  else
    pol2=x^2-x+divexact(1-j,4)
  end
  return number_field([pol1,pol2], cached = false, check = false)

end

function _C22_exts_abexts(bound::Int, only_real::Bool = false; unramified_outside::Vector{fmpz} = fmpz[])
  Qx, x = PolynomialRing(FlintZZ, "x")
  pairs = _find_pairs(bound, only_real, unramified_outside = unramified_outside)
  return (_ext_with_autos(Qx, x, i, j) for (i, j) in pairs)
end

function _ext_with_autos(Qx, x, i::Int, j::Int)
  first = i
  second = j
  g = gcd(i, j)
  if g != 1
    third = divexact(i*j, g^2)
    if gcd(first, third) == 1
      second = third
    elseif gcd(second, third) == 1
      first = third
    end
  end
  y1 = mod(first, 4)
  cp1 = Vector{fmpz}(undef, 3)
  cp1[3] = fmpz(1)
  if y1 != 1
    cp1[1] = fmpz(-first)
    cp1[2] = fmpz(0)
  else
    cp1[1] = fmpz(divexact(1-first,4))
    cp1[2] = fmpz(-1)
  end
  y2 = mod(second, 4)
  cp2 = Vector{fmpz}(undef, 3)
  cp2[3] = fmpz(1)
  if y2 != 1
    cp2[1] = fmpz(-second)
    cp2[2] = fmpz(0)
  else
    cp2[1] = fmpz(divexact(1-second, 4))
    cp2[2] = fmpz(-1)
  end
  return Qx(cp1), Qx(cp2)
end

function __get_term(a::fmpq_mpoly, exps::Vector{UInt})
   z = fmpq()
   ccall((:fmpq_mpoly_get_coeff_fmpq_ui, libflint), Nothing,
         (Ref{fmpq}, Ref{fmpq_mpoly}, Ptr{UInt}, Ref{FmpqMPolyRing}),
         z, a, exps, parent(a))
   return z
end

function _C22_with_max_ord(l)
  list = Vector{Tuple{AnticNumberField, Vector{NfToNfMor}, Vector{NfToNfMor}}}()
  Qx, x = PolynomialRing(FlintQQ, "x", cached = false)
  K = NumberField(x-1, cached = false)[1]
  for (p1, p2) in l
    Kns, g = number_field(fmpz_poly[p1, p2], check = false, cached = false)
    S, mS = simple_extension(Kns, check = false, cached = false, simplified = true)
    Hecke._assure_has_inverse_data(mS)
    gen1 = mS\(g[1])
    gen2 = mS\(g[2])
    d1 = discriminant(p1)
    d2 = discriminant(p2)
    cf = gcd(d1, d2)
    if isone(cf)
      B = Vector{nf_elem}(undef, 4)
      B[1] = S(1)
      B[2] = mS\(g[1])
      B[3] = mS\(g[2])
      B[4] = B[2] * B[3]
      M = basis_matrix(B, FakeFmpqMat)
      hnf_modular_eldiv!(M.num, M.den, :lowerleft)
      O = NfAbsOrd(S, FakeFmpqMat(M.num, M.den))
      O.disc = d1^2*d2^2
      d3 = numerator(discriminant(S))
      if d3 < 0
        O.disc = -O.disc
      end
      O.index = divexact(d3, O.disc)
      O.ismaximal = 1
      Hecke._set_maximal_order_of_nf(S, O)
    else
      maximal_order(S)
    end
    auts = small_generating_set(automorphisms(S, isabelian = true, copy = false))
    push!(list, (S, auts, NfToNfMor[hom(K, S, S(1), check = false)]))
  end
  return list
end

function _disc(a::Int, b::Int, c::Int, bound::Int)
  a1 = mod(a, 4)
  b1 = mod(b, 4)
  if a1 == 1 && b1 == 1
    return a*b*c <= bound
  end
  c1 = mod(c, 4)
  if a1 == 1 || b1 == 1 || c1 == 1
    return 16*a*b*c <= bound
  else
    return 64*a*b*c <= bound
  end
end

function _pairs_totally_real(pairs, ls, bound)
  b1=floor(Int, Base.sqrt(bound))
  pairs[1, 1:length(ls)] .= false
  pairs[1:length(ls), 1] .= false
  for j = 2:length(ls)
    for i = j:length(ls)
      pairs[i, j] = false
    end
  end
  for j = 2:length(ls)
    second = ls[j]
    for i = 2:j-1
      if pairs[i, j]
        first = ls[i]
        g = gcd(first, second)
        if isone(g)
          third = first*second
        else
          third = divexact(first*second, g^2)
        end
        if abs(third) > b1
          pairs[i, j] = false
          continue
        end
        k = searchsortedfirst(ls, third)
        if i < k
          pairs[i, k] = false
        else
          pairs[k, i] = false
        end
        if j < k
          pairs[j, k] = false
        else
          pairs[k, j] = false
        end
        if !_disc(first, second, third, bound)
          pairs[i, j] = false
        end
      end
    end
  end
  it = findall(pairs)
  totally_real_exts = Vector{Tuple{Int, Int}}(undef, length(it))
  ind = 1
  for I in it
    totally_real_exts[ind] = (ls[I[1]], ls[I[2]])
    ind += 1
  end
  return totally_real_exts

end

function _find_pairs(bound::Int, only_real::Bool = false; unramified_outside::Vector{fmpz} = fmpz[] )

  #first, we need the squarefree numbers
  b1=ceil(Int, Base.sqrt(bound))
  ls = squarefree_up_to(b1, prime_base = unramified_outside)
  #The first step is to enumerate all the totally real extensions.
  pairs = trues(length(ls), length(ls))
  real_exts = _pairs_totally_real(pairs, ls, bound)
  if only_real
    return real_exts
  end
  ls1 = -ls
  #Now, all the others.
  pairs[1:length(ls), 2:length(ls)] .= true
  for j = 2:length(ls)
    second = ls[j]
    for i = 1:length(ls)
      if pairs[i, j]
        first = ls1[i]
        g = gcd(first, second)
        if isone(g)
          third = first*second
        else
          third = divexact(first*second, g^2)
        end
        abt = -third
        if abt > b1
          pairs[i, j] = false
          continue
        end
        k = searchsortedfirst(ls, abt)
        pairs[k, j] = false
        if !_disc(first, second, third, bound)
          pairs[i, j] = false
        end
      end
    end
  end
  it = findall(pairs)
  ind = 1
  res = Vector{Tuple{Int, Int}}(undef, length(it))
  for I in it
    res[ind] = (ls1[I[1]], ls[I[2]])
    ind += 1
  end
  return vcat(res, real_exts)
end

function _from_relative_to_abs(L::NfRelNS{T}, auts::Vector{<: NumFieldMor{NfRelNS{T}, NfRelNS{T}}}) where T

  @vtime :AbExt 2 Ks, mKs = simplified_absolute_field(L, cached = false)
  #Now, we have to construct the maximal order of this field.
  #I am computing the preimages of mKs by hand, by inverting the matrix.
  #Now, the automorphisms.
  autos=Array{NfToNfMor, 1}(undef, length(auts))
  el = mKs(gen(Ks))
  for i = 1:length(auts)
    y = mKs\(auts[i](el))
    autos[i] = hom(Ks, Ks, y, check = false)
  end
  _set_automorphisms_nf(Ks, closure(autos, degree(Ks)))
  @vprint :AbExt 2 "Finished\n"
  return Ks, autos
end

#######################################################################################
#
#  relative discriminant for abelian extension function
#
#######################################################################################

function discriminant_conductor(C::ClassField{MapClassGrp, GrpAbFinGenMap}, bound::fmpz; lwp::Dict{Tuple{Int, Int}, Array{GrpAbFinGenElem, 1}} = Dict{Tuple{Int, Int}, Array{GrpAbFinGenElem, 1}}())
  return true
end

function discriminant_conductor(C::ClassField, bound::fmpz; lwp::Dict{Tuple{Int, Int}, Array{GrpAbFinGenElem, 1}} = Dict{Tuple{Int, Int}, Array{GrpAbFinGenElem, 1}}())

  mr = C.rayclassgroupmap
  O = base_ring(C)
  n = degree(C)
  e = Int(exponent(C))
  lp = mr.fact_mod
  abs_disc = factor(discriminant(O)^n).fac
  if isempty(lp)
    C.absolute_discriminant=abs_disc
    return true
  end
  K = nf(O)
  d = degree(K)
  discr = fmpz(1)
  mp = pseudo_inv(C.quotientmap) * mr
  R = domain(mp)
  a = minimum(defining_modulus(mr)[1])
  primes_done = fmpz[]
  if isprime(n)
    for (p, v) in lp
      if minimum(p, copy = false) in primes_done
        continue
      end
      push!(primes_done, minimum(p, copy = false))
      ap = n*v-v
      qw = divexact(d, p.splitting_type[1])*ap
      mul!(discr, discr, fmpz(p.minimum)^qw)
      if discr > bound
        @vprint :AbExt 2 "too large\n"
        return false
      else
        if haskey(abs_disc, p.minimum)
          abs_disc[p.minimum] += qw
        else
          abs_disc[p.minimum] = qw
        end
      end
    end
    return true
  end

  powers = mr.powers
  groups_and_maps = mr.groups_and_maps

  for i = 1:length(powers)
    p, q = powers[i]
    if p.minimum in primes_done
      continue
    end
    push!(primes_done, p.minimum)
    if p == q
      ap = n
      tmg = groups_and_maps[i][2].tame[p]
      el = C.quotientmap(tmg.disc_log)
      Q, mQ = quo(R, GrpAbFinGenElem[el], false)
      ap -= order(Q)
      qw = divexact(d, p.splitting_type[1])*ap
      mul!(discr, discr, fmpz(p.minimum)^qw)
      if discr > bound
        @vprint :AbExt 2 "too large\n"
        return false
      else
        if haskey(abs_disc, p.minimum)
          abs_disc[p.minimum] += qw
        else
          abs_disc[p.minimum] = qw
        end
      end
      continue
    end
    np = p.minimum^divexact(d, p.splitting_type[1])
    ap = n*lp[p]
    s = lp[p]
    @hassert :AbExt 1 s>=2
    els=GrpAbFinGenElem[]
    for k=2:lp[p]
      s = s-1
      pk = p^s
      pv = pk*p
      if haskey(lwp, (Int(p.minimum), s+1))
        gens = lwp[(Int(p.minimum), s+1)]
      else
        gens_els = _1pluspk_1pluspk1(K, p, pk, pv, powers, a, e)
        gens = Vector{GrpAbFinGenElem}(undef, length(gens_els))
        for i = 1:length(gens)
          gens[i] = mr\(ideal(O, gens_els[i]))
        end
        lwp[(Int(p.minimum), s+1)] = gens
      end
      for i = 1:length(gens)
        push!(els, C.quotientmap(gens[i]))
      end
      o = order(quo(R, els, false)[1])
      ap -= o
      tentative_ap = ap - (lp[p] - k + 1)*o
      tentative_discr = discr * (np^tentative_ap)
      if tentative_discr > bound
        return false
      end
      @hassert :AbExt 1 ap>0
    end
    if haskey(groups_and_maps[i][2].tame, p)
      v = groups_and_maps[i][2].tame[p]
      push!(els, C.quotientmap(v.disc_log))
    end
    ap -= order(quo(R, els, false)[1])
    @hassert :AbExt 1 ap>0
    np1 = np^ap
    mul!(discr, discr, np1)
    if discr > bound
      @vprint :AbExt 2 "too large\n"
      return false
    else
      if haskey(abs_disc, p.minimum)
        abs_disc[p.minimum] += ap*divexact(d, p.splitting_type[1])
      else
        abs_disc[p.minimum] = ap*divexact(d, p.splitting_type[1])
      end
    end
  end
  C.absolute_discriminant = abs_disc
  return true

end

#same function but for ray class groups over QQ

function discriminant_conductorQQ(O::NfOrd, C::ClassField, m::Int, bound::fmpz)

  n = degree(C)
  discr=fmpz(1)
  mp = pseudo_inv(C.quotientmap) * C.rayclassgroupmap
  G=domain(mp)

  cyc_prime= isprime(n)==true

  lp=factor(m).fac
  abs_disc=Dict{fmpz,Int}()

  R=ResidueRing(FlintZZ, m, cached=false)

  for (p,v) in lp
    if v==1
      ap=n
      if cyc_prime
        ap-=1
      else
        x=_unit_grp_residue_field_mod_n(Int(p),n)[1]
        s=divexact(m,Int(p))
        d,a,b=gcdx(s, Int(p))
        l=Int((R(x)*a*s+b*Int(p)).data)
        el=mp\ideal(O,l)
        q,mq=quo(G, GrpAbFinGenElem[el], false)
        ap-= order(q)
      end
      discr*=p^ap
      if discr>bound
        @vprint :AbExt 2 "too large\n"
        return false
      else
        abs_disc[p]=ap
      end
    else
      ap=n*v
      pow=Int(p)^Int(v)
      el = R(1)
      if cyc_prime
        ap-=v
      else
        if isodd(p)
          s=divexact(m,pow)
          d,a,b=gcdx(pow,s)
          s1=R(1+p)^(p-1)
          el=G[1]
          if v==2
            el=mp\ideal(O,Int((b*s*R(s1)+a*pow).data))
            ap-=order(quo(G,GrpAbFinGenElem[el], false)[1])
          else
            for k=0:v-2
              el=mp\ideal(O,Int((b*s*R(s1)^(p^k)+a*pow).data))
              ap-=order(quo(G, GrpAbFinGenElem[el], false)[1])
              @hassert :AbExt 1 ap>0
            end
          end
          if gcd(n,p-1)==1
            ap-=order(quo(G, GrpAbFinGenElem[mp\(ideal(O,fmpz((b*s*R(s1)+a*pow).data)))], false)[1])
          else
            x=_unit_grp_residue_field_mod_n(Int(p),n)[1]
            el1=mp\ideal(O,Int((R(x)*b*s+a*pow).data))
            ap-=order(quo(G, GrpAbFinGenElem[mp\(ideal(O,Int((b*s*R(s1)+a*pow).data))), el1], false)[1])
          end
        else
          s=divexact(m,2^v)
          d,a,b=gcdx(2^v,s)
          el=0*G[1]
          for k=v-3:-1:0
            el=mp\ideal(O,Int((R(5)^(2^k)*b*s+a*2^v).data))
            ap-=order(quo(G, GrpAbFinGenElem[el], false)[1])
          end
          el1=mp\ideal(O,Int((R(-1)*b*s+a*p^v).data))
          ap-=2*order(quo(G, GrpAbFinGenElem[el, el1], false)[1])
        end
      end
      discr*=p^ap
      if discr>bound
        @vprint :AbExt 2 "too large\n"
        return false
      else
        abs_disc[p]=ap
      end
    end
  end
  C.absolute_discriminant=abs_disc
  return true
end

function discriminantQQ(O::NfOrd, C::ClassField, m::Int)

  discr=fmpz(1)
  n = degree(C)
  mp = pseudo_inv(C.quotientmap) * C.rayclassgroupmap
  G = domain(mp)

  cyc_prime= isprime(n)==true

  lp=factor(m).fac
  abs_disc=Dict{fmpz,Int}()

  R=ResidueRing(FlintZZ, m, cached=false)

  for (p,v) in lp
    if v==1
      ap=n
      if cyc_prime
        ap-=1
      else
        x=_unit_grp_residue_field_mod_n(Int(p),n)[1]
        s=divexact(m,Int(p))
        d,a,b=gcdx(s, Int(p))
        l=Int((R(x)*a*s+b*Int(p)).data)
        el=mp\ideal(O,l)
        q,mq=quo(G, GrpAbFinGenElem[el], false)
        ap-= order(q)
      end
      discr*=p^ap
      abs_disc[p]=ap
    else
      ap=n*v
      pow=Int(p)^Int(v)
      el = R(1)
      if cyc_prime
        ap-=v
      else
        if isodd(p)
          s=divexact(m,pow)
          d,a,b=gcdx(pow,s)
          s1=R(1+p)^(p-1)
          el=G[1]
          if v==2
            el=mp\ideal(O,Int((b*s*R(s1)+a*pow).data))
            ap-=order(quo(G,GrpAbFinGenElem[el], false)[1])
          else
            for k=0:v-2
              el=mp\ideal(O,Int((b*s*R(s1)^(p^k)+a*pow).data))
              ap-=order(quo(G, GrpAbFinGenElem[el], false)[1])
              @hassert :AbExt 1 ap>0
            end
          end
          if gcd(n,p-1)==1
            ap-=order(quo(G, GrpAbFinGenElem[mp\(ideal(O,fmpz((b*s*R(s1)+a*pow).data)))], false)[1])
          else
            x=_unit_grp_residue_field_mod_n(Int(p),n)[1]
            el1=mp\ideal(O,Int((R(x)*b*s+a*pow).data))
            ap-=order(quo(G, GrpAbFinGenElem[mp\(ideal(O,Int((b*s*R(s1)+a*pow).data))), el1], false)[1])
          end
        else
          s=divexact(m,2^v)
          d,a,b=gcdx(2^v,s)
          el=0*G[1]
          for k=v-3:-1:0
            el=mp\ideal(O,Int((R(5)^(2^k)*b*s+a*2^v).data))
            ap-=order(quo(G, GrpAbFinGenElem[el], false)[1])
          end
          el1=mp\ideal(O,Int((R(-1)*b*s+a*p^v).data))
          ap-=2*order(quo(G, GrpAbFinGenElem[el, el1], false)[1])
        end
      end
      discr*=p^ap
      abs_disc[p]=ap
    end
  end
  C.absolute_discriminant=abs_disc
  return discr
end

###############################################################################
#
#  conductor function for abelian extension function
#
###############################################################################

#  For this function, we assume the base field to be normal over Q and the conductor of the extension we are considering to be invariant
#  Checks if the defining modulus is the conductor of C

function _is_conductor_min_normal(C::ClassField{MapClassGrp, GrpAbFinGenMap}; lwp::Dict{Int, Vector{GrpAbFinGenElem}} = Dict{Int, Vector{GrpAbFinGenElem}}())
  return true
end

function _is_conductor_min_normal(C::Hecke.ClassField; lwp::Dict{Int, Array{GrpAbFinGenElem, 1}} = Dict{Int, Array{GrpAbFinGenElem, 1}}())
  mr = C.rayclassgroupmap
  lp = mr.fact_mod
  if isempty(lp)
    return true
  end

  a = minimum(defining_modulus(mr)[1])
  mp = pseudo_inv(C.quotientmap) * mr
  R = domain(mp)
  e = Int(exponent(C))
  O = base_ring(C)
  K = nf(O)
  lp = mr.fact_mod
  powers = mr.powers
  groups_and_maps = mr.groups_and_maps
  #first, tame part
  primes_done = fmpz[]
  for i = 1:length(powers)
    p, q = powers[i]
    if p.minimum in primes_done
      continue
    end
    push!(primes_done, p.minimum)
    if p == q
      #The prime is tamely ramified
      v = groups_and_maps[i][2].tame[p]
      el = C.quotientmap(v.disc_log)
      if iszero(el)
        return false
      end
      continue
    end
    if haskey(lwp, Int(p.minimum))
      gens = lwp[Int(p.minimum)]
    else
      k = lp[p]-1
      pk = p^k
      pv = q
      gens_els = _1pluspk_1pluspk1(K, p, pk, pv, powers, a, e)
      gens = Vector{GrpAbFinGenElem}(undef, length(gens_els))
      for i = 1:length(gens)
        gens[i] = mr\(ideal(O, gens_els[i]))
      end
      lwp[Int(p.minimum)] = gens
    end
    iscond = false
    for i in 1:length(gens)
      if !iszero(C.quotientmap(gens[i]))
        iscond = true
        break
      end
    end
    if !iscond
      return false
    end
  end
  return true
end

#
#  Same function as above, but in the assumption that the map comes from a ray class group over QQ
#

function _is_conductor_minQQ(C::Hecke.ClassField, n::Int)

  mr = C.rayclassgroupmap
  mp = pseudo_inv(C.quotientmap) * mr
  m = defining_modulus(mr)[1]
  mm = Int(minimum(m))
  lp = factor(mm)

  O=order(m)
  K=nf(O)

  R=ResidueRing(FlintZZ, mm, cached=false)
  for (p,v) in lp.fac
    if isodd(p)
      if v==1
        x=_unit_grp_residue_field_mod_n(Int(p), n)[1]
        s=divexact(mm,Int(p))
        d,a,b=gcdx(s,Int(p))
        l=a*s*R(x)+p*b
        if iszero(mp\(ideal(O,Int(l.data))))
          return false
        end
      else
        s=divexact(mm,p^v)
        d,a,b=gcdx(s,p^v)
        x=a*s*R(1+p)^((p-1)*p^(v-2))+p^v*b
        if iszero(mp\(ideal(O,Int(x.data))))
          return false
        end
      end
    else
      if v==1
        return false
      end
      if v==2
        s=divexact(mm,4)
        d,a,b=gcdx(s,4)
        l=a*s*R(-1)+4*b
        if iszero(mp\(ideal(O,Int(l.data))))
          return false
        end
      else
        s=divexact(mm,2^v)
        d,a,b=gcdx(s, 2^v)
        l=a*s*R(5)^(2^(v-3))+2^v*b
        if iszero(mp\(ideal(O,Int(l.data))))
          return false
        end
      end
    end
  end
  return true

end

#Returns the cyclic extension of prime degree i with minimal discriminant
function minimal_prime_cyclic_extension(i::Int)
  k = 2
  while !isprime(k*i+1)
    k +=1
  end
  K, a = cyclotomic_field(k*i+1)
  auts = small_generating_set(automorphisms(K))
  auts[1] = auts[1]^i
  g = auts[1]
  el = g(a)
  for i = 2:k
    g *= auts[1]
    el += g(a)
  end
  f = minpoly(el)
  L = number_field(f, check = false, cached = false)[1]
  set_special(L, :isabelian => true)
  return simplify(L, cached = false)[1]
end


text:
using Test, ANSIColoredPrinters

@testset "isnormal" begin
    n = ANSIColoredPrinters.SGRColor()
    @test ANSIColoredPrinters.isnormal(n)

    c16 = ANSIColoredPrinters.SGRColor("30")
    @test !ANSIColoredPrinters.isnormal(c16)

    c256 = ANSIColoredPrinters.SGRColor("38_5", "000")
    @test !ANSIColoredPrinters.isnormal(c256)
end

@testset "is216color" begin
    c000 = ANSIColoredPrinters.SGRColor("38_5", "000")
    @test ANSIColoredPrinters.is216color(c000)

    c080808 = ANSIColoredPrinters.SGRColor("48_5", "080808")
    @test !ANSIColoredPrinters.is216color(c080808)

    cabc = ANSIColoredPrinters.SGRColor("38_2", "abc")
    @test !ANSIColoredPrinters.is216color(cabc)
end

@testset "set_16colors!" begin
    ctx = ANSIColoredPrinters.SGRContext()

    ANSIColoredPrinters.set_16colors!(ctx, "30")
    @test ctx.fg.class == "30"
    @test isempty(ctx.fg.hex)
    @test ANSIColoredPrinters.codes(ctx.fg) === (30,)

    ANSIColoredPrinters.set_16colors!(ctx, "47")
    @test ctx.bg.class == "47"
    @test isempty(ctx.bg.hex)
    @test ANSIColoredPrinters.codes(ctx.bg) === (47,)

    ANSIColoredPrinters.set_16colors!(ctx, "97")
    @test ctx.fg.class == "97"
    @test isempty(ctx.fg.hex)
    @test ANSIColoredPrinters.codes(ctx.fg) === (97,)

    ANSIColoredPrinters.set_16colors!(ctx, "100")
    @test ctx.bg.class == "100"
    @test isempty(ctx.bg.hex)
    @test ANSIColoredPrinters.codes(ctx.bg) === (100,)
end

@testset "set_256colors!" begin
    ctx = ANSIColoredPrinters.SGRContext()

    # 16 colors
    ANSIColoredPrinters.set_256colors!(ctx, "38", "1")
    @test ctx.fg.class == "31"
    @test isempty(ctx.fg.hex)
    @test ANSIColoredPrinters.codes(ctx.fg) === (31,)

    ANSIColoredPrinters.set_256colors!(ctx, "48", "6")
    @test ctx.bg.class == "46"
    @test isempty(ctx.bg.hex)
    @test ANSIColoredPrinters.codes(ctx.bg) === (46,)

    ANSIColoredPrinters.set_256colors!(ctx, "38", "15")
    @test ctx.fg.class == "97"
    @test isempty(ctx.fg.hex)
    @test ANSIColoredPrinters.codes(ctx.fg) === (97,)

    ANSIColoredPrinters.set_256colors!(ctx, "48", "8")
    @test ctx.bg.class == "100"
    @test isempty(ctx.bg.hex)
    @test ANSIColoredPrinters.codes(ctx.bg) === (100,)

    # 216 colors (6 * 6 * 6)
    ANSIColoredPrinters.set_256colors!(ctx, "38", "16")
    @test ctx.fg.class == "38_5"
    @test ctx.fg.hex == "000"
    @test ANSIColoredPrinters.codes(ctx.fg) === (38, 5, 16)

    ANSIColoredPrinters.set_256colors!(ctx, "48", "17")
    @test ctx.bg.class == "48_5"
    @test ctx.bg.hex == "00005f"
    @test ANSIColoredPrinters.codes(ctx.bg) === (48, 5, 17)

    ANSIColoredPrinters.set_256colors!(ctx, "38", "110")
    @test ctx.fg.class == "38_5"
    @test ctx.fg.hex == "87afd7"
    @test ANSIColoredPrinters.codes(ctx.fg) === (38, 5, 110)

    ANSIColoredPrinters.set_256colors!(ctx, "38", "230")
    @test ctx.fg.class == "38_5"
    @test ctx.fg.hex == "ffffd7"
    @test ANSIColoredPrinters.codes(ctx.fg) === (38, 5, 230)

    ANSIColoredPrinters.set_256colors!(ctx, "48", "231")
    @test ctx.bg.class == "48_5"
    @test ctx.bg.hex == "fff"
    @test ANSIColoredPrinters.codes(ctx.bg) === (48, 5, 231)

    # grays
    ANSIColoredPrinters.set_256colors!(ctx, "38", "232")
    @test ctx.fg.class == "38_5"
    @test ctx.fg.hex == "080808"
    @test ANSIColoredPrinters.codes(ctx.fg) === (38, 5, 232)

    ANSIColoredPrinters.set_256colors!(ctx, "48", "255")
    @test ctx.bg.class == "48_5"
    @test ctx.bg.hex == "eee"
    @test ANSIColoredPrinters.codes(ctx.bg) === (48, 5, 255)
end

@testset "set_24bitcolors" begin
    ctx = ANSIColoredPrinters.SGRContext()

    ANSIColoredPrinters.set_24bitcolors!(ctx, "38", "0", "128", "255")
    @test ctx.fg.class == "38_2"
    @test ctx.fg.hex == "0080ff"
    @test ANSIColoredPrinters.codes(ctx.fg) === (38, 2, 0, 128, 255)

    ANSIColoredPrinters.set_24bitcolors!(ctx, "48", "170", "187", "204")
    @test ctx.bg.class == "48_2"
    @test ctx.bg.hex == "abc"
    @test ANSIColoredPrinters.codes(ctx.bg) === (48, 2, 170, 187, 204)

    # 216 colors
    ANSIColoredPrinters.set_24bitcolors!(ctx, "38", "0", "0", "0")
    @test ctx.fg.class == "38_5"
    @test ctx.fg.hex == "000"
    @test ANSIColoredPrinters.codes(ctx.fg) === (38, 5, 16)

    ANSIColoredPrinters.set_24bitcolors!(ctx, "48", "0", "0", "95")
    @test ctx.bg.class == "48_5"
    @test ctx.bg.hex == "00005f"
    @test ANSIColoredPrinters.codes(ctx.bg) === (48, 5, 17)

    # grays
    ANSIColoredPrinters.set_24bitcolors!(ctx, "38", "8", "8", "8")
    @test ctx.fg.class == "38_5"
    @test ctx.fg.hex == "080808"
    @test ANSIColoredPrinters.codes(ctx.fg) === (38, 5, 232)

    ANSIColoredPrinters.set_24bitcolors!(ctx, "48", "238", "238", "238")
    @test ctx.bg.class == "48_5"
    @test ctx.bg.hex == "eee"
    @test ANSIColoredPrinters.codes(ctx.bg) === (48, 5, 255)
end


text:
using HiddenPoisson
using Parameters
using Base.Test

import HiddenPoisson: transitions, observation

"Count the (positive) integers in `itr`, returning a vector of counts."
function count_integers(itr, largest = maximum(itr))
    counts = zeros(Int, largest)
    for i in itr
        counts[i] += 1
    end
    counts
end

function test_competing_poisson(rates, N=10000)
    Tis = [HiddenPoisson.competing_poisson(rates) for _ in 1:N]
    @test isapprox(mean(first.(Tis)), 1/sum(rates); rtol = 0.05)
    @test isapprox(normalize(count_integers((last(Ti) for Ti in Tis), length(rates)), 1),
                   normalize(rates, 1); rtol = 0.1)
end

@testset "competing poisson" begin
    for _ in 1:1000
        N = rand(2:10)
        rates = abs.(randn(10))
        test_competing_poisson(rates)
    end
end

"""
Data generating process for the Diamond-Mortensen-Pissarides model in
continuous time.

# States

Employment (`:E`) and unemployment (`:U`).

# Shocks

employment: job finding (at rate λ)
unemployment: separation (at rate σ)

# Observations

Return :E and :U according to the state (deterministic).
"""
@with_kw immutable DMPmodel{T}
    λ::T                        # job finding rate
    σ::T                        # separation rate
end

transitions(m::DMPmodel, state) = Dict(:E => [m.σ => :U], :U => [m.λ => :E])[state]
observation(m::DMPmodel, state) = state

employment_rate(m::DMPmodel) = model.λ/(model.λ + model.σ)

model = DMPmodel(λ = 0.4, σ = 0.1)

@testset "hidden state steps" begin
    hs1 = HiddenState(model, :E, 2.5, :U)
    o1, hs2 = next_observation(hs1, 1)
    @test o1 == :E
    @test hs1.model ≡ hs2.model
    @test hs1.state ≡ hs2.state
    @test hs2.T == 1.5
    @test hs1.next_state == hs2.next_state
    o2, hs3 = next_observation(hs2, 1)
    @test o2 == :E
    @test hs1.model ≡ hs3.model
    @test hs1.state ≡ hs3.state
    @test hs3.T == 0.5
    @test hs1.next_state == hs3.next_state
end

@testset "hidden state constructor" begin
    @test isapprox(mean(HiddenState(model, :E).T for _ in 1:10000),
                   1/model.σ; rtol = 0.05)
end

@testset "hidden state simulation" begin
    hs = HiddenState(model, :E)
    E = 0
    N = 100000
    for _ in 1:N
        o, hs = next_observation(hs, 1/12)
        if o == :E
            E += 1
        end
    end
    @test isapprox(E/N, employment_rate(model); rtol = 0.03)
end

@testset "hidden state simulation weighted random duration" begin
    hs = HiddenState(model, :E)
    E = 0.0
    T = 0.0
    for _ in 1:10000
        t = abs(randn())
        o, hs = next_observation(hs, t)
        T += t
        if o == :E
            E += t
        end
    end
    @test isapprox(E/T, employment_rate(model); rtol = 0.03)
end

@testset "multiple observations" begin
    hs = HiddenState(model, :E)
    E = 0
    N = 0
    for _ in 1:10000
        ts = abs(randn(rand(1:10)))
        N += length(ts)
        os, hs = next_observations(hs, ts)
        @test length(os) == length(ts)
        for o in os
            if o == :E
                E += 1
            end
        end
    end
    @test isapprox(E/N, employment_rate(model); rtol = 0.03)
end

@testset "show" begin
    hs = HiddenState(nothing, 1, 2, 3)
    @test sprint(show, hs) == "HiddenState with model $(hs.model)\n    state $(hs.state)\n    time to next shock $(hs.T)\n    next state $(hs.next_state)"
end


text:
# This implements the variant when A is dense and m < n

using LinearAlgebra

struct IndGraphFat{T <: RealOrComplex} <: IndGraph
  m::Int
  n::Int
  A::Array{T,2}
  AA::Array{T,2}
  tmp::Array{T,1}
  # tmpx::Array{T,1}
  F::Cholesky{T, Array{T, 2}} # LL factorization
end

function IndGraphFat(A::Array{T,2}) where {T <: RealOrComplex}
  m, n = size(A)
  AA = A * A'

  F = cholesky(AA + I)
  #normrows = vec(sqrt.(sum(abs2.(A), 2)))

  IndGraphFat(m, n, A, AA, Array{T, 1}(undef, m), F)
end

function prox!(
    x::AbstractArray{T},
    y::AbstractArray{T},
    f::IndGraphFat,
    c::AbstractArray{T},
    d::AbstractArray{T},
    gamma=1.0
    ) where {T <: RealOrComplex}

  # y .= f.F \ (f.A * c + f.AA * d)
  mul!(f.tmp, f.A, c)
  mul!(y, f.AA, d)
  y .+= f.tmp
  ldiv!(f.F, y)

  # f.A' * (d - y) + c # note: for complex the complex conjugate is used
  copyto!(f.tmp, d)
  f.tmp .-= y
  mul!(x, adjoint(f.A), f.tmp)
  x .+= c
  return 0.0
end

function (f::IndGraphFat)(x::AbstractArray{T}, y::AbstractArray{T}) where
    {T <: RealOrComplex}
  # the tolerance in the following line should be customizable
  mul!(f.tmp, f.A, x)
  f.tmp .-= y
  if norm(f.tmp, Inf) <= 1e-10
    return 0.0
  end
  return +Inf
end

fun_name(f::IndGraphFat) = "Indicator of an operator graph defined by dense full row rank matrix"
# fun_dom(f::IndGraph) = "AbstractArray{Real,1}, AbstractArray{Complex,1}"
# fun_expr(f::IndGraph) = "x,y ↦ 0 if Ax = y, +∞ otherwise"
# fun_params(f::IndGraph) =
#   string( "A = ", typeof(f.A), " of size ", size(f.A))

function prox_naive(
    f::IndGraphFat,
    c::AbstractArray{T},
    d::AbstractArray{T},
    gamma=1.0
  ) where {T <: RealOrComplex}

  y = f.F \ (f.A * c + f.AA * d)
  return c + f.A' * (d - y), y, 0.0
end

function (f::IndGraphFat)(xy::AbstractVector{T}) where
  {T <: RealOrComplex}
  x, y = splitinput(f, xy)
  return f(x, y)
end

(f::IndGraphFat)(xy::Tuple{AbstractVector{T}, AbstractVector{T}}) where
  {T <: RealOrComplex} = f(xy[1], xy[2])


text:
using Images, ImageFiltering, FileIO

function path_weights(weights)
	ret = copy(weights)
	dim = size(weights)
	for i in 2 : dim[1]
		for j in 1 : dim[2]
			min = ret[i - 1, j]
			if j > 1 && ret[i - 1, j - 1] < min
				min = ret[i - 1, j - 1]
			end
			if j < dim[2] && ret[i - 1, j + 1] < min
				min = ret[i - 1, j + 1]
			end
			ret[i, j] += min
		end
	end
	return ret
end

function back_track(weights)
	dim = size(weights)

	buffer::Array{Int, 1} = ones(Int, dim[1])

	buffer[1] = 1
	for i in 2 : dim[2]
		if weights[dim[1], i] < weights[dim[1], buffer[1]]
			buffer[1] = i
		end
	end

	for i in dim[1] - 1 : -1 : 1
		k = dim[1] - i + 1
		j = buffer[k - 1]
		buffer[k] = j
		if j > 1 && weights[i, j - 1] < weights[i, buffer[k]]
			buffer[k] = j - 1
		end
		if j < dim[2] && weights[i, j + 1] < weights[i, buffer[k]]
			buffer[k] = j + 1
		end
	end

	ret::Array{Tuple{Int, Int}, 1} = Array{Tuple{Int, Int}, 1}(undef, dim[1])

	for i in 1 : dim[1]
		ret[i] = dim[1] - i + 1, buffer[i]
	end

	return ret
end

function gradiant(image)
	kernel = [1 0 -1; 2 0 -2;1 0 -1]
	sobel_x = imfilter(image, centered(kernel))
	grad = imfilter(sobel_x, centered(kernel'))
end

function absolute_gradiant(image)
	height, width = size(image)
	grad = gradiant(image)
	reshape([abs(pixel.r)+abs(pixel.g)+abs(pixel.b) for pixel in grad],height,width)
end

function reduce_width(image::Array{ColorTypes.RGBA{FixedPointNumbers.Normed{UInt8,8}},2}, reduction)
	# image = load(image_path)

	for i in 1:reduction
		h, w = size(image)

		#=
		Det finnes mange måter å gi hver piksel en vekt.
		Her brukes summen av absoluttverdiene avgradientene,
		men det kan endres om man vil.
		=#
		weights = absolute_gradiant(image)
		cumulative_weights = path_weights(weights)

		#=
		Her bruker vi stien vi finner med back_track til å lage en boolsk
		matrise der cellene vi ønsker å fjerne fra bildet har verdien false
		=#
		crop_matrix = ones(Bool,h,w)
		for (i,j) in back_track(cumulative_weights)
			crop_matrix[i,j] = false
		end

		#=
		Et boolsk oppslag i matrisen vår fjerner alle celler der den boolske
		matrisen er false. Fordi resultatet av et slikt oppslag er en 1-dimensjonell
		matrise så bruker vi reshape for å få originalformen tilbake, minus den ene
		pikselen vi har fjernet i bredden.
		=#
		image = reshape(image'[crop_matrix'],w-1,h)'
	end

	image
end


text:
create(c::Client, ::Type{Reaction}, m::Message, emoji::StringOrChar) = create_reaction(c, m.channel_id, m.id, emoji)
function create(c::Client, ::Type{Reaction}, m::Message, e::Emoji)
    emoji = e.id === nothing ? e.name : "$(e.name):$(e.id)"
    return create_reaction(c, m.channel_id, m.id, emoji)
end

retrieve(c::Client, ::Type{Reaction}, m::Message, emoji::StringOrChar) = get_reactions(c, m.channel_id, m.id, emoji)
retrieve(c::Client, ::Type{Reaction}, m::Message, e::Emoji) = get(Reaction, c, m, e.name)

delete(c::Client, ::Type{Reaction}, m::Message) = delete_all_reactions(c, m.channel_id, m.id)
function delete(c::Client, emoji::StringOrChar, m::Message, u::User)
    return if c.state.user !== nothing && u.id == c.state.user.id
        delete_own_reaction(c, m.channel_id, m.id, emoji)
    else
        delete_user_reaction(c, m.channel_id, m.id, emoji, u.id)
    end
end
function delete(c::Client, e::Emoji, m::Message, u::User)
    emoji = e.id === nothing ? e.name : "$(e.name):$(e.id)"
    delete(c, emoji, m, u)
end


text:
using Base.Test
using JointPosteriors

model = ( ProbabilityVector(3), )

struct BinaryClassificationData
  X::Array{Int64,1}
  freq::Array{Int64,1}
  NmX::Array{Int64, 1}
  α_m_m1::Float64
  β_m_m1::Float64
  α_p_m1::Float64
  β_p_m1::Float64
  α_τ_m1::Float64
  β_τ_m1::Float64
end

function log_density(data, p)
    @inbounds log_π = data.α_m_m1 * log(p[2]) + data.β_m_m1 * log(1 - p[2]) + data.α_p_m1*log(p[3]) + data.β_p_m1 * log(1 - p[3]) + data.α_τ_m1 * log(p[1]) + data.β_τ_m1 * log(1 - p[1])

    @inbounds @simd for i ∈ eachindex(data.X)
        log_π += data.freq[i] * log( p[1] * (1 - p[2])^data.X[i] * p[2]^(data.NmX[i]) + (1-p[1]) * p[3]^data.X[i] * (1-p[3])^data.NmX[i] )
    end
    log_π
end

function BinaryClassificationData(X::Array{Int, 1}, freq::Array{Int,1}, n::Int; αm::Real = 1, βm::Real = 1, αp::Real = 1, βp::Real = 1, ατ::Real = 1, βτ::Real = 1)
  BinaryClassificationData(X, freq, n .- X, αm - 1, βm - 1, αp - 1, βp - 1, ατ - 1, βτ - 1)
end
X = [0, 1, 2, 3, 4, 7, 8, 9];
freq = [10, 2, 2, 1, 2, 3, 2, 16];
data = BinaryClassificationData(X, freq, 9, βm = 2, βp = 2);

#mod = Model(model)
#md = LogDensities.ModelDiff(log_density, mod.diff_buffer, data);
#jp = fit(mod, data)

function run_tests(model, data, μ, σ, quants, f, tol = 1e-2)
    builds = (SparseQuadratureGrids.AdaptiveRaw, SparseQuadratureGrids.Adaptive, SparseQuadratureGrids.SmolyakRaw, SparseQuadratureGrids.Smolyak)
    q_rules = (SparseQuadratureGrids.GenzKeister, SparseQuadratureGrids.KronrodPatterson)
    @testset begin
        @testset "Build: $b, Q rule: $q." for b ∈ builds, q ∈ q_rules
            mod = Model(model, b{q})
            fit(mod, data)
            jp = fit(mod, data)
            m_grid = marginal(jp, f)
            m_norm = marginal(jp, f, Normal)
            @test isapprox(m_norm.μ, μ, rtol = tol)
            @test isapprox(m_norm.σ, σ, rtol = tol)
            @test isapprox(quantile(m_grid, .025), quants[1], rtol = tol)
            @test isapprox(quantile(m_grid, .25), quants[2], rtol = tol)
            @test isapprox(quantile(m_grid, .5), quants[3], rtol = tol)
            @test isapprox(quantile(m_grid, .75), quants[4], rtol = tol)
            @test isapprox(quantile(m_grid, .975), quants[5], rtol = tol)
            println("\n")
            println(quantile.(m_norm, [.025,.25,.5,.75,.975]))
            println("\n")
            @test isapprox(quantile(m_norm, .025), quants[1], rtol = tol)
            @test isapprox(quantile(m_norm, .25), quants[2], rtol = tol)
            @test isapprox(quantile(m_norm, .5), quants[3], rtol = tol)
            @test isapprox(quantile(m_norm, .75), quants[4], rtol = tol)
            @test isapprox(quantile(m_norm, .975), quants[5], rtol = tol)
        end
    end
end

quants = [0.391, 0.495, 0.55, 0.602, 0.696]
run_tests(model, data, 0.5504, 0.077, quants, p -> p[1], 10^-1.5)