text:
@enum SortOrder ascending descending ignore

mutable struct FieldType
    name::String
    doc::Union{String, Nothing}
    type::Schema
    default::Any
    order::SortOrder
    aliases::Union{Vector{String}, Nothing}
end

==(x::FieldType, y::FieldType) = x.name == y.name && x.type == y.type
FieldType(name::String, type::Schema, doc=nothing, default=nothing, order=ascending, aliases=nothing) =
    FieldType(name, doc, type, default, order, aliases)
FieldType() = FieldType("", nothing, "", nothing, ascending, nothing)
StructTypes.StructType(::Type{FieldType}) = StructTypes.Mutable()
StructTypes.omitempties(::Type{FieldType}) = true

mutable struct RecordType <: NamedType
    type::String
    namespace::String
    name::String
    doc::Union{String, Nothing}
    aliases::Union{Vector{String}, Nothing}
    fields::Vector{FieldType}
end

==(x::RecordType, y::RecordType) = x.name == y.name && x.fields == y.fields
RecordType(name, fields::Vector{FieldType}, namespace="", doc=nothing, aliases=nothing) =
    RecordType("record", namespace, name, doc, aliases, fields)
RecordType() = RecordType("record", "", "", nothing, nothing, FieldType[])
StructTypes.StructType(::Type{RecordType}) = StructTypes.Mutable()

struct Record{names, T, N} <: Tables.AbstractRow
    type::RecordType
    data::Vector{UInt8}
    values::NTuple{N, Int} # starting byte positions for each field in `data`
end

type(x::Record) = getfield(x, :type)
data(x::Record) = getfield(x, :data)
values(x::Record) = getfield(x, :values)
Base.NamedTuple(x::Record{names, types, N}) where {names, types, N} =
    NamedTuple{names, types}(
        ntuple(i -> readvalue(Binary(), type(x).fields[i].type, fieldtype(types, i), data(x), values(x)[i], length(data(x)), (;))[1], Val(N))
    )

juliatype(x::RecordType) = Record{symtup(map(f->f.name, x.fields)), Tuple{map(f->juliatype(f.type), x.fields)...}, length(x.fields)}

autoname(T::Type{Record{names, types, N}}) where {names, types, N} =
    Base.string("Record_", hash(T))

schematype(::StructTypes.CustomStruct, ::Type{Record{names, types, N}}, name=autoname(Record{names, types, N})) where {names, types, N} =
    schematype(Record{names, types})
schematype(::StructTypes.CustomStruct, ::Type{Record{names, types}}, name=autoname(Record{names, types, fieldcount(types)})) where {names, types} =
    RecordType(name, FieldType[FieldType(String(names[i]), schematype(fieldtype(types, i))) for i = 1:length(names)])
schematype(::Type{NamedTuple{names, types}}) where {names, types} =
    schematype(StructTypes.CustomStruct(), Record{names, types})
schematype(::StructTypes.DataType, ::Type{T}) where {T} =
    schematype(StructTypes.CustomStruct(), Record{fieldnames(T), Tuple{fieldtypes(T)...}}, Base.string(T))
schematype(::StructTypes.CustomStruct, ::Type{T}) where {T} =
    schematype(StructTypes.lowertype(T))
schematype(::StructTypes.ArrayType, ::Type{T}) where {T <: Tuple} =
    schematype(NamedTuple{symtup(1:fieldcount(T)), T})

StructTypes.StructType(::Type{<:Record}) = StructTypes.CustomStruct()
StructTypes.lower(x::Record) = NamedTuple(x)

Tables.columnnames(::Record{names}) where {names} = names
Tables.getcolumn(x::Record{names, T}, ::Type{S}, i::Int, nm::Symbol) where {names, T, S} =
    readvalue(Binary(), S, data(x), values(x)[i], length(data(x)), (;))[1]
Tables.getcolumn(x::Record{names, T}, i::Int) where {names, T} =
    readvalue(Binary(), fieldtype(T, i), data(x), values(x)[i], length(data(x)), (;))[1]
function Tables.getcolumn(x::Record{names, T}, nm::Symbol) where {names, T}
    i = Tables.columnindex(names, nm)
    return readvalue(Binary(), fieldtype(T, i), data(x), values(x)[i], length(data(x)), (;))[1]
end

mutable struct WriteClosure{B, T, KW}
    type::RecordType
    buf::T
    pos::Int
    len::Int
    opts::KW
end

function (f::WriteClosure{B, T, KW})(i, nm, typ, v) where {B, T, KW}
    @debug 2 "writing $typ at pos = $(f.pos)"
    f.pos = writevalue(B(), f.type.fields[i].type, v, f.buf, f.pos, f.len, f.opts)
end

function writevalue(B::Binary, T::RecordType, x, buf, pos, len, opts)
    c = WriteClosure{Binary, typeof(buf), typeof(opts)}(T, buf, pos, len, opts)
    StructTypes.foreachfield(c, x)
    return c.pos
end

mutable struct NBytesClosure
    RT::RecordType
    n::Int
end

function (f::NBytesClosure)(i, nm, typ, v)
    f.n += nbytes(f.RT.fields[i].type, v)
end

function nbytes(RT::RecordType, x)
    c = NBytesClosure(RT, 0)
    StructTypes.foreachfield(c, x)
    return c.n
end

mutable struct ReadClosure{B, T, KW}
    type::RecordType
    buf::T
    pos::Int
    len::Int
    opts::KW
end

function (f::ReadClosure{B, T, KW})(i, nm, ::Type{S}) where {B, T, KW, S}
    x_i, pos = readvalue(B(), f.type.fields[i].type, S, f.buf, f.pos, f.len, f.opts)
    f.pos = pos
    return x_i
end

function readvalue(B::Binary, RT::RecordType, ::Type{T}, buf, pos, len, opts) where {T}
    c = ReadClosure{Binary, typeof(buf), typeof(opts)}(RT, buf, pos, len, opts)
    x = StructTypes.construct(c, T)
    return x, c.pos
end

mutable struct ReadSkipClosure{B, T, KW}
    type::RecordType
    buf::T
    pos::Int
    len::Int
    opts::KW
end

function (f::ReadSkipClosure{B, T, KW})(i, nm, ::Type{S}) where {B, T, KW, S}
    f.pos = skipvalue(B(), f.type.fields[i].type, S, f.buf, f.pos, f.len, f.opts)
    return
end

function skipvalue(B::Binary, RT::RecordType, ::Type{T}, buf, pos, len, opts) where {T}
    c = ReadSkipClosure{Binary, typeof(buf), typeof(opts)}(RT, buf, pos, len, opts)
    x = StructTypes.foreachfield(c, T)
    return c.pos
end

mutable struct PosClosure
    pos::Int
end

function (f::PosClosure)(RT, buf, len, opts, i, ::Type{S}) where {S}
    pos = f.pos
    f.pos = skipvalue(Binary(), RT.fields[i].type, S, buf, pos, len, opts)
    return pos
end

readvalue(B::Binary, RT::RecordType, ::Type{Record{names, types}}, buf, pos, len, opts) where {names, types} =
    readvalue(B, RT, Record{names, types, fieldcount(types)}, buf, pos, len, opts)

function readvalue(B::Binary, RT::RecordType, ::Type{Record{names, types, N}}, buf, pos, len, opts) where {names, types, N}
    c = PosClosure(pos)
    values = ntuple(i -> c(RT, buf, len, opts, i, fieldtype(types, i)), Val(N))
    return Record{names, types, N}(RT, buf, values), c.pos
end

skipvalue(B::Binary, RT::RecordType, ::Type{Record{names, types}}, buf, pos, len, opts) where {names, types} =
    skipvalue(B, RT, Record{names, types, fieldcount(types)}, buf, pos, len, opts)

function skipvalue(B::Binary, RT::RecordType, ::Type{Record{names, types, N}}, buf, pos, len, opts) where {names, types, N}
    c = PosClosure(pos)
    values = ntuple(i -> c(RT, buf, len, opts, i, fieldtype(types, i)), Val(fieldcount(types)))
    return c.pos
end


text:
include("support.jl")

function return_kth_to_last(ll :: LinkedList, k :: Int)
  forward :: Node = ll.head
  backward :: Node = ll.head

  # Make the forward pointer move k - 1 steps
  for _ = 1:(k-1)
    @assert forward.next != nothing "[return_kth_to_last] LinkedList has fewer than k nodes"
    forward = forward.next
  end

  # Move until the forward pointer reaches the end
  while forward.next != nothing
    forward = forward.next
    backward = backward.next
  end

  return backward.item
end

# Time  : O(n)
# Space : O(1)

using Base: Test

ll = LinkedList()
append_to_tail!(ll, 1)
append_to_tail!(ll, 2)
append_to_tail!(ll, 3)
append_to_tail!(ll, 4)
append_to_tail!(ll, 5)
append_to_tail!(ll, 6)
@test return_kth_to_last(ll, 4) == 3

ext: jl

max_stars_repo_name: Kiandr/crackingcodinginterview

max_issues_repo_name: Kiandr/crackingcodinginterview

max_forks_repo_name: Kiandr/crackingcodinginterview


text:
# models/potentials
abstract type AbstractPotential end
abstract type AbstractMixedModel <: AbstractPotential end
abstract type AbstractSinglePotential <: AbstractPotential end
abstract type AbstractTrainablePotential <: AbstractSinglePotential end

# systems
abstract type AbstractSystem end
abstract type AbstractThermostattedSystem <: AbstractSystem end
abstract type AbstractActiveSystem <: AbstractThermostattedSystem end

# thermostats
abstract type AbstractThermostat end
abstract type AbstractLangevin <: AbstractThermostat end

# integrators
abstract type AbstractIntegrator end

# callbacks
abstract type AbstractCallback end

# reinforcement learning types
abstract type AbstractValueBaseline end

# optimizer types
abstract type AbstractOptimizer end


text:
################################################################################
#
#  (Generalized) norm relations
#
################################################################################

# Example
#
# julia> f = x^8+8*x^6+64*x^4-192*x^2+576
# x^8+8*x^6+64*x^4-192*x^2+576
#
# julia> K, a = number_field(f);
#
# julia> N  = Hecke._norm_relation_setup(K)
# Norm relation of
#   Number field over Rational Field with defining polynomial x^8+8*x^6+64*x^4-192*x^2+576
# of index 4 and the subfields
#   Number field over Rational Field with defining polynomial x^2+1
#   Number field over Rational Field with defining polynomial x^2-6
#   Number field over Rational Field with defining polynomial x^2-2
#   Number field over Rational Field with defining polynomial x^2+2
#   Number field over Rational Field with defining polynomial x^2+6
#   Number field over Rational Field with defining polynomial x^2-3
#   Number field over Rational Field with defining polynomial x^2-x+1

mutable struct NormRelation{T}
  K::AnticNumberField
  subfields::Vector{Tuple{AnticNumberField, NfToNfMor}}
  denominator::Int
  coefficients::Vector{Vector{Tuple{NfToNfMor, Int}}}
  ispure::Bool
  pure_coefficients::Vector{Int}
  isabelian::Bool
  isgeneric::Bool
  coefficients_gen::Vector{Vector{Tuple{Int, NfToNfMor, NfToNfMor}}}
  embed_cache::Dict{Tuple{Int, Int}, Dict{nf_elem, nf_elem}}
  mor_cache::Dict{NfToNfMor, Dict{nf_elem, nf_elem}}
  induced::Dict{NfToNfMor, perm{Int}}

  function NormRelation{T}() where {T}
    z = new{T}()
    z.ispure = false
    z.isabelian = false
    z.isgeneric = false
    z.denominator = 0
    z.embed_cache = Dict{Tuple{Int, Int}, Dict{nf_elem, nf_elem}}()
    z.mor_cache = Dict{NfToNfMor, Dict{nf_elem, nf_elem}}()
    z.induced = Dict{NfToNfMor, perm{Int}}()
    return z
  end
end

function Base.show(io::IO, N::NormRelation)
  print(io, "Norm relation of\n  ", N.K, "\nof index ", index(N), " and the subfields")
  for i in 1:length(N)
    print(io, "\n  ", subfields(N)[i][1])
  end
end

function Base.hash(x::AlgGrpElem, h::UInt)
  return Base.hash(x.coeffs, h)
end

function _norm_relation_setup_abelian(K::AnticNumberField; small_degree::Bool = true, pure::Bool = true, index::fmpz = zero(fmpz))
  G = automorphisms(K)
  A, GtoA, AtoG = find_isomorphism_with_abelian_group(G, *);
  if iszero(index)
    subs = [f for f in subgroups(A) if order(f[1]) > 1]
  else
    subs = [f for f in subgroups(A) if order(f[1]) > 1 && order(f[1]) == index || order(f[1]) == index^2]
  end

  b, den, ls = _has_norm_relation_abstract(A, subs, large_index = !small_degree, pure = pure)
  @assert b
  n = length(ls)

  z = NormRelation{Int}()
  z.K = K
  z.subfields = Vector{Tuple{AnticNumberField, NfToNfMor}}(undef, n)
  z.denominator = den
  z.ispure = false

  println("Computing the subfields ... ")

  for i in 1:n
    println("$i/$n")
    @time F, mF = fixed_field(K, NfToNfMor[AtoG[f] for f in ls[i][2]])
    println("I start simplifying")
    @time maximal_order(F)
    @time S, mS = simplify(F)
    println("DONE")
    L = S
    mL = mS * mF
    z.subfields[i] = L, mL
  end

  if pure
    z.ispure = true
    z.pure_coefficients = Vector{Int}(undef, n)
    for i in 1:n
      z.pure_coefficients[i] = ls[i][1]
    end
  else
    z.coefficients = Vector{Vector{Tuple{NfToNfMor, Int}}}(undef, n)
    for i in 1:n
      w = Vector{Tuple{NfToNfMor, Int}}(undef, length(ls[i][1]))
      for (j, (expo, auto)) in enumerate(ls[i][1])
        w[j] = (AtoG[auto], expo)
      end
      z.coefficients[i] = w
    end

  end

  z.isabelian = true

  return z
end

function _norm_relation_setup_generic(K::AnticNumberField; small_degree::Bool = true, target_den::fmpz = zero(fmpz))
  println("Compute automorphisms ...")
  A = automorphisms(K)
  println("Compute group table ...")
  G, AtoG, GtoA = generic_group(A, *)
  println("Compute the abstract norm relation")
  if iszero(target_den)
    b, den, ls = _has_norm_relation_abstract(G, [f for f in subgroups(G, conjugacy_classes = true) if order(f[1]) > 1])
  else
    b, den, ls = _has_norm_relation_abstract(G, [f for f in subgroups(G, conjugacy_classes = true) if order(f[1]) > 1], target_den = target_den)
  end
  println("done")
  @assert b
  n = length(ls)

  z = NormRelation{Int}()
  z.K = K
  z.subfields = Vector{Tuple{AnticNumberField, NfToNfMor}}(undef, n)
  z.denominator = den
  z.ispure = false

  for i in 1:n
    println("Computing fixed fields ...")
    F, mF = fixed_field(K, NfToNfMor[GtoA[f] for f in ls[i][1]])
    println("Simplifying ...")
    S, mS = simplify(F, cached = true)
    L = S
    println("done")
    mL = mS * mF
    z.subfields[i] = L, mL
  end

  z.coefficients_gen = Vector{Vector{Tuple{Int, NfToNfMor, NfToNfMor}}}(undef, n)

  for i in 1:n
    w = Vector{Tuple{Int, NfToNfMor, NfToNfMor}}(undef, length(ls[i][2]))
    for (j, (expo, auto_pre, auto_post)) in enumerate(ls[i][2])
      w[j] = (expo, GtoA[auto_pre], GtoA[auto_post])
    end
    z.coefficients_gen[i] = w
  end

  z.isabelian = false

  return z
end

Base.length(N::NormRelation) = length(N.subfields)

field(N::NormRelation) = N.K

subfields(N::NormRelation) = N.subfields

subfield(N::NormRelation, i::Int) = N.subfields[i]

embedding(N::NormRelation, i::Int) = N.subfields[i][2]

index(N::NormRelation) = N.denominator

norm(N::NormRelation, i::Int, a) = norm(N.embeddings[i], a)

ispure(N::NormRelation) = N.ispure

function image(f::NfToNfMor, a::FacElem{nf_elem, AnticNumberField})
  D = Dict{nf_elem, fmpz}(f(b) => e for (b, e) in a)
  return FacElem(D)
end

function check_relation(N::NormRelation, a::nf_elem)
  @assert !iszero(a)
  if ispure(N)
    z = one(field(N))
    for i in 1:length(N)
      z = z * N(norm(embedding(N, i), a), i)
    end
    return a^index(N) == z
  else
    # I applied the exponent when embedding
    z = one(field(N))
    for i in 1:length(N)
      z = z * N([ norm(embedding(N, i), auto(a)^exp) for  (exp, auto, _) in N.coefficients_gen[i]], i)
    end
    return a^index(N) == z
  end
end

function (N::NormRelation)(x::Union{nf_elem, FacElem{nf_elem, AnticNumberField}}, i::Int)
  z = one(N.K)
  _, mk = subfield(N, i)
  y = mk(x)
  if ispure(N)
    return y^N.pure_coefficients[i]
  end
end

function _apply_auto(N::NormRelation, auto::NfToNfMor, x::nf_elem)
  if haskey(N.mor_cache, auto)
    d = N.mor_cache[auto]
    if haskey(d, x)
      return d[x]
    else
      y = auto(x)
      d[x] = y
      return y
    end
  else
    d = Dict{nf_elem, nf_elem}()
    y = auto(x)
    d[x] = y
    N.mor_cache[auto] = d
    return y
  end
end

function (N::NormRelation)(x::Union{nf_elem, FacElem{nf_elem, AnticNumberField}}, i::Int, j::Int)
  @assert (1 <= j) && (j <= length(N.coefficients_gen[i]))
  if haskey(N.embed_cache, (i, j)) && haskey(N.embed_cache[i, j], x)
    return N.embed_cache[i, j][x]
  end

  _, mk = subfield(N, i)
  y = _apply_auto(N, mk, x)

  exp, _, auto = N.coefficients_gen[i][j]

  z = _apply_auto(N, auto, y)
  return z
end

function (N::NormRelation)(x::Vector{<:Union{nf_elem, FacElem{nf_elem, AnticNumberField}}}, i::Int)
  @assert length(x) == length(N.coefficients_gen[i])
  z = one(N.K)
  for j in 1:length(N.coefficients_gen[i])
    z = z * N(x[j], i, j)
  end
  return z
end

function induce_action(N::NormRelation, i, j, s, FB, cache)
  S = FB.ideals
  ZK = order(S[1])
  z = zero_matrix(SMat, FlintZZ, 0, length(S))
  mk = embedding(N, i)
  zk = order(s[1])

  _ , _, auto = N.coefficients_gen[i][j]
  #@show auto
  println("Call to induce(FB, auto)...")
  @time if haskey(N.induced, auto)
    p = N.induced[auto]
  else
    p = induce(FB, auto)
    N.induced[auto] = p
  end
  #@show p

  if length(cache) == 0
    cache = resize!(cache, length(s))
    cached = false
  else
    cached = true
  end

  @show cached

  @assert mod(degree(ZK), degree(zk)) == 0
  reldeg = divexact(degree(ZK), degree(zk))

  for l in 1:length(s)
    if cached
      v = cache[l]
    else
      v = Tuple{Int, fmpz}[]
      P = s[l]
      genofsl = elem_in_nf(_apply_auto(N, mk, P.gen_two.elem_in_nf))
      pmin = minimum(P, copy = false)
      # Compute the number of ideals
      # Assume that L/K and L/Q are normal
      rele = divexact(splitting_type((FB.fb[pmin].lp)[1][2])[1], splitting_type(P)[1])
      relf = divexact(splitting_type((FB.fb[pmin].lp)[1][2])[2], splitting_type(P)[2])
      @assert mod(reldeg, rele * relf) == 0
      numb_ideal = divexact(reldeg, rele * relf)
      found = 0
      for k in 1:length(S)
        Q = S[k]
        if minimum(Q, copy = false) == pmin
          ant = anti_uniformizer(Q)
          if (genofsl * ant) in ZK
            found += 1
            @assert mod(splitting_type(Q)[1], splitting_type(s[l])[1]) == 0
            push!(v, (k, divexact(splitting_type(Q)[1], splitting_type(s[l])[1])))
          end
        end
        if found == numb_ideal
          break
        end
      end
      cache[l] = v
    end

    # We still need to sort the positions of the non-zero entries
    #@show l, v
    sort!(v, by = x -> x[1])
    #@show l, v
    # Now apply the automorphism
    r = [(p[i], e) for (i, e) in v]
    sort!(r, lt = (a,b)->a[1]<b[1])
    #@show l, r
    push!(z, sparse_row(FlintZZ, r))
  end

  #for l in 1:length(s)
  #  @show l
  #  a = uniformizer(s[l])
  #  @show [ valuation(a, Q) for Q in s]
  #  @show [ valuation(auto(mk(a.elem_in_nf)), Q) for Q in FB.ideals]
  #  @show [ z.rows[l][m] for m in 1:length(FB.ideals) ]
  #  @assert [ valuation(auto(mk(a.elem_in_nf)), Q) for Q in FB.ideals] == [ z.rows[l][m] for m in 1:length(FB.ideals) ]
  #end
  @show z
  return z
end

function induce_action(N::NormRelation, i, s::Vector, S::Vector)
  if !ispure(N)
    throw(error("Not implemented yet"))
  end
  ZK = order(S[1])
  z = zero_matrix(SMat, FlintZZ, 0, length(S))
  mk = embedding(N, i)
  for j in 1:length(s)
    v = Tuple{Int, fmpz}[]
    for k in 1:length(S)
      Q = S[k]
      if minimum(Q, copy = false) == minimum(s[j], copy = false)
        ant = anti_uniformizer(Q)
        if (elem_in_nf(mk(s[j].gen_two.elem_in_nf)) * ant) in ZK
          push!(v, (k, divexact(splitting_type(Q)[1], splitting_type(s[j])[1])))
        end
      end
    end

    # We still need to sort the positions of the non-zero entries
    sort!(v, by = x -> x[1])

    push!(z, N.pure_coefficients[i] * sparse_row(FlintZZ, v))
  end
  return z
end

function _setup_for_norm_relation_fun(K, S = prime_ideals_up_to(maximal_order(K), factor_base_bound_grh(maximal_order(K))))
  ZK = order(S[1])
  FB = NfFactorBase(ZK, S)
  c = class_group_init(FB)
  UZK = Hecke.UnitGrpCtx{FacElem{nf_elem, AnticNumberField}}(ZK)
  return c, UZK
end

function _add_sunits_from_norm_relation!(c, N)
  cp = sort!(collect(Set(minimum(x) for x = c.FB.ideals)))
  K = N.K
  for i = 1:length(N)
    k, mk = subfield(N, i)
    println("Computing maximal order ...")
    @time zk = maximal_order(k)
    println("Computing lll basis ... ")
    @show k
    @time zk = lll(zk)
    print("Computing class group of $k... ")
    @time class_group(zk, redo = true, use_aut = true)
    println("done")
    lpk = NfOrdIdl[ P[1] for p in cp for P = prime_decomposition(zk, p)]
    println("Now computing the S-unit group for lp of length $(length(lpk))")
    @assert length(lpk) > 0
    @time Szk, mS = Hecke.sunit_mod_units_group_fac_elem(lpk)
    @show length(N.coefficients_gen[i])

    if ispure(N)
      @time z = induce_action(N, i, lpk, c.FB.ideals) # Careful, don't use S instead of FB.ideals

      # Embedding is expensive, so let us build a cache
      D = Dict{nf_elem, nf_elem}()
      function N_mk(x, D, i)
        if haskey(D, x)
          return D[x]
        else
          y = N(x, i)
          D[x] = y
          return y
        end
      end

      print("Feeding in the S-units of the small field ... ")

      for j=1:ngens(Szk)
        print(j, " ")
        u = mS(Szk[j])  #do compact rep here???
        #@show z.rows
        valofnewelement = mul(mS.valuations[j], z)
        #v = zeros(fmpz, length(FB.ideals))
        #v[valofnewelement.pos] = valofnewelement.values
        #@assert [valuation((FacElem(Dict((N_mk(x, D, i), v) for (x,v) = u.fac))), P) for P in FB.ideals] == v
        Hecke.class_group_add_relation(c, FacElem(Dict((N_mk(x, D, i), v) for (x,v) = u.fac)), valofnewelement)
      end
    else
      #D = Dict{nf_elem, nf_elem}()
      cc = Vector{Tuple{Int, fmpz}}[]
      for j in 1:length(N.coefficients_gen[i])
        println("Inducing the action ... ")
        @time z = induce_action(N, i, j, lpk, c.FB, cc)

        print("Feeding in the S-units of the small field ... ")

        for l=1:ngens(Szk)
          print(l, " ")
          u = mS(Szk[l])  #do compact rep here???
          #@show z.rows
          valofnewelement = mul(mS.valuations[l], z)
          #v = zeros(fmpz, length(c.FB.ideals))
          #@show valofnewelement
          #v[valofnewelement.pos] = valofnewelement.values
          #@assert [valuation((FacElem(Dict((N(x, i, j), v) for (x,v) = u.fac))), P) for P in c.FB.ideals] == v
          Hecke.class_group_add_relation(c, FacElem(Dict((N(x, i, j), v) for (x,v) = u.fac)), valofnewelement)
          #Hecke.class_group_add_relation(c, FacElem(Dict(K(1) => fmpz(1))), valofnewelement)
        end
        println("")
      end
    end

    println("done")

    # Skipping the units
    #U, mU = unit_group_fac_elem(zk)
    #for j=2:ngens(U) # I cannot add a torsion unit. He will hang forever.
    #  u = mU(U[j])
    #  Hecke._add_unit(UZK, FacElem(Dict((N_mk(x, D, i), v) for (x,v) = u.fac)))
    #end
    #UZK.units = Hecke.reduce(UZK.units, UZK.tors_prec)
  end

  #UZK.tentative_regulator = Hecke.regulator(UZK.units, 64)
end

function _compute_sunit_group_mod(K::AnticNumberField, N::NormRelation, S)
  ZK = order(S[1])
  FB = NfFactorBase(ZK, S)
  c = class_group_init(FB)
  cp = sort!(collect(Set(minimum(x) for x = c.FB.ideals)))

  UZK = Hecke.UnitGrpCtx{FacElem{nf_elem, AnticNumberField}}(ZK)

  for i = 1:length(N)
    k, mk = subfield(N, i)
    zk = maximal_order(k)
    print("Computing class group of $k... ")
    class_group(zk, use_aut = true)
    println("done")
    lpk = [ P[1] for p in cp for P = prime_decomposition(zk, p)]
    println("Now computing the S-unit group for lp of length $(length(lpk))")
    @assert length(lpk) > 0
    Szk, mS = Hecke.sunit_mod_units_group_fac_elem(lpk)
    z = induce_action(N, i, lpk, FB.ideals) # Careful, don't use S instead of FB.ideals

    # Embedding is expensive, so let us build a cache
    D = Dict{nf_elem, nf_elem}()
    function N_mk(x, D, i)
      if haskey(D, x)
        return D[x]
      else
        y = N(x, i)
        D[x] = y
        return y
      end
    end

    print("Feeding in the S-units of the small field ... ")

    for j=1:ngens(Szk)
      print(j, " ")
      u = mS(Szk[j])  #do compact rep here???
      #@show z.rows
      valofnewelement = mul(mS.valuations[j], z)
      #v = zeros(fmpz, length(FB.ideals))
      #v[valofnewelement.pos] = valofnewelement.values
      #@assert [valuation((FacElem(Dict((N_mk(x, D, i), v) for (x,v) = u.fac))), P) for P in FB.ideals] == v
      Hecke.class_group_add_relation(c, FacElem(Dict((N_mk(x, D, i), v) for (x,v) = u.fac)), valofnewelement)
    end

    println("done")

    U, mU = unit_group_fac_elem(zk)
    for j=2:ngens(U) # I cannot add a torsion unit. He will hang forever.
      u = mU(U[j])
      Hecke._add_unit(UZK, FacElem(Dict((N_mk(x, D, i), v) for (x,v) = u.fac)))
    end
    UZK.units = Hecke.reduce(UZK.units, UZK.tors_prec)
  end

  UZK.tentative_regulator = Hecke.regulator(UZK.units, 64)

  println("Now saturating ... ")

  for (p, e) in factor(index(N))
    @show p
    b = Hecke.saturate!(c, UZK, Int(p))
    while b
      @show "SATURATING AGAIN"
      b = Hecke.saturate!(c, UZK, Int(p))
    end
  end
  return c, UZK
end

one(T::FacElemMon{AnticNumberField}) = T()

function simplify(c::Hecke.ClassGrpCtx)
  d = Hecke.class_group_init(c.FB, SMat{fmpz}, add_rels = false)
  U = Hecke.UnitGrpCtx{FacElem{nf_elem, AnticNumberField}}(order(d))

  Hecke.module_trafo_assure(c.M)
  trafos = c.M.trafo

  for i=1:length(c.FB.ideals)
    x = zeros(fmpz, length(c.R_gen) + length(c.R_rel))
    x[i] = 1
    for j in length(trafos):-1:1
      Hecke.apply_right!(x, trafos[j])
    end
    y = FacElem(vcat(c.R_gen, c.R_rel), x)
    fl = Hecke.class_group_add_relation(d, y, deepcopy(c.M.basis.rows[i]))
    @assert fl
  end
  for i=1:nrows(c.M.rel_gens)
    if iszero(c.M.rel_gens.rows[i])
      Hecke._add_unit(U, c.R_rel[i])
    end
  end
  for i=1:length(U.units)
    Hecke.class_group_add_relation(d, U.units[i], SRow(FlintZZ))
  end
  return d, U
end

function units(c::Hecke.ClassGrpCtx)
  d = Hecke.class_group_init(c.FB, SMat{fmpz}, add_rels = false)
  U = Hecke.UnitGrpCtx{FacElem{nf_elem, AnticNumberField}}(order(d))

  Hecke.module_trafo_assure(c.M)
  trafos = c.M.trafo

  for i=1:nrows(c.M.rel_gens)
    if iszero(c.M.rel_gens.rows[i])
      Hecke._add_unit(U, c.R_rel[i])
    end
  end

  U.units = Hecke.reduce(U.units, U.tors_prec)
  U.tentative_regulator = Hecke.regulator(U.units, 64)

  return U
end

################################################################################
#
#  Code for abstract norm relations
#
################################################################################

# TODO:
# This is a mess.
# We first have to figure out the different relation types.
function _has_norm_relation_abstract(G::GrpAb,
                                     H::Vector{Tuple{GrpAbFinGen, GrpAbFinGenMap}};
                                     primitive::Bool = false,
                                     interactive::Bool = false,
                                     greedy::Bool = false,
                                     large_index::Bool = false,
                                     pure::Bool = false)
  QG = AlgGrp(FlintQQ, G)
  wd = decompose(QG)
  idempotents = [ f(one(A)) for (A, f) in wd]
  norms_rev = Dict{elem_type(QG), Int}()
  norms = Vector{elem_type(QG)}(undef, length(H))
  for i in 1:length(H)
    norms_rev[sum([QG(H[i][2](h)) for h in H[i][1]])] = i
    norms[i] = sum([QG(H[i][2](h)) for h in H[i][1]])
  end

  nonannihilating = Vector{Vector{Int}}(undef, length(idempotents))

  for j in 1:length(idempotents)
    nonannihilating_ej = Vector{Int}()
    for i in 1:length(norms)
      if !iszero(idempotents[j] * norms[i])
        push!(nonannihilating_ej, i)
      end
    end
    nonannihilating[j] = nonannihilating_ej
  end

  l = length(idempotents)

  subgroups_needed = Int[]

  if any(isempty, nonannihilating)
    return false, zero(FlintZZ), Vector{Tuple{Vector{Tuple{fmpz, GrpAbFinGenElem}}, Vector{GrpAbFinGenElem}}}()
  end

  for i in 1:length(nonannihilating)
    if greedy || !large_index
      sort!(nonannihilating[i], by = x -> order(H[x][1]), rev = true)
    elseif large_index
      sort!(nonannihilating[i], by = x -> order(H[x][1]))
    end
  end

  if interactive
    println("Which subgroups to choose?:")
    string_read = readline(stdin)
    subgroups_needed = eval(Meta.parse(string_read))
    if any(x -> length(intersect(subgroups_needed, x)) == 0, nonannihilating)
      error("They don't give a relation")
    end
  else
    k = nonannihilating[1][1]
    push!(subgroups_needed, k)
    i = 2
    for i in 2:l
      if length(intersect(subgroups_needed, nonannihilating[i])) > 0
        continue
      else
        push!(subgroups_needed, nonannihilating[i][1])
      end
    end
  end

  if primitive
    # Check if we can replace a subgroup by bigger subgroups
    for i in 1:length(subgroups_needed)
      Q, mQ = quo(G, H[subgroups_needed[i]][1])
      if length(Q) == 1
        continue
      end
      b, den, ls = _has_norm_relation_abstract(Q, [f for f in subgroups(Q) if order(f[1]) > 1])
      if b
        return _has_norm_relation_abstract(G, deleteat!(H, subgroups_needed[i]))
      end
    end
  end

  cyc = [h for h in H]# if iscyclic(h[1])]

  if pure
    n = Int(order(G))
    m = zero_matrix(FlintQQ, length(cyc), n)
    for i in 1:length(cyc)
      for j in 1:n
        m[i, j] = norms[i].coeffs[j]
      end
    end

    onee = matrix(FlintQQ, 1, n, coeffs(one(QG)))

    b, v, K = can_solve_with_kernel(m', onee')

    if !b
      return false, zero(FlintZZ), Vector{Tuple{Vector{Tuple{fmpz, GrpAbFinGenElem}}, Vector{GrpAbFinGenElem}}}()
    end

    subgroups_needed = Int[ i for i in 1:nrows(v) if !iszero(v[i, 1])]

    den = one(FlintZZ)

    for i in subgroups_needed
      den = lcm(den, denominator(v[i, 1]))
    end

    solutions_as_group_element_arrays = Vector{Tuple{fmpz, Vector{GrpAbFinGenElem}}}(undef, length(subgroups_needed))

    for i in 1:length(subgroups_needed)
    #  vv = Tuple{fmpz, GrpAbFinGenElem}[]
    #  for (j, g) in enumerate(group(QG))
    #    if !iszero(sol[i].coeffs[j])
    #      push!(vv, (FlintZZ(den*sol[i].coeffs[j]), g))
    #    end
    #  end
      solutions_as_group_element_arrays[i] = (FlintZZ(den * v[subgroups_needed[i], 1]), [H[subgroups_needed[i]][2](h) for h in H[subgroups_needed[i]][1]])
    end

    @assert b
    return true, den, solutions_as_group_element_arrays

  end

  # Now solve the equation to find a representation 1 = sum_{H} x_H N_H
  n = Int(order(G))

  onee = matrix(FlintQQ, 1, n, coeffs(one(QG)))

  m = zero_matrix(FlintQQ, length(subgroups_needed) * n, n)

  B = basis(QG)

  for i in 1:length(subgroups_needed)
    for j in 1:n
      for k in 1:n
        m[(i - 1)*n + k, j] = (B[k] * norms[subgroups_needed[i]]).coeffs[j]
      end
    end
  end

  b, v, K = can_solve_with_kernel(m', onee')

  # Now collect the coefficients again as elements of Q[G]
  sol = Vector{elem_type(QG)}(undef, length(subgroups_needed))
  for i in 1:length(subgroups_needed)
    sol[i] = QG([v[k, 1] for k in ((i - 1)*n + 1):i*n])
  end

  @assert one(QG) == sum(sol[i] * norms[subgroups_needed[i]] for i in 1:length(subgroups_needed))

  den = one(FlintZZ)

  for i in 1:length(subgroups_needed)
    for j in 1:n
      den = lcm(den, denominator(sol[i].coeffs[j]))
    end
  end

  solutions_as_group_element_arrays = Vector{Tuple{Vector{Tuple{fmpz, GrpAbFinGenElem}}, Vector{GrpAbFinGenElem}}}(undef, length(subgroups_needed))

  for i in 1:length(subgroups_needed)
    vv = Tuple{fmpz, GrpAbFinGenElem}[]
    for (j, g) in enumerate(group(QG))
      if !iszero(sol[i].coeffs[j])
        push!(vv, (FlintZZ(den*sol[i].coeffs[j]), g))
      end
    end
    solutions_as_group_element_arrays[i] = (vv, [H[subgroups_needed[i]][2](h) for h in H[subgroups_needed[i]][1]])
  end

  return true, den, solutions_as_group_element_arrays
end

# TODO: Make this great again
function _has_norm_relation_abstract(G::GrpGen, H::Vector{Tuple{GrpGen, GrpGenToGrpGenMor}};
                                     primitive::Bool = false,
                                     greedy::Bool = false,
                                     large_index::Bool = false,
                                     pure::Bool = false,
                                     target_den::fmpz = zero(fmpz))
  QG = AlgGrp(FlintQQ, G)
  wd = decompose(QG)
  idempotents = [ f(one(A)) for (A, f) in wd]
  norms_rev = Dict{elem_type(QG), Int}()
  norms = Vector{elem_type(QG)}(undef, length(H))
  for i in 1:length(H)
    norms_rev[sum([QG(H[i][2](h)) for h in H[i][1]])] = i
    norms[i] = sum([QG(H[i][2](h)) for h in H[i][1]])
  end

  n = Int(order(G))

  if pure
    m = zero_matrix(FlintQQ, length(H), n)
    for i in 1:length(H)
      for j in 1:n
        m[i, j] = norms[i].coeffs[j]
      end
    end

    onee = matrix(FlintQQ, 1, n, coeffs(one(QG)))

    b, v, K = can_solve_with_kernel(m', onee')

    return b
  end

  nonannihilating = Vector{Vector{Int}}(undef, length(idempotents))

  for j in 1:length(idempotents)
    nonannihilating_ej = Vector{Int}()
    for i in 1:length(norms)
      if !iszero(idempotents[j] * norms[i])
        push!(nonannihilating_ej, i)
      end
    end
    nonannihilating[j] = nonannihilating_ej
  end

  l = length(idempotents)

  subgroups_needed = Int[]

  if any(isempty, nonannihilating)
    return false, zero(FlintZZ), Vector{Tuple{Vector{Tuple{fmpz, GrpAbFinGenElem}}, Vector{GrpAbFinGenElem}}}()
  end

  subgroups_needed = Int[]

  for i in 1:length(nonannihilating)
    if greedy || !large_index
      sort!(nonannihilating[i], by = x -> order(H[x][1]), rev = true)
    elseif large_index
      sort!(nonannihilating[i], by = x -> order(H[x][1]))
    end
  end

  k = nonannihilating[1][1]
  push!(subgroups_needed, k)
  i = 2
  for i in 2:l
    if length(intersect(subgroups_needed, nonannihilating[i])) > 0
      continue
    else
      push!(subgroups_needed, nonannihilating[i][1])
    end
  end

  if !iszero(target_den)
    subgroups_needed = collect(1:length(H))
  end

  if primitive
    # Check if we can replace a subgroup by bigger subgroups
    for i in 1:length(subgroups_needed)
      Q, mQ = quo(G, H[subgroups_needed[i]][1])
      if length(Q) == 1
        continue
      end
      b, den, ls = _has_norm_relation_abstract(Q, [f for f in subgroups(Q) if order(f[1]) > 1])
      if b
        return _has_norm_relation_abstract(G, deleteat!(H, subgroups_needed[i]))
      end
    end
  end

  # Now solve the equation to find a representation 1 = sum_{H} x_H N_H y_H
  n = Int(order(G))

  # We use the following trick to reduce the number of equations
  # Instead of working with the Q-basis (g N_H h)_(H, g in G, h in H), we take
  # (g N_H h)_(H, g in G/H, h in G\N_G(H))

  left_cosets_for_sub = Vector{GrpGenElem}[]
  right_cosets_for_normalizer = Vector{GrpGenElem}[]

  for i in 1:length(subgroups_needed)
    lc = left_cosets(G, H[subgroups_needed[i]][2])
    push!(left_cosets_for_sub, lc)
    NGH = normalizer(G, H[subgroups_needed[i]][2])
    push!(right_cosets_for_normalizer, right_cosets(G, NGH[2]))
  end

  if iszero(target_den)
    onee = matrix(FlintQQ, 1, n, coeffs(one(QG)))

    m = zero_matrix(FlintQQ, dot(length.(left_cosets_for_sub), length.(right_cosets_for_normalizer)), n)

    B = basis(QG)

    help_with_indicies = Vector{Tuple{Int, Int, Int}}(undef, nrows(m))
    cc = 1
    for i in 1:length(subgroups_needed)
      lcosub = length(left_cosets_for_sub[i])
      lconor = length(right_cosets_for_normalizer[i])
      for (k, g) in enumerate(left_cosets_for_sub[i])
        for (l, h) in enumerate(right_cosets_for_normalizer[i])
          for j in 1:n
            m[cc, j] = (QG(g) * norms[subgroups_needed[i]] * QG(h)).coeffs[j]
          end
          help_with_indicies[cc] = (i, k, l)
          cc += 1
        end
      end
    end

    b, w, K = can_solve_with_kernel(m, onee, side = :left)
    v = w
  elseif true
    onee = matrix(FlintZZ, 1, n, coeffs(target_den * one(QG)))

    m = zero_matrix(FlintZZ, dot(length.(left_cosets_for_sub), length.(right_cosets_for_normalizer)), n)

    help_with_indicies = Vector{Tuple{Int, Int, Int}}(undef, nrows(m))
    cc = 1
    for i in 1:length(subgroups_needed)
      for (k, g) in enumerate(left_cosets_for_sub[i])
        for (l, h) in enumerate(right_cosets_for_normalizer[i])
          for j in 1:n
            m[cc, j] = FlintZZ((QG(g) * norms[subgroups_needed[i]] * QG(h)).coeffs[j])
          end
          help_with_indicies[cc] = (i, k, l)
          cc += 1
        end
      end
    end

    b, w, K = can_solve_with_kernel(m, onee, side = :left)

    v = 1//target_den * change_base_ring(w, FlintQQ)
  end

  @assert b

  @show v

  z = zero(QG)
  for cc in 1:ncols(v)
    if iszero(v[1, cc])
      continue
    else
      i, k, l = help_with_indicies[cc]
      z = z + v[1, cc] *  QG(left_cosets_for_sub[i][k]) * norms[subgroups_needed[i]] * QG(right_cosets_for_normalizer[i][l])
    end
  end

  @assert isone(z)

  den = one(FlintZZ)

  for i in 1:ncols(v)
    den = lcm(den, denominator(v[1, i]))
  end

  solutions = Vector{Tuple{Vector{GrpGenElem}, Vector{Tuple{fmpz, GrpGenElem, GrpGenElem}}}}()
  solutions_as_group_element_arrays = Vector{Tuple{Vector{Tuple{fmpz, GrpAbFinGenElem}}, Vector{GrpAbFinGenElem}}}(undef, length(subgroups_needed))

  vvv = Vector{Vector{Tuple{fmpz, GrpGenElem, GrpGenElem}}}(undef, length(subgroups_needed))
  for i in 1:length(vvv)
    vvv[i] = Vector{Tuple{fmpz, GrpGenElem, GrpGenElem}}()
  end

  for cc in 1:ncols(v)
    if iszero(v[1, cc])
      continue
    else
      i, k, l = help_with_indicies[cc]
      push!(vvv[i], (FlintZZ(den * v[1, cc]), left_cosets_for_sub[i][k], right_cosets_for_normalizer[i][l]))
    end
  end

  for i in 1:length(vvv)
    sgroup = [H[subgroups_needed[i]][2](h) for h in H[subgroups_needed[i]][1]]
    if length(vvv[i]) > 0
      push!(solutions, (sgroup, vvv[i]))
    end
  end

  #for i in 1:length(subgroups_needed)
  #  sgroup = [H[subgroups_needed[i]][2](h) for h in H[subgroups_needed[i]][1]]
  #  vv = Vector{Tuple{fmpz, GrpGenElem, GrpGenElem}}()
  #  for (k, g) in enumerate(group(QG))
  #    for (l, h) in enumerate(group(QG))
  #      if iszero(v[1, (i - 1)*n^2 + (k - 1) * n + l])
  #        continue
  #      else
  #        push!(vv, (FlintZZ(den * v[1, (i - 1)*n^2 + (k - 1) * n + l]), g, h))
  #      end
  #    end
  #  end
  #  if length(vv) == 0
  #    continue
  #  end
  #  push!(solutions, (sgroup, vv))
  #end

  if !iszero(target_den)
    @assert iszero(mod(target_den, den))
  end

  return true, den, solutions
end


text:
"""
    ZeroRange(n)

Defines an `AbstractUnitRange` corresponding to the half-open interval
[0,n), equivalent to `0:n-1` but with the lower bound guaranteed to be
zero by Julia's type system.
"""
immutable ZeroRange{T<:Integer} <: AbstractUnitRange{T}
    len::T
    ZeroRange(len) = new(max(zero(T), len))
end
ZeroRange{T<:Integer}(len::T) = ZeroRange{T}(len)

Base.length(r::ZeroRange) = r.len

Base.length{T<:Union{Int,Int64}}(r::ZeroRange{T}) = T(r.len)

let smallint = (Int === Int64 ?
                Union{Int8,UInt8,Int16,UInt16,Int32,UInt32} :
                Union{Int8,UInt8,Int16,UInt16})
    Base.length{T <: smallint}(r::ZeroRange{T}) = Int(r.len)
    Base.start{T<:smallint}(r::ZeroRange{T}) = 0
end

Base.first{T}(r::ZeroRange{T}) = zero(T)
Base.last{ T}(r::ZeroRange{T}) = r.len-one(T)

Base.start{T}(r::ZeroRange{T}) = oftype(one(T)+one(T), zero(T))
Base.done{T}(r::ZeroRange{T}, i) = i == oftype(i, r.len)

@inline function Base.getindex{T}(v::ZeroRange{T}, i::Integer)
    @boundscheck ((i > 0) & (i <= length(v))) || Base.throw_boundserror(v, i)
    convert(T, i-1)
end

@inline function Base.getindex{R,S<:Integer}(r::ZeroRange{R}, s::AbstractUnitRange{S})
    @boundscheck checkbounds(r, s)
    R(first(s)-1):R(last(s)-1)
end

Base.intersect(r::ZeroRange, s::ZeroRange) = ZeroRange(min(r.len,s.len))

Base.promote_rule{T1,T2}(::Type{ZeroRange{T1}},::Type{ZeroRange{T2}}) =
    ZeroRange{promote_type(T1,T2)}
Base.convert{T<:Real}(::Type{ZeroRange{T}}, r::ZeroRange{T}) = r
Base.convert{T<:Real}(::Type{ZeroRange{T}}, r::ZeroRange) = ZeroRange{T}(r.len)

Base.show(io::IO, r::ZeroRange) = print(io, typeof(r).name, "(", r.len, ")")


text:
abstract type AbstractTransform end
abstract type AbstractCompose <: AbstractTransform end

const EPS = Float32(1e-6)

function Base.show(io::IO, c::AbstractCompose)
    print(io, "$(typeof(c))(")
    join(io, c.transforms, ", ")
    print(io, ")")
end

Base.show(io::IO, t::AbstractTransform) = print(io, "$(typeof(t))(...)")


text:
module LazyArrays
using Base, Base.Broadcast, LinearAlgebra, FillArrays, StaticArrays
import LinearAlgebra.BLAS

import Base: ReinterpretArray, ReshapedArray, AbstractCartesianIndex, Slice,
             RangeIndex, BroadcastStyle, copyto!, length, broadcastable, axes,
             getindex, eltype, tail, IndexStyle, IndexLinear,
             *, +, -, /, \, ==, isinf, isfinite, sign, angle, show, isless,
         fld, cld, div, min, max, minimum, maximum, mod,
         <, ≤, >, ≥, promote_rule, convert, copy,
         size, step, isempty, length, first, last, ndims,
         getindex, setindex!, OneTo, intersect, @_inline_meta, inv,
         sort, sort!, issorted, sortperm, cumsum, sum, in, broadcast,
         eltype, parent, real, imag,
         conj, transpose, adjoint,
         exp, log, sqrt, cos, sin, tan, csc, sec, cot,
                   cosh, sinh, tanh, csch, sech, coth,
                   acos, asin, atan, acsc, asec, acot,
                   acosh, asinh, atanh, acsch, asech, acoth, (:),
         AbstractMatrix, AbstractArray, checkindex, unsafe_length, OneTo,
        to_shape, _sub2ind, print_matrix, print_matrix_row, print_matrix_vdots,
      checkindex, Slice, @propagate_inbounds, @_propagate_inbounds_meta,
      _in_range, _range, _rangestyle, Ordered,
      ArithmeticWraps, floatrange, reverse, unitrange_last,
      AbstractArray, AbstractVector, axes, (:), _sub2ind_recurse, broadcast, promote_eltypeof,
      similar, @_gc_preserve_end, @_gc_preserve_begin,
      has_offset_axes, @nexprs, @ncall, @ntuple

import Base.Broadcast: BroadcastStyle, AbstractArrayStyle, Broadcasted, broadcasted,
                        combine_eltypes, DefaultArrayStyle, instantiate, materialize,
                        materialize!

import LinearAlgebra: AbstractTriangular, checksquare, pinv

import LinearAlgebra.BLAS: BlasFloat, BlasReal, BlasComplex

import FillArrays: AbstractFill

import StaticArrays: StaticArrayStyle

export Mul, MulArray, MulVector, MulMatrix, InvMatrix, PInvMatrix,
        Hcat, Vcat, Kron, BroadcastArray, cache, Ldiv, Inv, PInv, Diff, Cumsum

include("memorylayout.jl")
include("cache.jl")
include("lazybroadcasting.jl")
include("lazyconcat.jl")
include("linalg/linalg.jl")
include("lazysetoperations.jl")
include("lazyoperations.jl")

end # module


text:
using StatsBase
using AdaptiveResonance
using ClusterValidityIndices
using Logging
using Plots
using Random

# Plotting options
theme(:dark)    # Plotting style
dpi = 300       # Plotting dots-per-inch
# gr()            # GR backend (default for Plots.jl)
unicodeplots()

# Parameters
n = 5                           # Window size

# Include the library definitions
# include(projectdir("julia/lib_sim.jl"))
include("test_utils.jl")

# Set the logging level to Info and standardize the random seed
LogLevel(Logging.Info)
Random.seed!(0)
# data_dir(args...) = projectdir("work/data/meta_icvi", args...)
# results_dir(args...) = projectdir("work/results/meta_icvi", args...)
data_dir(args...) = joinpath("../data/testing", args...)
results_dir(args...) = joinpath("../data/results", args...)

# Load the data and test across all supervised modules
data = load_iris(data_dir("Iris.csv"))
data.train_x = convert(Array{Float64, 2}, data.train_x)
data.train_y = convert(Array{Int64, 1}, data.train_y)

# data.train_x, data.train_y = sort_cvi_data(convert(Matrix{Real}, data.train_x), data.train_y)
data.test_x, data.test_y = sort_cvi_data(convert(Matrix{Real},data.test_x), data.test_y)
# data.test_x, data.test_y = sort_cvi_data(data.test_x, data.test_y)

# Initialize the ART module
art = DDVFA()
# Set up the data manually because the module can't infer from single samples
# data_setup!(art.config, data.train_x)

# Get the dimension and size of the data
dim, n_samples = AdaptiveResonance.get_data_shape(data.train_x)
y_hat_train = zeros(Int64, n_samples)
_, n_samples_test = AdaptiveResonance.get_data_shape(data.test_x)
y_hat = zeros(Int64, n_samples_test)

# Create the CVIs
cvis = [
    PS(),
    GD43()
]
cvis_test = deepcopy(cvis)
n_cvis = length(cvis)

# Train the DDVFA model
y_hat_train = train!(art, data.train_x, y=data.train_y)
# y_hat_train = train!(art, train_x)
y_hat = classify(art, data.test_x)

# Preallocate the criterion_values
criterion_values = zeros(n_samples, n_cvis)
criterion_values_test = zeros(n_samples_test, n_cvis)
labels_ordered = relabel_cvi_data(y_hat_train)
labels_ordered_test = relabel_cvi_data(y_hat)
for i = 1:n_samples
    sample = convert(Array{Float64, 1}, data.train_x[:, i])
    label = labels_ordered[i]
    for j = 1:n_cvis
        criterion_values[i, j] = get_icvi!(cvis[j], sample, label)
    end
end
for i = 1:n_samples_test
    sample = convert(Array{Float64, 1}, data.test_x[:, i])
    label = labels_ordered_test[i]
    for j = 1:n_cvis
        criterion_values_test[i, j] = get_icvi!(cvis_test[j], sample, label)
    end
end

# Get the train correlations
n_corr = n_samples - n + 1      # Number of correlations
corrs = zeros(n_corr)           # Preallocation for correlations
iter_corr = n:n_samples         # Iterator for the correlation generation
for i = iter_corr
    corrs[i-n+1] = corspearman(criterion_values[i-n+1:i, 1], criterion_values[i-n+1:i, 2])
end

# Get the test correlations
n_corr_test = n_samples_test - n + 1      # Number of correlations
corrs_test = zeros(n_corr_test)           # Preallocation for correlations
iter_corr_test = n:n_samples_test         # Iterator for the correlation generation
for i = iter_corr_test
    corrs_test[i-n+1] = corspearman(criterion_values_test[i-n+1:i, 1], criterion_values_test[i-n+1:i, 2])
end

# Plot the two incremental trends ("manual" and porcelain) atop one another
p = plot(dpi=dpi, legend=:bottomleft)
plot!(p, iter_corr, corrs, label="Spearman")
plot!(p, 1:n_samples, y_hat_train, label="Labels")
plot!(p, 1:n_samples, data.train_y, label="True Labels")
# plot!(p, 1:n_corr, criterion_values_p, label="Porcelain")
title!("DDVFA + Iris Train")
xlabel!("Sample Index")
ylabel!("Spearman")

# Plot the two incremental trends ("manual" and porcelain) atop one another
g = plot(dpi=dpi, legend=:bottomleft)
plot!(g, iter_corr_test, corrs_test, label="Spearman")
plot!(g, 1:n_samples_test, y_hat, label="Labels")
plot!(g, 1:n_samples_test, data.test_y, label="True Labels")
# plot!(p, 1:n_corr, criterion_values_p, label="Porcelain")
title!("DDVFA + Iris Test")
xlabel!("Sample Index")
ylabel!("Spearman")

# Show the plots together
l = @layout [a; b]
# pt = plot(p, g, layout = l, size=(500,500))
pt = plot(p, g, layout = l, size=(800,500))

# Save and show the plot
# png(pt, results_dir("1_correlation_iris"))
display(pt)

# xlims!(1, n_corr)
# ylims!(0, Inf)

# Calculate performance
perf_train = performance(y_hat_train, data.train_y)
perf_test = performance(y_hat, data.test_y)

@info "DDVFA Training Perf: $perf_train"
@info "DDVFA Testing Perf: $perf_test"


text:
# all test targets
available_targets = [
    "align",
    "phylo",
    "intervals",
    "seq",
    "var",
    "services",
    "structure",
    "tools",
    "util"
]

if isempty(ARGS)
    # run all available test targets
    targets = available_targets
else
    targets = ARGS
    invalids = setdiff(targets, available_targets)
    if !isempty(invalids)
        error("there are invalid test targets: ", join(invalids, ", "))
    end
end

include("TestFunctions.jl")

for target in targets
    include("$target/runtests.jl")
end


text:
using Images, CausticsEngineering

image = Images.load("./examples/cat_posing.jpg") # Check current working directory with pwd()
mesh, image_caustics = engineer_caustics(image);
nothing;

image_caustics'

max_stars_repo_name: Mehul-Bagaria/causticsEngineering

max_issues_repo_name: Mehul-Bagaria/causticsEngineering

max_forks_repo_name: Mehul-Bagaria/causticsEngineering


text:
using Test
using GeoParams

@testset "DiffusionCreepLaws" begin

# This tests the MaterialParameters structure
CharUnits_GEO   =   GEO_units(viscosity=1Pa*s, length=1m);

# Define a linear viscous creep law ---------------------------------
x1      =   DiffusionCreep()
@test x1.n.val == 1.0
@test x1.p.val == -3.0
@test x1.A.val == 1.5MPa^-1.0*s^-1*m^(3.0)

x2      =   DiffusionCreep(n=3, p=-3.0)
@test x2.A.val == 1.5MPa^-3.0*s^-1*m^(3.0)

# perform a computation with the dislocation creep laws
    # Calculate EpsII, using a set of pre-defined values
CharDim = GEO_units(length=1000km, viscosity=1e19Pa*s, stress=100MPa, temperature=1000C)
EpsII   = GeoUnit(1.0s^-1.0)
        Nondimensionalize!(EpsII,CharDim)
TauII   = GeoUnit(0.3MPa)
        Nondimensionalize!(TauII,CharDim)
P       = GeoUnit(1.0e9Pa)
        Nondimensionalize!(P,CharDim)
T       = GeoUnit(1400C)
        Nondimensionalize!(T,CharDim)
f       = GeoUnit(1000NoUnits)
        Nondimensionalize!(f,CharDim)
d       = GeoUnit(10mm)
        Nondimensionalize!(d,CharDim)
c       = CreepLawVariables(P=P,T=T,f=f,d=d)
Phase   = SetMaterialParams(Name="Viscous Matrix", Phase=1,
                                     Density   = ConstantDensity(),
                                     CreepLaws = DiffusionCreep(n=1.0NoUnits, r=1.0NoUnits, p=-3.0NoUnits, A=1.0e6MPa^(-1.0)*s^(-1.0)*µm^(3.0), E=335000.0J/mol, V=4.0e-6m^(3.0)/mol, Apparatus="Invariant"), CharDim = CharDim)
EpsII.val = ComputeCreepLaw_EpsII(TauII,Phase.CreepLaws[1],c)
# Tested value is Dimensional so need to non-dimensionalize it back to compare
@test EpsII.val  ≈ Nondimensionalize(7.823165072337786e-15s^(-1.0), CharDim) rtol = 1e-12
    # Check that once inverted, we get back the TauII that we used to calculate EpsII
NewTau = ComputeCreepLaw_TauII(EpsII,Phase.CreepLaws[1],c)
#@show NewTau
@test NewTau ≈ TauII.val

# Given stress
#@test ComputeCreepLaw_EpsII(1e6Pa, x1, CreepLawParams())==5e-13/s                # dimensional input

# Given strainrate
#@test ComputeCreepLaw_TauII(1e-13/s, x1, CreepLawParams())==1e18*2*1e-13Pa       # dimensional input
# -------------------------------------------------------------------

end


text:
#
# Functions defined for ITensors.jl objects
# that may be moved to ITensors.jl
#

#
# imports.jl
#

import Base: isapprox, eltype

using ITensors: AbstractMPS

######################################################
# TagSet
#

function random_tags()
  ts = TagSet()
  ntags = length(ts.data)
  tagtype = eltype(ts.data)
  for n in 1:length(ts.data)
    rand_tag = rand(tagtype)
    # Before ITensors v0.2, addtags expects a Tag.
    # Afterward, just needs an integer.
    rand_tag = ITensors.version() < v"0.2" ? ITensors.Tag(rand_tag) : rand_tag
    ts = addtags(ts, rand_tag)
  end
  return ts
end

######################################################
# IndexSet
#

function has_indpairs(is, plevs::Pair{Int,Int})
  return any(i -> hasplev(i, first(plevs)) && setprime(i, last(plevs)) in is, is)
end

######################################################
# ITensor
#

function sqrt_hermitian(ρ::ITensor; cutoff::Float64 = 1e-15)
  if !isapprox(swapprime(dag(ρ), 0 => 1), ρ)
    error("matrix is not hermitian")
  end
  D, U = eigen(ρ; ishermitian = true, cutoff = cutoff)
  sqrtD = D
  sqrtD .= sqrt.(D)
  return U' * sqrtD * dag(U)
end

######################################################
# MPS
#

## # For |ψ⟩ and |ϕ⟩, return |ψ⟩⊗⟨ϕ|
## function ITensors.outer(ψ::MPS, ϕ::MPS; kwargs...)
##   # XXX: implement by converting to MPOs and
##   # contracting the MPOs?
##   @assert ψ == ϕ'
##   return MPO(ϕ; kwargs...)
## end

eltype(ψ::MPS) = ITensor
eltype(ψ::MPO) = ITensor

function promote_leaf_eltypes(ψ::AbstractMPS)::Type{<:Number}
  eltypeψ = eltype(ψ[1])
  for n in 2:length(ψ)
    eltypeψ = promote_type(eltypeψ, eltype(ψ[n]))
  end
  return eltypeψ
end

function isapprox(
  x::AbstractMPS,
  y::AbstractMPS;
  atol::Real=0,
  rtol::Real=Base.rtoldefault(promote_leaf_eltypes(x), promote_leaf_eltypes(y), atol),
)
  #d = norm(x - y)
  normx² = inner(x, x)
  normy² = inner(y, y)
  d² = normx² + normy² - 2 * real(inner(x, y))
  @assert imag(d²) < 1e-15
  d = sqrt(abs(d²))
  normx = sqrt(abs(normx²))
  normy = sqrt(abs(normy²))
  return d <= max(atol, rtol * max(normx, normy))
end

isapprox(x::AbstractMPS, y::ITensor; kwargs...) = isapprox(prod(x), y; kwargs...)
isapprox(x::ITensor, y::AbstractMPS; kwargs...) = isapprox(y, x; kwargs...)


text:
################################################## Experiment Centric Functions

using mjWrap

## accepts model,
## current and next states,
## current and next observations,
## applied controls that took us from current to next state
function rewardfunction(x::mjWrap.mjSet,
                        s0::AbstractVector{Float64}, # pass in states
                        s1::AbstractVector{Float64},
                        o0::AbstractVector{Float64}, # observations
                        o1::AbstractVector{Float64},
                        ctrl::AbstractVector{Float64})::Float64

   pos0 = s0[1]
   pos1 = s1[1]
   height = s1[3]

   reward = (pos1 - pos0) / (x.dt*x.skip) # forward x direction
   #upright_bonus = 1.0
   #if height > 1.0
   #    reward += upright_bonus
   #end
   #reward -= 3.0 * (height - 0.65)^2 # mode == normal
   #reward = max(reward, -1.0)
   #reward -= norm(x.d.ctrl)^2

   return max(reward, -1.0)
end

function obsfunction(x::mjWrap.mjSet,
                     s::AbstractVector{Float64},
                     o::AbstractVector{Float64})
   o[1] = s[3] # copy z height and quat
   o[2] = s[4]
   o[3] = s[5]
   o[4] = s[6]
   o[5] = s[7]
end

## accepts model,
## default state,
## preallocated matrix of states (nqnv x K)
function initfunction!(x::mjWrap.mjSet,
                       s::mjWrap.TrajSamples)
   q0 = [0.0, 0.0, 0.65, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, -1.0, 0.0, -1.0, 0.0, 1.0]
   snew = view( s.state, :, 1, : )
   K = size(snew, 2)
   nq, nv, nu, ns = mjWrap.modelparams(x)
   s0 = s.s0
   angle = 2pi # diverse
   #angle = pi/20 # nondiverse
   conts = 0
   meanre = mean(sum(s.reward, 1))
   #stdre  = stdev(sum(s.reward, 1))
   for t=1:K
      re = sum(s.reward[:,t])
      val, idx = findmax(s.reward[:,t])
      #if re < meanre || s.state[3,end,t] < 0.8 # if not continuing
      if true #re < meanre # if not continuing
         snew[1:3,t]    = q0[1:3]    + rand() * 0.8 - 0.2  # init_qpos
         #snew[3] = 0.7
         snew[4:7,t]    = q0[4:7]

         zangle         = ( rand() * angle ) - angle/2.0
         snew[4,t]      = cos( zangle/2.0 ) # set W correctly
         snew[7,t]      = sin( zangle/2.0 ) # z rotation
         snew[8:x.nq,t] = q0[8:x.nq] + rand() * 1.0 - 0.5  # init_qpos

         snew[(1+x.nq):end,t] = rand(x.nv) * 0.1 - 0.05
      else
         conts += 1
         snew[:,t] = s.state[:,idx,t]
      end
   end
end

function ctrlfunction!(x::mjWrap.mjSet, ctrl::Array{Float64})
   randn!(ctrl)
end

# modelfunc should takes in time index (iteration index) to know change curric
function modelfunction!(x::mjWrap.mjSet, iter::Int=0)
end

## evaluation function: reward function shapes exploration, evaluation
## function defines success or failure
function evalfunction(x::mjWrap.mjSet, d::mjWrap.TrajSamples)
   return mean(@view(d.state[1,d.T,:])) # x position at time T, for all K
end

function modelfunction!(x::mjWrap.mjSet)
end

using ExpFunctions
myfuncs = FunctionSet(modelfunction!,
                      initfunction!,
                      ctrlfunction!,
                      obsfunction,
                      rewardfunction,
                      evalfunction)


text:
module Synchrony
using DSP

export frequencies

# Get next fast FFT size for a given signal length
nextfastfft(n) = nextprod([2, 3, 5, 7], n)

# Compute frequencies based on length of padded FFT
function frequencies(nfft::Int, fs::Real=1.0)
    freq = (0:div(nfft, 2))*(fs/nfft)
    (freq, 1:length(freq))
end

# Compute frequencies based on length of padded FFT, with limits
function frequencies(nfft::Int, fs::Real, fmin::Real, fmax::Real=Inf)
    (freq,) = frequencies(nfft, fs)
    freqrange = searchsortedfirst(freq, fmin):(fmax == Inf ? length(freq) : searchsortedlast(freq, fmax))
    (freq[freqrange], freqrange)
end

# Compute frequencies based on data length, assuming padding to nextfastfft
frequencies(A::Union{AbstractVector{T}, AbstractMatrix{T}, AbstractArray{T,3}},
            fs::Real=1.0, fmin::Real=0.0, fmax::Real=Inf) where {T<:Real} =
    frequencies(nextfastfft(size(A, 1)), fs, fmin, fmax)

# Get the equivalent complex type for a given type
complextype(::Type{T}) where {T<:Complex} = T
complextype(::Type{T}) where {T} = Complex{T}

# Get preferred output type for a given input type
outputtype(::Type{T}) where {T<:AbstractFloat} = T
outputtype(::Type{T}) where {T<:Real} = Float64

include("transform_stats.jl")
include("point_field.jl")
include("wavelet.jl")
end


text:
using BinDeps
using LinearAlgebra, Libdl

BinDeps.@setup

blas = library_dependency("libblas", alias=["libblas.dll"])
lapack = library_dependency("liblapack", alias=["liblapack.dll"])

const ENV_VAR      = "CSDP_USE_JULIA_LAPACK"
const JULIA_LAPACK = if haskey(ENV, ENV_VAR)
    value = ENV[ENV_VAR]
    if lowercase(value) in ["1", "true", "yes"]
        @info "Using the blas and lapack libraries shipped with Julia as the environment variable `$ENV_VAR` is set to `$value`."
        true
    elseif lowercase(value) in ["0", "false", "no"]
        @info "Using system blas and lapack libraries as the environment variable `$ENV_VAR` is set to `$value`."
        false
    else
        error("The environment variable `$ENV_VAR` is set to `$value`. Set it to `true` or `false` instead.")
    end
else
    if BinDeps.issatisfied(blas) && BinDeps.issatisfied(lapack)
        @info "Using system blas and lapack libraries. Set the environment variable `$ENV_VAR` to `true` to use the blas/lapack library shipped with Julia."
        false
    else
        @info "Using the blas and lapack libraries shipped with Julia as there is no system blas and lapack libraries installed."
        true
    end
end

include("constants.jl")
include("compile.jl")

# @info "libname = $libname"
const depends = if JULIA_LAPACK
    # Create a new `bindeps_context` global variable that does not
    # have `blas` and `lapack` in the list of dependencies.
    BinDeps.@setup
    []
else
    [blas, lapack]
end
depends = JULIA_LAPACK ? [] : [blas, lapack]

# LaPack/BLAS dependencies
if !JULIA_LAPACK
    @static if Sys.iswindows()
        # wheel = "numpy/windows-wheel-builder/raw/master/atlas-builds"
        # atlas = "https://github.com/$wheel"
        # atlasdll = "/atlas-3.11.38-sse2-64/lib/numpy-atlas.dll"
        # download("https://raw.githubusercontent.com/$wheel/$atlasdll"),
        #           "$libdir/libatlas.dll")
        ## at the end ...
        # push!(BinDeps.defaults, BuildProcess)
    end
end

csdp = library_dependency("csdp", aliases=[libname], depends=depends)

provides(Sources, URI(download_url), csdp, unpacked_dir=csdpversion)

provides(BuildProcess,
         (@build_steps begin
             GetSources(csdp)
             CreateDirectory(libdir)
             CreateDirectory(builddir)
             @build_steps begin
                  ChangeDirectory(srcdir)
                  patch_int
                  compile_objs
             end
         end),
         [csdp])

# Prebuilt DLLs for Windows
provides(Binaries,
   URI("https://github.com/EQt/winlapack/blob/49454aee32649dc52c5b64f408a17b5270bd30f4/win-csdp-$(Sys.WORD_SIZE).7z?raw=true"),
   [csdp, lapack, blas], unpacked_dir="usr", os = :Windows)

@BinDeps.install Dict(:csdp => :csdp)

open(joinpath(dirname(@__FILE__), "deps.jl"), write = true, append = true) do io
    print(io, "const CSDP_INT = ")
    if JULIA_LAPACK
        println(io, "Clong")
    else
        println(io, "Cint")
    end
end


text:
using Compat: Nothing
"""
    AbstractThings

Supertype for container of objects being observed, whether these are
species, sequences, tips of a phylogeny (which could be either), or
some other type of thing. This will contain the names of the things
being observed, and (optionally) metadata about them, such as a
phylogeny that connects them, taxonomic information, their sequences,
trait information, information on similarity between the different
things, etc.

"""
abstract type AbstractThings end

"""
    AbstractLocationData

Composed within AbstractPlaces in cases when geographic location data exists. It
can reference locations with some geographical component. This may be a
series of arbitrarily arranged points, a series of areas, or even grid
of regularly spaced quadrats (see subtype AbstractGrid).

"""
abstract type AbstractLocationData end

"""
    AbstractPlaces{LocationDataType <: Union{Nothing, AbstractLocationData}}

AbstractPlaces is the supertype for containers of the places where things are
found (see AbstractThings). This will contain names or a reference for the
places, and (optionally) metadata such as what kind of place these are.
AbstractPlaces is parameterised by the spatial location data type for the
places. This should be Nothing if the places have no associated spatial data, or
a subtype of AbstractLocationData if they have spatial data. Other metadata in
the AbstractPlaces subtype should be in the AbstractPlaces subtype.

"""
abstract type AbstractPlaces{LocationDataType <: Union{Nothing, AbstractLocationData}} end

"""
    AbstractPoints <: AbstractLocationData

Subtype of AbstractLocationData where locations are a series of points in space.
"""
abstract type AbstractPoints <: AbstractLocationData end

"""
    AbstractGrid <: AbstractLocationData

Subtype of AbstractLocationData where locations are a grid of regularly
spaced, identically shaped, locations.
"""
abstract type AbstractGrid <: AbstractLocationData end

"""
    AbstractAssemblage{D <: Real (e.g. Int, Float64, Bool),
                       T <: AbstractThings,
                       P <: AbstractPlaces}

An assemblage of things recorded as being present in one or more
places. These may, for instance, be species counts in quadrats over a
regular grid, relative abundance of viral sequences in a group of
individuals, or presence-absence of genera over multiple islands.

"""
abstract type AbstractAssemblage{D <: Real,
                                 T <: AbstractThings,
                                 P <: AbstractPlaces} end

# depend explicitly on AxisArrays (and/or sparse) in some implementation?


text:
struct Suppressed{T}
    item::T
end
Base.show(io::IO, x::Suppressed) = print(io, "<suppressed ", x.item, '>')

function print_var(io::IO, var::JuliaInterpreter.Variable)
    print(io, "  | ")
    T = typeof(var.value)
    local val
    try
        val = repr(var.value)
        if length(val) > 150
            val = Suppressed("$(length(val)) bytes of output")
        end
    catch
        val = Suppressed("printing error")
    end
    println(io, highlight_code(string(var.name, "::", T, " = ", val); context=io))
end

function print_locals(io::IO, frame::Frame)
    vars = JuliaInterpreter.locals(frame)
    for var in vars
        # Hide gensymmed variables
        startswith(string(var.name), "#") && continue
        print_var(io, var)
    end
end

function print_frame(io::IO, num::Integer, frame::Frame; current_line=false)
    print(io, "[$num] ")
    println(io, locdesc(frame; current_line=current_line))
    print_locals(io, frame)
end

function pattern_match_kw_call(expr)
    if isexpr(expr, :call)
        f = string(expr.args[1])
        is_kw = occursin("#kw#", f) || (startswith(f, "#") && endswith(f, "_kw"))
    else
        is_kw = false
    end
    is_kw || return expr
    args = length(expr.args) >= 4 ? expr.args[4:end] : []
    kws_nt = expr.args[2]
    kws = []
    for (k, w) in pairs(kws_nt)
        push!(kws, Expr(:kw, k, w))
    end
    f = expr.args[3]
    return :($f($(args...); $(kws...)))
end

@static if VERSION < v"1.3.0-DEV.179"
    const append_any = Base.append_any
else
    append_any(@nospecialize x...) = append!([], Core.svec((x...)...))
end

function pattern_match_apply_call(expr, frame)
    if !(isexpr(expr, :call) && expr.args[1] == Core._apply)
        return expr
    end
    args = [@lookup(frame, expr.args[i+2]) for i in 1:(length(expr.args)-2)]
    new_expr = Expr(:call, expr.args[2])
    argsflat = append_any(args...)
    for x in argsflat
        push!(new_expr.args, (isa(x, Symbol) || isa(x, Expr) || isa(x, QuoteNode)) ? QuoteNode(x) : x)
    end
    return new_expr
end

function print_next_expr(io::IO, frame::Frame)
    expr = pc_expr(frame)
    @assert expr !== nothing
    print(io, "About to run: ")
    isa(expr, Expr) && (expr = copy(expr))
    if isexpr(expr, :(=))
        expr = expr.args[2]
    end
    if isexpr(expr, :call) || isexpr(expr, :return)
        for i in 1:length(expr.args)
            val = try
                @lookup(frame, expr.args[i])
            catch err
                err isa UndefVarError || rethrow(err)
                expr.args[i]
            end
            expr.args[i] = maybe_quote(val)
        end
    end
    expr = pattern_match_kw_call(expr)
    expr = pattern_match_apply_call(expr, frame)
    if isa(expr, Expr)
        for (i, arg) in enumerate(expr.args)
            try
                nbytes = length(repr(arg))
                if nbytes > max(40, div(200, length(expr.args)))
                    expr.args[i] = Suppressed("$nbytes bytes of output")
                end
            catch
                expr.args[i] = Suppressed("printing error")
            end
        end
    end
    print(io, highlight_code(string(expr); context=io))
    println(io)
end

function breakpoint_linenumbers(frame::Frame; lowered=false)
    framecode = frame.framecode
    breakpoint_lines = Dict{Int, BreakpointState}()
    for stmtidx in 1:length(framecode.breakpoints)
        isassigned(framecode.breakpoints, stmtidx) || continue
        bp = framecode.breakpoints[stmtidx]
        line = lowered ? stmtidx : JuliaInterpreter.linenumber(frame, stmtidx)
        breakpoint_lines[line] = bp
    end
    return breakpoint_lines
end

function print_status(io::IO, frame::Frame; force_lowered=false)
    # Buffer to avoid flickering
    outbuf = IOContext(IOBuffer(), io)
    printstyled(outbuf, "In ", locdesc(frame), "\n")
    loc = locinfo(frame)

    if loc !== nothing && !force_lowered
        defline, current_line, body = loc
        breakpoint_lines = breakpoint_linenumbers(frame)
        ok = print_sourcecode(outbuf, body, current_line, defline, breakpoint_lines)
        if !ok
            printstyled(io, "failed to lookup source code, showing lowered code:\n"; color=Base.warn_color())
            print_codeinfo(outbuf, frame)
        end
    else
        print_codeinfo(outbuf, frame)
    end
    print_next_expr(outbuf, frame)
    print(io, String(take!(outbuf.io)))
end

const NUM_SOURCE_LINES_UP_DOWN = Ref(4)

function print_codeinfo(io::IO, frame::Frame)
    src = frame.framecode.src
    if isdefined(JuliaInterpreter, Symbol("replace_coretypes!"))
        src = JuliaInterpreter.copy_codeinfo(src)
        JuliaInterpreter.replace_coretypes!(src; rev=true)
    end
    code = JuliaInterpreter.framecode_lines(src)
    active_line = convert(Int, frame.pc[])
    startline = max(1, active_line - NUM_SOURCE_LINES_UP_DOWN[])
    endline = min(length(code), active_line + NUM_SOURCE_LINES_UP_DOWN[])
    code = code[startline:endline]
    breakpoint_lines = breakpoint_linenumbers(frame; lowered=true)
    print_lines(io, code, active_line, breakpoint_lines, startline)
end

"""
Determine the offsets in the source code to print, based on the offset of the
currently highlighted part of the code, and the start and stop line of the
entire function.
"""
function compute_source_offsets(code::AbstractString, current_offsetline::Integer, file::SourceFile)
    desired_startline = current_offsetline - NUM_SOURCE_LINES_UP_DOWN[]
    desired_stopline = current_offsetline + NUM_SOURCE_LINES_UP_DOWN[] + 1
    if desired_startline > length(file.offsets)
        return -1, -1
    end
    desired_startline = max(desired_startline, 1)
    startoffset = file.offsets[desired_startline]
    stopoffset = lastindex(code)-1
    if desired_stopline < lastindex(file.offsets)
        stopoffset = file.offsets[desired_stopline]
    end
    startoffset, stopoffset
end

@enum HighlightOption begin
    HIGHLIGHT_OFF
    HIGHLIGHT_SYSTEM_COLORS
    HIGHLIGHT_256_COLORS
    HIGHLIGHT_24_BIT
end

const _syntax_highlighting = Ref(Sys.iswindows() ? HIGHLIGHT_SYSTEM_COLORS : HIGHLIGHT_256_COLORS)
const _current_theme = Ref{Type{<:Highlights.AbstractTheme}}(Highlights.Themes.MonokaiMiniTheme)

set_theme(theme::Type{<:Highlights.AbstractTheme}) = _current_theme[] = theme
set_highlight(opt::HighlightOption) = _syntax_highlighting[] = opt

function Format.render(io::IO, ::MIME"text/ansi-debugger", tokens::Format.TokenIterator)
    for (str, id, style) in tokens
        fg = style.fg.active ? map(Int, (style.fg.r, style.fg.g, style.fg.b)) : :nothing
        bg = style.bg.active ? map(Int, (style.bg.r, style.bg.g, style.bg.b)) : :nothing
        crayon = Crayon(
            foreground = fg,
            background = bg,
            bold       = style.bold ? true : :nothing,
            italics    = style.italic ? true : :nothing,
            underline  = style.underline ? true : :nothing,
        )
        if _syntax_highlighting[] == HIGHLIGHT_256_COLORS
            crayon = Crayons.to_256_colors(crayon)
        elseif _syntax_highlighting[] == HIGHLIGHT_SYSTEM_COLORS
            crayon = Crayons.to_system_colors(crayon)
        end
        print(io, crayon, str, inv(crayon))
    end
end

function highlight_code(code; context=nothing)
    if _syntax_highlighting[] != HIGHLIGHT_OFF
        try
            sprint(highlight, MIME("text/ansi-debugger"), code, Lexers.JuliaLexer, _current_theme[]; context=context)
        catch e
            printstyled(stderr, "failed to highlight code, $e\n"; color=Base.warn_color())
            return code
        end
    else
        return code
    end
end

const RESET = Crayon(reset = true)
function breakpoint_char(bp::BreakpointState)
    if bp.isactive
        return bp.condition === JuliaInterpreter.truecondition ? '●' : '◐'
    end
    return bp.condition === JuliaInterpreter.falsecondition ? ' ' : '○'
end

function print_sourcecode(io::IO, code::AbstractString, current_line::Integer, defline::Integer, breakpoint_lines::Dict{Int, BreakpointState} = Dict{Int, BreakpointState}())
    code = highlight_code(code; context=io)
    file = SourceFile(code)
    current_offsetline = current_line - defline + 1
    checkbounds(Bool, file.offsets, current_offsetline) || return false

    startoffset, stopoffset = compute_source_offsets(code, current_offsetline, file)
    if startoffset == -1
        printstyled(io, "Line out of file range (bad debug info?)")
        return false
    end

    # Compute necessary data for line numbering
    startline = compute_line(file, startoffset)
    code = split(code[(startoffset+1):(stopoffset+1)], '\n')
    print_lines(io, code, current_line, breakpoint_lines, startline + defline - 1)
    return true
end

function print_lines(io, code, current_line, breakpoint_lines, startline)
    if !isempty(code) && all(isspace, code[end])
        pop!(code)
    end
    stopline = startline + length(code) - 1

    # Count indentation level (only count spaces for now)
    min_indentation = typemax(Int)
    for textline in code
        all(isspace, textline) && continue
        isempty(textline) && continue
        indent_line = 0
        for char in textline
            char != ' ' && break
            indent_line += 1
        end
        min_indentation = min(min_indentation, indent_line)
    end
    for i in 1:length(code)
        code[i] = code[i][min_indentation+1:end]
    end
    lineno = startline
    stoplinelength = ndigits(stopline)
    for textline in code
        break_on_line = haskey(breakpoint_lines, lineno)
        prefix = (" ", :normal)
        break_on_line           && (prefix = (breakpoint_char(breakpoint_lines[lineno]), :light_red))
        lineno == current_line  && (prefix = (">", :yellow))
        printstyled(io,
            string(prefix[1], lpad(lineno, stoplinelength), "  "),
            color = prefix[2])

        println(io, textline)
        lineno += 1
    end
    _syntax_highlighting[] == HIGHLIGHT_OFF || print(io, RESET)
    println(io)
end


text:
using Test
using CLIMAParameters: AbstractEarthParameterSet
using CLIMAParameters.Planet
using CLIMAParameters.Atmos.SubgridScale

@testset "Atmos SubgridScale" begin
    struct EarthParameterSet <: AbstractEarthParameterSet end

    ps = EarthParameterSet()
    # Atmos.SubgridScale
    @test c_1_KASM(ps) ≈ c_a_KASM(ps) * 0.76^2
    @test c_2_KASM(ps) ≈ c_e2_KASM(ps) + 2 * c_1_KASM(ps)
    @test c_3_KASM(ps) ≈ c_a_KASM(ps)^(3 / 2)

end

max_stars_repo_name: climate-machine/CLIMAParameters

max_issues_repo_name: climate-machine/CLIMAParameters

max_forks_repo_name: climate-machine/CLIMAParameters


text:
using CSV, DataFrames, MVDA, Latexify, Printf

function run(files)
    metrics = [:train, :validation, :test, :model]
    columns = Symbol[]
    for metric in metrics
        push!(columns, Symbol(metric, :_md))
        push!(columns, Symbol(metric, :_lo))
        push!(columns, Symbol(metric, :_hi))
    end
    columns = [:sparsity; columns]

    output = DataFrame()
    datasets = String[]
    for file in files
        # strip the dataset name
        push!(datasets, first(split(basename(file), '-')))

        df_full = CSV.read(file, DataFrame)
        df_reps = MVDA.cv_error(df_full)            # approximate CV error by averaging over folds
        df_ints = MVDA.credible_intervals(df_reps)  # create 95% equal-tailed credible intervals

        s_opt = df_ints.model_md[1]
        idx = searchsortedfirst(df_ints.sparsity, s_opt)
        row = df_ints[idx, columns]

        result = DataFrame(
            x1="$(round(row.model_md, digits=2))",
            x2="$(round(row.model_lo, digits=2))",
            x3="$(round(row.model_hi, digits=2))",
            x4="$(round(row.train_md, digits=2))",
            x5="$(round(row.train_lo, digits=2))",
            x6="$(round(row.train_hi, digits=2))",
            x7="$(round(row.validation_md, digits=2))",
            x8="$(round(row.validation_lo, digits=2))",
            x9="$(round(row.validation_hi, digits=2))",
            x10="$(round(row.test_md, digits=2))",
            x11="$(round(row.test_lo, digits=2))",
            x12="$(round(row.test_hi, digits=2))",
            x13=first(df_full.nfeatures)
        )

        output = vcat(output, result)
    end

    # sort by features and drop the last column
    idx = sortperm(output, [:x13])
    output = output[idx,1:end-1]
    datasets .= datasets[idx]

    return latexify(output,
        env=:table,
        booktabs=true,
        adjustment=:r,
        side=datasets,
        transpose=false,
        latex=false,
    )
end

dir = ARGS[1]
@info "Reading from directory: $(dir)"
files = readdir(dir, join=true)
filter!(contains(".dat"), files)
display(run(files))


text:
using CellMLModelRepository
using CellMLToolkit
using CSV, DataFrames
using Test
using Base.Threads

datadir = joinpath(@__DIR__, "../data/")

# curl
df = cellml_metadata()
@test df isa DataFrame
# display(df)

# cellml_models()
# fns = readdir(joinpath(@__DIR__, "../data/cellml_models/"); join=true)[1:10] # only do 10
# @show fns
# @test isfile(fns[1])

# clone
p = joinpath(datadir, "repos/")
mkpath(p)
reposdf = cellml_repo_table()
@test reposdf isa DataFrame

workspaces_df = cellml_workspaces(reposdf)
@test workspaces_df isa DataFrame

cleandf = unique(strip.(dropmissing(workspaces_df)), :repo)
CSV.write("$(datadir)workspaces.csv", cleandf)
clone_physiome(p, cleandf)

result_df = run_all_repos(p; limit=1)
display(result_df)
mkpath("logs")
CSV.write("logs/physiome.csv", result_df)


text:
using Documenter
import Pkg
using Trixi2Vtk
using Trixi2Img

# Get Trixi root directory
trixi_root_dir = dirname(@__DIR__)

# Install dependencies and import modules...
# ...Trixi
Pkg.activate(trixi_root_dir)
Pkg.instantiate()
import Trixi

# Define module-wide setups such that the respective modules are available in doctests
DocMeta.setdocmeta!(Trixi,
                    :DocTestSetup,
                    :(push!(LOAD_PATH, ".."); using Trixi);
                    recursive=true)
DocMeta.setdocmeta!(Trixi2Vtk, :DocTestSetup, :(using Trixi2Vtk); recursive=true)
DocMeta.setdocmeta!(Trixi2Img, :DocTestSetup, :(using Trixi2Img); recursive=true)

# Make documentation
makedocs(
    # Specify modules for which docstrings should be shown
    modules = [Trixi, Trixi2Vtk, Trixi2Img],
    # Set sitename to Trixi
    sitename="Trixi.jl",
    # Provide additional formatting options
    format = Documenter.HTML(
        # Disable pretty URLs during manual testing
        prettyurls = get(ENV, "CI", nothing) == "true",
        # Explicitly add favicon as asset
        assets = ["assets/favicon.ico"],
        # Set canonical URL to GitHub pages URL
        canonical = "https://trixi-framework.github.io/Trixi.jl/stable"
    ),
    # Explicitly specify documentation structure
    pages = [
        "Home" => "index.md",
        "Development" => "development.md",
        "Visualization" => "visualization.md",
        "Style guide" => "styleguide.md",
        "GitHub & Git" => "github-git.md",
        "Reference" => [
                        "Trixi.jl" => "reference-trixi.md",
                        "Trixi2Vtk.jl" => "reference-trixi2vtk.md",
                        "Trixi2Img.jl" => "reference-trixi2img.md",
                       ],
        "Authors" => "authors.md",
        "Contributing" => "contributing.md",
        "License" => "license.md"
    ]
)

deploydocs(
    repo = "github.com/trixi-framework/Trixi.jl",
)


text:
function coreMVI(logl::Function, gradlogl::Function, LAposteriors; seed = 1, S = 100, optimiser = Optim.LBFGS(), iterations = 1, numerical_verification = false, Stest=0, show_every=-1, inititerations=0)

    D = length(mean(LAposteriors[1]))

    @assert(D == size(cov(LAposteriors[1]), 1) == size(cov(LAposteriors[1]), 2))

    @printf("Running MVI with S=%d, D=%d for %d iterations\n", S, D, iterations)

    #----------------------------------------------------
    # generate latent variables
    #----------------------------------------------------

    Ztrain = generatelatentZ(S = S, D = D, seed = seed)

    Ztest  = generatelatentZ(S = Stest, D = D, seed = seed+1)

    #----------------------------------------------------
    function unpack(param)
    #----------------------------------------------------

        @assert(length(param) == 2*D)

        local μ     = @view param[1:D]

        local Esqrt = @view param[D+1:2*D]

        return μ, Esqrt

    end

    #----------------------------------------------------
    function minauxiliary(V, param)
    #----------------------------------------------------

        local μ, Esqrt = unpack(param)

        return -1.0 * elbo(μ, Esqrt, V, Ztrain)

    end

    #----------------------------------------------------
    function minauxiliary_grad(V, param)
    #----------------------------------------------------

        local μ, Esqrt = unpack(param)

        return -1.0 * elbo_grad(μ, Esqrt, V, Ztrain)

    end

    #----------------------------------------------------
    function getcov(V, Esqrt)
    #----------------------------------------------------

        local aux = V * Diagonal(Esqrt)

        return aux*aux'

    end

    #----------------------------------------------------
    function getcovroot(V, Esqrt)
    #----------------------------------------------------

        return V * Diagonal(Esqrt)

    end

    #----------------------------------------------------
    function elbo(μ, Esqrt, V, Z)
    #----------------------------------------------------

        local C = getcovroot(V, Esqrt)

        mean(map(z -> logl(μ .+ C*z), Z)) + ApproximateVI.entropy(Esqrt)

    end

    #----------------------------------------------------
    function elbo_grad(μ, Esqrt, V, Z)
    #----------------------------------------------------

        local grad_μ     = zeros(D)

        local grad_Esqrt = zeros(D)

        local C          = getcovroot(V, Esqrt)

        local S          = length(Z)

        for i=1:S

            local aux = gradlogl(μ .+ C * Z[i])

            grad_μ += aux / S

            grad_Esqrt .+= V' * aux .* Z[i] / S

        end

        # logdet contribution in entropy

        grad_Esqrt .+= 1.0./Esqrt

        return [grad_μ; grad_Esqrt]

    end

    gradhelper(V, storage, param) = copyto!(storage, minauxiliary_grad(V, param))

    #----------------------------------------------------
    # Numerically verify gradient
    #----------------------------------------------------

    if numerical_verification

        local V, Esqrt = eigendecomposition(cov(LAposteriors[1]))
        adgrad = ForwardDiff.gradient(x->minauxiliary(V, x), [mean(LAposteriors[1]); Esqrt])
        angrad = minauxiliary_grad(V, [mean(LAposteriors[1]); Esqrt])
        @printf("gradient from AD vs analytical gradient\n")
        display([vec(adgrad) vec(angrad)])
        @printf("maximum absolute difference is %f\n", maximum(abs.(vec(adgrad) - vec(angrad))))

    end

    #----------------------------------------------------
    # Evaluate initial solutions for few iterations
    #----------------------------------------------------

    function initoptimise(q)

        local V, Esqrt = eigendecomposition(cov(q))

        Optim.optimize(x->minauxiliary(V,x), (x1,x2)->gradhelper(V,x1,x2), [mean(q); Esqrt], optimiser, Optim.Options(iterations = inititerations))

    end

    results = if inititerations > 0
        @showprogress "Initial search with random start " map(initoptimise, LAposteriors)
    else
        map(initoptimise, LAposteriors)
    end

    bestindex = argmin(map(r -> r.minimum, results))

    bestinitialsolution = results[bestindex].minimizer

    V, Esqrt = eigendecomposition(cov(LAposteriors[bestindex]))

    #----------------------------------------------------
    # Call optimiser
    #----------------------------------------------------

    options = Optim.Options(extended_trace = false, store_trace = false, show_trace = true, show_every=show_every, iterations = iterations, g_tol = 1e-6)

    result  = Optim.optimize(x->minauxiliary(V, x), (x1,x2)->gradhelper(V,x1,x2), bestinitialsolution, optimiser, options)

    μopt, Esqrtopt = unpack(result.minimizer)

    #----------------------------------------------------
    # Return results
    #----------------------------------------------------

    Σopt = getcov(V, Esqrtopt)

    return MvNormal(μopt, Σopt), elbo(μopt, Esqrtopt, V, generatelatentZ(S = 10*S, D = D, seed = seed+2))

end


text:
# Autogenerated wrapper script for FLANN_jll for x86_64-linux-musl-cxx03
export libflann, libflann_cpp

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("FLANN")
JLLWrappers.@declare_library_product(libflann, "libflann.so.1.9")
JLLWrappers.@declare_library_product(libflann_cpp, "libflann_cpp.so.1.9")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libflann,
        "lib/libflann.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libflann_cpp,
        "lib/libflann_cpp.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module MethodOfLines
    using LinearAlgebra
    using SciMLBase
    using DiffEqBase
    using DiffEqOperators
    using ModelingToolkit
    using ModelingToolkit: operation, istree, arguments, variable
    using SymbolicUtils, Symbolics
    using SymbolicUtils: operation, arguments
    using IfElse
    import DomainSets

    include("discretization/fornberg.jl")

    include("grid_types.jl")
    include("MOLFiniteDifference.jl")

    include("discretization/discretize_vars.jl")
    include("MOL_utils.jl")
    include("interiormap.jl")

    include("discretization/differential_discretizer.jl")
    include("bcs/boundary_types.jl")

    include("periodic_map.jl")

    include("discretization/generate_finite_difference_rules.jl")

    include("bcs/generate_bc_eqs.jl")

    include("error_analysis.jl")
    include("discretization/MOL_discretization.jl")

    export MOLFiniteDifference, discretize, symbolic_discretize, ODEFunctionExpr, generate_code, grid_align, edge_align, center_align, get_discrete

end


text:
"""
    WeightedGraph

Utility graph type to store one extra property with vertices and edges.
"""
mutable struct WeightedGraph{T} <: AbstractGraph{T}
    graph::AbstractGraph{T}
    vertex_w::Dict{T,Any}
    edge_w::Dict{Pair{T,T},Any}
end

Base.convert(::Type{WeightedGraph}, g::SimpleGraph) = WeightedGraph(g)
Base.convert(::Type{WeightedGraph}, g::SimpleDiGraph) = WeightedGraph(g)

function WeightedGraph(graph::AbstractGraph{T}) where {T}
    return WeightedGraph(graph, Dict{T,Any}(), Dict{Pair{T,T},Any}())
end

for extend in [
    :is_directed,
    :edgetype,
    :ne,
    :nv,
    :vertices,
    :edges,
    :outneighbors,
    :inneighbors,
    :has_vertex,
    :has_edge,
    :add_vertex!,
    :add_edge!,
    :rem_vertex!,
    :rem_edge!
]
    eval(quote
        Graphs.$extend(wg::WeightedGraph, args...) = $extend(wg.graph, args...)
    end)
end

"""
    add_vertex!(wg::WeightedGraph)
    add_vertex!(wg::WeightedGraph, weight::Any)
"""
function Graphs.add_vertex!(wg::WeightedGraph, weight::Any)
    check = add_vertex!(wg.graph)
    if check
        wg.vertex_w[nv(wg.graph)] = weight
    end
    return check
end

"""
    rem_vertex!(wg::WeightedGraph, vertex)
"""
function Graphs.rem_vertex!(wg::WeightedGraph{T}, vertex::T) where {T}
    check = rem_vertex!(wg.graph, vertex)
    if check
        for i in vertex:nv(wg.graph)
            if i == nv(wg.graph)
                delete!(wg.vertex_w, i)
                break
            end
            wg.vertex_w[i] = wg.vertex_w[i + one(T)]
        end
    end
    return check
end

"""
    add_edge!(wg::WeightedGraph, from_vertex, to_vertex)
    add_edge!(wg::WeightedGraph, from_vertex, to_vertex, weight::Any)
"""
function Graphs.add_edge!(wg::WeightedGraph{T}, from_vertex::T, to_vertex::T, weight::Any) where {T}
    check = add_edge!(wg.graph, from_vertex, to_vertex)
    if check
        wg.edge_w[(from_vertex, to_vertex)] = weight
    end
    return check
end

"""
    rem_edge!(wg::WeightedGraph{T}, from_vertex, to_vertex)
"""
function Graphs.rem_edge!(wg::WeightedGraph{T}, from_vertex::T, to_vertex::T) where {T}
    check = rem_edge!(wg.graph, from_vertex, to_vertex)
    if check
        delete!(wg.edge_w, (from_vertex, to_vertex))
    end
    return check
end

"""
    get_prop(wg::WeightedGraph, vertex::Int)
    get_prop(wg::WeightedGraph, from_vertex::Int, to_vertex::Int)
"""
function get_prop(wg::WeightedGraph, vertex::Int)
    return get(wg.vertex_w, vertex, nothing)
end

function get_prop(wg::WeightedGraph, from_vertex::Int, to_vertex::Int)
    return get(wg.edge_w, from_vertex=>to_vertex, nothing)
end

"""
    set_prop!(wg::WeightedGraph, vertex::Int, prop)
    set_prop!(wg::WeightedGraph, from_vertex::Int, to_vertex::Int, prop)
"""
function set_prop!(wg::WeightedGraph, vertex::Int, prop::Any)
    wg.vertex_w[vertex] = prop
end

function set_prop!(wg::WeightedGraph, from_vertex::Int, to_vertex::Int, prop::Any)
    wg.edge_w[from_vertex=>to_vertex] = prop
end

"""
    weights(wg::WeightedGraph)
"""
weights(wg::WeightedGraph) = edge_props(wg)

"""
    vertex_props(wg::WeightedGraph)
"""
function vertex_props(wg::WeightedGraph)
    map((vertex) -> wg.vertex_w[vertex], vertices(wg))
end

"""
    edge_props(wg::WeightedGraph)
"""
function edge_props(wg::WeightedGraph)
    map((e) -> wg.edge_w[Pair(e)], edges(wg))
end


text:
using Documenter, ACE

makedocs(sitename="ACEsuit Documentation",
         pages = [
        "Home" => "index.md",
        "Introduction" => "intro.md",
        "Getting Started" => "gettingstarted.md",
        "Developer Docs" => "devel.md",
        "ED-Bonds" => "envpairbasis.md"
        # "Subsection" => [
        #     ...
        # ]
        ])

# deploydocs(
#     repo = "github.com/JuliaMolSim/ACE.jl.git",
# )

max_stars_repo_name: ACEsuit/ACEdocs

max_issues_repo_name: ACEsuit/ACEdocs

max_forks_repo_name: ACEsuit/ACEdocs


text:
include_path = joinpath(
    abspath(dirname(dirname(@__FILE__))),
    "src/FilePaths.jl"
)
include(include_path)

using FilePaths

import URIParser: URI

if VERSION >= v"0.5-"
    using Base.Test
else
    using BaseTestNext
    const Test = BaseTestNext
end

info("Beginning tests...")

include("mode.jl")
include("path.jl")

info("All tests passed.")


text:
function dlpRead_test()
  o=MetoceanTools.dlpRead("AnonymousDLP.dlp")
end

max_stars_repo_name: jadamc/MetoceanTools

max_issues_repo_name: jadamc/MetoceanTools

max_forks_repo_name: jadamc/MetoceanTools


text:
# Frank-Wolfe Methods
# CFW and BFW in Mitradijieva and Lindberg (2013)

# required packages: Graphs, Optim

include("misc.jl")

function BPR(x::Vector{Float64}, td::TA_Data)
    bpr = similar(x)
    for i in 1:length(bpr)
        bpr[i] = td.free_flow_time[i] * ( 1.0 + td.b[i] * (x[i]/td.capacity[i])^td.power[i] )
        bpr[i] += td.toll_factor * td.toll[i] + td.distance_factor * td.link_length[i]
    end
    return bpr
end

gradient = BPR

function objective(x::Vector{Float64}, td::TA_Data)
    # value = free_flow_time .* ( x + b.* ( x.^(power+1)) ./ (capacity.^power) ./ (power+1))
    # return sum(value)

    sum = 0.0
    for i in 1:length(x)
        sum += td.free_flow_time[i] * ( x[i] + td.b[i]* ( x[i]^(td.power[i]+1)) / (td.capacity[i]^td.power[i]) / (td.power[i]+1))
        sum += td.toll_factor * td.toll[i] + td.distance_factor * td.link_length[i]
    end
    return sum
end

function hessian(x::Vector{Float64}, td::TA_Data)
    no_arc = length(td.init_node)

    h = zeros(no_arc, no_arc)
    h_diag = hessian_diag(x, td)

    for i in 1:no_arc
        h[i,i] = h_diag[i]
    end

    return h

    #Link travel time = free flow time * ( 1 + b * (flow/capacity)^Power ).
end

function hessian_diag(x::Vector{Float64}, td::TA_Data)
    h_diag = Array{Float64}(undef, length(x))
    for i in 1:length(x)
        if td.power[i] >= 1.0
            h_diag[i] = td.free_flow_time[i] * td.b[i] * td.power[i] * (x[i]^(td.power[i]-1)) / (td.capacity[i]^td.power[i])
        else
            h_diag[i] = 0 # Some cases, power is zero.
        end
    end
    # h_diag = free_flow_time .* b .* power .* (x.^(power-1)) ./ (capacity.^power)

    return h_diag
    #Link travel time = free flow time * ( 1 + b * (flow/capacity)^Power ).
end

function all_or_nothing_single(travel_time::Vector{Float64}, td::TA_Data, graph, link_dic)
    local state::LightGraphs.DijkstraState{Float64, Int}
    x = zeros(size(td.init_node))

    for r in 1:size(td.travel_demand)[1]
        # for each origin node r, find shortest paths to all destination nodes
        state = TA_dijkstra_shortest_paths(graph, travel_time, r, td.init_node, td.term_node, td.first_thru_node)

        for s in 1:size(td.travel_demand)[2]
            # for each destination node s, find the shortest-path vector
            # load travel demand
            # x = x + travel_demand[r,s] * get_vector(state, r, s, link_dic)
            add_demand_vector!(x, td.travel_demand[r,s], state, r, s, link_dic)
        end
    end
    return x
end

# parallel computing version
# function all_or_nothing_parallel(travel_time::Vector{Float64}, td::TA_Data, graph, link_dic)
#     local state::LightGraphs.DijkstraState{Float64, Int}
#     vv = zeros(size(td.init_node))
#     x = zeros(size(td.init_node))

#     x = x .+ @distributed (+) for r in 1:size(td.travel_demand)[1]
#         # for each origin node r, find shortest paths to all destination nodes
#         # if there is any travel demand starting from node r.
#         vv = zeros(size(td.init_node))

#         if sum(td.travel_demand, dims=2)[r] > 0.0
#             state = TA_dijkstra_shortest_paths(graph, travel_time, r, td.init_node, td.term_node, td.first_thru_node)

#             for s=1:size(td.travel_demand)[2]
#                 # for each destination node s, find the shortest-path vector
#                 # v = get_vector(state, r, s, init_node, term_node)

#                 if td.travel_demand[r,s] > 0.0
#                     # load travel demand
#                     # vv = vv + travel_demand[r,s] * get_vector(state, r, s, link_dic)
#                     add_demand_vector!(vv, td.travel_demand[r,s], state, r, s, link_dic)
#                 end
#             end

#         end

#         vv
#     end

#     return x
# end

function all_or_nothing(travel_time::Vector{Float64}, td::TA_Data, graph, link_dic)
    # if nprocs() > 1 # if multiple CPU processes are available
    #     return all_or_nothing_parallel(travel_time, td, graph, link_dic)
    # else
    #     return all_or_nothing_single(travel_time, td, graph, link_dic)
    #     # when nprocs()==1, using @distributed just adds unnecessary setup time. I guess.
    # end

    return all_or_nothing_single(travel_time, td, graph, link_dic)
end

function ta_frank_wolfe(td::TA_Data; method=:bfw, max_iter_no=2000, step=:exact, log=:off, tol=1e-3)

    setup_time = time()

    if log==:on
        println("-------------------------------------")
        println("Network Name: $(ta_data.network_name)")
        println("Method: $method")
        println("Line Search Step: $step")
        println("Maximum Interation Number: $max_iter_no")
        println("Tolerance for AEC: $tol")
        println("Number of processors: ", nprocs())
    end

    # preparing a graph
    graph = create_graph(td.init_node, td.term_node)
    link_dic = sparse(td.init_node, td.term_node, collect(1:td.number_of_links))

    setup_time = time() - setup_time

    if log==:on
        println("Setup time = $setup_time seconds")
    end

    n_links = td.number_of_links

    iteration_time = time()

    # Finding a starting feasible solution
    travel_time = BPR(zeros(n_links), td)
    x0 = all_or_nothing(travel_time, td, graph, link_dic)

    # Initializing variables
    xk = x0
    tauk = 0.0
    yk_FW = x0
    sk_CFW = yk_FW
    Hk_diag = Vector{Float64}(undef, n_links)

    dk_FW = Vector{Float64}(undef, n_links)
    dk_bar = Vector{Float64}(undef, n_links)
    dk_CFW = Vector{Float64}(undef, n_links)
    dk = Vector{Float64}(undef, n_links)

    alphak = 0.0
    Nk = 0.0
    Dk = 0.0

    tauk = 0.0
    is_first_iteration = false
    is_second_iteration = false

    sk_BFW = yk_FW
    sk_BFW_old = yk_FW

    dk_bbar = Vector{Float64}(undef, n_links)
    muk = Vector{Float64}(undef, n_links)
    nuk = Vector{Float64}(undef, n_links)
    beta0 = 0.0
    beta1 = 0.0
    beta2 = 0.0

    for k in 1:max_iter_no
        # Finding yk
        travel_time = BPR(xk, td)
        yk_FW = all_or_nothing(travel_time, td, graph, link_dic)

        # Basic Frank-Wolfe Direction
        dk_FW = yk_FW - xk
        Hk_diag = hessian_diag(xk, td) # Hk_diag is a diagonal vector of matrix Hk

        # Finding a feasible direction
        if method == :fw # Original Frank-Wolfe
            dk = dk_FW
        elseif method == :cfw # Conjugate Direction F-W
            if k==1 || tauk > 0.999999 # If tauk=1, then start the process all over again.
                sk_CFW = yk_FW
                dk_CFW = sk_CFW - xk
            else
                dk_bar = sk_CFW - xk  # sk_CFW from the previous iteration k-1

                Nk = dot( dk_bar, Hk_diag .* dk_FW )
                Dk = dot( dk_bar, Hk_diag .* (dk_FW - dk_bar) )

                delta = 0.0001
                if Dk != 0.0 && 0.0 <= Nk/Dk <= 1.0 - delta
                    alphak = Nk/Dk
                elseif Dk != 0.0 && Nk/Dk > 1.0 - delta
                    alphak = 1.0 - delta
                else
                    alphak = 0.0
                end

                # Generating new sk_CFW and dk_CFW
                sk_CFW = alphak .* sk_CFW .+ (1.0 - alphak) .* yk_FW
                dk_CFW = sk_CFW .- xk
            end

            # Feasible Direction to Use for CFW
            dk = dk_CFW
        elseif method == :bfw # Bi-Conjugate Direction F-W

            if tauk > 0.999999
                is_first_iteration = true
                is_second_iteration = true
            end

            if k==1 || is_first_iteration       # First Iteration is like FW
                sk_BFW_old = yk_FW
                dk_BFW = dk_FW
                is_first_iteration = false
            elseif k==2 || is_second_iteration  # Second Iteration is like CFW
                dk_bar = sk_BFW_old - xk # sk_BFW_old from the previous iteration 1

                Nk = dot( dk_bar, Hk_diag .* dk_FW )
                Dk = dot( dk_bar, Hk_diag .* (dk_FW - dk_bar) )

                delta = 0.0001
                if Dk != 0.0 && 0.0 <= Nk/Dk <= 1.0 - delta
                    alphak = Nk/Dk
                elseif Dk != 0.0 && Nk/Dk > 1.0 - delta
                    alphak = 1.0 - delta
                else
                    alphak = 0.0
                end

                # Generating new sk_BFW and dk_BFW
                sk_BFW = alphak .* sk_BFW_old .+ (1-alphak) .* yk_FW
                dk_BFW = sk_BFW .- xk

                is_second_iteration = false
            else
                # println("over there $tauk")
                # sk_BFW, tauk is from iteration k-1
                # sk_BFW_old is from iteration k-2

                dk_bar  = sk_BFW - xk
                dk_bbar = tauk * sk_BFW - xk + (1.0 - tauk) * sk_BFW_old

                muk = - dot( dk_bbar, Hk_diag .* dk_FW ) / dot( dk_bbar, Hk_diag .* (sk_BFW_old - sk_BFW) )
                nuk = - dot( dk_bar, Hk_diag .* dk_FW ) / dot( dk_bar, Hk_diag .* dk_bar) + muk*tauk/(1-tauk)

                muk = max(0.0, muk)
                nuk = max(0.0, nuk)

                beta0 = 1.0 / ( 1.0 + muk + nuk )
                beta1 = nuk * beta0
                beta2 = muk * beta0

                sk_BFW_new = beta0 * yk_FW + beta1 * sk_BFW + beta2 * sk_BFW_old
                dk_BFW = sk_BFW_new - xk

                sk_BFW_old = sk_BFW
                sk_BFW = sk_BFW_new
            end

            # Feasible Direction to Use for BFW
            dk = dk_BFW
        else
            error("The type of Frank-Wolfe method is specified incorrectly. Use :fw, :cfw, or :bfw.")
        end
        # dk is now identified.

        if step == :exact
            # Line Search from xk in the direction dk
            optk = optimize(tau -> objective(xk+tau*dk, td), 0.0, 1.0, GoldenSection())
            tauk = optk.minimizer
        elseif step == :newton
            # Newton step
            tauk = - dot( gradient(xk, td), dk ) / dot( dk, Hk_diag.*dk )
            tauk = max(0.0, min(1.0, tauk))
        end

        # Average Excess Cost
        average_excess_cost = ( xk' * travel_time - yk_FW' * travel_time ) / sum(td.travel_demand)
        if log==:on
            @printf("k=%4d, tauk=%15.10f, objective=%15f, aec=%15.10f\n", k, tauk, objective(xk, td), average_excess_cost)
        end

        # rel_gap = ( objective(xk) - best_objective ) / best_objective

        # Convergence Test
        if average_excess_cost < tol
        # if rel_gap < tol
            break
        end

        # Update x
        new_x = xk + tauk * dk
        xk = new_x

        @assert minimum(xk) >= 0

    end

    iteration_time = time() - iteration_time

    if log==:on
        println("Iteration time = $iteration_time seconds")
    end

    return xk, travel_time, objective(xk, td)

end

#