text:
export group_algebra, galois_module, group

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

base_ring(A::AlgGrp{T}) where {T} = A.base_ring::parent_type(T)

Generic.dim(A::AlgGrp) = size(multiplication_table(A, copy = false), 1)

elem_type(::Type{AlgGrp{T, S, R}}) where {T, S, R} = AlgGrpElem{T, AlgGrp{T, S, R}}

order_type(::AlgGrp{fmpq, S, R}) where { S, R } = AlgAssAbsOrd{AlgGrp{fmpq, S, R}, elem_type(AlgGrp{fmpq, S, R})}

order_type(::Type{AlgGrp{fmpq, S, R}}) where { S, R } = AlgAssAbsOrd{AlgGrp{fmpq, S, R}, elem_type(AlgGrp{fmpq, S, R})}

order_type(::AlgGrp{T, S, R}) where { T <: NumFieldElem, S, R } = AlgAssRelOrd{T, fractional_ideal_type(order_type(parent_type(T)))}
order_type(::Type{AlgGrp{T, S, R}}) where { T <: NumFieldElem, S, R } = AlgAssRelOrd{T, fractional_ideal_type(order_type(parent_type(T)))}

@doc Markdown.doc"""
    group(A::AlgGrp) -> Group

Returns the group defining $A$.
"""
group(A::AlgGrp) = A.group

has_one(A::AlgGrp) = true

function (A::AlgGrp{T, S, R})(c::Array{T, 1}) where {T, S, R}
  length(c) != dim(A) && error("Dimensions don't match.")
  return AlgGrpElem{T, typeof(A)}(A, c)
end

@doc Markdown.doc"""
    multiplication_table(A::AlgGrp; copy::Bool = true) -> Array{RingElem, 2}

Given a group algebra $A$ this function returns the multiplication table of
$A$: If the function returns $M$ and the basis of $A$ is $g_1,\dots, g_n$ then
it holds $g_i \cdot g_j = g_{M[i, j]}$.
"""
function multiplication_table(A::AlgGrp; copy::Bool = true)
  if copy
    return deepcopy(A.mult_table)
  else
    return A.mult_table
  end
end

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

@doc Markdown.doc"""
    group_algebra(K::Ring, G; op = *) -> AlgGrp

Returns the group ring $K[G]$.
$G$ may be any set and `op` a group operation on $G$.
"""
group_algebra(K::Ring, G; op = *) = AlgGrp(K, G, op = op)

group_algebra(K::Ring, G::GrpAbFinGen) = AlgGrp(K, G)

function group_algebra(K::Field, G; op = *)
  A = AlgGrp(K, G, op = op)
  if iszero(characteristic(K))
    A.issemisimple = 1
  else
    A.issemisimple = isone(gcd(characteristic(K), order(G))) ? 1 : 2
  end
  return A
end

function group_algebra(K::Field, G::GrpAbFinGen)
  A = group_algebra(K, G, op = +)
  A.iscommutative = true
  return A
end

@doc Markdown.doc"""
    (K::Ring)[G::Group] -> AlgGrp
    (K::Ring)[G::GrpAbFinGen] -> AlgGrp

Returns the group ring $K[G]$.
"""
getindex(K::Ring, G::Group) = group_algebra(K, G)
getindex(K::Ring, G::GrpAbFinGen) = group_algebra(K, G)

################################################################################
#
#  Commutativity
#
################################################################################

iscommutative_known(A::AlgGrp) = (A.iscommutative != 0)

@doc Markdown.doc"""
    iscommutative(A::AlgGrp) -> Bool

Returns `true` if $A$ is a commutative ring and `false` otherwise.
"""
function iscommutative(A::AlgGrp)
  if iscommutative_known(A)
    return A.iscommutative == 1
  end
  for i in 1:dim(A)
    for j in 1:dim(A)
      if multiplication_table(A, copy = false)[i, j] != multiplication_table(A, copy = false)[j, i]
        A.iscommutative = 2
        return false
      end
    end
  end
  A.iscommutative = 1
  return true
end

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

function show(io::IO, A::AlgGrp)
  compact = get(io, :compact, false)
  if compact
    print(io, "Group algebra of dimension ", dim(A), " over ", base_ring(A))
  else
    print(io, "Group algebra of group\n", group(A), "\nover\n", base_ring(A))
  end
end

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

# TODO: This is broken. I have to copy everything carefully by hand.
#function Base.deepcopy_internal(A::AlgGrp, dict::IdDict)
#  B = AlgGrp(base_ring(A), group(A))
#  for x in fieldnames(typeof(A))
#    if x != :base_ring && x != :group && isdefined(A, x)
#      setfield!(B, x, Base.deepcopy_internal(getfield(A, x), dict))
#    end
#  end
#  B.base_ring = A.base_ring
#  B.group = A.group
#  return B
#end

###############################################################################
#
#  Trace basis
#
###############################################################################

function _assure_trace_basis(A::AlgGrp{T}) where {T}
  if isdefined(A, :trace_basis_elem)
    return nothing
  end

  A.trace_basis_elem = Vector{T}(undef, dim(A))
  A.trace_basis_elem[1] = base_ring(A)(dim(A))
  for i = 2:dim(A)
    A.trace_basis_elem[i] = zero(base_ring(A))
  end
  return nothing
end

################################################################################
#
#  Center
#
################################################################################

@doc Markdown.doc"""
    center(A::AlgGrp) -> AlgAss, AbsAlgAssMor

Returns the center $C$ of $A$ and the inclusion $C \to A$.
"""
function center(A::AlgGrp{T}) where {T}
  if isdefined(A, :center)
    return A.center::Tuple{AlgAss{T}, morphism_type(AlgAss{T}, typeof(A))}
  end

  # Unlike for AlgAss, we should cache the centre even if A is commutative
  # since it is of a different type, so A !== center(A)[1].
  # Otherwise center(A)[1] !== center(A)[1] which is really annoying.
  B, mB = AlgAss(A)
  C, mC = center(B)
  mD = compose_and_squash(mB, mC)
  A.center = C, mD
  return (C, mD)::Tuple{AlgAss{T}, morphism_type(AlgAss{T}, typeof(A))}
end

################################################################################
#
#  Conversion to AlgAss
#
################################################################################

function AlgAss(A::AlgGrp{T, S, R}) where {T, S, R}
  K = base_ring(A)
  mult = Array{T, 3}(undef, dim(A), dim(A), dim(A))
  B = basis(A)
  d = dim(A)
  for i in 1:d
    for j in 1:d
      l = multiplication_table(A, copy = false)[i, j]
      for k in 1:d
        if k == l
          mult[i, j, k] = one(K)
        else
          mult[i, j, k] = zero(K)
        end
      end
    end
  end
  B = AlgAss(K, mult, one(A).coeffs)
  B.iscommutative = A.iscommutative
  B.issimple = A.issimple
  B.issemisimple = A.issemisimple
  AtoB = hom(A, B, identity_matrix(K, dim(A)), identity_matrix(K, dim(A)))
  if isdefined(A, :center)
    Z, ZtoA = center(A)
    B.center = (Z, compose_and_squash(AtoB, ZtoA))
  end
  if isdefined(A, :decomposition)
    dec = Tuple{AlgAss{T}, morphism_type(AlgAss{T}, typeof(B))}[]
    for (C, CtoA) in A.decomposition::Vector{Tuple{AlgAss{T}, morphism_type(AlgAss{T}, typeof(A))}}
      CtoB = compose_and_squash(AtoB, CtoA)
      push!(dec, (C, CtoB))
    end
    B.decomposition = dec
  end
  if isdefined(A, :maps_to_numberfields)
    NF = A.maps_to_numberfields::Vector{Tuple{_ext_type(T),_abs_alg_ass_to_nf_abs_mor_type(A)}}
    fields_and_maps = Tuple{AnticNumberField,_abs_alg_ass_to_nf_abs_mor_type(B)}[]
    for (K, AtoK) in NF
      BtoK = AbsAlgAssToNfAbsMor(B, K, AtoK.mat, AtoK.imat)
      push!(fields_and_maps, (K, BtoK))
    end
    B.maps_to_numberfields = fields_and_maps
  end
  return B, hom(B, A, identity_matrix(K, dim(A)), identity_matrix(K, dim(A)))
end

################################################################################
#
#  Misc
#
################################################################################

Base.enumerate(G::Generic.SymmetricGroup{T}) where T <: Integer = enumerate(Generic.AllPerms(G.n))

################################################################################
#
#  Generators
#
################################################################################

# Helper function for gens, changes mid in place
function _merge_elts_in_gens!(left::Vector{Tuple{Int, Int}}, mid::Vector{Tuple{Int, Int}}, right::Vector{Tuple{Int, Int}})
  nl = length(left)
  if length(mid) == 0
    mid = deepcopy(left)
  elseif nl != 0
    if left[nl][1] == mid[1][1]
      t = popfirst!(mid)
      prepend!(mid, left)
      tt = (t[1], mid[nl][2] + t[2])
      mid[nl] = tt
    else
      prepend!(mid, left)
    end
  end
  nm = length(mid)
  if nm == 0
    return deepcopy(right)
  end
  if length(right) == 0
    return mid
  end
  if mid[nm][1] == right[1][1]
    t = pop!(mid)
    append!(mid, right)
    tt = (t[1], mid[nm][2] + t[2])
    mid[nm] = tt
  else
    append!(mid, right)
  end
  return mid
end

@doc Markdown.doc"""
    gens(A::AlgGrp, return_full_basis::Type{Val{T}} = Val{false})
      -> Vector{AlgGrpElem}

Returns a subset of `basis(A)`, which generates $A$ as an algebra over
`base_ring(A)`.
If `return_full_basis` is set to `Val{true}`, the function also returns a
`Vector{AbsAlgAssElem}` containing a full basis consisting of monomials in
the generators and a `Vector{Vector{Tuple{Int, Int}}}` containing the
information on how these monomials are built. E. g.: If the function returns
`g`, `full_basis` and `v`, then we have
`full_basis[i] = prod( g[j]^k for (j, k) in v[i] )`.
"""
function gens(A::AlgGrp, return_full_basis::Type{Val{T}} = Val{false}) where T
  G = group(A)
  group_gens = gens(G)

  return_full_basis == Val{true} ? nothing : return map(A, group_gens)

  full_group = elem_type(G)[ id(G) ]
  elts_in_gens = Vector{Tuple{Int, Int}}[ Tuple{Int, Int}[] ]
  constructed_elements = falses(BigInt(order(G)))
  constructed_elements[A.group_to_base[id(G)]] = true
  new_elements = Set{Int}()
  for i = 1:length(group_gens)
    g = group_gens[i]
    push!(full_group, g)
    push!(elts_in_gens, Tuple{Int, Int}[ (i, 1) ])
    constructed_elements[A.group_to_base[g]] = true
    push!(new_elements, length(full_group))
  end
  k = 1 + length(group_gens) # == number of constructed elements, i. e. #{ i | constructed_elements[i] == true }

  while k != dim(A) || !isempty(new_elements)
    i = pop!(new_elements)
    g = full_group[i]

    n = length(full_group)
    for r = 1:n
      s = op(g, full_group[r])
      for l = 1:n
        if !iscommutative(A)
          t = op(full_group[l], s)
        else
          t = s
        end
        if constructed_elements[A.group_to_base[t]]
          continue
        end
        constructed_elements[A.group_to_base[t]] = true
        k += 1
        push!(full_group, t)
        coord = _merge_elts_in_gens!(elts_in_gens[l], deepcopy(elts_in_gens[i]), elts_in_gens[r])
        push!(elts_in_gens, coord)
        push!(new_elements, length(full_group))
        if iscommutative(A)
          break
        end
        k == dim(A) ? break : nothing
      end
      k == dim(A) ? break : nothing
    end
  end

  return map(A, group_gens), map(A, full_group), elts_in_gens
end

################################################################################
#
#  Isomorphisms to number fields
#
################################################################################

# Assumes that Gal(K/k) == group(A), where k = base_field(K) and that group(A) is
# abelian.
# Returns a k-linear map from K to A and one from A to K
function _find_isomorphism(K::Union{ AnticNumberField, NfRel{nf_elem} }, A::AlgGrp)
  G = group(A)
  aut = automorphisms(K)

  aut_dict = Dict{elem_type(K), Int}()
  n = length(aut)
  identity = 0
  for i = 1:n
    b = image_primitive_element(aut[i])
    aut_dict[b] = i
    if b == gen(K)
      identity = i
    end
  end

  op_array = zeros(Int, n, n)
  for i = 1:n
    for j = i:n
      if i == identity
        k = j
      elseif j == identity
        k = i
      else
        b = aut[i](image_primitive_element(aut[j]))
        k = aut_dict[b]
      end
      op_array[i, j] = k
      # It is assumed, that the group is abelian
      op_array[j, i] = k
    end
  end

  local op
  let op_array = op_array
    function op(iiii::Int, jjjj::Int)
      return op_array[iiii, jjjj]
    end
  end

  auttoG, Gtoaut = find_isomorphism(Int[ i for i in 1:n ], op, G)

  alpha = normal_basis(K)
  basis_alpha = Vector{elem_type(K)}(undef, dim(A))
  for (i, g) in enumerate(G)
    f = aut[Gtoaut[g]]
    basis_alpha[A.group_to_base[g]] = f(alpha)
  end

  M = zero_matrix(base_field(K), degree(K), degree(K))
  for i = 1:degree(K)
    a = basis_alpha[i]
    for j = 1:degree(K)
      M[i, j] = coeff(a, j - 1)
    end
  end

  invM = inv(M)

  local KtoA
  let K = K, invM = invM, A = A
    function KtoA(x::Union{ nf_elem, NfRelElem })
      t = zero_matrix(base_field(K), 1, degree(K))
      for i = 1:degree(K)
        t[1, i] = coeff(x, i - 1)
      end
      y = t*invM
      return A([ y[1, i] for i = 1:degree(K) ])
    end
  end

  local AtoK
  let K = K, M = M
    function AtoK(x::AlgGrpElem)
      t = matrix(base_field(K), 1, degree(K), coefficients(x))
      y = t*M
      return K(parent(K.pol)([ y[1, i] for i = 1:degree(K) ]))
    end
  end

  return KtoA, AtoK
end

mutable struct NfToAlgGrpMor{S, T, U} <: Map{AnticNumberField, AlgGrp{S, T, U}, HeckeMap, AbsAlgAssMor}
  K::AnticNumberField
  mG::GrpGenToNfMorSet{AnticNumberField}
  A::AlgGrp{S, T, U}
  M::fmpq_mat
  Minv::fmpq_mat

  function NfToAlgGrpMor{S, T, U}() where {S, T, U}
    return new{S, T, U}()
  end
end

function show(io::IO, f::NfToAlgGrpMor)
  print(io, "Galois module map from\n")
  print(io, f.K)
  print(io, "\nto\n")
  print(io, f.A)
end

function (f::NfToAlgGrpMor)(O::NfAbsOrd)
  A = codomain(f)
  B = basis(O)
  G = group(A)
  ZG = Order(A, collect(G))
  return ideal_from_lattice_gens(A, ZG, [f(elem_in_nf(b)) for b in B], :right)
end

automorphism_map(f::NfToAlgGrpMor) = f.mG

function galois_module(K::AnticNumberField, aut::Map = automorphism_group(K)[2]; normal_basis_generator = normal_basis(K))
  G = domain(aut)
  A = FlintQQ[G]
  return _galois_module(K, A, aut, normal_basis_generator = normal_basis_generator)
end

function _galois_module(K::AnticNumberField, A, aut::Map = automorphism_group(K)[2]; normal_basis_generator = normal_basis(K))
  G = domain(aut)
  alpha = normal_basis_generator

  basis_alpha = Vector{elem_type(K)}(undef, dim(A))
  for (i, g) in enumerate(G)
    f = aut(g)
    basis_alpha[A.group_to_base[g]] = f(alpha)
  end

  M = zero_matrix(base_field(K), degree(K), degree(K))
  for i = 1:degree(K)
    a = basis_alpha[i]
    for j = 1:degree(K)
      M[i, j] = coeff(a, j - 1)
    end
  end

  invM = inv(M)

  z = NfToAlgGrpMor{fmpq, GrpGen, GrpGenElem}()
  z.K = K
  z.mG = aut
  z.A = A
  z.M = M
  z.Minv = invM

  return A, z
end

function galois_module(K::AnticNumberField, A::AlgGrp; normal_basis_generator = normal_basis(K))
  G = group(A)
  Au, mAu = automorphism_group(K)
  fl, f = isisomorphic(G, Au)
  @assert fl
  aut = Vector{NfToNfMor}(undef, order(G))
  for g in G
    aut[g[]] = mAu(f(g))
  end
  h = GrpGenToNfMorSet(G, aut, K)

  return _galois_module(K, A, h, normal_basis_generator = normal_basis(K))
end

domain(f::NfToAlgGrpMor) = f.K

codomain(f::NfToAlgGrpMor) = f.A

function image(f::NfToAlgGrpMor, x::nf_elem)
  K = domain(f)
  @assert parent(x) === K
  A = codomain(f)

  t = zero_matrix(base_field(K), 1, degree(K))
  for i = 1:degree(K)
    t[1, i] = coeff(x, i - 1)
  end
  y = t*f.Minv
  return A([ y[1, i] for i = 1:degree(K) ])
end

function preimage(f::NfToAlgGrpMor, x::AlgGrpElem)
  K = domain(f)
  t = matrix(base_field(K), 1, degree(K), coefficients(x))
  y = t*f.M
  v = fmpq[ y[1, i] for i = 1:degree(K) ]
  return K(v)
end

# Returns the group algebra Q[G] where G = Gal(K/Q) and a Q-linear map from K
# to Q[G] and one from Q[G] to K
function _galois_module(K::AnticNumberField, to_automorphisms::Map = automorphism_group(K)[2]; normal_basis_generator = normal_basis(K))
  G = domain(to_automorphisms)
  A = FlintQQ[G]
  alpha = normal_basis_generator

  basis_alpha = Vector{elem_type(K)}(undef, dim(A))
  for (i, g) in enumerate(G)
    f = to_automorphisms(g)
    basis_alpha[A.group_to_base[g]] = f(alpha)
  end

  M = zero_matrix(base_field(K), degree(K), degree(K))
  for i = 1:degree(K)
    a = basis_alpha[i]
    for j = 1:degree(K)
      M[i, j] = coeff(a, j - 1)
    end
  end

  invM = inv(M)

  function KtoA(x::nf_elem)
    t = zero_matrix(base_field(K), 1, degree(K))
    for i = 1:degree(K)
      t[1, i] = coeff(x, i - 1)
    end
    y = t*invM
    return A([ y[1, i] for i = 1:degree(K) ])
  end

  function AtoK(x::AlgGrpElem)
    t = matrix(base_field(K), 1, degree(K), coefficients(x))
    y = t*M
    return K(parent(K.pol)([ y[1, i] for i = 1:degree(K) ]))
  end

  return A, KtoA, AtoK
end

const _reps = [(i=24,j=12,n=5,dims=(1,1,2,3,3),
                reps=Vector{Vector{Rational{BigInt}}}[[[1],[1],[1],[1]],
                  [[-1],[1],[1],[1]],
                  [[-1,0,-1,1],[0,-1,1,-1],[1,0,0,1],[1,0,0,1]],
                  [[0,1,0,1,0,0,0,0,1],[0,0,1,1,0,0,0,1,0],[-1,0,0,0,1,0,0,0,-1],[1,0,0,0,-1,0,0,0,-1]],
                  [[0,-1,0,-1,0,0,0,0,-1],[0,0,1,1,0,0,0,1,0],[-1,0,0,0,1,0,0,0,-1],[1,0,0,0,-1,0,0,0,-1]]]),
               (i=48,j=48,n=10,dims=Int[1,1,1,1,2,2,3,3,3,3],
                reps=Vector{Vector{Rational{BigInt}}}[[[1],[1],[1],[1],[1]],[[-1],[1],[1],[1],[1]],[[1],[-1],[1],[1],[1]],[[-1],[-1],[1],[1],[1]],[[-1,0,-1,1],[-1,0,0,-1],[0,-1,1,-1],[1,0,0,1],[1,0,0,1]],[[1,3//2,0,-1],[1,0,0,1],[1,3//2,-2,-2],[1,0,0,1],[1,0,0,1]],[[0,1,0,1,0,0,0,0,1],[1,0,0,0,1,0,0,0,1],[0,0,1,1,0,0,0,1,0],[-1,0,0,0,1,0,0,0,-1],[1,0,0,0,-1,0,0,0,-1]],[[0,-1,0,-1,0,0,0,0,-1],[1,0,0,0,1,0,0,0,1],[0,0,1,1,0,0,0,1,0],[-1,0,0,0,1,0,0,0,-1],[1,0,0,0,-1,0,0,0,-1]],[[0,1,0,1,0,0,0,0,1],[-1,0,0,0,-1,0,0,0,-1],[0,0,1,1,0,0,0,1,0],[-1,0,0,0,1,0,0,0,-1],[1,0,0,0,-1,0,0,0,-1]],[[0,-1,0,-1,0,0,0,0,-1],[-1,0,0,0,-1,0,0,0,-1],[0,0,1,1,0,0,0,1,0],[-1,0,0,0,1,0,0,0,-1],[1,0,0,0,-1,0,0,0,-1]]])
              ]

################################################################################
#
#  Wedderburn decomposition
#
################################################################################

mutable struct AbsAlgAssMorGen{S, T, U, V} <: Map{S, T, HeckeMap, AbsAlgAssMorGen}
  domain::S
  codomain::T
  M::U
  Minv::V
end

#function AbsAlgAssMorGen(A::S, B::T, M::U, N::V) where {S, T, U, V}
#  return AbsAlgAssMorGen{S, T, U, V}(A, B, M, N)
#end

domain(f::AbsAlgAssMorGen) = f.domain

codomain(f::AbsAlgAssMorGen) = f.codomain

function Base.show(io::IO, f::AbsAlgAssMorGen)
  print(io, "Algebra morphism from \n$(domain(f))\n to\n$(codomain(f))\n")
end

function image(f::AbsAlgAssMorGen, z)
  @assert parent(z) == domain(f)
  v = matrix(base_ring(codomain(f)), 1, dim(domain(f)), coefficients(z))
  return codomain(f)(_eltseq(v * f.M))
end

(f::AbsAlgAssMorGen)(z::AbsAlgAssElem) = image(f, z)

function preimage(f::AbsAlgAssMorGen, z)
  @assert parent(z) == codomain(f)
  v = matrix(FlintQQ, 1, dim(domain(f)), _coefficients_of_restricted_scalars(z)) * f.Minv
  return domain(f)(_eltseq(v))
end

# Write M_n(K) as M_n(Q) if [K : Q] = 1
# We use the "restricted scalar map" to modell M_n(Q) -> M_n(K)
function _as_full_matrix_algebra_over_Q(A::AlgMat{nf_elem})
  K = base_ring(A)
  @assert isabsolute(K) && degree(K) == 1
  B = matrix_algebra(FlintQQ, degree(A))

  M = identity_matrix(K, dim(B))
  Minv = identity_matrix(FlintQQ, dim(B))

  return B, AbsAlgAssMorGen(B, A, M, Minv)
end

################################################################################
#
#  Central primitive idempotents
#
################################################################################

# This is Corollary 2.1 of
# Eric Jespers, Guilherme Leal and Antonio Paques
# Central idempotents in the rational group algebra of a finite nilpotent group
# https://www.worldscientific.com/doi/10.1142/S0219498803000398
function _central_primitive_idempotents_abelian(A::AlgGrp)
  G = group(A)
  @assert base_ring(A) isa FlintRationalField
  @assert isabelian(G)
  S = subgroups(G, fun = (x, m) -> sub(x, m, false))
  o = one(A)
  idem = elem_type(A)[]
  push!(idem, 1//order(G) * sum(basis(A)))
  for (s, ms) in S
    Q, mQ = quo(G, ms, false)
    if !iscyclic(Q)
      continue
    end
    e = 1//(order(s)) * sum([A(ms(x)) for x in s])
    M = minimal_subgroups(Q, false)
    for (H, mH) in M
      U, mU = sub(G, append!([mQ\(mH(x)) for x in gens(H)], [ms(x) for x in gens(s)]))
      uhat = 1//(order(U)) * sum([A(mU(x)) for x in U])
      e = e * (o - uhat)
    end
    push!(idem, e)
  end
  return idem
end

function __decompose_abelian_group_algebra(A::AlgGrp)
  T = elem_type(base_ring(A))
  idems = _central_primitive_idempotents_abelian(A)
  res = Vector{Tuple{AlgAss{T}, morphism_type(AlgAss{T}, typeof(A))}}()
  for idem in idems
    S, StoA = subalgebra(A, idem, true)
    S.issimple = 1
    push!(res, (S, StoA))
  end
  return res
end

function decompose(A::AlgGrp)
  if isdefined(A, :decomposition)
    return A.decomposition
  end
  if group(A) isa GrpAbFinGen && (base_ring(A) isa FlintRationalField)
    res = __decompose_abelian_group_algebra(A)
    A.decomposition = res
    return res
  end
  G = group(A)
  res = __decompose(A)

  #if !isdefined(res[1][1], :isomorphic_full_matrix_algebra)
  #  if order(G) == 24 && find_small_group(G)[1] == (24, 12) &&
  #      base_ring(A) isa FlintRationalField
  #    @assert G.isfromdb
  #    _compute_matrix_algebras_from_reps(A, res, _reps[1])
  #  end
  #
  #  if order(G) == 48 && find_small_group(G)[1] == (48, 48) &&
  #      base_ring(A) isa FlintRationalField
  #    @assert G.isfromdb
  #    _compute_matrix_algebras_from_reps(A, res, _reps[2])
  #  end
  #end

  return res
end

function _compute_matrix_algebras_from_reps2(A, res)
  G = group(A)
  smallid, H, HtoG = find_small_group(G)
  idempotents = elem_type(A)[r[2](one(r[1])) for r in res]
  data = DefaultSmallGroupDB().db[smallid[1]][smallid[2]]
  Qx = Globals.Qx
  for j in data.galrep
    #@show j, data.schur[j]
    if data.schur[j] != 1
      continue
    end
    field, _ = number_field(Qx(data.fields[j]), "a", cached = false)
    d = data.dims[j]
    mats = dense_matrix_type(nf_elem)[ matrix(field, d, d, map(field, data.mod[j][k])) for k in 1:length(data.mod[j])]
    D = Tuple{GrpGenElem, dense_matrix_type(nf_elem)}[(H[H.gens[i]], mats[i]) for i in 1:length(H.gens)]
    op = (x, y) -> (x[1] * y[1], x[2] * y[2])
    id = (Hecke.id(H), identity_matrix(field, d))
    cl = closure(D, op, id)

    k0 = 0
    for k in 1:length(idempotents)
      c = coefficients(idempotents[k])
      z = zero_matrix(field, d, d)
      for (h, M) in cl
        i = A.group_to_base[HtoG(h)]
        z += field(c[i]) * M
      end
      if isone(z)
        k0 = k
        break
      end
    end

    #@show k0

    B, mB = res[k0]
    basisB = basis(B)

    MB = matrix_algebra(field, d)

    forward_matrix = zero_matrix(field, dim(B), dim(MB))

    back_matrix = zero_matrix(FlintQQ, dim(B), dim(B))

    BinMB = elem_type(MB)[]

    for i in 1:dim(B)
      img = MB(_evaluate_rep(mB(basisB[i]), d, cl, HtoG))
      push!(BinMB, img)
      elem_to_mat_row!(forward_matrix, i, img)
      v = _coefficients_of_restricted_scalars(img)
      for j in 1:length(v)
        back_matrix[i, j] = v[j]
      end
    end

    back_matrix = inv(back_matrix)

    # now comes the horror
    #
    #
    #@show back_matrix

    #v = matrix(FlintQQ, 1, dim(B), coefficients(rand(B, -10:10)))
    #@show v
    #@show matrix(FlintQQ, 1, dim(B), _coefficients_of_restricted_scalars(MB(collect(change_base_ring(field, v) * forward_matrix)))) * back_matrix

    #BtoMB = function(z)
    #  v = matrix(base_ring(MB), 1, dim(B), coefficients(z))
    #  return MB(collect(v * forward_matrix))
    #end

    #MBtoB = function(z)
    #  v = matrix(FlintQQ, 1, dim(B), _coefficients_of_restricted_scalars(z)) * back_matrix
    #  return B(collect(v))
    #end

    #println("Adding to $k0: $MB")
    f = AbsAlgAssMorGen(B, MB, forward_matrix, back_matrix)
    B.isomorphic_full_matrix_algebra = MB, f
  end
end

function _compute_matrix_algebras_from_reps(A, res, reps)
  G = group(A)
  idempotents = elem_type(A)[r[2](one(r[1])) for r in res]
  for j in 1:reps.n
    d = reps.dims[j]
    @assert length(reps.reps[j]) == length(G.gens)
    mats = fmpq_mat[ matrix(FlintQQ, d, d, reps.reps[j][k]) for k in 1:length(reps.reps[j])]
    D = Tuple{GrpGenElem, fmpq_mat}[(G[G.gens[i]], mats[i]) for i in 1:length(G.gens)]
    op = (x, y) -> (x[1] * y[1], x[2] * y[2])
    id = (Hecke.id(G), identity_matrix(FlintQQ, d))
    cl = closure(D, op, id)
    @assert length(cl) == order(G)
    k0 = 0
    for k in 1:length(idempotents)
      e = idempotents[k]
      c = coefficients(e)
      z = _evaluate_rep(e, d, cl)
      if isone(z)
        k0 = k
        break
      end
    end

    @assert k0 != 0

    B, mB = res[k0]

    @assert dim(B) == d^2

    basisB = basis(B)

    MB = matrix_algebra(FlintQQ, d)

    h = zero_matrix(FlintQQ, d^2, d^2)

    for i in 1:dim(B)
      img = MB(_evaluate_rep(mB(basisB[i]), d, cl))
      elem_to_mat_row!(h, i, img)
    end
    B.isomorphic_full_matrix_algebra = (MB, hom(B, MB, h, inv(h)))
  end
end

function _coefficients_of_restricted_scalars(x)
  A = parent(x)
  K = base_ring(A)
  m = dim(A)
  n = degree(K)
  nm = n * m
  y = Vector{fmpq}(undef, nm)
  yy = coefficients(x, copy = false)
  k = 1
  for i = 1:m
    for j = 1:n
      y[k] = coeff(yy[i], j - 1)
      k += 1
    end
  end
  return y
end

function _absolute_basis(A)
  m = dim(A)
  K = base_ring(A)
  n = degree(K)
  B = Vector{elem_type(A)}()
  bK = basis(K)
  for i in 1:n
    for j in 1:m
      v = Vector{elem_type(K)}(undef, m)
      for k in 1:m
        v[i] = zero(K)
      end
      v[j] = bK[j]
      push!(B, A(v))
    end
  end
  return B
end

function _evaluate_rep(el, d, rep)
  c = coefficients(el)
  A = parent(el)
  z = zero_matrix(FlintQQ, d, d)
  for (g, M) in rep
    i = A.group_to_base[g]
    z += c[i] * M
  end
  return z
end

function _evaluate_rep(el, d, rep, f)
  c = coefficients(el)
  A = parent(el)
  z = zero_matrix(base_ring(rep[1][2]), d, d)
  for (g, M) in rep
    i = A.group_to_base[f(g)]
    z += base_ring(M)(c[i]) * M
  end
  return z
end

elem_type(::Type{NfMorSet{AnticNumberField}}) = NfToNfMor

################################################################################
#
#  Opposite algebra
#
################################################################################

function opposite_algebra(A::AlgGrp{S, GrpAbFinGen, GrpAbFinGenElem}) where S
  G = group(A)
  z = zero_matrix(base_ring(A), dim(A), dim(A))
  for g in G
    i = A.group_to_base[g]
    elem_to_mat_row!(z, i, A(-g))
  end
  return A, hom(A, A, z, z)
end

function opposite_algebra(A::AlgGrp{S, GrpGen, GrpGenElem}) where S
  G = group(A)
  z = zero_matrix(base_ring(A), dim(A), dim(A))
  for g in G
    i = A.group_to_base[g]
    elem_to_mat_row!(z, i, A(inv(g)))
  end
  return A, hom(A, A, z, z)
end

################################################################################
#
#  Ad hoc methods
#
################################################################################

function isfree_s4_fabi(K::AnticNumberField)
  if istamely_ramified(K, fmpz(2))
    println("fabi 1")
    return true
  end

  P = prime_decomposition(maximal_order(K), 2)[1][1]

  D = decomposition_group(P)

  if length(D) == 24 && isweakly_ramified(K, P)
    println("fabi 2")
    return true
  end

  G, _, _ = generic_group(D, *, false)
  id, _, _ = find_small_group(G)
  if id == (12, 3) # A4
    println("fabi 3")
    return true
  end

  if id == (8, 3) # D4
    if isweakly_ramified(K, P)
      A, mA = automorphism_group(K)
      I, mI = inertia_subgroup(K, P, mA)
      fl = _isnormal([mI(i) for i in I])
      if order(I) == 4 && fl
        println("fabi 4")
        return true
      end
    end
  end

  return false
end

function isfree_a4_fabi(K::AnticNumberField)
  if istamely_ramified(K, fmpz(2))
    println("fabi 1")
    return true
  end

  P = prime_decomposition(maximal_order(K), 2)[1][1]

  D = decomposition_group(P)

  if length(D) == 12
    println("fabi 2")
    return true
  end

  return false
end

function isfree_a5_fabi(K::AnticNumberField)
  if !istamely_ramified(K, fmpz(2))
    println("fabi 1")
    return false
  end

  if !(istamely_ramified(K, fmpz(3)) || !isalmost_maximally_ramified(K, fmpz(3)))
    println("fabi 2")
    return false
  end

  if !(istamely_ramified(K, fmpz(5)) || !isalmost_maximally_ramified(K, fmpz(5)))
    println("fabi 3")
    return false
  end

  return true
end

function isalmost_maximally_ramified(K::AnticNumberField, p::fmpz)
  P = prime_decomposition(maximal_order(K), p)[1][1]
  G, mG = automorphism_group(K)
  D, mD = decomposition_group(K, P, mG) # this is the local Galois group
  ram_groups = []
  t = 1
  while true
    R, mR = ramification_group(K, P, t, mG)
    push!(ram_groups, [mR(r) for r in R])
    if order(R) == 1
      break
    end
    t += 1
  end

  QG = group_algebra(FlintQQ, G)
  A, KtoA = galois_module(K)
  I = KtoA(maximal_order(K))
  assOrd = right_order(I)

  S = subgroups(D)
  for (H, mH) in S
    HinG = [ mD(mH(h)) for h in H]
    for t in 1:length(ram_groups)
      if all(h in ram_groups[t] for h in HinG) && all(h in HinG for h in ram_groups[t + 1])
        eH = sum(QG(h) for h in HinG) * QG(1//Order(H))
        if !(eH in assOrd)
          return false
        end
      end
    end
  end

  return true
end

function hom(KG::AlgGrp, KH::AlgGrp, m::GrpGenToGrpGenMor)
  @assert base_ring(KG) === base_ring(KH)
  K = base_ring(KG)
  M = zero_matrix(K, dim(KG), dim(KH))
  for i in 1:dim(KG)
    j = KH.group_to_base[m(KG.base_to_group[i])]
    M[i, j] = one(K)
  end
  return hom(KG, KH, M)
end


text:
# Routines having to do with the normalizing function values on a grid with respect
# to a given measure.

normalizing_weights(grid, measure) = sqrt.(quadweights(grid, measure))

"""
Return a normalization for a sampling operator on a grid. The normalization ensures
a norm equivalence between the continuous norm of the function space associated
with the measure, and the discrete norm of the sampling vector.

In practice, this amounts to a diagonal scaling with the square roots of the
weights of a quadrature rule on the grid that is exact on the span of the dictionary
in the space.
"""
sampling_normalization(gb::GridBasis{T}, measure) where T =
	DiagonalOperator(gb, gb, sqrt.(quadweights(grid(gb), measure)))

ext: jl


text:
include("linear_steps.jl")

"""
Iterate through the layers in the net to get a prediction.
"""
function forward_propagate_model_weights(DMatrix, parameters)
    master_cache = []
    A = DMatrix
    L = trunc(Int, length(parameters) / 2)

    # Forawrd propagate until the last layer
    for l = 1:(L-1)
        A_prev = A
        A, cache = linear_forward_activate(A_prev, parameters[string("W_", (l))], parameters[string("b_", (l))], activation_function = "tanh")
        # Debugging println statements
        #==if l == L-1
            println(string(l, " : ", A[1:2]))
        end # if==#
        push!(master_cache, cache)
    end # for

    # For the final layer (outputs)
    Ŷ, cache = linear_forward_activate(A, parameters[string("W_", (L))], parameters[string("b_", (L))], activation_function = "softmax")
    # Make sure the outputs are within bounds too
    push!(master_cache, cache)

    return Ŷ, master_cache
end # function

max_stars_repo_name: Squalm/NeuralNetworks

max_issues_repo_name: Squalm/NeuralNetworks

max_forks_repo_name: Squalm/NeuralNetworks


text:
module StructuredDistances

using Distances

# Write your package code here.

include("utils.jl")
include("threaded_pairwise.jl")
include("distances/helpers.jl")

# Path distances
include("distances/paths.jl")

# Multiset distances
include("distances/multisets/matching_distances.jl")
include("distances/multisets/matching_distances_generalised.jl")
include("distances/multisets/matching_distances_eval.jl")
include("distances/multisets/earth_movers_distance.jl")

# Sequence distances
include("distances/sequences/edit_distance.jl")
include("distances/sequences/dtw_distance.jl")

# Basic set distances (Hamming + Jaccard)
include("distances/set_distances.jl")

# Graph distances
include("distances/graph_distances.jl")

# Normalised distances (normalises any metric)
include("distances/normalised.jl")

# Plotting recipies
include("plot_recipes/seq_plot.jl")
include("plot_recipes/pathseq_plot.jl")

end


text:
export iauP06e
"""
Precession angles, IAU 2006, equinox based.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  canonical models.

Given:
   date1,date2   double   TT as a 2-part Julian Date (Note 1)

Returned (see Note 2):
   eps0          double   epsilon_0
   psia          double   psi_A
   oma           double   omega_A
   bpa           double   P_A
   bqa           double   Q_A
   pia           double   pi_A
   bpia          double   Pi_A
   epsa          double   obliquity epsilon_A
   chia          double   chi_A
   za            double   z_A
   zetaa         double   zeta_A
   thetaa        double   theta_A
   pa            double   p_A
   gam           double   F-W angle gamma_J2000
   phi           double   F-W angle phi_J2000
   psi           double   F-W angle psi_J2000

Notes:

   1. The TT date date1+date2 is a Julian Date, apportioned in any
   convenient way between the two arguments.  For example,
   JD(TT)=2450123.7 could be expressed in any of these ways,
   among others:

         date1          date2

      2450123.7           0.0       (JD method)
      2451545.0       -1421.3       (J2000 method)
      2400000.5       50123.2       (MJD method)
      2450123.5           0.2       (date & time method)

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable.  The J2000 method is best matched to the way
   the argument is handled internally and will deliver the
   optimum resolution.  The MJD method and the date & time methods
   are both good compromises between resolution and convenience.

   2. This function returns the set of equinox based angles for the
   Capitaine et al. "P03" precession theory, adopted by the IAU in
   2006.  The angles are set out in Table 1 of Hilton et al. (2006):

   eps0   epsilon_0   obliquity at J2000.0
   psia   psi_A       luni-solar precession
   oma    omega_A     inclination of equator wrt J2000.0 ecliptic
   bpa    P_A         ecliptic pole x, J2000.0 ecliptic triad
   bqa    Q_A         ecliptic pole -y, J2000.0 ecliptic triad
   pia    pi_A        angle between moving and J2000.0 ecliptics
   bpia   Pi_A        longitude of ascending node of the ecliptic
   epsa   epsilon_A   obliquity of the ecliptic
   chia   chi_A       planetary precession
   za     z_A         equatorial precession: -3rd 323 Euler angle
   zetaa  zeta_A      equatorial precession: -1st 323 Euler angle
   thetaa theta_A     equatorial precession: 2nd 323 Euler angle
   pa     p_A         general precession
   gam    gamma_J2000 J2000.0 RA difference of ecliptic poles
   phi    phi_J2000   J2000.0 codeclination of ecliptic pole
   psi    psi_J2000   longitude difference of equator poles, J2000.0

   The returned values are all radians.

   3. Hilton et al. (2006) Table 1 also contains angles that depend on
   models distinct from the P03 precession theory itself, namely the
   IAU 2000A frame bias and nutation.  The quoted polynomials are
   used in other SOFA functions:

   . iauXy06  contains the polynomial parts of the X and Y series.

   . iauS06  contains the polynomial part of the s+XY/2 series.

   . iauPfw06  implements the series for the Fukushima-Williams
      angles that are with respect to the GCRS pole (i.e. the variants
      that include frame bias).

   4. The IAU resolution stipulated that the choice of parameterization
   was left to the user, and so an IAU compliant precession
   implementation can be constructed using various combinations of
   the angles returned by the present function.

   5. The parameterization used by SOFA is the version of the Fukushima-
   Williams angles that refers directly to the GCRS pole.  These
   angles may be calculated by calling the function iauPfw06.  SOFA
   also supports the direct computation of the CIP GCRS X,Y by
   series, available by calling iauXy06.

   6. The agreement between the different parameterizations is at the
   1 microarcsecond level in the present era.

   7> When constructing a precession formulation that refers to the GCRS
   pole rather than the dynamical pole, it may (depending on the
   choice of angles) be necessary to introduce the frame bias
   explicitly.

   8. It is permissible to re-use the same variable in the returned
   arguments.  The quantities are stored in the stated order.

Reference:

   Hilton, J. et al., 2006, Celest.Mech.Dyn.Astron. 94, 351

Called:
   iauObl06     mean obliquity, IAU 2006

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauP06e(double date1, double date2,
#              double *eps0, double *psia, double *oma, double *bpa,
#              double *bqa, double *pia, double *bpia,
#              double *epsa, double *chia, double *za, double *zetaa,
#              double *thetaa, double *pa,
#              double *gam, double *phi, double *psi)

function iauP06e(date1::Real, date2::Real)
   ref_eps0   = Ref{Float64}(0.0)
   ref_psia   = Ref{Float64}(0.0)
   ref_oma    = Ref{Float64}(0.0)
   ref_bpa    = Ref{Float64}(0.0)
   ref_bqa    = Ref{Float64}(0.0)
   ref_pia    = Ref{Float64}(0.0)
   ref_bpia   = Ref{Float64}(0.0)
   ref_epsa   = Ref{Float64}(0.0)
   ref_chia   = Ref{Float64}(0.0)
   ref_za     = Ref{Float64}(0.0)
   ref_zetaa  = Ref{Float64}(0.0)
   ref_thetaa = Ref{Float64}(0.0)
   ref_pa     = Ref{Float64}(0.0)
   ref_gam    = Ref{Float64}(0.0)
   ref_phi    = Ref{Float64}(0.0)
   ref_psi    = Ref{Float64}(0.0)

   ccall((:iauP06e, libsofa_c), Cvoid,
       (Cdouble, Cdouble,
       Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble},
       Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble},
       Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble},
       Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}),
       convert(Float64, date1), convert(Float64, date2),
       ref_eps0, ref_psia, ref_oma, ref_bpa, ref_bqa, ref_pia,
       ref_bpia, ref_epsa, ref_chia, ref_za, ref_zetaa, ref_thetaa,
       ref_pa, ref_gam, ref_phi, ref_psi)

   return ref_eps0[], ref_psia[], ref_oma[], ref_bpa[], ref_bqa[],
          ref_pia[], ref_bpia[], ref_epsa[], ref_chia[], ref_za[],
          ref_zetaa[], ref_thetaa[], ref_pa[], ref_gam[], ref_phi[],
          ref_psi[]
end


text:
using Dierckx
using Plots
using RollingFunctions
using LinearAlgebra

#using PyPlot
include("heaviside.jl")
include("delta.jl")
include("angle_of_incidence.jl")
include("kpp_psim.jl")
include("kpp_psis.jl")
include("bulk.jl")
include("holtslag.jl")
include("holtslag_psim.jl")
include("holtslag_psis.jl")
include("large.jl")

cpl=1;
moist=1;

dt=60; #seconds
ndays=60; # duration of simulation in days

n=Int(round(ndays*24*60*60/dt)); # number of time steps

H=50; # atmospheric layer thickness
dz=1; # oceanic layer thickness
vlevo=50; # number of oceanic vertical levels
vleva=50; # number of atmospheric vertical levels
rhow0=1029; # kg/m^3
cp=4184; #J/K/kg
cpa=1005; #J/K/kg
sb = 5.670e-8;
d2k=273.15;
Av=2500000;
gamma_blk=0.01;
cvapor_fac=640380;
cvapor_exp=5107.4;
cvapor_exp_ice=5897.8;
saltsat=0.980;
p0=101300; #reference pressure
gravity_mks=9.81;
humid_fac=0.606;
tref=25; #reference temperature
sref=25;
SO=sref; #salinity not implemented
alpha=2e-4; #expansion coefficients for temperature
beta=7.4e-4; #expansion coefficients for salt
emissivity=0.97; # ocean emissivity
kappa=1e-5; # m^2/kg air absorption coefficient
Rgas=287.05;
karman=0.4;

ZA=H*[1:vleva;];
ZO=zeros(vlevo+1);
ZO[2:vlevo+1]=-[1:vlevo;]*dz;
kT0=ones(vlevo)*1e-4; # SST vertical mixing
kT=kT0; # SST vertical mixing

Aimp=0.5; # Atmospheric implicitness
Oimp=0.5; # Oceanic implicitness

R=0.62; g1=0.6; g2=20;

f=R*exp.(ZO/g1)+(1-R)*exp.(ZO/g2);
f[ZO.<-200].=0;
f=-diff(f,dims=1);

albedo=0.06;
lat=30;
lon=0;
t=[1:n;]*dt/3600;
jstart=70;# julian day, borial summer
jadd=round.((t.-t[1])/24);
jd=jstart.+jadd;

AOI=angle_of_incidence.(lat,lon,jd,t); AOI[AOI.<0].=0;
SW=(1362*(1-albedo).*f'.*AOI);
#mean(SW(:,1))
#SW=(repmat((sin(2*pi/n*(1:n)*ndays+pi/2)+1)/2*1362*(1-albedo),vlev,1).*
#repmat(f',1,n))'.*cos(lat/180*pi);

TO=zeros(n+1,vlevo);
UO=zeros(n+1,vlevo);
SST=zeros(n+1);

UA=zeros(n+1,vleva);
qA=zeros(n+1,vleva);
THETA=zeros(n+1,vleva);
rhoa=zeros(n+1,vleva);
TA=zeros(n+1,vleva)

UOst=zeros(vlevo); TOst=zeros(vlevo);
UAst=zeros(vleva); qAst=zeros(vleva); THETAst=zeros(vleva);

############################################################################
# Initial conditions
############################################################################
#U0=ones(n,1)*5;
#THETA0=(sin(2*pi/n*(1:n)*ndays/5-pi/2)+21);
TA[1,:]=25 .-0.007.*H.*[1:vleva;];
P=zeros(vleva);
P[1]=101300;
P[2:vleva]=P[1].*exp.(-cumsum(H./(Rgas.*(TA[1,2:end].+d2k)./gravity_mks)));
rhoa[1,:]=P./(Rgas.*(TA[1,:].+d2k));
THETA[1,:]=(TA[1,:].+d2k).*(P[1]./P).^(2/7).-d2k;
THETA[1,1:20].=THETA[1,1];
THETA0=ones(n,1)*THETA[1,end].+1;#+gamma_blk*H+0.01;#+0.51;#+SW*kappa/cpa*dt;
T0=ones(n+1,vlevo).*26;#repmat((16+(vlev-1:-1:0)/vlev*8),n+1,1);
TO[:,:].=T0;
TO[:,1]=T0[:,1]; TO[1,:]=T0[1,1].+15 .*([vlevo:-1:1;]/vlevo.-1);
SST.=TO[:,1];
ssq = saltsat*cvapor_fac*exp(-cvapor_exp./(TA[1]+d2k))./rhoa[1,1];
RH=0.8;
q0=0.5*ssq;
#Tv=(THETA(i-1,:)+d2k).*(1+humid_fac*q(i-1));
qA[1,:]=saltsat*cvapor_fac*exp.(-cvapor_exp./(TA[1,:].+d2k))./rhoa[1,:].*RH;
UA[1,:]=(0.2/karman).*(log.((ZA.-H/2)./1e-4));
#U(1,:)=U(1,end);
U0=(0.2/karman).*(log((ZA[end]+H/2)/1e-4)).*ones(1,n+1);
atemp=THETA[1,1]+d2k; aqh=qA[1,1];
speed=abs(UA[1,1]); sst=SST[1];
tmp0=bulk(atemp,aqh,speed,sst,H/2,H/2,H/2,H/2,rhoa[1,1]);
huol0=tmp0[8];
TAup=270;
#U0=(sin(2*pi/n*(1:n)*ndays/5+2*pi/n*3600/dt*sf)*8+U(1,end));
#UO(:)=0.;
#U0(3600*24*5/dt:end)=U(1,end); one pertubation
############################################################################

LH=zeros(1,n); SH=zeros(1,n); TAU=zeros(1,n); TAUO=zeros(1,n);
LW=zeros(1,n); E=zeros(1,n); ssq=zeros(1,n);
xi=zeros(1,n); ce=zeros(1,n); ch=zeros(1,n);
ustar=zeros(1,n); qstar=zeros(1,n); tstar=zeros(1,n);
psimh=zeros(1,n); psixh=zeros(1,n); Ri=zeros(n,vleva);
rd=zeros(1,n); rh=zeros(1,n); re=zeros(1,n);
Qnet=zeros(1,n);
KAm=zeros(n,vleva);KAt=zeros(n,vleva);
KOm=zeros(n,vlevo);KOt=zeros(n,vlevo);

for i=2:n+1
    println(i)
    if rem(i,1440*60/dt)==0
        println(string("Day  ",round(i/(1440*60/dt),digits=5)))
    end

    ########################################################################
    # surface
    ########################################################################
    hl,hs,_,ce[i-1],ch[i-1],_,ssq[i-1],xi[i-1],rd[i-1],
    rh[i-1],re[i-1],ustar[i-1],qstar[i-1],tstar[i-1],psimh[i-1],psixh[i-1]=
    bulk(TA[i-1,1]+d2k,qA[i-1,1],abs(UA[i-1,1]-UO[i-1,1]),SST[i-1,1],
        H/2,H/2,H/2,10,rhoa[i-1,1]);
    LH[i-1]=-hl;
    SH[i-1]=-hs;
    E[i-1]=-hl/Av;
    ϵa=0.725+0.17*log10(sum(qA[1,:]/1000*1.2*H)*100)
    LW[i-1]=emissivity*sb*((SST[i-1]+d2k).^4 .- ϵa.*(TAup).^4);
    TAU[i-1]=rhoa[i-1,1].*ustar[i-1].^2 .*sign(ustar[i-1]);
    Qnet[i-1]=SW[i-1,1]-LH[i-1]-SH[i-1]-LW[i-1];
    ########################################################################
    # Atmospheric vertical diffusion
    ########################################################################

    Tv=(THETA[i-1,:].+d2k).*(1 .+humid_fac*qA[i-1,:]);
    KAm[i-1,:], KAt[i-1,:]= holtslag(Tv,qA[i-1,1],UA[i-1,:],
                 UA[i-1,:].*0,LH[i-1],SH[i-1],ustar[i-1],[0; ZA]);
    kam=KAm[i-1,:];
    kat=KAt[i-1,:];

    ########################################################################
    # Oceanic vertical diffusion
    ########################################################################

    rhow=rhow0.*(1 .-alpha.*(TO[i-1,:].-tref).+beta.*(SO.-sref));
    B=gravity_mks.*(rhow0.-rhow)./rhow0;
    Uo=(vlevo:-1:1)/vlevo/10; Vo=zeros(size(Uo)); FWflux=zeros(size(Qnet));
    KOm[i-1,:],KOt[i-1,:],_=
        large(TO[i-1,:],TO[i-1,:].*0,UO[i,:],Vo,rhow,Qnet[i-1]
             ,0,ustar[i-1],ZO,lat);
    kot=KOt[i-1,:];
    kom=KOm[i-1,:];

    ########################################################################
    # Atmospheric explicit time stepping
    ########################################################################

    TAUS=sign(UA[i-1,1]-UO[i-1,1]);

    UAst[1]=UA[i-1,1]-TAU[i-1]*TAUS/(H*rhoa[i-1,1])*dt-(1-Aimp)*(
        kam[1].*(UA[i-1,1]-UA[i-1,2])/H^2*dt);
    UAst[end]=UA[i-1,end]-(1-Aimp)*(
        kam[end]*(UA[i-1,end]-U0[i-1])/H^2*dt);
    UAst[2:end-1]=UA[i-1,2:end-1]+(1-Aimp)*(
        kam[1:end-2].*(UA[i-1,1:end-2]-UA[i-1,2:end-1])/H^2*dt-
        kam[2:end-1].*(UA[i-1,2:end-1]-UA[i-1,3:end  ])/H^2*dt);

    qAst[1]=qA[i-1,1]-E[i-1]/(H*rhow[1])*dt-
        (1-Aimp)*(kat[1]*(qA[i-1,1]-qA[i-1,2])/H^2*dt);
    qAst[end]=qA[i-1,end]-(1-Aimp)*(kat[end]*(qA[i-1,end]-q0)/H^2*dt);
    qAst[2:end-1]=qA[i-1,2:end-1]+(1-Aimp)*(
        kat[1:end-2].*(qA[i-1,1:end-2]-qA[i-1,2:end-1])/H^2*dt-
        kat[2:end-1].*(qA[i-1,2:end-1]-qA[i-1,3:end  ])/H^2*dt);

    THETAst[1]=THETA[i-1,1]+(SH[i-1])/cpa/(H*rhoa[i-1,1])*dt-
         (1-Aimp)*(kat[1]*(THETA[i-1,1]-THETA[i-1,2])/H^2*dt);
    THETAst[end]=THETA[i-1,end]-
         (1-Aimp)*(kat[end]*(THETA[i-1,end]-THETA0[i-1])/H^2*dt);
    THETAst[2:end-1]=THETA[i-1,2:end-1]+(1-Aimp)*(
         kat[1:end-2].*(THETA[i-1,1:end-2]-THETA[i-1,2:end-1])/H^2*dt-
         kat[2:end-1].*(THETA[i-1,2:end-1]-THETA[i-1,3:end  ])/H^2*dt);

    ########################################################################
    # Oceanic explicit time stepping
    ########################################################################

    if cpl==1

        TOst[1]=TO[i-1,1]+(SW[i-1,1]-LH[i-1]-SH[i-1]-LW[i-1])/cp/(dz*rhow[1])*dt-
            (1-Oimp)*(kot[1]*(TO[i-1,1]-TO[i-1,2])/dz^2*dt);
        TOst[end]=TO[i-1,end]+SW[i-1,end]/cp/(dz*rhow[end])*dt+
            (1-Oimp)*(kot[end-1].*(TO[i-1,end-1]-TO[i-1,end])/dz^2*dt);
        TOst[2:end-1]=TO[i-1,2:end-1]+(1-Oimp)*(
            kot[1:end-2].*(TO[i-1,1:end-2]-TO[i-1,2:end-1])/dz^2*dt-
            kot[2:end-1].*(TO[i-1,2:end-1]-TO[i-1,3:end  ])/dz^2*dt)+
            SW[i-1,2:end-1]./cp./(dz*rhow[2:end-1])*dt;

        UOst[1]=UO[i-1,1]+
             (TAU[i-1])*TAUS/(rhow[1])/dz*dt-
             (1-Oimp)*(kom[1]*(UO[i-1,1]-UO[i-1,2])/dz^2*dt);
        UOst[end]=UO[i-1,end]+
             (1-Oimp)*(kom[end-1].*(UO[i-1,end-1]-UO[i-1,end])/dz^2*dt);
        UOst[2:end-1]=UO[i-1,2:end-1]+(1-Oimp)*(
             kom[1:end-2].*(UO[i-1,1:end-2]-UO[i-1,2:end-1])/dz^2*dt-
             kom[2:end-1].*(UO[i-1,2:end-1]-UO[i-1,3:end  ])/dz^2*dt);

    end

    if Aimp==0
        UA[i,:]=UAst; qA[i,:]=qAst; THETA[i,:]=THETAst;
    else
        du=-Aimp*kam[1:end-1].*dt./H^2; dl=-Aimp*kam[2:end].*dt./H^2;
        d=(1 .+Aimp.*dt./H^2 .*([kam[1:end-1]; 0]+[0; kam[2:end]]));
        A=inv(Tridiagonal(dl,d,du));
        UA[i,:]=A*UAst;
        du=-Aimp*kat[1:end-1].*dt./H^2; dl=-Aimp*kat[2:end].*dt./H^2;
        d=(1 .+Aimp.*dt./H^2 .*([kat[1:end-1]; 0]+[0; kat[2:end]]));
        A=inv(Tridiagonal(dl,d,du));
        qA[i,:]=A*qAst; THETA[i,:]=A*THETAst;
    end

    if Oimp==0
        UO[i,:]=UOst; TO[i,:]=TOst;
    else
        du=-Oimp*kom[1:end-1].*dt./dz^2; dl=-Oimp*kom[2:end].*dt./dz^2;
        d=(1 .+Oimp.*dt./dz^2 .*([kom[1:end-1]; 0]+[0; kom[2:end]]));
        A=inv(Tridiagonal(dl,d,du));
        UO[i,:]=A*UOst;
        du=-Oimp*kot[1:end-1].*dt./dz^2; dl=-Oimp*kot[2:end].*dt./dz^2;
        d=(1 .+Oimp.*dt./dz^2 .*([kot[1:end-1]; 0]+[0; kot[2:end]]));
        A=inv(Tridiagonal(dl,d,du));
        TO[i,:]=A*TOst;
    end
    SST[i]=TO[i,1]

    TA[i,:]=(THETA[i,:].+d2k).*(P./P[1]).^(2/7).-d2k;
    rhoa[i,:]=P./(Rgas.*(TA[i,:].+d2k));

    if moist==1
        qsat=saltsat*cvapor_fac*exp.(-cvapor_exp./(TA[i,:].+d2k))./rhoa[i-1,:];
        dq=(qA[i,:].-qsat).*((qA[i,:].-qsat).>0);
        dq=min.(dq,qA[i,:]);
        dTHETA=Av.*dq./cp; #[J/kg][j-1 kg K]
        #println([i maximum(dTHETA) maximum(dq)])
        ind=qA[i,:].>qsat;
        qA[i,ind]=qA[i,ind]-dq[ind];
        #THETA[i,:]=THETA[i,: ]+dTHETA;
    end

    if rem(i,3600*6/dt)==0;#3600*6/dt
      println(i)
      f1=plot(-ZO[1:end-1],TO[i,:],ylabel=("TO"));#,ylims=(19,21)
      f2=plot(-ZO[1:end-1],UO[i,:],ylabel=("UO"));
      f3=plot(-ZO[1:end-1],kom,ylabel=("KOm"));
      f4=plot(ZA,UA[i,:],ylabel=("UA"));
      f5=plot(ZA,Tv.-273.15,ylabel=("Tv")); #ylim([19 21]);...
      f6=plot(ZA,qA[i,:],ylabel=("q"));
      f7=plot(ZA,kam,ylabel=("KAm"));
      plt=plot(f1,f2,f3,f4,f5,f6,f7,layout=(7,1),legend=false,titlefontsize=6)
      display(plt)
    end

end

k1=plot(rollmean(U0[1,:],Int(24*60*60/dt)),
              xticks=(0:24*60*60/dt:n,string.(0:ndays)),
              ylabel="U0");
k2=plot(rollmean(UA[:,end],Int(24*60*60/dt)),
              xticks=(0:24*60*60/dt:n,string.(0:ndays)),
              ylabel="Utop");
k3=plot(rollmean(UA[:,end]-U0[1,:],Int(24*60*60/dt)),
              xticks=(0:24*60*60/dt:n,string.(0:ndays)),
              ylabel="Utop-U0");
k4=plot(rollmean(UA[:,1],Int(24*60*60/dt)),
              xticks=(0:24*60*60/dt:n,string.(0:ndays)),
              ylabel="U");
k5=plot(rollmean(THETA[:,2],Int(24*60*60/dt)),
              xticks=(0:24*60*60/dt:n,string.(0:ndays)),
              ylabel="THETA");
k6=plot(rollmean(SST,Int(24*60*60/dt)),
              xticks=(0:24*60*60/dt:n,string.(0:ndays)),
              ylabel="SST");
k7=plot(rollmean(TO[:,10],Int(24*60*60/dt)),
              xticks=(0:24*60*60/dt:n,string.(0:ndays)),
              ylabel="T");
#tt=string("DZ=",dz)
plt=plot(k1,k2,k3,k4,k5,k6,k7,layout=(7,1),legend=false)
display(plt)

#PyPlot.suptitle(tt)

k1=plot(rollmean(xi[1,:],Int(24*60*60/dt)),
              xticks=(0:24*60*60/dt:n,string.(0:ndays)),
              title="xi");
k2=plot(rollmean(ustar[1,:],Int(24*60*60/dt)),
              xticks=(0:24*60*60/dt:n,string.(0:ndays)),
              title="ustar");
k3=plot(rollmean(psimh[1,:],Int(24*60*60/dt)),
              xticks=(0:24*60*60/dt:n,string.(0:ndays)),
              title="psimh");
k4=plot(rollmean(psixh[1,:],Int(24*60*60/dt)),
              xticks=(0:24*60*60/dt:n,string.(0:ndays)),
              title="psixh");
k5=plot(rollmean((rd[1,:].^2),Int(24*60*60/dt)),
              xticks=(0:24*60*60/dt:n,string.(0:ndays)),
              title="CD");
k6=plot(rollmean((rh[1,:].*rd[1,:]),Int(24*60*60/dt)),
              xticks=(0:24*60*60/dt:n,string.(0:ndays)),
              title="CH");
k7=plot(rollmean((re[1,:].*rd[1,:]),Int(24*60*60/dt)),
              xticks=(0:24*60*60/dt:n,string.(0:ndays)),
              title="CE");
#tt=string("DZ=",dz)
plt=plot(k1,k2,k3,k4,k5,k6,k7,layout=(7,1),legend=false)
display(plt)

#PyPlot.suptitle(tt)

k1=plot(rollmean(SW[:,1]/cp/(dz*rhow0)*dt,Int(24*60*60/dt)),
              xticks=(0:24*60*60/dt:n,string.(0:ndays)),
              title="SW");
k2=plot(rollmean(LW[1,:]/cp/(dz*rhow0)*dt,Int(24*60*60/dt)),
              xticks=(0:24*60*60/dt:n,string.(0:ndays)),
              title="LW");
k3=plot(rollmean(SH[1,:]/cp/(dz*rhow0)*dt,Int(24*60*60/dt)),
              xticks=(0:24*60*60/dt:n,string.(0:ndays)),
              title="SH");
k4=plot(rollmean(LH[1,:]/cp/(dz*rhow0)*dt,Int(24*60*60/dt)),
              xticks=(0:24*60*60/dt:n,string.(0:ndays)),
              title="LH");
k5=plot(rollmean((SW[:,1]-LW[1,:]-SH[1,:]-LH[1,:])/cp/(dz*rhow0)*dt,Int(24*60*60/dt)),
              xticks=(0:24*60*60/dt:n,string.(0:ndays)),
              title="QNET");
plt=plot(k1,k2,k3,k4,k5,layout=(5,1),legend=false);
display(plt)

#PyPlot.suptitle(tt)


text:
# guide src: https://www.kaggle.com/c/street-view-getting-started-with-julia/overview/knn-tutorial

# Main idea: implement k-Nearest Neighbor algorithm with Leave-One-Out-Fold
    # cross validation.

# Leave-One-Out-Fold-Cross-Validation (LOOF-CV) is similar to k-fold CV,
    # but k is set to be equal to the number of training points.
    # The model is tested on each individual data point after being
    # trained with the remaining points in the data.
    # Because LOOF-CV uses all but one of the data points for training,
    # it is not biased to any particular point or split, and could arguably give
    # a better estimation for the performance of the model.

# todo: @time comparison:
#      transposed and NON transposed data
#      vectorized and for loop for euclidian distance

# Parallelization use
using Distributed

processes_cnt = 7
addprocs(processes_cnt)  # 8 cores

# Macro @everywhere - make function arguments, the definition of the functions
#    available to each process.
# Macro @distributed [reducer] - make a parallel for-loop

;

# Load packages
using Images
using DataFrames
using CSV
using DecisionTree(vcat)
using Statistics

# Helping function
showln(x) = (show(x); println())

showln("Libraries loaded")

;

# Load datasets (see julia_guide_kaggle.ipynb for details)
function read_image(dataset_type, labels_info, image_size, dataset_dir_path)
    # Initialize X matrix
    X = zeros(size(labels_info, 1), image_size)
    for (idx, image_id) in enumerate(labels_info[:ID])
        # Load image
        image_full_path = "$(dataset_dir_path)/$(dataset_type)Resized/$(image_id).Bmp"
        image = load(image_full_path)
        image_greyscale = Gray.(image)
        # Transform image matrix to a vector and store it in data matrix
        X[idx, :] = reshape(image_greyscale, 1, image_size)
    end
    return X
end

image_size = 20 * 20
dataset_dir_path = "/home/max/Documents/ai/julia/julia_kaggle/data"

train_labels_info = CSV.read("$(dataset_dir_path)/trainLabels.csv")
test_labels_info = CSV.read("$(dataset_dir_path)/sampleSubmission.csv")
X_train = read_image("train", train_labels_info, image_size, dataset_dir_path)
X_test = read_image("test", test_labels_info, image_size, dataset_dir_path)
y_train = map(col_value -> col_value[1], train_labels_info[:Class])
y_train = convert(Array{Int64, 1}, y_train)

showln((size(X_train), size(X_test), size(y_train)))

showln("Loaded X_train, X_test, y_train")

;

# Prepare data for kNN algorithm

# Transpose X_train and X_test matrices
# In Julia, iteration over columns is faster than iteration over rows.
# Now, each column == one image file
X_train_tp = X_train'
X_test_tp = X_test'

showln((size(X_train_tp), size(X_test_tp)))

;

# Define distance function for kNN algorithm between 2 data points == 2 images

# NOTE: In Julia, 'for' loops can be faster than vectorized operations

@everywhere function euclidean_distance_vectorized(a, b)
    return dot(a-b, a-b)
end

@everywhere function euclidean_distance_forloop(a, b)
    distance = 0.0
    for i in 1:size(a, 1)
        distance += (a[i] - b[i]) * (a[i] - b[i])
    end
    return distance
end

showln("Created vectorized and forloop functions to calculate euclidean distance")

;

# Define kNN function

# NOTE: in Julia, creating an empty vector and filling it with each
# element at a time is faster than copying the entire vector at once.

# NOTE: Since the code calculates the distance between the i data point
# and all the points in the training data, the closest point to
# the i point is itself with a distance of zero. Hence, we EXCLUDE it and
# select the next k points.

# Find k nearest neighbors of data point i
@everywhere function get_k_nearest_neighbors(X, i, k)
    n_rows, n_cols = size(X)
    ith_image = Array{Float32}(undef, n_rows)
    for i_idx in 1:n_rows
        ith_image[i_idx] = X[i_idx, i]
    end
    distances = Array{Float32}(undef, n_cols) # empty vector that will contain distances between i data point and each data point in the X matrix.
    jth_image = Array{Float32}(undef, n_rows)  # empty vector that will contain the j data point.
    for j_idx in 1:n_cols
        for i_idx in 1:n_rows
            jth_image[i_idx] = X[i_idx, j_idx]
        end
        # Calculate the distance and save the result
        distances[j_idx] = euclidean_distance_forloop(ith_image, jth_image)
    end
    sorted_neighbors = sortperm(distances)
    nearest_neighbors = sorted_neighbors[2:k+1]  # note: start with the 2nd one
    return nearest_neighbors
end

showln("Created function to get k nearest neighbors")

;

# Assign label to the ith point according to the labels of the k nearest neighbors

# training data: X matrix
# labels: y vector
@everywhere function assign_label(X, y, k, i)
    nearest_neighbors = get_k_nearest_neighbors(X, i, k)
    labels_counts = Dict{Int, Int}()
    highest_count = 0
    most_popular_label = 0
    for neighbor in nearest_neighbors
        neighbor_label = y[neighbor]
        if !haskey(labels_counts, neighbor_label)
            labels_counts[neighbor_label] = 0
        end
        labels_counts[neighbor_label] += 1
        if labels_counts[neighbor_label] > highest_count
            highest_count = labels_counts[neighbor_label]
            most_popular_label = neighbor_label
        end
    end
    return most_popular_label
end

showln("Created function to assign label using kNN")

;

# Assign label to each point in the training data
k = 3
y_predictions = @distributed (vcat) for i in 1:size(y_train, 1)
    assign_label(X_train_tp, y_train, k, i)
end

showln(size(y_predictions))

display(size(y_predictions))
display(size(y_train))

;

# Measure accuracy of the model by comparing our predictions with the true labels
loof_cv_accuracy = mean(y_predictions .== y_train)
showln("LOOF-CV accuracyy of 1-NN is $(loof_cv_accuracy)")  # k=1: 0.444; k=3: 0.445; k=5: 0.428

;

max_stars_repo_name: honchardev/Fun

max_issues_repo_name: honchardev/Fun

max_forks_repo_name: honchardev/Fun


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

"""
    treewalk(f::Function, tree::GitTree, payload=Any[], post::Bool=false)

Traverse the entries in `tree` and its subtrees in post or pre order. Preorder
means beginning at the root and then traversing the leftmost subtree (and
recursively on down through that subtree's leftmost subtrees) and moving right
through the subtrees. Postorder means beginning at the bottom of the leftmost
subtree, traversing upwards through it, then traversing the next right subtree
(again beginning at the bottom) and finally visiting the tree root last of all.

The function parameter `f` should have following signature:

    (Cstring, Ptr{Void}, Ptr{Void}) -> Cint

A negative value returned from `f` stops the tree walk. A positive value means
that the entry will be skipped if `post` is `false`.
"""
function treewalk(f::Function, tree::GitTree, payload=Any[], post::Bool = false)
    cbf = cfunction(f, Cint, Tuple{Cstring, Ptr{Void}, Ptr{Void}})
    cbf_payload = Ref{typeof(payload)}(payload)
    @check ccall((:git_tree_walk, :libgit2), Cint,
                  (Ptr{Void}, Cint, Ptr{Void}, Ptr{Void}),
                   tree.ptr, post, cbf, cbf_payload)
    return cbf_payload
end

repository(tree::GitTree) = tree.owner
repository(te::GitTreeEntry) = repository(te.owner)

"""
    filename(te::GitTreeEntry)

Return the filename of the object on disk to which `te` refers.
"""
function filename(te::GitTreeEntry)
    str = ccall((:git_tree_entry_name, :libgit2), Cstring, (Ptr{Void},), te.ptr)
    str != C_NULL && return unsafe_string(str)
    return nothing
end

"""
    filemode(te::GitTreeEntry) -> Cint

Return the UNIX filemode of the object on disk to which `te` refers as an integer.
"""
function filemode(te::GitTreeEntry)
    return ccall((:git_tree_entry_filemode, :libgit2), Cint, (Ptr{Void},), te.ptr)
end

"""
    entrytype(te::GitTreeEntry)

Return the type of the object to which `te` refers. The result will be
one of the types which [`objtype`](@ref) returns, e.g. a `GitTree` or `GitBlob`.
"""
function entrytype(te::GitTreeEntry)
    otype = ccall((:git_tree_entry_type, :libgit2), Cint, (Ptr{Void},), te.ptr)
    return objtype(Consts.OBJECT(otype))
end

"""
    entryid(te::GitTreeEntry)

Return the [`GitHash`](@ref) of the object to which `te` refers.
"""
function entryid(te::GitTreeEntry)
    oid_ptr = ccall((:git_tree_entry_id, :libgit2), Ptr{UInt8}, (Ptr{Void},), te.ptr)
    return GitHash(oid_ptr)
end

function Base.count(tree::GitTree)
    return ccall((:git_tree_entrycount, :libgit2), Csize_t, (Ptr{Void},), tree.ptr)
end

function Base.getindex(tree::GitTree, i::Integer)
    if i < 1 || i > count(tree)
        throw(BoundsError(tree, i))
    end
    te_ptr = ccall((:git_tree_entry_byindex, :libgit2),
                   Ptr{Void},
                   (Ptr{Void}, Csize_t), tree.ptr, i-1)
    return GitTreeEntry(tree, te_ptr, false)
end

"""
    (::Type{T})(te::GitTreeEntry) where T<:GitObject

Get the git object to which `te` refers and return it as its actual type (the type
[`entrytype`](@ref) would show), for instance a `GitBlob` or `GitTag`.

# Examples
```julia
tree = LibGit2.GitTree(repo, "HEAD^{tree}")
tree_entry = tree[1]
blob = LibGit2.GitBlob(tree_entry)
```
"""
function GitObject(e::GitTreeEntry) end
function (::Type{T})(te::GitTreeEntry) where T<:GitObject
    repo = repository(te)
    obj_ptr_ptr = Ref{Ptr{Void}}(C_NULL)
    @check ccall((:git_tree_entry_to_object, :libgit2), Cint,
                  (Ptr{Ptr{Void}}, Ptr{Void}, Ref{Void}),
                   obj_ptr_ptr, repo.ptr, te.ptr)
    return T(repo, obj_ptr_ptr[])
end

function Base.show(io::IO, te::GitTreeEntry)
    println(io, "GitTreeEntry:")
    println(io, "Entry name: ", filename(te))
    println(io, "Entry type: ", entrytype(te))
    println(io, "Entry OID: ", entryid(te))
end

function Base.show(io::IO, tree::GitTree)
    println(io, "GitTree:")
    println(io, "Owner: ", repository(tree))
    println(io, "Number of entries: ", count(tree))
end

max_stars_repo_name: xorJane/julia

max_issues_repo_name: xorJane/julia

max_forks_repo_name: xorJane/julia


text:
export writeplist_string, writeplist

writeplist_string(obj) = encode(obj)

function writeplist(stream::IO, obj)
    print(stream, encode(obj))
end

function writeplist(filename::AbstractString, obj)
    open(filename, "w") do stream
        writeplist(stream, obj)
    end
end

encode(key::AbstractString, obj) = string(key, " = ", encode(obj), ";")

function encode(dict::Dict{K, V}) where {K, V}
    string("{", join([encode(key, value) for (key, value) in dict]), "}")
end

function encode(array::Vector{T}) where T
    string( "(", join(map(encode, array), ", "), ")")
end

# Quotes strings which are not simple alphanumeric.
function encode(s::AbstractString)
    if all(isalnum, s) && length(s) > 0 && isletter(s[1])
        s
    else
        "\"$s\""
    end
end

function encode(binary::Vector{UInt8})
    string("<", bytes2hex(binary), ">")
end

encode(num::Real) = string(num)

function encode(obj)
    error("Can't encode $obj because it is of type $(typeof(obj))")
end


text:
using Test
using LinearAlgebra

using GAFramework
using GAFramework.CoordinateGA

function test1()
    model = CoordinateModel(x -> x[1]==0.0 ? 0.0 : x[1] * sin(1/x[1]), [-1.0], [1.0])
    state = GAState(model, ngen=50, npop=300, elite_fraction=0.01,
                    mutation_params=Dict(:rate=>0.1),
                    print_fitness_iter=10)
    best = ga(state)
    x = best.value[1]
    y = best.objvalue
    println("$best $x $y")
    abs(abs(x) - 0.223126) < 0.1 && abs(y - (-0.217219)) < 0.1
end

function test2()
    model = CoordinateModel(x -> any(z -> z==0.0, x) ? 0.0 : dot(x, sin.(1 ./x)),
                            -ones(Float64,15), ones(Float64,15))
    # do simulated annealing when mutating
    state = GAState(model, ngen=500, npop=300, elite_fraction=0.01,
                    mutation_params=Dict(:rate=>0.1,:sa_rate=>0.1,:k=>1,
                                         :lambda=>1/1000,:maxiter=>1000), print_fitness_iter=50)
    best = ga(state)
    x = best.value
    y = best.objvalue
    println("$best $x $y")
    all(abs.(abs.(x) .- 0.222549) .< 0.1) && abs(y - (-0.21723*15)) < 0.1
end

function test3()
    model = CoordinateModel(x -> 10length(x) + sum(z -> z^2 - 10cos(2pi*z), x),
                            -5.12 .* ones(15), 5.12 .* ones(15))
    state = GAState(model, ngen=500, npop=300, elite_fraction=0.01,
                    mutation_params=Dict(:rate=>0.1,:sa_rate=>0.1,:k=>1,
                                         :lambda=>1/1000,:maxiter=>1000), print_fitness_iter=50)
    best = ga(state)
    x = best.value
    y = best.objvalue
    println("$best $x $y")
    all(abs.(x) .< 0.1) && abs(y) < 0.1
end

@test test1()
@test test2()
@test test3()


text:
function get_right_endpoint(
    theta_init::Vector{Float64}, # initial point of parameters
    theta_num::Int, # number of parameter to scan
    loss_func::Function, # lambda(theta) - labmbda_min - delta_lambda
    method::Val{:QUADR_EXTRAPOL}; # function works only for method QUADR_EXTRAPOL

    theta_bounds::Vector{Tuple{Float64,Float64}} = fill(
        (-Inf, Inf), length(theta_init)
        ),
    scan_bound::Float64 = 9.0,
    scan_tol::Float64 = 1e-3,
    loss_tol::Float64 = 0., # 1e-3,
    # method args
    scan_hini = 1.,
    scan_hmax = Inf,
    # local alg args
    local_alg::Symbol = :LN_NELDERMEAD,
    max_iter::Int = 10^5,
    ftol_abs::Float64 = 1e-3,
    kwargs... # options for local fitter
    )
    # dim of the theta vector
    n_theta = length(theta_init)

    # checking arguments
    # methods which are not supported
    if local_alg in [:LN_BOBYQA, :LN_SBPLX, :LN_NEWUOA]
        @warn "Using local_alg = :"*String(local_alg)*" may result in wrong output."
    end
    # when using :LN_NELDERMEAD initial parameters should not be zero
    if local_alg == :LN_NELDERMEAD
        zeroParameter = [ isapprox(theta_init[i], 0., atol=1e-2) for i in 1:n_theta]
        if any(zeroParameter)
            @warn "Close-to-zero parameters found when using :LN_NELDERMEAD."
            show(findall(zeroParameter))
        end
    end

    # to count loss function calls inside this function, accumulation
    accum_counter::Int = 0
    # empty container
    pps = ProfilePoint[]

    prof = profile(
        theta_init,
        theta_num,
        loss_func;
        theta_bounds = theta_bounds,
        local_alg = local_alg,
        ftol_abs = loss_tol
    )

    # first iteration
    x_3 = theta_init[theta_num]
    theta_init_3 = theta_init
    iteration_count = 0

    # other iterations
    while true
        # preparation
        global x_2 # not initialized for the first iteration
        global x_1 # not initialized for the first iteration
        global point_2 # not initialized for the first iteration
        global point_1 # not initialized for the first iteration
        iteration_count += 1 # to understand if this is a first iteration

        # get profile point
        point_3 = prof(
            x_3;
            theta_init_i = theta_init_3, # hypothetically this makes optimization more effective
            maxeval = max_iter - accum_counter # how many calls left
            )
        push!(pps, point_3)
        accum_counter += point_3.counter # update counter

        if point_3.ret == :MAXEVAL_REACHED
            return (nothing, pps, :MAX_ITER_STOP) # break
        elseif point_3.ret == :FORCED_STOP
            return (nothing, pps, :LOSS_ERROR_STOP)
        elseif x_3 >= scan_bound && point_3.loss < 0.
            return (nothing, pps, :SCAN_BOUND_REACHED) # break
        # no checking for the first and second iterations
        # elseif iteration_count>2 && isapprox(x_3, x_2, atol = scan_tol)
        elseif iteration_count>1 && isapprox((x_3 - x_2) * point_3.loss / (point_3.loss - point_2.loss), 0., atol = scan_tol)
            return (x_3, pps, :BORDER_FOUND_BY_SCAN_TOL) # break
        elseif isapprox(point_3.loss, 0., atol = loss_tol)
            return (x_3, pps, :BORDER_FOUND_BY_LOSS_TOL) # break
        else
        end

        # next step
        if iteration_count==1
            x_4_extrapol = x_3 + scan_hini
            x_4 = minimum([x_3+scan_hmax, x_4_extrapol, theta_bounds[theta_num][2]])
            (point_2, x_2, x_3, theta_init_3) = (point_3, x_3, x_4, point_3.params)
        elseif iteration_count==2
            D = [
                x_3^2 x_3 1;
                x_2^2 x_2 1;
                2*x_2 1 0
                ]
            Da = [
                point_3.loss x_3 1;
                point_2.loss x_2 1;
                0 1 0
                ]
            Db = [
                x_3^2 point_3.loss 1;
                x_2^2 point_2.loss 1;
                2*x_2 0 0
                ]
            Dc = [
                x_3^2 x_3 point_3.loss;
                x_2^2 x_2 point_2.loss;
                2*x_2 1 0
                ]
            a = det(Da) / det(D)
            b = det(Db) / det(D)
            c = det(Dc) / det(D)
            if a <= 0.
                x_4_extrapol = x_3 + scan_hini
            else
                x_4_extrapol = (-b + sqrt(b^2-4*a*c))/2/a
            end
            x_4 = minimum([x_3+scan_hmax, x_4_extrapol, theta_bounds[theta_num][2]])
            (x_1, x_2, x_3, point_1, point_2, theta_init_3) = (x_2, x_3, x_4, point_2, point_3, point_3.params)
        else
            D = [
                x_3^2 x_3 1;
                x_2^2 x_2 1;
                x_1^2 x_1 1
                ]
            Da = [
                point_3.loss x_3 1;
                point_2.loss x_2 1;
                point_1.loss x_1 1
                ]
            Db = [
                x_3^2 point_3.loss 1;
                x_2^2 point_2.loss 1;
                x_1^2 point_1.loss 1
                ]
            Dc = [
                x_3^2 x_3 point_3.loss;
                x_2^2 x_2 point_2.loss;
                x_1^2 x_1 point_1.loss
                ]
            a = det(Da) / det(D)
            b = det(Db) / det(D)
            c = det(Dc) / det(D)
            if a <= 0. || b^2-4*a*c < 0.
                x_4_extrapol = x_3 + scan_hini
            else
                x_4_extrapol = (-b + sqrt(b^2-4*a*c))/2/a
            end
            x_4 = minimum([x_3+scan_hmax, x_4_extrapol, theta_bounds[theta_num][2]])
            (x_1, x_2, x_3, point_1, point_2, theta_init_3) = (x_2, x_3, x_4, point_2, point_3, point_3.params)
        end
    end
end


text:
module BitBlit2D

include("geometry/Geometry.jl")
include("graphics/Graphics.jl")

end


text:
export inverselines
"""
    inverselines(nbit::Int; n_reg::Int=nbit) -> ChainBlock

inverse first `n_reg` lines

TODO:
deprecate this function, it is not used.
"""
function inverselines(nbit::Int; n_reg::Int=nbit)
    c = chain(nbit)
    for i = 1:(n_reg ÷ 2)
        push!(c, swap(i,(n_reg-i+1)))
    end
    c
end

max_stars_repo_name: QuantumBFS/YaoExtensions

max_issues_repo_name: QuantumBFS/YaoExtensions

max_forks_repo_name: QuantumBFS/YaoExtensions


text:
abstract Driver{Name}

type Device{C <: Class, D <: Driver}
    io::C
    commands::Set{ASCIIString}

    function Device(brick, path)
        io = C(brick, path)
        commands = Set{ASCIIString}(io.commands())
        new(io, commands)
    end
end

class_path{C <: Class, D <: Driver}(::Type{Device{C, D}}) = class_path(C)
name{N}(::Type{Driver{N}}) = N
driver_name{C, D <: Driver}(::Type{Device{C, D}}) = replace(string(name(D)), '_', '-')

function find_devices{T <: Device}(::Type{T}, brick::Brick)
    path = joinpath(brick.root_path, class_path(T))
    devices = Vector{T}()
    driver = driver_name(T)
    for dir in readdir(path)
        device_driver = open(joinpath(path, dir, "driver_name")) do f
            chomp(strip(readline(f)))
        end
        if device_driver == driver
            push!(devices, T(brick, joinpath(path, dir)))
        end
    end
    devices
end

function find_device_at_address{T <: Device}(::Type{T}, brick::Brick, address::AbstractString)
    for dev in find_devices(T, brick)
        if dev.io.address() == address
            return dev
        end
    end
    nothing
end

# Motors
typealias TachoMotor{Driver} Device{TachoMotorClass, Driver}

typealias LargeMotor TachoMotor{Driver{:lego_ev3_l_motor}}
typealias MediumMotor TachoMotor{Driver{:lego_ev3_m_motor}}

# Sensors
typealias LegoSensor{Driver} Device{LegoSensorClass, Driver}

typealias UltrasoundSensor LegoSensor{Driver{:lego_ev3_us}}
typealias ColorSensor LegoSensor{Driver{:lego_ev3_color}}
typealias TouchSensor LegoSensor{Driver{:lego_ev3_touch}}
typealias GyroSensor LegoSensor{Driver{:lego_ev3_gyro}}


text:
baremodule UnixIORPiHeaders

import Base
using Base: @ccall, Cstring, Cint, Sys, Vector, signed, Dict,
            zero, hasmethod, fieldtypes,
            (=>), (&), (+), (>>), (>), (==)

const __spawn_action = Cvoid

ioctl(fd, cmd, arg) = @ccall ioctl(fd::Cint, cmd::Cint, arg::Ptr{Cint})::Cint

Base.include(UnixIORPiHeaders, "generated.jl")

# Need Ptr{Ptr{UInt8}} not Ptr{String} for NULL-terminated string vectors:
execv(path, args::Vector{Ptr{UInt8}}) =
    @ccall execv(path::Ptr{UInt8}, args::Ptr{Ptr{UInt8}})::Cint

posix_spawn(pid, path, file_actions, attrp, argv::Vector{Ptr{UInt8}},
                                            envp::Vector{Ptr{UInt8}}) =
    @ccall posix_spawn(pid::Ptr{pid_t},
                       path::Cstring,
                       file_actions::Ptr{posix_spawn_file_actions_t},
                       attrp::Ptr{posix_spawnattr_t},
                       argv::Ptr{Ptr{UInt8}},
                       envp::Ptr{Ptr{UInt8}})::Cint

const SIG_DFL = sig_t(0)

# Need multiple methods for `open` and `fcntl`:
open(pathname::AbstractString, flags) =
    @ccall open(pathname::Cstring, flags::Cint)::Cint

open(pathname::AbstractString, flags, mode) =
    @ccall open(pathname::Cstring, flags::Cint, mode::Cint)::Cint

fcntl(fd, cmd) = @ccall fcntl(fd::Cint, cmd::Cint)::Cint
fcntl(fd, cmd, arg) = @ccall fcntl(fd::Cint, cmd::Cint, arg::Cint)::Cint

# Need variants for _m struct
tcsetattr_m(fd, action, p) =
    @ccall tcsetattr(fd::Cint, action::Cint, p::Ptr{termios_m})::Cint

tcgetattr_m(fd, p) =
    @ccall tcgetattr(fd::Cint, p::Ptr{termios_m})::Cint

cfsetspeed_m(p, speed) =
    @ccall cfsetspeed(p::Ptr{termios_m}, speed::speed_t)::Cint

# Not yet in glibc.
const SYS_pidfd_open=434 # https://git.io/J4j1A
pifd_open(pid, flags) = @ccall syscall(SYS_pidfd_open::Cint,
                                       pid::pid_t, flags::Cint)::Cint

const P_PIDFD = 3
waitid(idtype, id, infop, options) =
    @ccall syscall(SYS_waitid::Cint,
                   idtype::Cint,
                   id::id_t,
                   infop::Ptr{siginfo_t},
                   options::Cint,
                   Base.C_NULL::Ptr{Cvoid})::Cint

# Function-like mactos not yet wrapped.
WIFSIGNALED(x) = (((x & 0x7f) + 1) >> 1) > 0
WTERMSIG(x) = (x & 0x7f)
WIFEXITED(x) = WTERMSIG(x) == 0
WEXITSTATUS(x) = signed(UInt8((x >> 8) & 0xff))
WIFSTOPPED(x) = (x & 0xff) == 0x7f
WSTOPSIG(x) = WEXITSTATUS(x)
WIFCONTINUED(x) = x == 0xffff
WCOREDUMP(x) = x & 0x80

end # module


text:
import FEM.stiffness
import FEM.GaussPoint2
import FEM.Point2

facts("FEM.AbstractMaterial") do

context("FEM.Material.LinearIsotropic") do

    mat = FEM.linearisotropicmod().LinearIsotropic(200e9, 0.3)
    gp = GaussPoint2(Point2(0.0, 0.0), 0.0)
    context("FEM.Material.LinearIsotropic.stiffness") do

        # Stiffness
        #########################################################################################
        stiff = stiffness(mat, gp)
        # Below calculated using Calfem with command hooke(4, 200e9, 0.3)
        stiff_calfem = [
        [269230769230.769    115384615384.615    115384615384.615                  0];
        [115384615384.615    269230769230.769    115384615384.615                  0];
        [115384615384.615    115384615384.615    269230769230.769                  0];
        [               0                   0                   0   76923076923.0769]]

        @fact norm(stiff - stiff_calfem) / norm(stiff) => roughly(0.0, atol=10.0^(-14))
    end

end # context

end # facts


text:
### Random Variable subtypes

## Abstract parameter

abstract type Parameter{S<:ValueSupport, F<:VariateForm} <: Variable{Random} end

DiscreteParameter{F<:VariateForm} = Parameter{Discrete, F}
ContinuousParameter{F<:VariateForm} = Parameter{Continuous, F}
UnivariateParameter{S<:ValueSupport} = Parameter{S, Univariate}
MultivariateParameter{S<:ValueSupport} = Parameter{S, Multivariate}
MatrixvariateParameter{S<:ValueSupport} = Parameter{S, Matrixvariate}

ParameterVector{P<:Parameter} = Vector{P}

setpdf!(parameter::Parameter{S, F}, state::ParameterState{S, F}) where {S<:ValueSupport, F<:VariateForm} =
  parameter.setpdf(state, parameter.states)

setprior!(parameter::Parameter{S, F}, state::ParameterState{S, F}) where {S<:ValueSupport, F<:VariateForm} =
  parameter.setprior(state, parameter.states)

value_support{S<:ValueSupport, F<:VariateForm}(::Type{Parameter{S, F}}) = S
variate_form{S<:ValueSupport, F<:VariateForm}(::Type{Parameter{S, F}}) = F

function check_support(parameter::Parameter{S, F}, state::ParameterState{S, F}) where {S<:ValueSupport, F<:VariateForm}
  if value_support(parameter) != value_support(state)
    warn("Value support of parameter ($(value_support(parameter))) and of ($(value_support(state))) not in agreement")
  end

  if variate_form(parameter) != variate_form(state)
    error("Variate form of parameter ($(variate_form(parameter))) and of ($(variate_form(state))) not in agreement")
  end
end

dotshape(variable::Parameter) = "circle"


text:
module MyFiles

using SpecialFunctions, Dierckx, SymEngine, DelimitedFiles, ProgressMeter, DSP, FFTW, Roots, Parameters, JSON, NumericalIntegration, Dates
import DifferentialEquations: ODEProblem, solve, Tsit5

export  gnlse, Fiber, createTimeGrid, Pulse, gaussianPulse, sechPulse,
        pulseFromFile, FWHM, Pulse, getElectricField, spectrogram, wbeta2,
        lbeta2, lbeta3, sorder, dlength, L_operator, trapz, beta0, beta1,
        beta2, lfiss, χ3, bound, adk, geisslerPlasma, plasmaTerms, gnlsetest,
        wbeta0, wbeta1, uppe, getElectricFieldE, mfft, mifft, lbeta0,
        uppe2, LE_operator, βnl, dis_wave, zero_disp, lbeta2,
        soliton_order, getGamma, Waist, systemInfo,  soliton_period, smoothstep,
        flipsmoothstep, diffsion, test, diffsion0, rhs2, testrhs, ionRateAdk,
        geisslerPlasmaUPPE, cumtrapz, arb_f, chirpSechPulse, System,
        c, ϵ0, rindex, χ30,refractive_index, O3_disp, dispersion_coeffecient,
        gas, Tk_0, grid, Realχ, splbeta, Realχ, testdis, beta

function __init__()
    global x = symbols(:x)
end

include("constants.jl")
include("grid.jl")
include("raman.jl")
include("functions.jl")

include("gas.jl")
include("plasma.jl")
include("pulse.jl")
include("fiber.jl")
include("dispersion.jl")
include("gnlse.jl")
include("system.jl")

include("uppe.jl")
include("test.jl")

end # module

max_stars_repo_name: SabbahMohammed/Dev

max_issues_repo_name: SabbahMohammed/Dev

max_forks_repo_name: SabbahMohammed/Dev


text:
using Revise
using RemoteGtkIDE
using Base.Test

#value  = String(read(`julia repl.jl`))

port, server = RemoteGtkIDE.start_server()

client = connect(port)
serialize(client,("test",))
deserialize(client)

# @async while true
#     response = deserialize(client)
#     println(response)
# end

@test length( RemoteGtkIDE.remote_eval(client, Main, :(x=rand(3))) ) == 3

@test RemoteGtkIDE.remotecall_fetch(sin, client, 0) ≈ 0.0

@test_throws MethodError throw(RemoteGtkIDE.remotecall_fetch(sin, client, "pi"))


text:
export Parameters
using YAML

struct Parameters{T<:AbstractFloat,S<:Flag,U<:Bool,V<:Int,W<:String}
    # T float, S flag, U bool, V int, W string
    M::T
    mdot::T
    spin::T
    mu_nucleon::T
    reprocessing::U
    hard_xray_fraction::T
    corona_electron_energy::T
    warm_electron_energy::T
    warm_photon_index::T
    reflection_albedo::T
    min_energy::T
    max_energy::T
    uv_min_energy::T
    uv_max_energy::T
    xray_min_energy::T
    xray_max_energy::T
    function Parameters(;
        M=1e8,
        mdot=0.5,
        spin=0.0,
        mu_nucleon=0.61,
        reprocessing=true,
        hard_xray_fraction=0.02,
        corona_electron_energy=100.0,
        warm_electron_energy=0.2,
        warm_photon_index=2.5,
        reflection_albedo=0.3,
        min_energy=1e-4, # kev
        max_energy=200.0,
        uv_min_energy=0.00387,
        uv_max_energy=0.06,
        xray_min_energy=0.06,
        xray_max_energy=1.0,
    )
        return new{Float64, Flag, Bool, Int, String}(
            M,
            mdot,
            spin,
            mu_nucleon,
            reprocessing,
            hard_xray_fraction,
            corona_electron_energy,
            warm_electron_energy,
            warm_photon_index,
            reflection_albedo,
            min_energy,
            max_energy,
            uv_min_energy,
            uv_max_energy,
            xray_min_energy,
            xray_max_energy,
        )
    end
end

function Parameters(config::Dict)
    return Parameters(;config...)
end
function Parameters(file_path::String)
    config = YAML.load_file(file_path, dicttype=Dict{Symbol, Any})
    return Parameters(config)
end


text:
using StackViews, OffsetArrays
using Test
using Aqua, Documenter

@testset "Project meta quality checks" begin
    # Not checking compat section for test-only dependencies
    Aqua.test_all(StackViews; project_extras=true, deps_compat=true, stale_deps=true, project_toml_formatting=true)
    if VERSION >= v"1.2"
        doctest(StackViews, manual = false)
    end
end

@testset "StackViews.jl" begin
    A = reshape(collect(1:12), 3, 4)
    B = reshape(collect(13:24), 3, 4)

    sv = @inferred StackView([A, B])
    @test size(sv) == (3, 4, 2)
    @test collect(sv) == cat(A, B; dims=3)

    sv = if VERSION >= v"1.2.0"
        @inferred StackView((A, B))
    else
        StackView((A, B))
    end
    @test size(sv) == (3, 4, 2)
    @test collect(sv) == cat(A, B; dims=3)

    sv = @inferred StackView([A, B], Val(3))
    @test size(sv) == (3, 4, 2)
    @test collect(sv) == cat(A, B; dims=3)

    sv = @inferred StackView([A, B], Val(1))
    @test size(sv) == (2, 3, 4)
    @test sv[1, :, :] == A
    @test sv[2, :, :] == B

    sv = @inferred StackView([A, B], Val(2))
    @test size(sv) == (3, 2, 4)
    @test sv[:, 1, :] == A
    @test sv[:, 2, :] == B

    sv = @inferred StackView([A, B], Val(3))
    @test size(sv) == (3, 4, 2)
    @test sv[:, :, 1] == A
    @test sv[:, :, 2] == B

    sv = @inferred StackView([A, B], Val(4))
    @test size(sv) == (3, 4, 1, 2)
    @test sv[:, :, 1, 1] == A
    @test sv[:, :, 1, 2] == B

    sv = @inferred StackView([A, B], Val(5))
    @test size(sv) == (3, 4, 1, 1, 2)
    @test sv[:, :, 1, 1, 1] == A
    @test sv[:, :, 1, 1, 2] == B

    @test StackView([A, B], 1) == StackView([A, B], Val(1))
    @test StackView(A) == StackView([A, ]) == StackView((A, ))

    @test_throws ArgumentError StackView([A, B[:]])
    @test_throws ArgumentError StackView([A, B], 0)
    @test_throws ArgumentError StackView([A, B], Val(0))

    # repeat the array multiple times is allowed
    A = [1, 2, 3, 4]
    @test StackView([A for i in 1:3]) == repeat(A, 1, 3)

    @testset "heterogeneous arrays" begin
        A = Int[1, 2, 3, 4]
        B = Union{Missing, Float64}[missing, 6, 7, 8]
        sv = @inferred StackView([A, B])
        @test eltype(sv) == Union{Missing, Float64}
        T = Union{Missing, Float32}
        sv = @inferred StackView{T}([A, B])
        @test eltype(sv) == T

        A = Int[1, 2, 3, 4]
        B = [[1, ], [2, ], [3, ], [4, ]]
        @test_throws ArgumentError StackView(A, B)
        # Perhaps suprisingly, with explicit type T it bypasses the type check
        # It's an undefined behavior, just put it here for the record
        sv = @test_nowarn StackView{Int}(A, B)
        @test_throws MethodError collect(sv)
    end

    @testset "axes" begin
        # axes are unified to 1-based
        A = OffsetArray(collect(reshape(1:8, 2, 4)), -1, 1)
        B = OffsetArray(collect(reshape(9:16, 2, 4)), 1, -1)
        sv = StackView([A, B], 1)
        @inferred axes(sv)
        @test axes(sv) === (Base.OneTo(2), Base.OneTo(2), Base.OneTo(4))

        A = collect(reshape(1:8, 2, 4))
        B = collect(reshape(9:16, 2, 4))
        sv = StackView([A, B], 1)
        @inferred axes(sv)
        @test axes(sv) === (Base.OneTo(2), Base.OneTo(2), Base.OneTo(4))
    end

    @testset "setindex!" begin
        # setindex! modifies original arrays
        for dim in 1:3
            A = [0, 1, 2, 3]
            B = [4, 5, 6, 7]
            sv = StackView([A, B], dim)
            fill!(sv, -1)
            @test A == B == fill(-1, 4)
        end

        A = [0, 1, 2, 3]
        B = [4, 5, 6, 7]
        sv = StackView(A, B)
        sv[4, 2] = -1
        @test sv[4, 2] == B[4] == -1
    end
end


text:
module Photometry

using Reexport

include("aperture/Aperture.jl")
include("background/Background.jl")
include("detection/Detection.jl")

@reexport using .Aperture
@reexport using .Background
@reexport using .Detection

end


text:
# ---
# title: 1203. Sort Items by Groups Respecting Dependencies
# id: problem1203
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Depth-first Search, Graph, Topological Sort
# link: <https://leetcode.com/problems/sort-items-by-groups-respecting-dependencies/description/>
# hidden: true
# ---
#
# There are `n` items each belonging to zero or one of `m` groups where
# `group[i]` is the group that the `i`-th item belongs to and it's equal to `-1`
# if the `i`-th item belongs to no group. The items and the groups are zero
# indexed. A group can have no item belonging to it.
#
# Return a sorted list of the items such that:
#
#   * The items that belong to the same group are next to each other in the sorted list.
#   * There are some relations between these items where `beforeItems[i]` is a list containing all the items that should come before the `i`-th item in the sorted array (to the left of the `i`-th item).
#
# Return any solution if there is more than one solution and return an **empty
# list**  if there is no solution.
#
#
#
# **Example 1:**
#
# **![](https://assets.leetcode.com/uploads/2019/09/11/1359_ex1.png)**
#
#
#
#     Input: n = 8, m = 2, group = [-1,-1,1,0,0,1,0,-1], beforeItems = [[],[6],[5],[6],[3,6],[],[],[]]
#     Output: [6,3,4,1,5,2,0,7]
#
#
# **Example 2:**
#
#
#
#     Input: n = 8, m = 2, group = [-1,-1,1,0,0,1,0,-1], beforeItems = [[],[6],[5],[6],[3],[],[4],[]]
#     Output: []
#     Explanation:  This is the same as example 1 except that 4 needs to be before 6 in the sorted list.
#
#
#
#
# **Constraints:**
#
#   * `1 <= m <= n <= 3*10^4`
#   * `group.length == beforeItems.length == n`
#   * `-1 <= group[i] <= m-1`
#   * `0 <= beforeItems[i].length <= n-1`
#   * `0 <= beforeItems[i][j] <= n-1`
#   * `i != beforeItems[i][j]`
#   * `beforeItems[i] `does not contain duplicates elements.
#
#
## @lc code=start
using LeetCode

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


text:
using LaTeXTabulars

# for testing
using LaTeXTabulars: latex_cell

using Test
using LaTeXStrings

"Normalize whitespace, for more convenient testing."
squash_whitespace(string) = strip(replace(string, r"[ \n\t]+" => " "))

@test squash_whitespace(" something  \n with line  breaks \n  and   stuff \n") ==
    "something with line breaks and stuff"

"Comparison using normalized whitespace. For testing."
≅(a, b) = squash_whitespace(a) == squash_whitespace(b)

@testset "tabular" begin
    tb = Tabular("lcl")
    tlines = [Rule(:top),
              [L"\alpha", L"\beta", "sum"],
              Rule(:mid),
              [1, 2, 3],
              Rule(),           # a nice \hline to make it ugly
              [4.0 "5" "six";   # a matrix
               7 8 9],
              (CMidRule(1, 2), CMidRule("lr", 1, 1)), # just to test tuples
              [MultiColumn(2, :c, "centered")],       # ragged!
              Rule(:bottom)]
    tlatex = raw"\begin{tabular}{lcl}
                 \toprule
                 $\alpha$ & $\beta$ & sum \\
                 \midrule
                 1 & 2 & 3 \\
                 \hline
                 4.0 & 5 & six \\
                 7 & 8 & 9 \\ \cmidrule{1-2} \cmidrule(lr){1-1}
                 \multicolumn{2}{c}{centered} \\
                 \bottomrule
                 \end{tabular}"

    tlatex = replace(tlatex, "\r\n"=>"\n")
    @test latex_tabular(String, tb, tlines) ≅ tlatex
    tmp = tempname()
    latex_tabular(tmp, tb, tlines)
    @test isfile(tmp) && read(tmp, String) ≅ tlatex
    @test read(tmp, String) ≅ tlatex
end

@test_throws ArgumentError latex_cell(stdout, MultiColumn(2, :BAD, ""))
@test_throws ArgumentError CMidRule(3, 1)     # not ≤
@test_throws MethodError latex_cell(stdout, ("un", "supported"))
@test_throws MethodError CMidRule(1, 1, 1, 2) # invalid types


text:
using SumOfSquaresOptimization
using Base.Test

objtol = 1e-4

function equals_upto_perm(a,b)
    if length(a) != length(b)
        return false
    end
    mapreduce(p -> b[p] == a, Base.OrFun(), permutations(1:length(b)))
end
@test equals_upto_perm([1,2],[2,1])

# UNIT TESTS

# monomials & polynomials
a = Dict( :x1 => 3, :x2 => 2 )
b = Dict( :x1 => 1, :x3 => 5 )
c = Dict( :x1 => 4, :x2 => 2, :x3 => 5 )
@test a * b == c
@test SumOfSquaresOptimization.deg(a) == 5
p = Dict(a=>1.0, b=>2.0)
@test Dict(a=>1.0) + Dict(b=>2.0) == p
@test SumOfSquaresOptimization.deg(p) == 6
@test equals_upto_perm(SumOfSquaresOptimization.monoms([:x1,:x2],2), [ Dict(:x1=>2), Dict(:x1=>1, :x2=>1), Dict(:x2=>2) ])
@test SumOfSquaresOptimization.monoms(0) == [SumOfSquaresOptimization.SoSMonom()]
@test SumOfSquaresOptimization.monoms(3) == Array{SumOfSquaresOptimization.SoSMonom,1}()
@test equals_upto_perm(SumOfSquaresOptimization.decomp(Dict(:x1=>4),3), [ (Dict(:x1=>3),Dict(:x1=>1)), (Dict(:x1=>2),Dict(:x1=>2)), (Dict(:x1=>1),Dict(:x1=>3)) ])

# parsing
@test SumOfSquaresOptimization.parsepoly(Program(), parse("x1^3 * x2^2")) == Dict(a => 1.0)
@test SumOfSquaresOptimization.parsepoly(Program(), parse("x1^3 * x2^2 + 2.0 * x1 * x3^5")) == p

# END-TO-END TESTS

# unconstrained minimization: a quadratic bowl
prog = Program(minimize=true)
objective(prog, "x^2 + y^2 - 1")
sol = sossolve(prog,2)
@test status(sol) == :Normal
@test abs(primalobj(sol) + 1) < objtol
@test_throws ArgumentError moment(sol,"x^4")

# unconstrained, unbounded maximization
prog = Program(maximize=true)
objective(prog, "x^2 + y^2 - 1")
sol = sossolve(prog,2)
@test status(sol) == :Unbounded
@test isinf(primalobj(sol))
@test !signbit(primalobj(sol))
@test_throws ArgumentError moment(sol,"x")
@test_throws ArgumentError dumpsol(sol)

# min bisection in the 6-cycle
prog = Program(minimize=true)
for i in 1:6
    for j in [-1,1]
        jp = (i + j + 2*6 - 1) % 6 + 1
        @partobjective(prog, "x%d * (1-x%d)", i, jp)
    end
    @constraint(prog, "x%d^2 == x%d", i, i)
    @partconstraint(prog, :bisect, "x%d", i)
end
@partconstraint(prog, :bisect, "-%d", 6/2)

sol = sossolve(prog,2)
@test status(sol) == :Normal
@test abs(primalobj(sol) - 1.5) < objtol

# test dihedral symmetry on the same program
symmetrize_dihedral!(prog, [symbol("x$i") for i in 1:6])
sol = sossolve(prog,2)
@test status(sol) == :Normal
@test abs(primalobj(sol) - 1.5) < objtol
sol = sossolve(prog,4)
@test status(sol) == :Normal
@test abs(primalobj(sol) - 2.0) < objtol

max_stars_repo_name: willperry/SumOfSquaresOptimization


text:
include("help.jl")
function thecallofthemountain(r::Robot)
    steps_up = moves!(r,Nord) #вверх
    steps_left = moves!(r,West) #влево
    steps_down = moves!(r,Sud) #вниз

    #Робот - в юго-западном углу и в num_steps - закодирован пройденный путь
    side = Nord
    snake!(r,side) #Двигаемся змейкой

   moves!(r,West)
   moves!(r,Sud)

   movements!(r,Nord, steps_down)
   movements!(r,Ost,steps_left)
   movements!(r,Sud,steps_up)
   # Робот - в исходном положении
end

function snake!(r::Robot, side::HorizonSide)
    fl = false
    while (isborder(r,side) == false) #Если сверху или снизу нет стенки идем туда
        if (isborder(r,Ost)==true) #Если справа стенка
            fl = true
            break
        else #Если стенки нет
            move!(r,side)
        end
    end
    if (fl == false)
        move!(r,Ost)
        side = invers(side)
        snake!(r,side) #Рекурсия
    else

        side = Ost
        for way in (Nord, Ost, Sud, West)
            while isborder(r, side) == true && ismarker(r) == false
                putmarker!(r)
                move!(r,way)
            end
            putmarker!(r)
            move!(r,side)
            side = HorizonSide(mod(Int(side) - 1,4))
        end
        putmarker!(r)
    end
end

max_stars_repo_name: MCPavlo/mirea-progs

max_issues_repo_name: MCPavlo/mirea-progs

max_forks_repo_name: MCPavlo/mirea-progs


text:
using Visualize, GeometryTypes, ModernGL
using Visualize.JLRasterization
using Visualize: Debugging
using Visualize.JLRasterization: JLWindow
using ModernGL, FileIO
using Visualize: GLRasterization, Resolution, View, ProjectionView, EyePosition, Projection
using Visualize: VertexN, Light, MeshUniforms, Shading, vert_mesh, frag_mesh, GLCanvas
using Visualize: Model, Color, ShadingFunction
using Colors

window = JLWindow((Area => (800, 600),))

canvas = window[Canvas];
for field in FieldTraits.Fields(eltype(canvas))
    if haskey(cam, field)
        FieldTraits.link!(field, cam => canvas)
    end
end

light = Light(
    Vec3f0(10),
    Vec3f0(0.0),
    Vec3f0(0.9),
    1.0f0,
    Vec3f0(0.2),
    1.0f0
)
shading = Shading(
    Vec3f0(0.0),
    Vec3f0(0.3),
    8.0f0
)
# Create a functor
catmesh = normalmesh(load(Pkg.dir("GLVisualize", "assets", "cat.obj")))
catmesh = Base.view(
    reinterpret(Visualize.VertexN, catmesh.parent), catmesh.indexes[1]
)
show!(window)

vbo = VertexArray(catmesh);

args = MeshUniforms((
    ShadingFunction => ((V, N, L, color, shading, light) -> return L)
,))

uniforms = (canvas, map(UniformBuffer, (light, shading, args))...)

draw_cat = GLRasterizer(
    vbo,
    uniforms,
    vert_mesh, Visualize.frag_mesh
)

glEnable(GL_DEPTH_TEST)
glDepthMask(GL_TRUE)
glDepthFunc(GL_LEQUAL)
glDisable(GL_CULL_FACE)
GLAbstraction.enabletransparency()
glClearColor(1, 1, 1, 1)

@async begin
    while isopen(window)
        GLFW.PollEvents()
        glViewport(0, 0, widths(window[Area])...)
        glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
        canvas = GLCanvas(
            (
                (elem => cam[elem] for elem in (Projection, View, ProjectionView, EyePosition))...,
                Resolution => (800, 600),
            )
        )
        uniforms[1][1] = canvas # update
        draw_cat(vbo, uniforms)
        swapbuffers!(window)
        sleep(0.001)
        yield()
    end
    destroy!(window)
end


text:
module MeshCat

using GeometryBasics: GeometryBasics,
  AbstractGeometry,
  AbstractMesh,
  AbstractNgonFace,
  AbstractFace,
  Cylinder,
  GLTriangleFace,
  GeometryPrimitive,
  HyperRectangle,
  HyperSphere,
  NgonFace,
  OffsetInteger,
  Point,
  Point3f0,
  Polytope,
  SimplexFace,
  Vec,
  decompose,
  origin,
  radius,
  texturecoordinates,
  raw,
  widths

using CoordinateTransformations
using Rotations: rotation_between, Rotation, RotMatrix, UnitQuaternion
using Colors: Color, Colorant, RGB, RGBA, alpha, hex, red, green, blue
using StaticArrays: StaticVector, SVector, SDiagonal, SMatrix
using Parameters: @with_kw
using DocStringExtensions: SIGNATURES, TYPEDSIGNATURES
using Requires: @require
using Base.Filesystem: rm
using BinDeps: download_cmd, unpack_cmd
using UUIDs: UUID, uuid1
using LinearAlgebra: UniformScaling, Diagonal, norm
using Sockets: listen, @ip_str, IPAddr, IPv4, IPv6
using Base64: base64encode
using MsgPack: MsgPack, pack
import Mux
import Logging
import Mux.WebSockets
import Cassette
import FFMPEG

import Base: delete!

export Visualizer,
       IJuliaCell,
       render,
       render_static,
       static_html

export setobject!,
       settransform!,
       delete!,
       setprop!,
       setanimation!,
       setvisible!

export AbstractVisualizer,
       AbstractMaterial,
       AbstractObject,
       GeometryLike

export Object,
       HyperEllipsoid,
       PointCloud,
       Cone,
       Triad,
       MeshObject,
       MeshFileGeometry,
       MeshFileObject,
       Points,
       Line,
       LineLoop,
       LineSegments

export PointsMaterial,
       MeshLambertMaterial,
       MeshBasicMaterial,
       MeshPhongMaterial,
       LineBasicMaterial,
       Texture,
       PngImage

export Animation,
       atframe,
       @animation

export ArrowVisualizer

abstract type AbstractObject end
abstract type AbstractMaterial end

include("util.jl")
include("trees.jl")
using .SceneTrees
include("mesh_files.jl")
include("geometry.jl")
include("objects.jl")
include("animations.jl")
include("commands.jl")
include("abstract_visualizer.jl")
include("lowering.jl")
include("msgpack.jl")
include("visualizer.jl")
include("atframe.jl")
include("arrow_visualizer.jl")
include("render.jl")
include("servers.jl")
include("assets.jl")
include("integrations.jl")

const VIEWER_ROOT = joinpath(@__DIR__, "..", "assets", "meshcat", "dist")

function __init__()
    main_js = abspath(joinpath(VIEWER_ROOT, "main.min.js"))
    if !isfile(main_js)
        error("""
        main.min.js not found at $main_js.
        Please build MeshCat using `import Pkg; Pkg.build("MeshCat")`""")
    end
    setup_integrations()
end

end


text:
abstract type DiseaseStateSequence end
abstract type SEIR <: DiseaseStateSequence end
abstract type SIR <: DiseaseStateSequence end
abstract type SEI <: DiseaseStateSequence end
abstract type SI <: DiseaseStateSequence end


text:
function estimate(setup::Setup)
    y, m, P, S, β, δ = setup.y, setup.m, setup.P, setup.S, setup.β, setup.δ

    T, p = size(y)
    d = get_d(setup)
    a = d*p + 1
    F = [1; zeros(a - 1)]
    k = get_k(p, β)
    ν = get_ν(β)

    out = Estimation(y,               # y
                     zeros(T, a, p),  # m
                     zeros(T, a, a),  # P
                     zeros(T, p, p),  # S
                     zeros(T, p),     # μ
                     zeros(T, p, p),  # Σ
                     zeros(T, p),     # e
                     zeros(T, p),     # u
                     ν,
                     β,
                     δ)

    for t = 1:T
        if d > 1 && t > d
            F[2:end] = vec(y[t-1:-1:t-d, :])
        end

        # Predict parameters at time t|t-1
        P, S, Q = predict_parameters(P, S, F, δ, k)

        # Predict at time t|t-1
        μ, Σ = predict(m, Q, S, F, β)

        # Prediction error
        e = error(y[t, :], μ)
        u = standardised_error(e, Σ)

        # Storage
        out.m[t, :, :] = m
        out.P[t, :, :] = P
        out.S[t, :, :] = S
        out.μ[t, :]    = μ
        out.Σ[t, :, :] = Σ
        out.e[t, :]    = e
        out.u[t, :]    = u

        # Update at time t|t
        if t > d
            m, P, S = update_parameters(m, P, Q, S, F, e)
        end

    end
    return out
end

function predict_parameters(P, S, F, δ, k)
    # note: m is not updated
    P = P / δ
    Q = F' * P * F + 1.0
    S = S / k
    return (P, S, Q)
end

function predict(m, Q, S, F, β)
    μ = m' * F
    Σ = Q * (1 - β) / (3β - 2) * S
    return (μ, Σ)
end

function error(y, μ)
    return y - μ
end

function standardised_error(e, Σ)
    return inv(cholesky(Σ).L) * e
end

function update_parameters(m, P, Q, S, F, e)
    K = P * F / Q
    m = m + K * e'
    P = P - K * K' * Q
    S = S + e*e'/Q
    return (m, P, S)
end