text:
__precompile__()

module MixedModels

using ArgCheck, BlockArrays, CategoricalArrays, DataFrames, Distributions, GLM,
    LinearAlgebra, NLopt, Random, ProgressMeter, Showoff, SparseArrays, StaticArrays,
    Statistics, StatsBase, StatsModels

using LinearAlgebra: BlasFloat, BlasReal, HermOrSym, PosDefException, copytri!
using NamedArrays: NamedArray, setnames!
using Printf: @printf, @sprintf

using StatsFuns: log2π

import Base: *
import NLopt: Opt

export
       @formula,
       AbstractFactorReTerm,
       AbstractReTerm,
       Bernoulli,
       Binomial,
       Block,
       BlockedSparse,
       Gamma,
       LogitLink,
       LogLink,
       InverseGaussian,
       InverseLink,
       GeneralizedLinearMixedModel,
       LinearMixedModel,
       MatrixTerm,
       MixedModel,
       OptSummary,
       Poisson,
       RaggedArray,
       RepeatedBlockDiagonal,
       ScalarFactorReTerm,
       UniformBlockDiagonal,
       VarCorr,
       VectorFactorReTerm,

       bootstrap,
       bootstrap!,
       coef,
       coeftable,
       cond,
       describeblocks,
       condVar,
       deviance,
       dof,
       fit,
       fit!,
       fitted,
       fixef,      # extract the fixed-effects parameter estimates
       fnames,
       getΛ,
       getθ,
       GHnorm,
       Λ,
       Lambda,
       loglikelihood,
       lowerbd,    # lower bounds on the covariance parameters
       model_response,
       nblocks,
       nobs,
       objective,  # the objective function in fitting a model
       pwrss,      # penalized, weighted residual sum-of-squares
       pirls!,     # use Penalized Iteratively Reweighted Least Squares to obtain conditional modes of random effects
       predict,
       ranef,      # extract the conditional modes of the random effects
       refit!,     # install a response and refit the model
       residuals,
       sdest,      # the estimate of the standard deviation of the per-observation noise
       setθ!,
       simulate!,  # simulate a new response in place
       std,
       updateL!,   # update the lower-triangular, blocked matrix L to a new θ
       varest,     # estimate of the residual variance
       vcov

import Base: ==, *

include("types.jl")
include("gausshermite.jl")
include("modelterms.jl")
include("linalg/cholUnblocked.jl")
include("linalg/rankUpdate.jl")
include("linalg/scaleInflate.jl")
include("linalg/logdet.jl")
include("linalg.jl")
include("pls.jl")
include("simulate.jl")
include("PIRLS.jl")
include("mixedmodel.jl")
include("deprecates.jl")

end # module


text:
using StructJuMP, JuMP
using StructJuMPSolverInterface

include("select_solver.jl")

scen = 2
m = StructuredModel(num_scenarios=scen)
@variable(m, x[1:2])
@NLconstraint(m, x[1] + x[2] == 100)
@NLobjective(m, Min, x[1]^2 + x[2]^2)

for i in 1:scen
    bl = StructuredModel(parent=m)
    @variable(bl, y[1:2])
    @NLconstraint(bl, x[1] + y[1]+y[2] ≥  0)
    @NLconstraint(bl, x[1] + y[1]+y[2]  ≤ 50)
    @NLobjective(bl, Min, y[1]^2 + y[2]^2)
end

structJuMPSolve(m)

getVarValue(m)

ext: jl


text:
module MeasureIR

using Compat: @warn
using Unitful: s, Hz, Time, Frequency, uconvert, NoUnits
using SampledSignals: SampleBuf, samplerate
using DSP: DSP, gaussian, hanning, hilbert
using DSP: FIRFilter, FIRWindow, filt, db2amp, pow2db, xcorr, nextfastfft
using DSP: digitalfilter, Bandpass, Lowpass
using DSP: resample
using FFTW: rfft, irfft
using LinearAlgebra: dot
using Statistics: mean, quantile
using Optim: NewtonTrustRegion, optimize, converged, minimizer
using Printf

export stimulus, analyze, noisefloor, prepadding, snr
export expsweep, golay, mls, rpms, impulse
export pilot, findpilot, pilotsync, getperiod

"""
Abstract supertype for impulse response measurements. Subtypes should define a
constructor that allows the user to specify their system's expected decay time.

Subtypes should also define a `stimulus(sig::MyMeasurementType)` function that
returns a vector suitable for playing directly into the system, as well as a
`analyze(sig::MyMeasurementType, response)` that takes the system's response
and generates an impulse response.
"""
abstract type IRMeasurement end

"""
    prepadding(meas::IRMeasurement)

Get the number of samples of prepadding silence for the given measurement.
"""
function prepadding end

"""
    stimulus(sig)

Generate the stimulus signal for the given measurement process. This is the
signal that should be run through your test system to measure the impulse
response. After recording the response you can use the `analyze` function to
generate the impulse response.
"""
function stimulus end

"""
    analyze(sig::IRMeasurement, response::AbstractArray)

Convert a given stimulus response into an impulse response. Generally you will
first create the test signal with e.g. `golay`, then use `stimulus` to give you
the actual stimulus, then pass the original measurement object and the measures
response to `analyze`.

If `response` is a 2D array it is treated as an NxC multichannel response with
`N` frames and `C` channels.
"""
function analyze end

"""
    noisefloor(sig::IRMeasurement, response::AbstractArray)

Estimate the noise floor in the given signal using the prepadded silence
specified by the measurement. Throws an error if the measurement has no
prepadding.

For a multichannel response returns a row vector with power per channel, or
returns a single number if the input is a Vector.
"""
function noisefloor(sig::IRMeasurement, response::AbstractArray)
    N = prepadding(sig)
    N > 0 || throw(ArgumentError("Can't estimate noise without prepadding"))
    _noisefloor(response, N)
end

_noisefloor(response::AbstractMatrix, N) = sum(response[1:N, :].^2, dims=1) / N
_noisefloor(response::AbstractVector, N) = sum(response[1:N].^2) / N

include("options.jl")
include("util.jl")
include("golay.jl")
include("impulse.jl")
include("expsweep.jl")
include("mls.jl")
include("rpms.jl")
include("pilot.jl")
include("period.jl")

# workaround needed because DSP.jl doesn't handle SampleBufs
# and Float32s well. We dispatch to an internal _analyze to avoid method
# ambiguity issues if each measurement type defined its own `analyze` method.
analyze(sig::IRMeasurement, response::AbstractArray; kwargs...) =
        _analyze(sig, response; kwargs...)
analyze(sig::IRMeasurement, response::SampleBuf; kwargs...) =
        SampleBuf(_analyze(sig, Float64.(response.data); kwargs...), samplerate(response))

end # module


text:
#=

  Bus scheduling in Julia ConstraintSolver.jl

  Problem from Taha "Introduction to Operations Research", page 58.
  Scheduling of buses during a day.

  This is a slightly more general model than Taha's.

  Model created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
const CS = ConstraintSolver
include("constraints_utils.jl")

function bus_schedule(demands,print_solutions=true,all_solutions=true)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   # "all_solutions"=> all_solutions,
                                                            "all_optimal_solutions"=>all_solutions,
                                                            "logging"=>[],

                                                            "traverse_strategy"=>:BFS,
                                                            # "traverse_strategy"=>:DFS,
                                                            # "traverse_strategy"=>:DBFS,

                                                            # "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            "branch_split"=>:InHalf,

                                                            # https://wikunia.github.io/ConstraintSolver.jl/stable/options/#branch_strategy-(:Auto)
                                                            "branch_strategy" => :IMPS, # default
                                                            # "branch_strategy" => :ABS, # Activity Based Search
                                                            # "activity.decay" => 0.999, # default 0.999
                                                            # "activity.max_probes" => 10, # default, 10
                                                            # "activity.max_confidence_deviation" => 20, # default 20

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>6,

                                                            # "backtrack" => false, # default true
                                                            # "backtrack_sorting" => false, # default true

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            # "lp_optimizer" => glpk_optimizer,
                                                            # "lp_optimizer" => ipopt_optimizer,
                                        ))
    timeslots = length(demands)

    # result: how many buses start the schedule at time slot t?
    @variable(model, 0 <= x[1:timeslots] <= sum(demands), Int)
    @variable(model, 0 <= num_buses <= sum(demands), Int)

    # meet the demands for this and the next time slot
    for i in 1:timeslots-1
        @constraint(model, x[i]+x[i+1] >= demands[i])
    end
    # demand "around the clock"
    @constraint(model, x[timeslots]+x[1] >= demands[timeslots])

    @constraint(model, num_buses == sum(x))

    # Symmetry breaking
    my_min(model, x, x[1])

    @objective(model,Min, num_buses)

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                println("solution #$sol")
                num_buses_val = convert.(Integer,JuMP.value.(num_buses; result=sol))
                x_val = convert.(Integer,JuMP.value.(x; result=sol))
                println("x:$x_val num_buses:$num_buses_val")

            end
        end
    else
        println("status:$status")
    end

    return status
end

# demand: minimum number of buses at time t
demands = [8, 10, 7, 12, 4, 4]
@time bus_schedule(demands,true,true)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
module TKTDsimulations

    import DiffEqBase
    import DifferentialEquations

    include("model.jl")

    export runTK, runTK_MCMC
    export runSD, runSD_MCMC
    export runIT, runIT_MCMC

end


text:
# ================
# Homogeneous case
# ================

# X is the universal set => it is ignored
# in this case order reduction only has to be applied to the
# initial set because the loop does not create new generators
# it is assumed that order(Ω0) <= max_order
function reach_homog_GLGM06!(F::Vector{ReachSet{N, Zonotope{N, VN, MN}}},
                             Ω0::Zonotope{N, VN, MN},
                             Φ::AbstractMatrix,
                             NSTEPS::Integer,
                             δ::Float64,
                             max_order::Integer,
                             X::Universe,
                             preallocate::Val{false},
                             Δt0::TimeInterval,
                             disjointness_method::AbstractDisjointnessMethod) where {N, VN, MN}

    # initial reach set
    Δt = (zero(N) .. δ) + Δt0
    @inbounds F[1] = ReachSet(Ω0, Δt)

    k = 2
    @inbounds while k <= NSTEPS
        Rₖ = linear_map(Φ, set(F[k-1]))
        Δt += δ
        F[k] = ReachSet(Rₖ, Δt)
        k += 1
    end
    return F
end

# version that preallocates the output zonotopes
function reach_homog_GLGM06!(F::Vector{ReachSet{N, Zonotope{N, Vector{N}, Matrix{N}}}},
                             Ω0::Zonotope{N, Vector{N}, Matrix{N}},
                             Φ::AbstractMatrix,
                             NSTEPS::Integer,
                             δ::Float64,
                             max_order::Integer,
                             X::Universe,
                             preallocate::Val{true},
                             Δt0::TimeInterval,
                             disjointness_method::AbstractDisjointnessMethod) where {N}

    # initial reach set
    Δt = (zero(N) .. δ) + Δt0
    @inbounds F[1] = ReachSet(Ω0, Δt)

    n, p = size(Ω0.generators)

    # preallocate output
    Zout = Vector{Zonotope{N, Vector{N}, Matrix{N}}}(undef, NSTEPS)

    if p == 0
        Zout[1] = Zonotope(Ω0.center, zeros(N, n, 1))
        p = 1
    else
        Zout[1] = Ω0
    end

    @inbounds for i in 2:NSTEPS
        c = Vector{N}(undef, n)
        G = Matrix{N}(undef, n, p)
        Zout[i] = Zonotope(c, G)
    end

    k = 2
    @inbounds while k <= NSTEPS
        linear_map!(Zout[k], Φ, Zout[k-1])
        Δt += δ
        F[k] = ReachSet(Zout[k], Δt)
        k += 1
    end
    return F
end

# check interection with invariant on the loop
function reach_homog_GLGM06!(F::Vector{ReachSet{N, Zonotope{N, VN, MN}}},
                             Ω0::Zonotope{N, VN, MN},
                             Φ::AbstractMatrix,
                             NSTEPS::Integer,
                             δ::Float64,
                             max_order::Integer,
                             X::LazySet,
                             preallocate::Val{false},
                             Δt0::TimeInterval,
                             disjointness_method::AbstractDisjointnessMethod) where {N, VN, MN}
    # initial reach set
    Δt = (zero(N) .. δ) + Δt0
    @inbounds F[1] = ReachSet(Ω0, Δt)

    k = 2
    while k <= NSTEPS
        Rₖ = linear_map(Φ, set(F[k-1]))
        _is_intersection_empty(X, Rₖ, disjointness_method) && break
        Δt += δ
        F[k] = ReachSet(Rₖ, Δt)
        k += 1
    end
    if k < NSTEPS + 1
        resize!(F, k-1)
    end
    return F
end

# check interection with invariant on the loop, implementation with zonotope preallocation
function reach_homog_GLGM06!(F::Vector{ReachSet{N, Zonotope{N, Vector{N}, Matrix{N}}}},
                             Ω0::Zonotope{N, Vector{N}, Matrix{N}},
                             Φ::AbstractMatrix,
                             NSTEPS::Integer,
                             δ::Float64,
                             max_order::Integer,
                             X::LazySet,
                             preallocate::Val{true},
                             Δt0::TimeInterval,
                             disjointness_method::AbstractDisjointnessMethod) where {N, VN, MN}
    # initial reach set
    Δt = (zero(N) .. δ) + Δt0
    @inbounds F[1] = ReachSet(Ω0, Δt)

    n, p = size(Ω0.generators)

    # preallocate output
    Zout = Vector{Zonotope{N, Vector{N}, Matrix{N}}}(undef, NSTEPS)

    if p == 0
        Zout[1] = Zonotope(Ω0.center, zeros(N, n, 1))
        p = 1
    else
        Zout[1] = Ω0
    end

    @inbounds for i in 2:NSTEPS
        c = Vector{N}(undef, n)
        G = Matrix{N}(undef, n, p)
        Zout[i] = Zonotope(c, G)
    end

    k = 2
    while k <= NSTEPS
        linear_map!(Zout[k], Φ, Zout[k-1])
        _is_intersection_empty(X, Zout[k], disjointness_method) && break
        Δt += δ
        F[k] = ReachSet(Zout[k], Δt)
        k += 1
    end
    if k < NSTEPS + 1
        resize!(F, k-1)
    end
    return F
end

#= O L D
# homogeneous case using StaticArrays
# it is assumed that the order of Ω0 is at most max_order
function reach_homog_GLGM06!(F::Vector{ReachSet{N, Zonotope{N, VN, MN}}},
                             Ω0::Zonotope{N, VN, MN},
                             Φ::SMatrix,
                             NSTEPS::Integer,
                             δ::Float64,
                             max_order::Integer,
                             X::Universe) where {N, VN<:SVector, MN<:SMatrix}

    # initial reach set
    Δt = zero(N) .. δ
    F[1] = ReachSet(Ω0, Δt)

    G0 = genmat(Ω0)
    c0 = center(Ω0)
    n, p = size(G0)
    np = n * p

    # cache for powers of Φ
    Φ_power_k = MMatrix{n, n, N, n * n}(Φ)
    Φ_power_k_cache = similar(Φ)

    k = 1
    @inbounds while k <= NSTEPS - 1
        ck = Φ_power_k * c0
        Gk = SMatrix{n, p, N, np}(Φ_power_k * G0)
        Zk = Zonotope(ck, Gk)
        Δt += δ
        k += 1
        F[k] = ReachSet(Zk, Δt)

        mul!(Φ_power_k_cache, Φ_power_k, Φ)
        copyto!(Φ_power_k, Φ_power_k_cache)
    end
    return F
end
=#

#=
# follows the implementation for the general case
function reach_homog_GLGM06!(F::Vector{ReachSet{N, Zonotope{N, VN, MN}}},
                             Ω0::Zonotope{N, VN, MN},
                             Φ::SMatrix,
                             NSTEPS::Integer,
                             δ::Float64,
                             max_order::Integer,
                             X::Universe) where {N, VN<:SVector, MN<:SMatrix}

    # initial reach set
    Δt = zero(N) .. δ
    F[1] = ReachSet(Ω0, Δt)

    k = 2
    @inbounds while k <= NSTEPS
        # TODO: benchmark wrt other method
        # .. substitute with doing Rₖ = linear_map(Φ, set(F[k-1]))
        X = set(F[k-1])
        ck = Φ * X.center
        Gk = Φ * X.generators

        Δt += δ
        Zk = Zonotope(ck, Gk)
        F[k] = ReachSet(Zk, Δt)
        k += 1
    end
    return F
end
=#

# OLD
#=
# in-place linear map
function reach_homog_GLGM06_inplace!(F::Vector{ReachSet{N, Zonotope{N, VN, MN}}},
                                     Ω0::Zonotope{N, VN, MN},
                                     Φ::AbstractMatrix,
                                     NSTEPS::Integer,
                                     δ::Float64,
                                     max_order::Integer,
                                     X::Universe,
                                     reduction_method::AbstractReductionMethod) where {N, VN, MN}
    # initial reach set
    Δt = zero(N) .. δ
    Ω0red = reduce_order(Ω0, max_order, reduction_method)
    F[1] = ReachSet(Ω0red, Δt)
    c_Ω0 = center(Ω0red)
    G_Ω0 = genmat(Ω0red)
    n, p = size(G_Ω0)

    # cache for powers of Φ
    Φ_power_k = copy(Φ)
    Φ_power_k_cache = similar(Φ)

    # vector of generator matrices
    cvec = Vector{VN}(undef, NSTEPS-1)
    Gvec = Vector{MN}(undef, NSTEPS-1)
    k = 1
    @inbounds while k <= NSTEPS-1
        cvec[k] = Vector{N}(undef, n)
        Gvec[k] = Matrix{N}(undef, n, p)
        mul!(cvec[k], Φ_power_k, c_Ω0)
        mul!(Gvec[k], Φ_power_k, G_Ω0)
        Rₖ = Zonotope(cvec[k], Gvec[k])
        Δt += δ
        mul!(Φ_power_k_cache, Φ_power_k, Φ)
        copyto!(Φ_power_k, Φ_power_k_cache)
        k += 1
        F[k] = ReachSet(Rₖ, Δt)
    end
    return F
end
=#

#=
# in-place linear map
function reach_homog_GLGM06_inplace_v2!(F::Vector{ReachSet{N, Zonotope{N, VN, MN}}},
                                     Ω0::Zonotope{N, VN, MN},
                                     Φ::AbstractMatrix,
                                     NSTEPS::Integer,
                                     δ::Float64,
                                     max_order::Integer,
                                     X::Universe,
                                     reduction_method::AbstractReductionMethod) where {N, VN, MN}
    # initial reach set
    Δt = zero(N) .. δ
    Ω0red = reduce_order(Ω0, max_order, reduction_method) # may not preserve SArrays

    c_Ω0 = center(Ω0red)
    G_Ω0 = genmat(Ω0red)
    n, p = size(G_Ω0)

    c_Ω0 = MVector{n}(c_Ω0)
    G_Ω0 = MMatrix{n, p}(G_Ω0)

    F[1] = ReachSet(Ω0red, Δt)

    k = 1
    @inbounds while k <= NSTEPS - 1
        #c_Rₖ = Φ * F[k].X.center
        mul!(F[k+1].X.center, Φ, F[k].X.center)
        mul!(F[k+1].X.generators, Φ, F[k].X.generators)
        #F[k+1].Δt
        #G_Rₖ = Φ * F[k].X.generators
        #Rₖ = Zonotope(c_Rₖ, G_Rₖ)
        Δt += δ
        k += 1
        F[k] = ReachSet(Rₖ, Δt)
    end
    return F
end
=#


text:
using PyPlot

"""
`histplot(x)` plot a histogram of the values in `x` and
`histplot(x,n)` gives a plot with `n` bins.
"""
histplot = plt[:hist]


text:
# ---------------------------------------------------------------------------- #
#
#   pgauss.jl
#
#   Abstract type for polynomial order
#
#   λυτέος
#   Fall 2017
#
#   Max Opgenoord
#
# ---------------------------------------------------------------------------- #

"""
    PGauss

PGauss abstract type:
Overarching abstract type for polynomial order types.
"""
abstract type PGauss

end

### PG = 1
struct PG1 <: PGauss
  p::Int64
end
function PG1(  )
  PG1(1)
end
function comporder( P::PG1 )
  return ( PG1(), PG1(), PG1() )
end

### PG = 2
struct PG2 <: PGauss
  p::Int64
end
function PG2(  )
  PG2(2)
end
function comporder( P::PG2 )
  return ( PG3(), PG1(), PG1() )
end

### PG = 3
struct PG3 <: PGauss
  p::Int64
end
function PG3(  )
  PG3(3)
end
function comporder( P::PG3 )
  return ( PG3(), PG2(), PG2() )
end

### PG = 4
struct PG4 <: PGauss
  p::Int64
end
function PG4(  )
  PG4(4)
end
function comporder( P::PG4 )
  return ( PG5(), PG3(), PG3() )
end

### PG = 5
struct PG5 <: PGauss
  p::Int64
end
function PG5(  )
  PG5(5)
end
function comporder( P::PG5 )
  return ( PG5(), PG4(), PG3() )
end

### PG = 6
struct PG6 <: PGauss
  p::Int64
end
function PG6(  )
  PG6(6)
end
function comporder( P::PG6 )
  return ( PG7(), PG5(), PG4() )
end

### PG = 7
struct PG7 <: PGauss
  p::Int64
end
function PG7(  )
  PG7(7)
end
function comporder( P::PG7 )
  return ( PG7(), PG5(), PG5() )
end

### PG = 8
struct PG8 <: PGauss
  p::Int64
end
function PG8(  )
  PG8(8)
end
function comporder( P::PG8 )
  return ( PG8(), PG5(), PG6() )
end

### PG = 9
struct PG9 <: PGauss
  p::Int64
end
function PG9(  )
  PG9(9)
end
function comporder( P::PG9 )
  return ( PG9(), PG8(), PG7() )
end

### PG = 10
struct PG10 <: PGauss
  p::Int64
end
function PG10(  )
  PG10(10)
end
function comporder( P::PG10 )
  return ( PG11(), PG8(), PG8() )
end

### PG = 11
struct PG11 <: PGauss
  p::Int64
end
function PG11(  )
  PG11(11)
end
function comporder( P::PG11 )
  return ( PG11(), PG8(), PG9() )
end

### PG = 12
struct PG12 <: PGauss
  p::Int64
end
function PG12(  )
  PG12(12)
end
function comporder( P::PG12 )
  return ( PG12(), PG8(), PG10() )
end

### PG = 13
struct PG13 <: PGauss
  p::Int64
end
function PG13(  )
  PG13(13)
end
function comporder( P::PG13 )
  return ( PG13(), PG9(), PG11() )
end


text:
export HalfSpace

struct HalfSpace{K} <: AbstractSimpleShape{K}
    x₀::SFloat{K}  # point on boundary
    n̂::SFloat{K}  # outward direction normal
    HalfSpace{K}(x₀,n̂) where {K} = new(x₀,n̂)  # suppress default outer constructor
end

HalfSpace(x₀::AbsVecReal, n::AbsVecReal) = HalfSpace{length(x₀)}(x₀, normalize(n))

Base.:(==)(s1::HalfSpace, s2::HalfSpace) = s1.x₀==s2.x₀ && s1.n̂==s2.n̂
Base.isapprox(s1::HalfSpace, s2::HalfSpace) = s1.x₀≈s2.x₀ && s1.n̂≈s2.n̂
Base.hash(s::HalfSpace, h::UInt) = hash(s.x₀, hash(s.n̂, hash(:HalfSpace, h)))

level(x::SReal{K}, s::HalfSpace{K}) where {K} = s.n̂ ⋅ (x-s.x₀)

# The following functions cannot be implemented for HalfSpace.
# center()
# bounds()

max_stars_repo_name: wsshin/LevelSetShapes

max_issues_repo_name: wsshin/LevelSetShapes

max_forks_repo_name: wsshin/LevelSetShapes


text:
# Julia translation of http://nbviewer.jupyter.org/github/barbagroup/AeroPython/blob/master/lessons/04_Lesson04_vortex.ipynb
# Lesson 3 doublet

using PyPlot
using Distributions

close("all")
meshgrid(x,y) = (repmat(x',length(y),1),repmat(y,1,length(x)))

N = 50                                # number of points in each direction
x_start, x_end = -2.0, 2.0            # boundaries in the x-direction
y_start, y_end = -1.0, 1.0            # boundaries in the y-direction
x = linspace(x_start, x_end, N)    # creates a 1D-array with the x-coordinates
y = linspace(y_start, y_end, N)    # creates a 1D-array with the y-coordinates

X,Y=meshgrid(x,y)

gamma = 5.0

x_vortex, y_vortex = 0.0, 0.0   # location of the source

function get_velocity_vortex(strength, xv, yv, X, Y):
      """Returns the velocity field generated by a vortex.

      Arguments
      ---------
      strength -- strength of the vortex.
      xv, yv -- coordinates of the vortex.
      X, Y -- mesh grid.
      """
      u = + strength/(2*pi).*(Y-yv)./((X-xv).^2+(Y-yv).^2)
      v = - strength/(2*pi).*(X-xv)./((X-xv).^2+(Y-yv).^2)

      return u, v
    end

function get_stream_function_vortex(strength, xv, yv, X, Y):
    """Returns the stream-function generated by a vortex.

    Arguments
    ---------
    strength -- strength of the vortex.
    xv, yv -- coordinates of the vortex.
    X, Y -- mesh grid.
    """
    psi = strength/(4*pi).*log((X-xv).^2+(Y-yv).^2)

    return psi
  end

  # computes the velocity field on the mesh grid
  u_vortex, v_vortex = get_velocity_vortex(gamma, x_vortex, y_vortex, X, Y)

  # computes the stream-function on the mesh grid
  psi_vortex = get_stream_function_vortex(gamma, x_vortex, y_vortex, X, Y)

  # plots the streamlines
  size = 10
  PyPlot.figure(figsize=(size, (y_end-y_start)/(x_end-x_start)*size))
  PyPlot.grid(true)
  PyPlot.xlabel("x", fontsize=16)
  PyPlot.ylabel("y", fontsize=16)
  PyPlot.xlim(x_start, x_end)
  PyPlot.ylim(y_start, y_end)
  PyPlot.streamplot(X, Y, u_vortex, v_vortex, density=2, linewidth=1, arrowsize=1, arrowstyle="->")
  PyPlot.scatter(x_vortex, y_vortex, color="#CD2305", s=80, marker="o")

strength_sink = -1.0            # strength of the sink
x_sink, y_sink = 0.0, 0.0       # location of the sink

function get_velocity_sink(strength, xs, ys, X, Y):
    """Returns the velocity field generated by a sink.

    Arguments
    ---------
    strength -- strength of the sink.
    xs, ys -- coordinates of the sink.
    X, Y -- mesh grid.
    """
    u = strength/(2*pi).*(X-xs)./((X-xs).^2+(Y-ys).^2)
    v = strength/(2*pi).*(Y-ys)./((X-xs).^2+(Y-ys).^2)

    return u, v
  end

function get_stream_function_sink(strength, xs, ys, X, Y):
      """Returns the stream-function generated by a sink.

      Arguments
      ---------
      strength -- strength of the sink.
      xs, ys -- coordinates of the sink.
      X, Y -- mesh grid.
      """
      psi = strength/(2*pi)*atan2((Y-ys), (X-xs))

      return psi
    end

    # computes the velocity field on the mesh grid
u_sink, v_sink = get_velocity_sink(strength_sink, x_sink, y_sink, X, Y)

# computes the stream-function on the mesh grid
psi_sink = get_stream_function_sink(strength_sink, x_sink, y_sink, X, Y)

# superposition of the sink and the vortex
u = u_vortex + u_sink
v = v_vortex + v_sink
psi = psi_vortex + psi_sink

# plots the streamlines
size = 10
PyPlot.figure(figsize=(size, (y_end-y_start)/(x_end-x_start)*size))
PyPlot.xlabel("x", fontsize=16)
PyPlot.ylabel("y", fontsize=16)
PyPlot.xlim(x_start, x_end)
PyPlot.ylim(y_start, y_end)
PyPlot.streamplot(X, Y, u, v, density=2, linewidth=1, arrowsize=1, arrowstyle="->")
PyPlot.scatter(x_vortex, y_vortex, color="#CD2305", s=80, marker="o");


text:
using Documenter, NoncommutativeGraphs

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

makedocs(
    sitename="NoncommutativeGraphs.jl",
    modules=[NoncommutativeGraphs],
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"
    ),
    pages = [
        "Home" => "index.md",
        "Reference" => "reference.md",
    ],
)

deploydocs(
    repo = "github.com/dstahlke/NoncommutativeGraphs.jl.git",
    devbranch = "main",
    branch = "gh-pages",
)

max_stars_repo_name: dstahlke/NoncommutativeGraphs

max_issues_repo_name: dstahlke/NoncommutativeGraphs

max_forks_repo_name: dstahlke/NoncommutativeGraphs


text:
type WEdge
  direction::Symbol
  weight::Union{Number,Void}
  WEdge(direction::Symbol) = new(direction, nothing)
  WEdge(direction::Symbol, weight) = new(direction, weight)
end

function Base.show(io::IO, e::WEdge)
  print(io, string(isa(e.weight, Void) ? "" : e.weight, e.direction))
end

# Undirected
adjacency(e::Symbol, w::Number, b::Symbol) = adjacency(e, w, string(b))
adjacency(e::Symbol, w::Number, b::Label) = [WEdge(e, w); b]
adjacency(e::Symbol, a::Label, w::Number) = [a; WEdge(e, w)]
adjacency(e::Symbol, a::Label, b::Label) = [a; WEdge(e); b]
adjacency(e::Symbol, ex::Expr, b::Label) = [adjacency(ex.args...); WEdge(e); b]
adjacency(e::Symbol, ex::Expr, b::Symbol) = adjacency(e, ex, string(b))
adjacency(e::Symbol, a::Symbol, b::Symbol) = [string(a); WEdge(e); string(b)]

# Directed
function adjacency(w::Number, ex::Expr)
  if ex.head == :block
    for el in ex.args
      if isa(el, Expr)
        if el.head == :->
          return [WEdge(:->, w); adjacency(el.args...)]
        end
      elseif isa(el, Symbol)
        return [WEdge(:->, w); string(el)]
      else
        return [WEdge(:->, w); el]
      end
    end
  end
end

adjacency(sym::Symbol, ex::Expr) = adjacency(string(sym), ex)
function adjacency(a::Label, ex::Expr)
  if ex.head == :block
    for el in ex.args
      if isa(el, Expr)
        if el.head == :->
          return [a; WEdge(:->, nothing); adjacency(el.args...)]
        end
      elseif isa(el, Symbol)
        return [a; WEdge(:->, nothing); string(el)]
      else
        return [a; WEdge(:->, nothing); el]
      end
    end
  end
end

function parse(list::AbstractArray)
  range = 1:2:length(list)
  is_directed = range.stop>1 && list[2].direction==:->
  g = inclist(ExVertex, ExEdge{ExVertex}, is_directed=is_directed)
  dv = Dict()
  cnt = 0
  if range.stop==1
    v = first(list)
    vtx = ExVertex(1, v)
    vtx.attributes["label"] = v
    add_vertex!(g, vtx)
  else
    for idx in range[1:end-1]
      from,edge,to = list[idx:idx+2]
      for v in [from to]
        if !haskey(dv, v)
          cnt += 1
          vtx = ExVertex(cnt, v)
          vtx.attributes["label"] = v
          dv[v] = add_vertex!(g, vtx)
        end
      end
      attrs = AttributeDict()
      if !isa(edge.weight, Void)
        attrs["label"] = edge.weight
        attrs["weight"] = edge.weight
      end
      e = ExEdge{ExVertex}(num_edges(g), dv[from], dv[to], attrs)
      add_edge!(g, e)
    end
  end
  GraphStore(g)
end


text:
################################################################################
#
#  Contained in the alternating group
#
################################################################################

function is_discriminant_square(IdG::Main.ForeignGAP.MPtr)
  G = GAP.Globals.SmallGroup(IdG)
  mp = GAP.Globals.RegularActionHomomorphism(G)
  S = GAP.Globals.ImagesSource(mp)
  lg = GAP.Globals.GeneratorsOfGroup(S)
  for i = 1:length(lg)
    s = GAP.Globals.SignPerm(lg[i])
    if !isone(s)
      return false
    end
  end
  return true
end

################################################################################
#
#  Ramification at infinite places
#
################################################################################

function _real_level(L::Main.ForeignGAP.MPtr)
  #I compute the index of the first derived subgroup containing all the elements of order 2
  #First, I compute the subgroup generated by them.
  lElem = Main.ForeignGAP.MPtr[]
  G = L[1]
  ElemsG = GAP.Globals.Elements(G)
  for i = 1:length(ElemsG)
    g = ElemsG[i]
    if GAP.Globals.Order(g) == 2
      push!(lElem, g)
    end
  end
  S = GAP.Globals.Subgroup(G, GAP.julia_to_gap(lElem))
  #Now, I check containment.
  k = 0
  for i = 2:length(L)
    k += 1
    if !GAP.Globals.IsSubgroup(L[i], S)
      break
    end
  end
  return k
end

###############################################################################
#
#  Conductors
#
###############################################################################

function conductors_after_Brauer(F::FieldsTower, st::Vector{Int}, l_cond::Vector)
  lp = ramified_primes(F)
  autsH = F.auts_for_conductors
  auts = automorphisms(F.field, copy = false)
  imgH = F.imgs_autos
  proj = F.proj_ext
  n = prod(st)
  Hn = GAP.Globals.ImagesSource(proj)
  G = GAP.Globals.Source(proj)
  O = maximal_order(F)
  for p in lp
    lP = prime_decomposition(O, p)
    F, mF = Hecke.ResidueFieldSmall(O, lP[1][1])
    Gp = decomposition_group(auts, lP[1][1])
    Hp = inertia_subgroup(F, mF, Gp)
    gHp = small_generating_set(Hp)
    gensind = Int[]
    for x in gHp
      i = 1
      while x != auts[i]
        i += 1
      end
      push!(gensind, i)
    end
    els = Main.ForeignGAP.MPtr[imgH[i] for i in gensind]
    sub = GAP.Globals.Subgroup(Hn, GAP.julia_to_gap(els))
    ord = GAP.Globals.Size(sub)
    sizes_preimages = Int[]
    for aut in autsH
      subgs = Vector{Vector{Main.ForeignGAP.MPtr}}(undef, n)
      for i = 1:n
        subgs[i] = Vector{Main.ForeignGAP.MPtr}(undef, length(els))
      end
      for j = 1:length(els)
        el2 = GAP.GAP.Globals.Images(aut, els[j])
        pels = GAP.GAP.Globals.List(GAP.GAP.Globals.PreImages(proj, el2))
        for i = 1:length(pels)
          subgs[i][j] = pels[i]
        end
      end
      for lelem in subgs
        sub = GAP.Globals.Subgroup(G, GAP.julia_to_gap(lelem))
        onew = GAP.Globals.Size(sub)
        push!(sizes_preimages, onew)
      end
    end
    if minimum(sizes_preimages) != ord
      #The prime must ramify!
      l1 = Vector{Tuple{Int, Dict{NfOrdIdl, Int}}}()
      if !divisible(fmpz(n), p)
        for x in l_cond
          if divisible(fmpz(x[1]), p)
            push!(l1, x)
          end
        end
      else
        for x in l_cond
          if !isempty(x[2])
            lI = keys(x[2])
            found = false
            for k in lI
              if minimum(k) == p
                found = true
                break
              end
            end
            if found
              push!(l1, x)
            end
          end
        end
      end
      l_cond = l1
    elseif maximum(sizes_preimages) == ord && !divisible(fmpz(n), p)
      #The prime must be unramified!
      l1 = Vector{Tuple{Int, Dict{NfOrdIdl, Int}}}()
      if !divisible(fmpz(n), p)
        for x in l_cond
          if !divisible(fmpz(x[1]), p)
            push!(l1, x)
          end
        end
      else
        for x in l_cond
          if !isempty(x[2])
            lI = keys(x[2])
            found = false
            for k in lI
              if minimum(k) == p
                found = true
                break
              end
            end
            if !found
              push!(l1, x)
            end
          end
        end
      end
      l_cond = l1
    end
  end
  return l_cond

end

function conductor_general_case(F::FieldsTower, st::Vector{Int}, IdG::Main.ForeignGAP.MPtr, bound::fmpz, l_cond::Vector)
  n = prod(st)
  O = maximal_order(F)
  auts = automorphisms(nf(O), copy = false)
  Hperm = _from_autos_to_perm(auts)
  H = _perm_to_gap_grp(Hperm)
  G = GAP.Globals.SmallGroup(IdG)
  D = GAP.Globals.DerivedSeries(G)
  proj = GAP.Globals.NaturalHomomorphismByNormalSubgroup(G, D[end-1])
  Hn = GAP.Globals.ImagesSource(proj)
  iso = GAP.Globals.IsomorphismGroups(Hn, H)
  AutHn = GAP.Globals.AutomorphismGroup(Hn)
  autH = GAP.Globals.Elements(AutHn)

  lp = Hecke.ramified_primes(O)
  for p in lp
    lP = prime_decomposition(O, p)
    F, mF = Hecke.ResidueFieldSmall(O, lP[1][1])
    Gp = decomposition_group(auts, lP[1][1])
    Hp = inertia_subgroup(F, mF, Gp)
    gHp = small_generating_set(Hp)
    gensind = Int[]
    for x in gHp
      i = 1
      while x != auts[i]
        i += 1
      end
      push!(gensind, i)
    end
    els = Main.ForeignGAP.MPtr[GAP.Globals.PreImagesRepresentative(iso, _perm_to_gap_perm(Hperm[i])) for i in gensind]
    sub = GAP.Globals.Subgroup(Hn, GAP.julia_to_gap(els))
    ord = GAP.Globals.Size(sub)
    sizes_preimages = Int[]
    for s = 1:length(autH)
      aut = autH[s]
      subgs = Vector{Vector{Main.ForeignGAP.MPtr}}(undef, n)
      for i = 1:n
        subgs[i] = Vector{Main.ForeignGAP.MPtr}(undef, length(els))
      end
      for j = 1:length(els)
        el2 = GAP.Globals.Images(aut, els[j])
        pels = GAP.Globals.List(GAP.Globals.PreImages(proj, el2))
        for i = 1:length(pels)
          subgs[i][j] = pels[i]
        end
      end
      for lelem in subgs
        sub = GAP.Globals.Subgroup(G, GAP.julia_to_gap(lelem))
        onew = GAP.Globals.Size(sub)
        push!(sizes_preimages, onew)
      end
    end
    if minimum(sizes_preimages) != ord
      #The prime must ramify!
      l1 = Vector{Tuple{Int, Dict{NfOrdIdl, Int}}}()
      if !divisible(fmpz(n), p)
        for x in l_cond
          if divisible(fmpz(x[1]), p)
            push!(l1, x)
          end
        end
      else
        for x in l_cond
          if !isempty(x[2])
            lI = keys(x[2])
            found = false
            for k in lI
              if minimum(k) == p
                found = true
                break
              end
            end
            if found
              push!(l1, x)
            end
          end
        end
      end
      l_cond = l1
    elseif maximum(sizes_preimages) == ord && !divisible(fmpz(n), p)
      #The prime must be unramified!
      l1 = Vector{Tuple{Int, Dict{NfOrdIdl, Int}}}()
      if !divisible(fmpz(n), p)
        for x in l_cond
          if !divisible(fmpz(x[1]), p)
            push!(l1, x)
          end
        end
      else
        for x in l_cond
          if !isempty(x[2])
            lI = keys(x[2])
            found = false
            for k in lI
              if minimum(k) == p
                found = true
                break
              end
            end
            if !found
              push!(l1, x)
            end
          end
        end
      end
      l_cond = l1
    end
  end
  return l_cond

end

function conductors_with_restrictions(F::FieldsTower, st::Vector{Int}, IdG::Main.ForeignGAP.MPtr, bound::fmpz)

  O = maximal_order(F)
  l_cond = Hecke.conductors(O, st, bound)
  if isdefined(F, :auts_for_conductors)
    new_conds = conductors_after_Brauer(F, st, l_cond)
  else
    new_conds = conductor_general_case(F, st, IdG, bound, l_cond)
  end
  if length(st) != 1 || !isprime(st[1]) || isempty(new_conds)
    return new_conds
  end
  #If the extension is cyclic, I take care of the discriminant being a square or not for the wild ramification
  issquare = is_discriminant_square(IdG)
  p = st[1]
  v = valuation(discriminant(O), p)
  is_square_disc_base_field = iszero(mod(v*p, 2))
  td = prime_decomposition_type(O, p)
  if iszero(mod(length(td) * td[1][1]*(p-1), 2))
    #Regardless of the exponents, the norm of the discriminant will be a square
    if issquare && is_square_disc_base_field
      return new_conds
    elseif issquare
      return typeof(new_conds)()
    else
      return new_conds
    end
  end
  #Now, p must be 2.
  if issquare && is_square_disc_base_field
    #Only the even exponents are allowed!
    newer_conds = typeof(new_conds)()
    for i = 1:length(new_conds)
      if isempty(new_conds[i][2])
        push!(newer_conds, new_conds[i])
        continue
      end
      if iszero(mod(first(values(new_conds[i][2])), 2))
        push!(newer_conds, new_conds[i])
      end
    end
  elseif issquare
    #Only the odd exponents are allowed!
    newer_conds = typeof(new_conds)()
    for i = 1:length(new_conds)
      if !isempty(new_conds[i][2]) && !iszero(mod(first(values(new_conds[i][2])), 2))
        push!(newer_conds, new_conds[i])
      end
    end
  else
    newer_conds = new_conds
  end

  #Now, tame ramification.
  list_tame = Int[x[1] for x in newer_conds]
  list_tame = coprime_base(list_tame)
  l = length(list_tame)
  for i = 1:length(list_tame)
    x = list_tame[i]
    if !isone(x) && !isprime(x)
      append!(list_tame, Hecke.divisors(x))
    end
  end
  list_tame = coprime_base(list_tame)
  for q in list_tame
    q == 1 && continue
    @assert isprime(q)
    v = valuation(discriminant(O), q)
    is_square_disc_base_field = iszero(mod(v*p, 2))
    td = prime_decomposition_type(O, q)
    if iszero(mod(length(td) * td[1][1] * (p-1), 2))
      #Regardless of the exponents, the norm of the discriminant will be a square
      if issquare && is_square_disc_base_field
        continue
      elseif issquare || is_square_disc_base_field
        return typeof(new_conds)()
      else
        continue
      end
    end
    #Now, p must be 2.
    if issquare && is_square_disc_base_field
      #Only the even exponents are allowed!
      #Therefore the prime can't ramify
      newest_conds = typeof(new_conds)()
      for i = 1:length(newer_conds)
        if !iszero(mod(newer_conds[i][1], q))
          push!(newer_conds, newer_conds[i])
        end
      end
    elseif issquare
      #Only the odd exponents are allowed!
      #Therefore the prime must ramify
      newest_conds = typeof(new_conds)()
      for i = 1:length(newer_conds)
        if iszero(mod(newer_conds[i][1], q))
          push!(newest_conds, newer_conds[i])
        end
      end
    else
      newest_conds = newer_conds
    end
    newer_conds = newest_conds
  end
  return newer_conds
end

###############################################################################
#
#  Abelian extensions of QQ
#
###############################################################################

function abelian_extensionsQQ(gtype::Array{Int, 1}, bound::fmpz, only_real::Bool = false)

  gtype = map(Int, snf(DiagonalGroup(gtype))[1].snf)
  #Easy case: quadratic and biquadratic extensions
  if gtype == [2]
    lq = Hecke._quad_ext(Int(bound), only_real)
    res = Vector{FieldsTower}(undef, length(lq))
    for i = 1:length(lq)
      x = lq[i]
      E = EquationOrder(x[1])
      E.ismaximal = 1
      E.index = fmpz(1)
      E.gen_index = fmpq(1)
      E.disc = discriminant(E)
      Hecke._set_maximal_order(x[1], E)
      auts = Vector{NfToNfMor}(undef, 2)
      auts[1] = NfToNfMor(x[1], x[1], gen(x[1]))
      auts[2] = x[2][1]
      Hecke._set_automorphisms_nf(x[1], auts)
      res[i] = FieldsTower(x[1], x[2], x[3])
    end
    return res
  end
  if gtype == [2,2]
    l = Hecke._C22_exts_abexts(Int(bound), only_real)
    @vprint :Fields 1 "Computing maximal orders\n"
    @vprint :FieldsNonFancy 1 "Computing maximal orders\n"
    lf = Hecke._C22_with_max_ord(l)
    res1 = Vector{FieldsTower}(undef, length(lf))
    for i = 1:length(lf)
      res1[i] = FieldsTower(lf[i][1], lf[i][2], lf[i][3])
    end
    return res1
  end
  l1 = _abelian_extensionsQQ(gtype, bound, only_real)
  list1 = Vector{FieldsTower}(undef, length(l1))
  @vprint :Fields 1 "Computing maximal orders\n\n"
  @vprint :FieldsNonFancy 1 "Computing maximal orders\n"
  for i = 1:length(l1)
    @vprint :Fields 1 "\e[1FComputing maximal order $(i) /$(length(l1)) \n"
    x = l1[i]
    K, auts = Hecke._from_relative_to_absQQ(x[1], x[2])
    if length(gtype) == 1
      #If the group is cyclic, I prefer to have a generator!
      new_auts = Vector{NfToNfMor}(undef, 1)
      new_auts[1] = auts[1]
      for i = 2:length(auts)
        new_auts[1] *= auts[i]
      end
      list1[i] = FieldsTower(K, new_auts, [Hecke.NfToNfMor(base_field(x[1]), K, K(1))])
    else
      list1[i] = FieldsTower(K, auts, [Hecke.NfToNfMor(base_field(x[1]), K, K(1))])
    end
  end
  @vprint :Fields 1 "\e[1F$(Hecke.set_cursor_col())$(Hecke.clear_to_eol())"
  return list1

end

function _abelian_extensionsQQ(gtype::Array{Int,1}, absolute_discriminant_bound::fmpz, only_real::Bool = false)

  Qx, x = PolynomialRing(FlintQQ, "x", cached = false)
  K, _ = NumberField(x-1, "a", cached = false)
  O = maximal_order(K)
  n = prod(gtype)
  expo = lcm(gtype)

  #Getting conductors
  l_conductors = Hecke.conductorsQQ(O, gtype, absolute_discriminant_bound)
  sort!(l_conductors, rev = true)
  len = length(l_conductors)
  @vprint :Fields 1 "Number of conductors: $(len) \n\n"
  @vprint :FieldsNonFancy 1 "Number of conductors: $(len) \n"

  complex = iseven(expo) && !only_real

  #Now, the big loop
  class_fields = Vector{Hecke.ClassField}()
  for (i, k) in enumerate(l_conductors)
    if iszero(mod(i, 1000))
      pt = len - i
      @vprint :Fields 1 "\e[1F$(Hecke.set_cursor_col())$(Hecke.clear_to_eol()) Conductors to test: $(pt)\n"
    end
    r, mr = Hecke.ray_class_groupQQ(O, k, complex, expo)
    if !Hecke._are_there_subs(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, n)
        push!(class_fields, C)
      end
    end
  end
  fields = Vector{Tuple{Hecke.NfRel_ns{nf_elem}, Array{Hecke.NfRel_nsToNfRel_nsMor{nf_elem},1}}}(undef, length(class_fields))
  for i = 1:length(class_fields)
    @vprint :Fields 1 "\e[1FComputing class field $(i) /$(length(class_fields)) \n"
    C = class_fields[i]
    fields[i] = (number_field(C), Hecke.automorphism_groupQQ(C))
  end
  @vprint :Fields 1 "\e[1F$(Hecke.set_cursor_col())$(Hecke.clear_to_eol())Number of fields found: $(length(fields)) \n"
  @vprint :FieldsNonFancy 1 "Number of fields found: $(length(fields)) \n"
  return fields

end

###############################################################################
#
#  General abelian layer
#
###############################################################################

function check_ext(C::Hecke.ClassField, bound::fmpz, Dcond::Dict, Ddisc::Dict)

  @vtime :Fields 3 fl2 = Hecke._is_conductor_min_normal(C, lwp = Dcond)
  if !fl2
    return false
  end
  @vtime :Fields 3 fl3 = Hecke.discriminant_conductor(C, bound, lwp = Ddisc)
  if !fl3
    return false
  end
  return true

end

function _construct_grp(IdH::Main.ForeignGAP.MPtr, uncom::Int)
  G = GAP.Globals.SmallGroup(IdH)
  L = GAP.Globals.DerivedSeries(G)
  mH1 = GAP.Globals.NaturalHomomorphismByNormalSubgroup(G, L[end-1])
  K = GAP.Globals.Kernel(mH1)
  gens = GAP.gap_to_julia(GAP.Globals.MinimalGeneratingSet(K))
  for i = 1:length(gens)
    o = GAP.Globals.Order(gens[i])
    ex = ppio(o, uncom)[1]
    gens[i] = gens[i]^ex
  end
  S = GAP.Globals.Subgroup(G, GAP.julia_to_gap(gens))
  Q = GAP.Globals.FactorGroup(G, S)
  IdCheck = GAP.Globals.IdGroup(Q)
  return IdCheck
end

function max_ramified_prime(O::NfOrd, gtype::Vector{Int}, bound::fmpz)
  n = prod(gtype)
  fac = factor(n)
  m = Int(maximum(keys(fac.fac)))
  k = divexact(n, m)
  b1 = Int(root(bound, degree(O)*(m-1)*k))
  return b1
end

function _abelian_normal_extensions(F::FieldsTower, gtype::Array{Int, 1}, absbound::fmpz, IdCheck::Main.ForeignGAP.MPtr, only_real::Bool, IdG::Main.ForeignGAP.MPtr)
  K = F.field
  O = maximal_order(K)
  n = prod(gtype)
  inf_plc = Array{InfPlc, 1}()
  if abs(discriminant(O))^n > absbound
    return Vector{Hecke.ClassField{Hecke.MapRayClassGrp, GrpAbFinGenMap}}[]
  end
  bound = div(absbound, abs(discriminant(O))^n)
  l_conductors = conductors_with_restrictions(F, gtype, IdG, bound)
  @vprint :Fields 1 "   Number of conductors: $(length(l_conductors))\n"
  @vprint :FieldsNonFancy 1 "Number of conductors: $(length(l_conductors))\n"
  if length(l_conductors) == 0
    @vprint :Fields 1 "\n"
    return Vector{Hecke.ClassField{Hecke.MapRayClassGrp, GrpAbFinGenMap}}[]
  end
  @vprint :Fields 2 "\n"
  @vtime :Fields 2 Cl, mCl = class_group(O, use_aut = true)
  if mod(n, 2) == 0 && !only_real
    inf_plc = real_places(K)
  end
  expo = lcm(gtype)
  if length(l_conductors) == 1 && isone(l_conductors[1][1]) && isempty(l_conductors[1][2]) && !divisible(order(Cl)* (2^length(inf_plc)), fmpz(n))
    @vprint :Fields 1 "\n"
    return Vector{Hecke.ClassField{Hecke.MapRayClassGrp, GrpAbFinGenMap}}[]
  end
  Hecke.allow_cache!(mCl)
  @vtime :Fields 3 rcg_ctx = Hecke.rayclassgrp_ctx(O, expo)
  @vprint :Fields 1 "\n"
  j = -1
  first_group = true
  autos = F.generators_of_automorphisms
  class_fields_with_act = Tuple{Hecke.ClassField{Hecke.MapRayClassGrp, GrpAbFinGenMap}, Vector{GrpAbFinGenMap}}[]
  for k in l_conductors
    j += 1
    @vprint :Fields 1 "\e[1FConductors to test: $(length(l_conductors)-j) \n"
    @vprint :Fields 3 "\n\n"
    @vtime :Fields 3 r, mr = Hecke.ray_class_group_quo(O, k[1], k[2], inf_plc, rcg_ctx)
    if !Hecke._are_there_subs(r, gtype)
      continue
    end
    if first_group
      B1 = maximum([next_prime(max_ramified_prime(O, gtype, bound)), 211])
      lP = Hecke.find_gens(pseudo_inv(rcg_ctx.mC), PrimesSet(B1, -1))[1]
      rcg_ctx.mC.small_gens = lP
      first_group = false
    end
    mr.clgrpmap.small_gens = rcg_ctx.mC.small_gens
    @vtime :Fields 3 act = Hecke.induce_action_new(mr, autos)
    @vtime :Fields 3 ls = stable_subgroups(r, act, op = (x, y) -> quo(x, y, false)[2], quotype = gtype)
    Hecke.totally_positive_generators(mr)
    Dcond = Dict{Int, Array{GrpAbFinGenElem, 1}}()
    Ddisc = Dict{Tuple{Int, Int}, Array{GrpAbFinGenElem, 1}}()
    for s in ls
      @hassert :Fields 1 order(codomain(s)) == n
      C = ray_class_field(mr, s)
      fl = check_ext(C, bound, Dcond, Ddisc)
      if fl
        res_act = _action(s, act)
        @vtime :Fields 3 fl1 = check_group_extension(IdCheck, autos, res_act)
        if fl1
          push!(class_fields_with_act, (C, res_act))
        end
      end
    end
    @vprint :Fields 3 "\n\n"
  end
  if isempty(class_fields_with_act)
    return Vector{Hecke.ClassField{Hecke.MapRayClassGrp, GrpAbFinGenMap}}[]
  end
  emb_sub = F.subfields[end]
  @vprint :Fields 1 "\e[1F$(Hecke.clear_to_eol())Sieving $(length(class_fields_with_act)) abelian extensions\n"
  class_fields = Hecke.check_abelian_extensions(class_fields_with_act, autos, emb_sub)
  return class_fields
end

################################################################################
#
#  From the class fields to the computation of the defining polynomials
#
################################################################################

function from_class_fields_to_fields(class_fields::Vector{Hecke.ClassField{Hecke.MapRayClassGrp, GrpAbFinGenMap}}, autos::Vector{NfToNfMor}, grp_to_be_checked::Dict{Int, Main.ForeignGAP.MPtr})

  if isempty(class_fields)
    @vprint :Fields 1 "\e[1F$(Hecke.set_cursor_col())$(Hecke.clear_to_eol())"
    return Tuple{Hecke.NfRel_ns{nf_elem}, Vector{Hecke.NfRel_nsToNfRel_nsMor{nf_elem}}}[]
  end
  K = base_ring(class_fields[1])
  divisors_of_n = collect(keys(grp_to_be_checked))
  sort!(divisors_of_n)
  pclassfields = Vector{Vector{typeof(class_fields[1])}}(undef, length(divisors_of_n))
  right_grp = trues(length(class_fields))
  ind = 1
  for p in divisors_of_n
    it = findall(right_grp)
    if iszero(length(it))
      break
    end
    cfieldsp = Vector{typeof(class_fields[1])}(undef, length(class_fields))
    for i in it
      cfieldsp[i] = Hecke.maximal_p_subfield(class_fields[i], Int(p))
    end
    idE = grp_to_be_checked[p]
    if mod(order(torsion_unit_group(maximal_order(K))[1]), exponent(cfieldsp[1])) == 0
      compute_fields(cfieldsp, autos, idE, right_grp)
      pclassfields[ind] = cfieldsp
      ind += 1
      continue
    end
    E = GAP.Globals.SmallGroup(idE)
    S, H = max_ab_norm_sub_containing(E)
    if S == H
      compute_fields(cfieldsp, autos, idE, right_grp)
      pclassfields[ind] = cfieldsp
      ind += 1
      continue
    end
    @vprint :Fields 1 "\e[1F$(Hecke.set_cursor_col())$(Hecke.clear_to_eol())"
    # I can compute the fields over a subfield.
    #First, I need the subfields.
    K = base_field(class_fields[1])
    subfields = compute_subfields(K, E, H, S)
    computing_over_subfields(cfieldsp, subfields, idE, autos, right_grp)
    pclassfields[ind] = cfieldsp
    ind += 1
  end
  it = findall(right_grp)
  fields = Vector{Tuple{NfRel_ns{nf_elem}, Vector{Hecke.NfRel_nsToNfRel_nsMor{nf_elem}}}}(undef, length(it))
  ind = 1
  for i in it
    res = Vector{typeof(class_fields[1])}(undef, length(divisors_of_n))
    for j = 1:length(divisors_of_n)
      res[j] = pclassfields[j][i]
    end
    fields[ind] = _ext_and_autos(res, autos)
    ind += 1
  end
  return fields

end

function compute_fields(class_fields::Vector{Hecke.ClassField{Hecke.MapRayClassGrp, GrpAbFinGenMap}}, autos::Vector{NfToNfMor}, grp_to_be_checked::Main.ForeignGAP.MPtr, right_grp)
  K = base_field(class_fields[1])
  fields = Tuple{Hecke.NfRel_ns{nf_elem}, Vector{Hecke.NfRel_nsToNfRel_nsMor{nf_elem}}}[]
  expo = Int(exponent(codomain(class_fields[1].quotientmap)))
  #Since I want to compute as few Frobenius as possible, I want to first compute the extensions
  #whose set of divisors is maximal
  it = findall(right_grp)
  mins = Vector{fmpz}(undef, length(it))
  for i = 1:length(mins)
    mins[i] = minimum(defining_modulus(class_fields[it[i]])[1])
  end
  ismax = trues(length(mins))
  for i = 1:length(ismax)
    for j = i+1:length(ismax)
      if ismax[j]
        i2 = ppio(mins[i], mins[j])[2]
        if isone(i2)
          ismax[i] = false
          break
        end
        i3 = ppio(mins[j], mins[i])[2]
        if isone(i3)
          ismax[j] = false
        end
      end
    end
  end
  ord_class_fields = Vector{Int}(undef, length(ismax))
  j1 = 1
  j2 = length(ismax)
  for i = 1:length(ismax)
    if ismax[i]
      ord_class_fields[j1] = it[i]
      j1 += 1
    else
      ord_class_fields[j2] = it[i]
      j2 -= 1
    end
  end
  set_up_cycl_ext(K, expo, autos)

  for i in it
    C = class_fields[ord_class_fields[i]]
    L = number_field(C)
    autL = Hecke.absolute_automorphism_group(C, autos)
    Cpperm = perm_grp(autL)
    if !isone(gcd(degree(K), expo)) && GAP.Globals.IdGroup(Cpperm) != grp_to_be_checked
      right_grp[ord_class_fields[i]] = false
    end
  end
  return right_grp
end

function _ext_and_autos(resul::Vector{Hecke.ClassField{S, T}}, autos::Vector{NfToNfMor}) where S where T
  if length(resul) == 1
    return resul[1].A, resul[1].AbsAutGrpA
  end
  K = domain(autos[1])
  pols = Generic.Poly{nf_elem}[]
  for i = 1:length(resul)
    append!(pols, [Hecke.isunivariate(resul[i].A.pol[w])[2] for w = 1:length(resul[i].A.pol)])
  end
  L, gL = number_field(pols, cached = false, check = false)
  autL = Vector{Hecke.NfRel_nsToNfRel_nsMor{nf_elem}}()
  imgs_auts_base_field = Vector{Vector{Hecke.NfRel_nsElem{nf_elem}}}(undef, length(autos))
  for i = 1:length(autos)
    imgs_auts_base_field[i] = gens(L)
  end
  #Now, the automorphisms...
  for i = 1:length(resul)
    Cp = resul[i]
    p = ispower(degree(Cp))[2]
    w = 1
    while mod(degree(pols[w]), p) != 0
      w += 1
    end
    auts = Cp.AbsAutGrpA
    for phi in auts
      if phi.coeff_aut.prim_img == gen(K)
        imgsphi = gens(L)
        for ind = w:(w+length(Cp.A.pol)-1)
          imgsphi[ind] = evaluate(Hecke.isunivariate(phi.emb[ind-w+1].data)[2], imgsphi[ind])
        end
        push!(autL, Hecke.NfRel_nsToNfRel_nsMor(L, L, imgsphi))
      else
        ind_aut = 1
        while autos[ind_aut] != phi.coeff_aut
          ind_aut += 1
        end
        for ind = w:(w+length(Cp.A.pol)-1)
          imgs_auts_base_field[ind_aut][ind] = Hecke.msubst(phi.emb[ind-w+1].data, gL[w:(w+length(Cp.A.pol)-1)])
        end
      end
    end
  end
  for i = 1:length(imgs_auts_base_field)
    push!(autL, Hecke.NfRel_nsToNfRel_nsMor(L, L, autos[i], imgs_auts_base_field[i]))
  end
  return L, autL

end

function set_up_cycl_ext(K::AnticNumberField, n::Int, autK::Array{NfToNfMor, 1})
  fac = factor(n)
  for (p, v) in fac
    e = Int(p)^v
    if e == 2
      continue
    end
    C = cyclotomic_extension(K, e)
    if degree(C.Kr) == 1
      continue
    end
    auts = automorphisms(C, gens = autK, copy = false)
    @vprint :Fields 1 ": computing class group of cyclotomic extension of order $e"
    @vprint :FieldsNonFancy 1 "computing class group of cyclotomic extension of order $e\n"
    Cl, mCl = class_group(maximal_order(C.Ka), use_aut = true)
    Hecke.allow_cache!(mCl)
  end
  return nothing

end

function _action(t::Hecke.GrpAbFinGenMap, act::Array{Hecke.GrpAbFinGenMap,1})

  T = codomain(t)
  S, mS= snf(T)
  new_act = Array{Hecke.GrpAbFinGenMap, 1}(undef, length(act))
  for i = 1:length(act)
    res = mS.map*act[i].map*mS.imap
    new_act[i] = Hecke.GrpAbFinGenMap(S, S, res)
  end
  return new_act

end

################################################################################
#
#  Maximal subgroups which is in direct product with the derived subgroup
#
################################################################################

function max_ab_norm_sub_containing(G::Main.ForeignGAP.MPtr)

  D = GAP.Globals.DerivedSeries(G)
  H = D[end-1]
  G1 = GAP.Globals.Centralizer(G, H)
  if G1 == H
    return H, H
  end
  #First, I check if the centralizer split directly as a direct product.
  sc = GAP.Globals.ComplementClassesRepresentatives(G1, H)
  if !GAP.Globals.IsEmpty(sc)
    return sc[1], H
  end
  lS = GAP.gap_to_julia(GAP.Globals.ConjugacyClassesSubgroups(G1))
  #TODO: Subgroups in the quotient by H and not in the full group
  candidate = H
  sizecandidate = GAP.Globals.Size(H)
  for S in lS
    s = GAP.Globals.Representative(S)
    if !GAP.Globals.IsSubset(s, H)
      continue
    end
    sc = GAP.Globals.ComplementClassesRepresentatives(s, H)
    if GAP.Globals.IsEmpty(sc)
      continue
    end
    sncandidate = GAP.Globals.Size(s)
    if sncandidate > sizecandidate
      candidate = s
      sizecandidate = sncandidate
    end
  end
  return candidate, H
end

################################################################################
#
#  Computing the extension over subfields
#
################################################################################

function computing_over_subfields(class_fields, subfields, idE, autos, right_grp)
  it = findall(right_grp)
  new_class_fields, subs = translate_class_field_down(subfields, class_fields, it)
  translate_fields_up(class_fields, new_class_fields, subs, it)
  #Now, finally, the automorphisms computation and the isomorphism check
  for i in it
    C = class_fields[i]
    C1 = new_class_fields[i]
    mL = subs[1]
    indsubf = 1
    while domain(mL) != base_field(C1)
      indsubf += 1
      mL = subs[indsubf]
    end
    maprel = Hecke.NfRel_nsToNfRel_nsMor(C1.A, C.A, mL, gens(C.A))
    autsrelC1 = Hecke.rel_auto(C1)
    autsrelC = Vector{Hecke.NfRel_nsToNfRel_nsMor{nf_elem}}(undef, length(autsrelC1))
    for s = 1:length(autsrelC1)
      el = autsrelC1[s]
      autsrelC[s] = Hecke.NfRel_nsToNfRel_nsMor(C.A, C.A, [maprel(x) for x in el.emb])
    end
    rel_extend = Hecke.new_extend_aut(C, autos)
    autsA = vcat(rel_extend, autsrelC)
    C.AbsAutGrpA = autsA
    Cpperm = perm_grp(autsA)
    if GAP.Globals.IdGroup(Cpperm) != idE
      right_grp[i] = false
    end
  end
  it = findall(right_grp)
  fields = Vector{Tuple{Hecke.NfRel_ns{nf_elem}, Vector{Hecke.NfRel_nsToNfRel_nsMor{nf_elem}}}}(undef, length(it))
  ind = 1
  for i in it
    C = class_fields[i]
    fields[ind] = (C.A, C.AbsAutGrpA)
    ind += 1
  end
  return fields
end

function translate_extensions(mL::NfToNfMor, class_fields, new_class_fields, ctxK, it)

  to_be_done = Int[]
  L = domain(mL)
  OL = maximal_order(L)
  K = codomain(mL)
  OK = maximal_order(K)
  n = Int(exponent(codomain(class_fields[1].quotientmap)))
  ordext = Int(order(codomain(class_fields[1].quotientmap)))
  ctx = Hecke.rayclassgrp_ctx(OL, n)
  d = divexact(discriminant(maximal_order(K)), discriminant(OL)^(divexact(degree(K), degree(L))))
  f = factor(ideal(OL, d))
  F = factor(ideal(OK, d))
  for indclf in it
    if isassigned(new_class_fields, indclf)
      continue
    end
    C = class_fields[indclf]
    #First, I need a modulus.
    #I take the intersection of the modulus of C with L
    mR = C.rayclassgroupmap
    fM0 = mR.fact_mod
    fm0 = Dict{NfOrdIdl, Int}()
    for (p, v) in fM0
      p1 = Hecke.intersect_prime(mL, p)
      if !haskey(fm0, p1)
        fm0[p1] = v
      end
    end
    #Now, I have problems, so I need to add the ramification of the other extension.
    for (p, v) in f
      if !haskey(fm0, p)
        if isone(gcd(minimum(p), n))
          fm0[p] = 1
        else
          fm0[p] = v
        end
      else
        if !isone(gcd(minimum(p), n))
          fm0[p] = max(v, fm0[p])
        end
      end
    end
    infplc = InfPlc[]
    if iszero(mod(n, 2))
      infplc = real_places(L)
    end
    r, mr = Hecke.ray_class_group_quo(OL, fm0, infplc, ctx, check = false)
    if order(r) < degree(C)
      push!(to_be_done, indclf)
      continue
    end
    #Now, the norm group of K over L
    ngL, mngL = Hecke.norm_group(mL, mr)
    if !divisible(order(ngL), degree(C))
      push!(to_be_done, indclf)
      continue
    end
    #Finally, the norm group of C over L
    #I use the usual strategy, as in check_abelian_extension
    fM0 = merge(F, fM0)
    RM, mRM = Hecke.ray_class_group_quo(OK, fM0, InfPlc[], ctxK, check = false)
    lP, gS = Hecke.find_gens(mRM)
    listn = NfOrdIdl[norm(mL, x) for x in lP]
    # Create the map between R and r by taking norms
    preimgs = Vector{GrpAbFinGenElem}(undef, length(listn))
    for i = 1:length(preimgs)
      preimgs[i] = mr\listn[i]
    end
    proj = hom(gS, preimgs)
    #compute the norm group of C in RM
    prms = Vector{GrpAbFinGenElem}(undef, length(lP))
    for i = 1:length(lP)
      prms[i] = C.quotientmap(mR\lP[i])
    end
    RMtoR = hom(gS, prms)
    k, mk = kernel(RMtoR)
    mp = mk*proj
    ck, mck = cokernel(mp)
    expected_order = divexact(order(r), order(ngL)) * degree(C)
    if order(ck) != expected_order
      push!(to_be_done, indclf)
      continue
    end
    @vprint :Fields 1 "$(Hecke.set_cursor_col())$(Hecke.clear_to_eol())Doing field $(indclf) / $(length(class_fields))"
    s, ms = sub(ck, GrpAbFinGenElem[mck(mngL(x)) for x in gens(ngL)])
    mq1 = find_complement(ms)
    mqq = mck * mq1
    @assert domain(mqq) == r
    C1 = ray_class_field(mr, mqq)
    number_field(C1)
    @assert degree(C1) == degree(C)
    new_class_fields[indclf] = C1
  end
  @vprint :Fields 1 "$(Hecke.set_cursor_col())$(Hecke.clear_to_eol())"
  return to_be_done

end

#m is the injection representing a subgroup
#we want a complement of m
function find_complement(m::GrpAbFinGenMap)
  q, mq = cokernel(m)
  s, ms = snf(q)
  cmpl = GrpAbFinGenElem[]
  for x in gens(s)
    push!(cmpl, mq\(ms(x)))
  end
  return quo(codomain(m), cmpl, false)[2]
end

mutable struct BadField <: Exception
  K
end

function Base.show(io::IO, E::BadField)
  if isdefined(E, :K)
    println("Bad field $(K) encountered")
  else
    println("Bad field encountered")
  end
end

function create_sub(ss, iso, PermGAP, auts, K)
  lS = GAP.Globals.MinimalGeneratingSet(ss)
  lS1 = Main.ForeignGAP.MPtr[GAP.Globals.Image(iso, lS[i]) for i in 1:length(lS)]
  inds = Vector{Int}(undef, length(lS1))
  for j = 1:length(inds)
    el = lS1[j]
    i = 1
    while el != PermGAP[i]
      i += 1
    end
    inds[j] = i
  end
  mL = Hecke.fixed_field(K, [auts[i] for i in inds])[2]
  return mL
end

function compute_subfields(K::AnticNumberField, E, H, S)

  proj = GAP.Globals.NaturalHomomorphismByNormalSubgroup(E, H)
  Hn = GAP.Globals.ImagesSource(proj)
  imgsS = GAP.Globals.Image(proj, S)
  autsHn = GAP.Globals.Elements(GAP.Globals.AutomorphismGroup(Hn))
  orbitS = [imgsS]
  for i = 1:length(autsHn)
    varphi = autsHn[i]
    Sn = GAP.Globals.Image(varphi, imgsS)
    if !(Sn in orbitS)
      push!(orbitS, Sn)
    end
  end

  auts = automorphisms(K, copy = false)
  Hperm = _from_autos_to_perm(auts)
  Hauts = _perm_to_gap_grp(Hperm)
  PermGAP = Vector{Main.ForeignGAP.MPtr}(undef, length(Hperm))
  for w = 1:length(Hperm)
    PermGAP[w] = _perm_to_gap_perm(Hperm[w])
  end
  iso = GAP.Globals.IsomorphismGroups(Hn, Hauts)
  return (create_sub(ss, iso, PermGAP, auts, K) for ss in orbitS)

end

global deb = []

function translate_class_field_down(subfields, class_fields, it)

  new_class_fields = similar(class_fields)
  #Now, I translate the fields over the subfields.
  to_be_done = Int[i for i in it]
  created_subfields = NfToNfMor[]
  K = base_field(class_fields[1])
  OK = maximal_order(K)
  ctxK = Hecke.rayclassgrp_ctx(OK, exponent(class_fields[1]))
  for mL in subfields
    push!(created_subfields, mL)
    to_be_done_new = translate_extensions(mL, class_fields, new_class_fields, ctxK, it)
    if length(to_be_done_new) == 0
      return new_class_fields, created_subfields
    end
    to_be_done = to_be_done_new
  end
  push!(deb, (subfields, class_fields, it))
  error("Something went wrong!")
end

function translate_fields_up(class_fields, new_class_fields, subfields, it)
  K = base_field(class_fields[1])
  Ky = PolynomialRing(K, "y", cached = false)[1]
  for i in it
    C = class_fields[i]
    C1 = new_class_fields[i]
    mL = subfields[1]
    indsubf = 1
    while domain(mL) != base_field(C1)
      indsubf += 1
      mL = subfields[indsubf]
    end
    L = domain(mL)
    D = Dict{Int, NfToNfMor}()
    for i = 1:length(new_class_fields[i].cyc)
      d = degree(new_class_fields[1].cyc[i])
      if !haskey(D, d)
        CEK = cyclotomic_extension(K, d)
        CEL = cyclotomic_extension(L, d)
        mrel = Hecke.NfRelToNfRelMor(CEL.Kr, CEK.Kr, mL, gen(CEK.Kr))
        g = mrel(CEL.mp[1](gen(CEL.Ka)))
        mp = hom(CEL.Ka, CEK.Ka, CEK.mp[1]\(g))
        D[d] = mp
      end
    end
    cyc = Vector{Hecke.ClassField_pp{Hecke.MapRayClassGrp, GrpAbFinGenMap}}(undef, length(C1.cyc))
    for j = 1:length(cyc)
      Ccyc = C1.cyc[j]
      d = degree(Ccyc)
      #First, easy: the degree
      Cpp = Hecke.ClassField_pp{Hecke.MapRayClassGrp, GrpAbFinGenMap}()
      Cpp.rayclassgroupmap = C.rayclassgroupmap
      Cpp.degree = d
      #Then, the fac elem corresponding to the generator of the Kummer Extension
      Cpp.a = FacElem(Dict{nf_elem, fmpz}(D[d](x) => v for (x, v) in Ccyc.a))
      #Now, the Kummer extension
      Lzeta = codomain(D[d])
      Lt = PolynomialRing(Lzeta, "t", cached = false)[1]
      d1 = degree(Ccyc.K)
      coeffs = Vector{nf_elem}(undef, d1 + 1)
      coeffs[1] = D[d](coeff(Ccyc.K.pol, 0))
      for s = 2:length(coeffs)-1
        coeffs[s] = zero(Lzeta)
      end
      coeffs[end] = one(Lzeta)
      Cpp.K = number_field(Lt(coeffs), cached = false, check = false)[1]
      #The target extension
      fdef = Ccyc.A.pol
      coeffs1 = Vector{nf_elem}(undef, degree(fdef)+1)
      for s = 1:length(coeffs1)
        coeffs1[s] = mL(coeff(fdef, s-1))
      end
      Cpp.A = number_field(Ky(coeffs1), cached = false, check = false)[1]
      #Now, the primitive element of the target extension seen in Cpp.K
      mrel2 = Hecke.NfRelToNfRelMor(Ccyc.K, Cpp.K, D[d], gen(Cpp.K))
      Cpp.pe = mrel2(Ccyc.pe)
      Cpp.o = d1
      cyc[j] = Cpp
    end
    C.cyc = cyc
    C.A = number_field([cyc[j].A.pol for j = 1:length(cyc)], cached = false, check = false)[1]
  end
  return nothing

end


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

bat_default(::typeof(bat_valgrad), ::Val{:algorithm}, ::Function) = ForwardDiffAD()


text:
#==
\/\/\Parts/\/\/
==#
include("HControl.jl")
include("HDraw.jl")
using Compose
using DataFrames

function Scatter(x::Array, y::Array; shape = Circle(.5,.5,25),
        grid = Grid(4),
        features = [Axis(:X), Axis(:Y)],
        frame=Frame(1280,720,0mm,0mm),
        )
    points = Points(x, y, frame, shape)
    frame.add(points)
    glabels = GridLabels(x,y, grid)
    frame.add(grid)
    frame.add(glabels)
    for feature in features
        frame.add(feature);
    end
    return(frame)
end
function Scatter(x::DataFrame,
     y::Symbol, shape::Array, axiscolor=:lightblue, grid=Grid(3),
    frame=Frame(Frame(1280, 720, 0mm, 0mm))
    )
    topy = maximum(x[!, y])
    axisx = Line([(-1,-1), (-1,1), (1,1)],:blue)
   axisx_tag = axisx.update([(-1,-1), (-1,1), (1,1)])
    axisy = Line([(0,0), (0,1), (0,1)],:blue)
    axisy_tag = axisy.update([(0,0), (0,1), (0,1)])
    expression = string("compose(context()", ",", axisx_tag, axisy_tag)
    # Coordinate parsing -------
    y = x[!, y]
    counter = 0
    for col in eachcol(x, false)
        counter = counter + 1
        topx = maximum(col)
       if col !=  y
            for (i, w) in zip(col, y)
                expression = string(expression, "(context(), ")
                current_shape = shape[counter]
                inputx = (i / topx) * frame.width
                inputy = (w / topy) * frame.height
                exp = current_shape.update(inputx,inputy)
                expression = string(expression,string(exp),"),")
            end
        end
    end
    expression = string(expression,")")
    tt = transfertype(expression)
    frame.add(expression)
    return(frame)
end

mutable struct transfertype
   tag
end
@doc """Grid\n
The Grid function constructs a Grid type using percentages based on canvas sizes and divisions.\n
-------------------------------\n
======== PARAMETERS ======\n
(divisions,xlen=1280,ylen=720,colorx=:lightblue,colory=:lightblue,thickness=.2)\n
divisions:: An Int64 that determines the number of lines the grid will contain.\n
xlen:: An Int64 or Float64 representing the width of the canvas or area.\n
ylen:: An Int64 or Float64 representing the height of the canvas or area.\n
colorx:: A Hone Color symbol that will determine the color of the lines drawn on
the x axis. For more information on Hone Colors, refer to ?(HoneColors).\n
colory:: A Hone Color symbol that will determine the color of the lines drawn on
the y axis. For more information on Hone Colors, refer to ?(HoneColors).\n
thickness:: Will determine the width of each line in the grid.
--------------------------------\n
========= METHODS ========\n
Grid.show() - Shows the Grid's composition.\n
Grid.update() - Returns the Grid's Meta-tag.\n
Grid.save(URI) - Saves Grid as a Scalable Vector Graphic (SVG) File in the provided URI"""
function Grid(divisions,frame=Frame(1280,720,0mm,0mm,0mm,0mm),colorx=:lightblue,colory=:lightblue,thickness=.2)
    xlen = frame.width
    ylen = frame.height
    division_amountx = ylen / divisions
    division_amounty = xlen / divisions
    total = 0
    Xexpression = "(context(), "
    while total < xlen
        total = total + division_amountx
        linedraw = Line([(0,total),(xlen,total)],:lightblue,thickness)
        exp = linedraw.update(:This_symbol_means_nothing)
        Xexpression = string(Xexpression,string(exp))
    end
    Xexpression = string(Xexpression, "),")
    total = 0
    Yexpression = "(context(),"
    while total < ylen
        total = total + division_amounty
        println(total)
        linedraw = Line([(total,0),(total,ylen)],:lightblue,thickness)
        exp = linedraw.update(:hi)
        Yexpression = string(Yexpression,string(exp))
    end
    composexp = string(string(Xexpression), string(Yexpression),"),")
    vv = Meta.parse(string("compose(context(),",composexp,")"))
    composition = eval(vv)
    save(name) = draw(SVG(name), composition);
    update() = string(composexp)
    show() = composition
    tag = composexp
    (var)->(update;composexp;show;save;tag;division_amountx;division_amounty;frame;divisions)
end
#--------
function _arrayline(x,y,color=:orange,weight=2,axiscolor=:lightblue,
    grid=Grid(3), frame=Frame(1280,720,0mm,0mm,0mm,0mm))
    pairs = []
    topy = maximum(y)
    topx = maximum(x)
    for (i,w) in zip(x,y)
        x = (i / topx * frame.width)
        y = (w / topy * frame.height)
        pair = Tuple([x,y])
        push!(pairs,pair)
    end
    pairs = Array{Tuple{Float64,Real},1}(pairs)
    lin = Line(pairs,color,weight)
    axisx = Line([(0,frame.height), (frame.width,frame.height)],axiscolor)
    axisy = Line([(0,0), (0,frame.height)],axiscolor)
    expression = string("(context(),",lin.update(:foo),"), (context(), ",grid.update(),
    "), (context(), ",axisx.update(:foo),axisy.update(:foo),")")
    tt = transfertype(expression)
    frame.add(tt)
    show() = frame.show()
    tree() = introspect(composition)
    save(name) = draw(SVG(name), composition);
    get_frame() = frame
    (var)->(show;composition;tree;save;get_frame)
end
@doc """Line\n
The Line function takes either a DataFrame and Symbol OR two Arrays and returns
     a Hone line plot object.\n
-------------------------------\n
======== PARAMETERS ======\n
(x,y,shape=Circle(.5,.5,25),axiscolor=:lightblue,debug=false, grid=Grid(3), custom="", frame=Frame(1280,720,0mm,0mm,0mm,0mm)\n
x:: An Array of X coordinates OR a DataFrame containing both X's and Y's to be plotted.\n
y:: An array of corresponding Y coordinates OR a Symbol representing which DataFrame column to use as Y.\n
shape:: Takes a Hone.jl Shape object OR an iterable list of shapes with the same length as the number of X's
in the provided DataFrame, or X\n
axiscolor:: Symbol representing the color of both the X and Y axis's. For more information on available colors, please
use ?(HoneColors)\n
debug:: Boolean that determines whether or not the function will print the meta expression on creation.\n
grid:: Takes a Hone Grid object. For more information on grids, please use ?(Grid)\n
custom:: Custom takes a meta expression (as a an unparsed string), and can be used to add customizable features
 to a plot.\n
 frame:: Takes a Hone Frame object. For more information on Frames, please use ?(Frame)\n
--------------------------------\n
========= METHODS ========\n
Line.show() - Displays the frame which contains the plot's composition.\n
Line.tree() - Shows an introspection of contexts contained in the plot's composition.\n
Line.save(URI) - Saves the plot as a Scalable Vector Graphic (SVG) file at the given URI.\n
Line.get_frame() - Returns the frame which contains the plot as a child object."""
Linear(x::Array,y::Array,shape) = _arrayline(x,y,shape)
@doc """Colors\n
You can use HoneColors(color) to view a color.\n
color:: Symbol used to represent the desired color.\n
----------------------------------------------\n
========== COLORS ============\n
:blue, :green, :lightblue, :pink, :orange, :red, :black, :gray, :white, :yellow,
 :purple, :magenta, :lime
"""
function HoneColors(color)
    d = Circle(.5,.5,.5,color)
    d.show()
end
@doc """GridLabels\n
GridLabels provide labels to individual lines on the grid of your plot.\n
-------------------------------\n
======== PARAMETERS ======\n
(x,y,grid)\n
x:: An Array of X coordinates OR a DataFrame containing both X's and Y's to be plotted.\n
y:: An array of corresponding Y coordinates OR a Symbol representing which DataFrame column to use as Y.\n
grid:: Takes a Hone Grid object. For more information on grids, please use ?(Grid)\n
--------------------------------\n
"""
function GridLabels(x,y,grid)
    frame = grid.frame
    divamounty = grid.division_amounty
    divamountx = grid.division_amountx
    total = divamountx
    topx = maximum(x)
    topy = maximum(y)
    xlabels = []
    while total < (divamountx * grid.divisions)
        percentage = total / frame.height
        curr_label = topy * percentage
        push!(xlabels,(curr_label, frame.height - total))
        total += divamountx
    end
    xtags = ""
    for (key,data) in xlabels
        textlabel = Label(string(round(key)), 40 , data, "", 3)
        xtags = string(xtags, textlabel.tag)
    end
    # Ys
    total = divamounty
    ylabels = []
    while total < (divamounty * grid.divisions)
        percentage = total / frame.width
        curr_label = topx * percentage
        push!(ylabels,(curr_label, frame.width - total))
        total += divamounty
    end
    ytags = ""
    for (key,data) in ylabels
        textlabel = Label(string(round(key)), frame.width - data ,frame.height - 40 , "", 3)
        ytags = string(ytags, textlabel.tag)
    end
    tag = string(xtags,ytags)
    ()->(tag;xtags;ytags)
end
function Points(x, y, frame=Frame(1280,720,0mm,0mm,0mm,0mm), shape=Circle(.5, .5, 25))
    fheight = frame.height
    fwidth = frame.width
    topx = maximum(x)
    topy = maximum(y)
   express = string("")
    # Coordinate parsing -------
    for (i, w) in zip(x, y)
        inputx = (i / topx) * fwidth
        inputy = (w / topy) * fheight
        inputy = fheight - inputy
        println(inputx)
        println(inputy)
        exp = shape.update(inputx,inputy)
        express = string(express,string(exp))
    end
    tag = express
    show() = eval(Meta.parse(string("compose(context(), ", tag,")")))
    (var)->(tag;show)
end
function Axis(orientation=:X, axiscolor = :gray, frame=Frame(1280,720,0mm,0mm,0mm,0mm),)
    if orientation == :X
        pairs = [(0,frame.height), (frame.width,frame.height)]
    else orientation == :Y
        pairs = [(0,0),(0, frame.height)]
    end
    axis = Line(pairs,axiscolor)
    tag = axis.update([pairs])
    (var)->(tag)
end
#---------------------------


text:
# Wrapper for current backend as `Symbol`
mutable struct CurrentADBackend
    backend::Symbol
end

const _current_adbackend = CurrentADBackend(:none)
const _adbackends = Symbol[]

"""
    adbackend(backend::Symbol = :default) -> Symbol

Check if autodiff `backend` is valid (see [`adbackends`](@ref)) and return
formatted. Note that `backend=:default` returns the current backend.
"""
adbackend() = _current_adbackend.backend
adbackend(backend::Symbol) = adbackend(Val(backend))
adbackend(backend::Val{:default}) = adbackend()

function adbackend(backend::Val{T}) where {T}
    T in _adbackends && return T
    throw(ArgumentError("Invalid AD backend $(T). Valid options are $(adbackends())."))
end

"""
    adbackend!(backend::Symbol)

Set current backend for autodiff to `backend`.
"""
function adbackend!(backend)
    backend = adbackend(backend)
    _current_adbackend.backend = backend
    return backend
end

"""
    adbackends() -> Vector{Symbol}

Get vector of currently valid AD backends.
"""
adbackends() = _adbackends

"""
    _gradient(f, p::Number, backend = adbackend()) -> Number
    _gradient(f, p::Array, backend = adbackend()) -> Array

Compute gradient of function `f` at `x` using AD backend `backend`.
"""
_gradient(f, p) = _gradient(f, p, Val(adbackend()))
_gradient(f, p, backend::Symbol) = _gradient(f, p, Val(adbackend(backend)))

"""
    _jacobian(f, p, backend = adbackend()) -> Array

Compute Jacobian matrix of function `f` at `p` using AD backend `backend`.
Inputs and outputs of `f` are vectorized.
"""
_jacobian(f, p) = _jacobian(f, p, Val(adbackend()))
_jacobian(f, p, backend::Symbol) = _jacobian(f, p, Val(adbackend(backend)))

# Finite differences

push!(_adbackends, :finitedifferences)
adbackend!(:finitedifferences)

_default_fdm() = central_fdm(5, 1)

_gradient(f, p, backend::Val{:finitedifferences}) = _gradient(f, p, _default_fdm())

function _gradient(f, p, fdm::FiniteDifferences.FiniteDifferenceMethod)
    return FiniteDifferences.grad(fdm, f, p)[1]
end

_jacobian(f, p, backend::Val{:finitedifferences}) = _jacobian(f, p, _default_fdm())

function _jacobian(f, p, fdm::FiniteDifferences.FiniteDifferenceMethod)
    return FiniteDifferences.jacobian(fdm, f, p)[1]
end


text:
"""
    rk4_step(f_xy, h, x, y)
    -> (x1, y1)
Returns coordinates (x1, y1) after one Runge-Kutta 4th order step.
        f is a tuple-valued gradient in two dimensions
        (x, y) is the starting point
        h is a step-size quantity. If the integration variable is time,
        it defines the length of time betweeen points. It can be negative.
"""
function rk4_step(f_xy, h, x, y)
    fx0, fy0 = f_xy(x, y)
    x1 = x + fx0 * h * 0.5
    y1 = y + fy0 * h * 0.5
    fx1, fy1 = f_xy(x1, y1)
    x2 = x + fx1 * h * 0.5
    y2 = y + fy1 * h * 0.5
    fx2, fy2  = f_xy(x2, y2)
    x3 = x + fx2 * h
    y3 = y + fy2 * h
    fx3, fy3 = f_xy(x3, y3)
    x +=  h * ( fx0  + 2∙fx1 + 2∙fx2 + fx3 ) / 6
    y +=  h * ( fy0  + 2∙fy1 + 2∙fy2 + fy3 ) / 6
    x, y
end

"""
    rk4_step!(f, vx, vy, h, n)

Coordinates (vx[n + 1], vy[n + 1]) are updated with the estimated position,
along function f.
        f is a tuple-valued gradient in two dimensions
        vx and vy are vectors representing coordinates like (vx[n], vy[n])
        h is a step-size quantity. If the integration variable is time,
        it defines the length of time betweeen points. It can be negative.
"""
function rk4_step!(f, vx, vy, h, n)
    fx0, fy0 = f(vx[n], vy[n])
    x1 = vx[n] + fx0 * h * 0.5
    y1 = vy[n] + fy0 * h * 0.5
    fx1, fy1 = f(x1, y1)
    x2 = vx[n] + fx1 * h * 0.5
    y2 = vy[n] + fy1 * h * 0.5
    fx2, fy2  = f(x2, y2)
    x3 = vx[n] + fx2 * h
    y3 = vy[n] + fy2 * h
    fx3, fy3 = f(x3, y3)
    vx[n + 1] = vx[n] + 1/6 * h * ( fx0  + 2∙fx1 + 2∙fx2 + fx3 )
    vy[n + 1] = vy[n] + 1/6 * h * ( fy0  + 2∙fy1 + 2∙fy2 + fy3 )
end

"""
    rk4_steps!(f, vx, vy, h)

Coordinates (vx[n > 1], vy[n > 1]) are updated with the estimated positions,
along function f.
     f(x,y) returns a tuple-valued gradient in two dimensions
     vx and vy are vectors representing coordinates like (vx[n], vy[n])
     h is a step-size quantity. If the integration variable is time,
     it defines the length of time betweeen points. It can be negative
"""
function rk4_steps!(f, vx, vy, h)
    @assert length(vx) == length(vy)
    if !isnan(f(vx[1], vy[1])[1])
        for n in 1:(length(vx) - 1)
            rk4_step!(f, vx, vy, h, n)
        end
    else
        fill!(vx, NaN * vx[1])
        fill!(vy, NaN * vy[1])
    end
    nothing
end


text:
N = 64; G = path_graph(N)
L = Matrix(laplacian_matrix(G))
𝛌, 𝚽 = eigen(L); 𝚽 = 𝚽 .* sign.(𝚽[1,:])'
Q = incidence_matrix(G; oriented = true)
∇𝚽 = Q' * 𝚽


text:
# # Vibration of a cube of nearly incompressible material

# Source code: [`unit_cube_modes_tut.jl`](unit_cube_modes_tut.jl)

# ## Description

# Compute the free-vibration spectrum of a unit cube of nearly
# incompressible isotropic material, E = 1, ν = 0.499, and ρ = 1 (refer to [1]).

# The solution with the `FinEtools` package is compared with a commercial
# software  solution, and hence we also export the model to Abaqus.

# ## References

# [1] Puso MA, Solberg J (2006) A stabilized nodally integrated tetrahedral.
# International Journal for Numerical Methods in Engineering 67: 841-867.
# [2] P. Krysl, Mean-strain 8-node hexahedron with optimized energy-sampling
# stabilization, Finite Elements in Analysis and Design 108 (2016) 41–53.

# ![](unit_cube-mode7.png)

# ## Goals

# - Show how to generate simple mesh.
# - Show how to set up the discrete model for a free vibration problem.
# - Show how to export the model to Abaqus.

##
# ## Definitions

# This is the finite element toolkit itself.
using FinEtools

# The linear stress analysis application is implemented in this package.
using FinEtoolsDeforLinear

# Convenience import.
using FinEtools.MeshExportModule

# The eigenvalue problem is solved with the Lanczos algorithm from this package.
using Arpack

# The material properties and dimensions are defined with physical units.
E = 1*phun("PA");
nu = 0.499;
rho = 1*phun("KG/M^3");
a = 1*phun("M"); # length of the side of the cube

# We generate a mesh of  5 x 5 x 5 serendipity 20-node hexahedral elements in a
# regular grid.
fens, fes  = H20block(a, a, a, 5, 5, 5);

# The problem is solved in three dimensions and hence we create the displacement
# field as three-dimensional with three displacement components per node. The
# degrees of freedom are then numbered  (note that no essential boundary
# conditions are applied, since the cube is free-floating).
geom = NodalField(fens.xyz)
u = NodalField(zeros(size(fens.xyz,1),3)) # displacement field
numberdofs!(u);

# The model is fully three-dimensional, and hence the material model  and the
# FEMM created below need to refer to an appropriate model-reduction scheme.
MR = DeforModelRed3D
material = MatDeforElastIso(MR, rho, E, nu, 0.0);

# Note that we compute the stiffness  and the mass matrix using different FEMMs.
# The difference  is only the quadrature rule chosen: in order to make the mass
# matrix  non-singular, an accurate  Gauss rule  needs to be used, whereas for
# the stiffness matrix we want to avoid the excessive stiffness  and therefore
# the reduced Gauss rule is used.
femm = FEMMDeforLinear(MR, IntegDomain(fes, GaussRule(3,2)), material);
K = stiffness(femm, geom, u)
femm = FEMMDeforLinear(MR, IntegDomain(fes, GaussRule(3,3)), material)
M = mass(femm, geom, u);

# The free vibration problem  can now be solved.   In order for the eigenvalue
# solver  to work well, we apply mass-shifting (otherwise the first matrix
# given to the solver – stiffness – would be singular). We specify the number
# of eigenvalues to solve for, and we  guess the frequency  with which to shift
# as 0.01 Hz.
neigvs = 20 # how many eigenvalues
OmegaShift = (0.01*2*pi)^2; # The frequency with which to shift

# The `eigs` routine can now be invoked to solve for a given number of
# frequencies from the smallest-magnitude end of the spectrum. Note that the
# mass shifting  needs to be undone when the solution is obtained.
evals, evecs, nconv = eigs(K+OmegaShift*M, M; nev=neigvs, which=:SM)
@show nconv == neigvs
evals = evals .- OmegaShift;
fs = real(sqrt.(complex(evals)))/(2*pi)
sigdig(n) = round(n * 10000) / 10000
println("Eigenvalues: $(sigdig.(fs)) [Hz]")
# The first nonzero frequency, frequency 7, should be around .263 Hz.

# The computed mode can be visualized in Paraview. Use the  "Animation view" to
# produce moving pictures for the mode.
mode = 7
scattersysvec!(u, evecs[:,mode])
File =  "unit_cube_modes.vtk"
vtkexportmesh(File, fens, fes; vectors=[("mode$mode", u.values)])
@async run(`"paraview.exe" $File`);

# Finally  we export the model to Abaqus.  Note that we specify the mass
# density property (necessary for dynamics).
AE = AbaqusExporter("unit_cube_modes_h20");
HEADING(AE, "Vibration modes of unit cube  of almost incompressible material.");
COMMENT(AE, "The  first six frequencies are rigid body modes.");
COMMENT(AE, "The  first nonzero frequency (7) should be around 0.26 Hz");
PART(AE, "part1");
END_PART(AE);
ASSEMBLY(AE, "ASSEM1");
INSTANCE(AE, "INSTNC1", "PART1");
NODE(AE, fens.xyz);
COMMENT(AE, "The hybrid form of the serendipity hexahedron is chosen because");
COMMENT(AE, "the material is  nearly incompressible.");
ELEMENT(AE, "C3D20RH", "AllElements", 1, connasarray(fes))
ORIENTATION(AE, "GlobalOrientation", vec([1. 0 0]), vec([0 1. 0]));
SOLID_SECTION(AE, "elasticity", "GlobalOrientation", "AllElements");
END_INSTANCE(AE);
END_ASSEMBLY(AE);
MATERIAL(AE, "elasticity")
ELASTIC(AE, E, nu)
DENSITY(AE, rho)
STEP_FREQUENCY(AE, neigvs)
END_STEP(AE)
close(AE)

# It remains is to load the model into Abaqus and execute it as a job.
# Alternatively Abaqus can be called on the input file to carry out the
# analysis at the command line as

# ```
# abaqus job=unit_cube_modes_h20.inp
# ```

# The output database `unit_cube_modes_h20.odb` can then be loaded for
# postprocessing, for instance from the command line as

# ```
# abaqus viewer database=unit_cube_modes_h20.odb
# ```

# Don't forget to compare the computed frequencies and the mode shapes.  For
# instance, the first six frequencies should be nearly 0, and the seventh
# frequency should be approximately  0.262 Hz. There may be  very minor
# differences due to the fact that  the FinEtools formulation is purely
# displacement-based, whereas the Abaqus model is hybrid (displacement plus
# pressure).


text:
import Distributions: Dirichlet

# draw random stochastic matrix
rsm(dx::Int,dy::Int)=mapslices(x->rand(Dirichlet(x)),ones(dx,dy),2)
rsm(k::Int)=rsm(k,k)

# draw a sparse random stochastic matrix
function rssm(dim::Int,density=.2)
	m=Array(Float64,(dim,dim))
	for i=1:dim
		m[i,:]=Base.sprand(1,dim,density)
		m[i,rand(1:dim)]+=rand()
	end
	sparse(m./sum(m,2))
end

#stochastic matrix normalization, ie each row is normalized to sum to 1
function nsm!(m::Array{Float64,2})
	for i in 1:size(m)[1]
		m[i,:]/=sum(m[i,:])
	end
end

function nsm(m::Array{Float64,2})
	nsm!(m)
	m
end

# REPARAMETRIZATION
# frequently useful when optimizing over a stochastic matrix

#assumes q is square
z2q(z::Array{Float64,1})=z2q(z,Int(.5*(1+sqrt(1+4*length(z)))))

function z2q(z::Array{Float64,1},dx::Int)
	#add some size checks
	dy=Int(length(z)/dx)+1
	q=zeros(dx,dy)
	mz=reshape(z,dx,dy-1)
	for ix=1:dx
		ez=[exp(slice(mz,ix,:));1]
		q[ix,:]=ez/sum(ez)
	end
	q
end

function q2z(q::Array{Float64,2})
	dx=size(q,1)
	z=zeros(dx,dx-1)
	for ix=1:dx
		z[ix,:]=log(sub(q,ix,1:dx-1)/q[ix,dx])
	end
	vec(z)
end

# jacobian d(vec(q))/d(vec(z))
function z2qjac(z::Array{Float64,1})
	dz=length(z)  #dz=dx^2-dx
	dx=round(Int,.5*(1+sqrt(1+4*dz)))
	# println(dx)
	jac=zeros(dx,dx,dx,dx-1)
	mz=reshape(z,dx,dx-1)
	ez=zeros(dx-1)
	for ix=1:dx
		ez[:]=exp(mz[ix,:])
		sez=(1+sum(ez))
		sez2=sez^2
		for jx=1:dx-1
			for jz=1:dx-1
				jac[ix,jx,ix,jz]=((jx==jz)*ez[jz]*sez-ez[jz]*ez[jx])/sez2
			end
		end
		#jx=dx:
		for jz=1:dx-1
			jac[ix,dx,ix,jz]=-ez[jz]/sez2
		end
	end
	reshape(jac,dx*dx,dz)
end


text:
include("bonds.jl")
include("tensorkit_timers.jl")

import LinearAlgebra
LinearAlgebra.BLAS.set_num_threads(1)

import TensorOperations
TensorOperations.enable_cache(; maxrelsize = 0.7)
import TensorKit
using TensorKit: ℂ, Z2Space, U1Space
using DelimitedFiles

K = 100 # number of repititions; scale up for more reliable testing.

# MPO Contractions

# TRIVIAL = no symmetry
N = length(mpodims_triv[1]) # number of different cases
mpo_triv_times = zeros(Float64, (K, N))
mpo_triv_times_gc = zeros(Float64, (K, N))
for i = 1:N
    Dmps = mpodims_triv[1][i]
    Dmpo = mpodims_triv[2][i]
    Dphys = mpodims_triv[3][i]

    Vmps = ℂ^Dmps
    Vmpo = ℂ^Dmpo
    Vphys = ℂ^Dphys

    mpo_timer = TensorKitTimers.mpo_timer(; Vmps = Vmps, Vmpo = Vmpo, Vphys = Vphys)

    times, times_gc = mpo_timer(; outer = K)
    mpo_triv_times[:, i] = times
    mpo_triv_times_gc[:, i] = times_gc
    empty!(TensorOperations.cache)

    tavg = sum(times)/K

    println("MPO Contraction, trivial symmetry: case $i out of $N finished at average time $tavg")
end

writedlm("results/tensorkit_mpo_triv.txt", mpo_triv_times)
writedlm("results/tensorkit_mpo_triv_gc.txt", mpo_triv_times_gc)

# Z2 Symmetry
N = length(mpodims_z2[1]) # number of different cases
mpo_z2_times = zeros(Float64, (K, N))
mpo_z2_times_gc = zeros(Float64, (K, N))
for i = 1:N
    Dmps = distributeZ2(mpodims_z2[1][i])
    Dmpo = distributeZ2(mpodims_z2[2][i])
    Dphys = distributeZ2(mpodims_z2[3][i])

    Vmps = Z2Space((a=>b for (a,b) in Dmps)...)
    Vmpo = Z2Space((a=>b for (a,b) in Dmpo)...)
    Vphys = Z2Space((a=>b for (a,b) in Dphys)...)

    mpo_timer = TensorKitTimers.mpo_timer(; Vmps = Vmps, Vmpo = Vmpo, Vphys = Vphys)

    times, times_gc = mpo_timer(; outer = K)
    mpo_z2_times[:, i] = times
    mpo_z2_times_gc[:, i] = times_gc
    empty!(TensorOperations.cache)

    tavg = sum(times)/K

    println("MPO Contraction, Z2 symmetry: case $i out of $N finished at average time $tavg")
end

writedlm("results/tensorkit_mpo_z2.txt", mpo_z2_times)
writedlm("results/tensorkit_mpo_z2_gc.txt", mpo_z2_times_gc)

# U1 Symmetry
N = length(mpodims_u1[1]) # number of different cases
mpo_u1_times = zeros(Float64, (K, N))
mpo_u1_times_gc = zeros(Float64, (K, N))
for i = 1:N
    Dmps = distributeU1(mpodims_u1[1][i])
    Dmpo = distributeU1(mpodims_u1[2][i])
    Dphys = distributeU1(mpodims_u1[3][i])

    Vmps = U1Space((a=>b for (a,b) in Dmps)...)
    Vmpo = U1Space((a=>b for (a,b) in Dmpo)...)
    Vphys = U1Space((a=>b for (a,b) in Dphys)...)

    mpo_timer = TensorKitTimers.mpo_timer(; Vmps = Vmps, Vmpo = Vmpo, Vphys = Vphys)

    times, times_gc = mpo_timer(; outer = K)
    mpo_u1_times[:, i] = times
    mpo_u1_times_gc[:, i] = times_gc
    empty!(TensorOperations.cache)

    tavg = sum(times)/K

    println("MPO Contraction, U1 symmetry: case $i out of $N finished at average time $tavg")
end

writedlm("results/tensorkit_mpo_u1.txt", mpo_u1_times)
writedlm("results/tensorkit_mpo_u1_gc.txt", mpo_u1_times_gc)

## PEPO Contractions

# TRIVIAL = no symmetry
N = length(pepodims_triv[1]) # number of different cases
pepo_triv_times = zeros(Float64, (K, N))
pepo_triv_times_gc = zeros(Float64, (K, N))
for i = 1:N
    Dpeps = pepodims_triv[1][i]
    Dpepo = pepodims_triv[2][i]
    Dphys = pepodims_triv[3][i]
    Denv = pepodims_triv[4][i]

    Vpeps = ℂ^Dpeps
    Vpepo = ℂ^Dpepo
    Vphys = ℂ^Dphys
    Venv = ℂ^Denv

    pepo_timer = TensorKitTimers.pepo_timer(; Vpeps = Vpeps, Vpepo = Vpepo, Vphys = Vphys, Venv = Venv)

    times, times_gc = pepo_timer(; outer = K)
    pepo_triv_times[:, i] = times
    pepo_triv_times_gc[:, i] = times_gc
    empty!(TensorOperations.cache)

    tavg = sum(times)/K

    println("PEPO Contraction, trivial symmetry: case $i out of $N finished at average time $tavg")
end

writedlm("results/tensorkit_pepo_triv.txt", pepo_triv_times)
writedlm("results/tensorkit_pepo_triv_gc.txt", pepo_triv_times_gc)

# Z2 Symmetry
N = length(pepodims_z2[1]) # number of different cases
pepo_z2_times = zeros(Float64, (K, N))
pepo_z2_times_gc = zeros(Float64, (K, N))
for i = 1:N
    Dpeps = distributeZ2(pepodims_z2[1][i])
    Dpepo = distributeZ2(pepodims_z2[2][i])
    Dphys = distributeZ2(pepodims_z2[3][i])
    Denv = distributeZ2(pepodims_z2[4][i])

    Vpeps = Z2Space((a=>b for (a,b) in Dpeps)...)
    Vpepo = Z2Space((a=>b for (a,b) in Dpepo)...)
    Vphys = Z2Space((a=>b for (a,b) in Dphys)...)
    Venv = Z2Space((a=>b for (a,b) in Denv)...)

    pepo_timer = TensorKitTimers.pepo_timer(; Vpeps = Vpeps, Vpepo = Vpepo, Vphys = Vphys, Venv = Venv)

    times, times_gc = pepo_timer(; outer = K)
    pepo_z2_times[:, i] = times
    pepo_z2_times_gc[:, i] = times_gc
    empty!(TensorOperations.cache)

    tavg = sum(times)/K

    println("PEPO Contraction, Z2 symmetry: case $i out of $N finished at average time $tavg")
end

writedlm("results/tensorkit_pepo_z2.txt", pepo_z2_times)
writedlm("results/tensorkit_pepo_z2_gc.txt", pepo_z2_times_gc)

# U1 Symmetry
N = length(pepodims_u1[1]) # number of different cases
pepo_u1_times = zeros(Float64, (K, N))
pepo_u1_times_gc = zeros(Float64, (K, N))
for i = 1:N
    Dpeps = distributeU1(pepodims_u1[1][i])
    Dpepo = distributeU1(pepodims_u1[2][i])
    Dphys = distributeU1(pepodims_u1[3][i])
    Denv = distributeU1(pepodims_u1[4][i])

    Vpeps = U1Space((a=>b for (a,b) in Dpeps)...)
    Vpepo = U1Space((a=>b for (a,b) in Dpepo)...)
    Vphys = U1Space((a=>b for (a,b) in Dphys)...)
    Venv = U1Space((a=>b for (a,b) in Denv)...)

    pepo_timer = TensorKitTimers.pepo_timer(; Vpeps = Vpeps, Vpepo = Vpepo, Vphys = Vphys, Venv = Venv)

    times, times_gc = pepo_timer(; outer = K)
    pepo_u1_times[:, i] = times
    pepo_u1_times_gc[:, i] = times_gc
    empty!(TensorOperations.cache)

    tavg = sum(times)/K

    println("PEPO Contraction, U1 symmetry: case $i out of $N finished at average time $tavg")
end

writedlm("results/tensorkit_pepo_u1.txt", pepo_u1_times)
writedlm("results/tensorkit_pepo_u1_gc.txt", pepo_u1_times_gc)

# ## MERA Contractions

# TRIVIAL = no symmetry
N = length(meradims_triv[1]) # number of different cases
mera_triv_times = zeros(Float64, (K, N))
mera_triv_times_gc = zeros(Float64, (K, N))
for i = 1:N
    Dmera = meradims_triv[1][i]

    Vmera = ℂ^Dmera

    mera_timer = TensorKitTimers.mera_timer(; Vmera = Vmera)

    times, times_gc = mera_timer(; outer = K)

    mera_triv_times[:, i] = times
    mera_triv_times_gc[:, i] = times_gc
    empty!(TensorOperations.cache)

    tavg = sum(times)/K

    println("MERA Contraction, trivial symmetry: case $i out of $N finished at average time $tavg")
end

writedlm("results/tensorkit_mera_triv.txt", mera_triv_times)
writedlm("results/tensorkit_mera_triv_gc.txt", mera_triv_times_gc)

# Z2 Symmetry
N = length(meradims_z2[1]) # number of different cases
mera_z2_times = zeros(Float64, (K, N))
mera_z2_times_gc = zeros(Float64, (K, N))
for i = 1:N
    Dmera = distributeZ2(meradims_z2[1][i])

    Vmera = Z2Space((a=>b for (a,b) in Dmera)...)

    mera_timer = TensorKitTimers.mera_timer(; Vmera = Vmera)

    times, times_gc = mera_timer(; outer = K)

    mera_z2_times[:, i] = times
    mera_z2_times_gc[:, i] = times_gc
    empty!(TensorOperations.cache)

    tavg = sum(times)/K

    println("MERA Contraction, Z2 symmetry: case $i out of $N finished at average time $tavg")
end

writedlm("results/tensorkit_mera_z2.txt", mera_z2_times)
writedlm("results/tensorkit_mera_z2_gc.txt", mera_z2_times_gc)

# U1 Symmetry
N = length(meradims_u1[1]) # number of different cases
mera_u1_times = zeros(Float64, (K, N))
mera_u1_times_gc = zeros(Float64, (K, N))
for i = 1:N
    Dmera = distributeU1(meradims_u1[1][i])

    Vmera = U1Space((a=>b for (a,b) in Dmera)...)

    mera_timer = TensorKitTimers.mera_timer(; Vmera = Vmera)

    times, times_gc = mera_timer(; outer = K)

    mera_u1_times[:, i] = times
    mera_u1_times_gc[:, i] = times_gc
    empty!(TensorOperations.cache)

    tavg = sum(times)/K

    println("MERA Contraction, U1 symmetry: case $i out of $N finished at average time $tavg")
end

writedlm("results/tensorkit_mera_u1.txt", mera_u1_times)
writedlm("results/tensorkit_mera_u1_gc.txt", mera_u1_times_gc)


text:
cd(dirname(@__FILE__))
lines = open("seven.txt", "r") do f
    readlines(f)
end

parse_bag(c::String) =
    c[findfirst(isequal(' '), c) + 1:findlast(isequal(' '), c) - 1]

parse_line(l::String) =
    let (k,v) = split(l, "contain")
        bags = filter(!isempty, [c == " no other bags." ? "" : parse_bag(string(strip(c))) for c in split(v, ",")])
        (string(strip(k)[1:end-5]), Set(bags))
    end

d = Dict(map(parse_line, lines))

function recurse_bags(d, s::Set{String})
    if s == Set(["shiny gold"]) || isempty(s)
        s
    else
        recurse_bags(d, union((e == "shiny gold" ? Set([e]) : d[e] for e in s if !isempty(e))...))
    end
end

n = count(x -> x == Set(["shiny gold"]), map(x -> recurse_bags(d, x), values(d)))

println(n)