text:
export ambient_space, rank, gram_matrix, inner_product, involution,
       islocal_square, isequivalent, isrationally_equivalent, quadratic_space,
       hermitian_space, diagonal, invariants, hasse_invariant, witt_invariant, orthogonal_basis, fixed_field

################################################################################
#
#  Types and constructors
#
################################################################################

abstract type AbsSpace{S} end

abstract type AbsLat{S} end

mutable struct QuadSpace{S, T} <: AbsSpace{S}
  K::S
  gram::T
  @declare_other

  function QuadSpace(K::S, G::T) where {S, T}
    # I also need to check if the gram matrix is Hermitian
    if dense_matrix_type(elem_type(S)) === T
      z = new{S, T}(K, G)
    else
      try
        Gc = change_base_ring(K, G)
        if typeof(Gc) !== dense_matrix_type(elem_type(S))
          throw(error("Cannot convert entries of the matrix to the number field"))
        end
        z = new{S, dense_matrix_type(elem_type(S))}(K, Gc)
        return z
      catch e
        rethrow(e)
        throw(error("Cannot convert entries of the matrix to the number field"))
      end
    end
  end
end

@doc Markdown.doc"""
    quadratic_space(K::NumField, n::Int) -> QuadSpace

Create the quadratic space over `K` with dimension `n` and Gram matrix
equal to the identity matrix.
"""
function quadratic_space(K::NumField, n::Int)
  G = identity_matrix(K, n)
  return QuadSpace(K, G)
end

@doc Markdown.doc"""
    quadratic_space(K::NumField, G::Int) -> QuadSpace

Create the quadratic space over `K` with Gram matrix `G`.
The matrix `G` must be square and symmetric.
"""
function quadratic_space(K::NumField, G::MatElem)
  return QuadSpace(K, G)
end

function quadratic_space(K::Field, G::MatElem)
  return QuadSpace(K, G)
end

mutable struct HermSpace{S, T, U, W} <: AbsSpace{S}
  E::S
  K::T
  gram::U
  involution::W
  @declare_other

  function HermSpace(E::S, gram::U) where {S, U}
    # I also need to check if the gram matrix is Hermitian
    if dense_matrix_type(elem_type(S)) === U
      gramc = gram
    else
      try
        gramc = change_base_ring(E, gram)
        if typeof(gramc) !== dense_matrix_type(elem_type(S))
          throw(error("Cannot convert entries of the matrix to the number field"))
        end
      catch e
        throw(error("Cannot convert entries of the matrix to the number field"))
      end
    end

    @assert degree(E) == 2
    A = automorphisms(E)
    a = gen(E)
    if A[1](a) == a
      involution = A[2]
    else
      involution = A[1]
    end

    K = base_field(E)

    z = new{S, typeof(K), dense_matrix_type(elem_type(S)), typeof(involution)}(E, K, gramc, involution)
    return z
  end
end

@doc Markdown.doc"""
    hermitian_space(K::NumField, n::Int) -> HermSpace

Create the Hermitian space over `K` with dimension `n` and Gram matrix equal to
the identity matrix. The number field `K` must be a quadratic extension, that
is, `degree(K) == 2` must hold.
"""
function hermitian_space(K::NumField, n::Int)
  G = identity_matrix(K, n)
  return HermSpace(K, G)
end

@doc Markdown.doc"""
    hermitian_space(K::NumField, G::MatElem) -> HermSpace

Create the Hermitian space over `K` with Gram matrix equal to the identity
matrix. The matrix `G` must be square and Hermitian with respect to the non-trivial
automorphism of `K`. The number field `K` must be a quadratic extension, that
is, `degree(K) == 2` must hold.
"""
function hermitian_space(K::NumField, G::MatElem)
  return HermSpace(K, G)
end

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

function Base.show(io::IO, V::QuadSpace)
  print(io, "Quadratic space over\n")
  println(io, base_ring(V))
  println(io, "with Gram matrix")
  print(io, gram_matrix(V))
end

function Base.show(io::IO, V::HermSpace)
  print(io, "Hermitian space over\n")
  println(io, base_ring(V))
  println(io, "with Gram matrix")
  print(io, gram_matrix(V))
end

################################################################################
#
#  Basic invariants
#
################################################################################

@doc Markdown.doc"""
    isquadratic(V::AbsSpace) -> Bool

Returns whether $V$ is quadratic.
"""
isquadratic(V::AbsSpace)

isquadratic(V::QuadSpace) = true

isquadratic(V::HermSpace) = false

@doc Markdown.doc"""
    ishermitian(V::AbsSpace) -> Bool

Returns whether $V$ is hermitian.
"""
ishermitian(V::AbsSpace)

ishermitian(V::QuadSpace) = true

ishermitian(V::HermSpace) = false

@doc Markdown.doc"""
    rank(V::AbsSpace) -> Int

Return the rank of the space `V`.
"""
rank(L::AbsSpace) = nrows(L.gram)

@doc Markdown.doc"""
    dim(V::AbsSpace) -> Int

Return the dimension of the space `V`.
"""
dim(V::AbsSpace) = rank(V)

@doc Markdown.doc"""
    gram_matrix(V::AbsSpace) -> MatElem

Return the Gram matrix of `V`.
"""
gram_matrix(V::AbsSpace) = V.gram

# Once we have quaternion spaces the following makes more sense
_base_algebra(V::QuadSpace) = V.K

_base_algebra(V::HermSpace) = V.E

@doc Markdown.doc"""
    base_ring(V::AbsSpace) -> NumField

Return the base field of `V`.
"""
base_ring(V::AbsSpace) = _base_algebra(V)

involution(V::QuadSpace) = identity

involution(V::HermSpace) = V.involution

@doc Markdown.doc"""
    fixed_field(V::AbsSpace) -> NumField

Return the fixed field of `V`.
"""
fixed_field(::AbsSpace)

fixed_field(V::QuadSpace) = base_ring(V)

fixed_field(V::HermSpace) = V.K

@doc Markdown.doc"""
    involution(V::AbsSpace) -> NumField

Return the involution of `V`.
"""
involution(V::AbsSpace)

# Maybe cache this

@doc Markdown.doc"""
    isregular(V::AbsSpace) -> Bool

Return whether `V` is regular, that is, if the Gram matrix
has full rank.
"""
function isregular(V::AbsSpace)
  G = gram_matrix(V)
  return rank(G) == nrows(G)
end

################################################################################
#
#  Determinant and discriminant
#
################################################################################

function det(V::QuadSpace)
  return det(gram_matrix(V))
end

function det(V::HermSpace)
  d = det(gram_matrix(V))
  @assert all(iszero(coeff(d, i)) for i in 1:degree(base_ring(V)) - 1)
  return coeff(d, 0)
end

@doc Markdown.doc"""
    det(V::AbsSpace) -> FieldElem

Returns the determinant of the space `V`.

In case `V` is Hermitian, the result is an element of the fixed field.
"""
det(::AbsSpace)

@doc Markdown.doc"""
    discriminant(V::AbsSpace) -> FieldElem

Returns the discriminant of the space `V`.

In case `V` is Hermitian, the result is an element of the "smaller field".
"""
function discriminant(V::AbsSpace)
  d = det(V)
  n = mod(rank(V), 4)
  if n == 0 || n == 1
    return d
  else
    return -d
  end
end

function _discriminant(G)
  d = det(G)
  if n == 0 || n == 1
    return d
  else
    return -d
  end
end

################################################################################
#
#  Computing inner products
#
################################################################################

@doc Markdown.doc"""
    gram_matrix(V::AbsSpace, M::MatElem) -> MatElem

Returns the gram matrix of the rows of `M`.
"""
function gram_matrix(V::AbsSpace{T}, M::MatElem{S}) where {S, T}
  if S === elem_type(T)
    return M * gram_matrix(V) * transpose(_map(M, involution(V)))
  else
    Mc = change_base_ring(base_ring(V), M)
    return Mc * gram_matrix(V) * transpose(_map(Mc, involution(V)))
  end
end

@doc Markdown.doc"""
    gram_matrix(V::AbsSpace, S::Vector{Vector}) -> MatElem

Returns the gram matrix of the sequence `S`.
"""
function gram_matrix(V::AbsSpace{T}, S::Vector{Vector{U}}) where {T, U}
  m = zero_matrix(base_ring(V), length(S), rank(V))
  for i in 1:length(S)
    if length(S[i]) != rank(V)
      throw(error("Vectors must be of length $(rank(V))"))
    end
    for j in 1:rank(V)
      m[i, j] = S[i][j]
    end
  end
  return gram_matrix(V, m)
end

function _inner_product(V, v, w)
  mv = matrix(base_ring(V), 1, nrows(V), v)
  mw = matrix(base_ring(V), ncols(V), 1, w)
  return (mv * V * mw)[1, 1]
end

function _inner_product(G, v, w, involution)
  return _inner_product(G, v, [involution(x) for x in w])
end

inner_product(V::QuadSpace, v::Vector, w::Vector) = _inner_product(gram_matrix(V), v, w)

inner_product(V::HermSpace, v::Vector, w::Vector) = _inner_product(gram_matrix(V), v, w, involution(V))

@doc Markdown.doc"""
    inner_product(V::AbsSpace, v::Vector, w::Vector) -> FieldElem

Returns the inner product of `v` and `w` with respect to `V`.
"""
inner_product(V::AbsSpace, v::Vector, w::Vector)

################################################################################
#
#  Diagonalization
#
################################################################################

function diagonal(V::QuadSpace)
  D, _ = _gram_schmidt(gram_matrix(V), involution(V))
  return diagonal(D)
end

function diagonal(V::HermSpace)
  D, _ = _gram_schmidt(gram_matrix(V), involution(V))
  return [ coeff(d, 0) for d in diagonal(D) ]
end

@doc Markdown.doc"""
    orthogonal_basis(V::AbsSpace) -> MatElem

Returns a matrix `M`, such that the rows of `M` form an orthgonal basis of `V`.
"""
function orthogonal_basis(V::AbsSpace)
  _, B = _gram_schmidt(gram_matrix(V), involution(V))
  return B
end

@doc Markdown.doc"""
    diagonal(V::AbsSpace) -> Vector{FieldElem}

Returns a vector of elements $a_1,\dotsc,a_n$ such that `V` is isometric to
the diagonal space $\langle a_1,\dotsc,a_n \rangle$.

The elements will be contained in the fixed field of `V`.
"""
diagonal(V::AbsSpace)

################################################################################
#
#  Gram-Schmidt
#
################################################################################

# Clean this up
function _gram_schmidt(M::MatElem, a)
  F = deepcopy(M)
  K = base_ring(F)
  n = nrows(F)
  S = identity_matrix(K, n)
  okk = isdiagonal(F)
  if !okk
    for i in 1:n
      if iszero(F[i,i])
        T = identity_matrix(K, n)
        ok = 0
        for j in (i + 1):n
          if !iszero(F[j, j])
            ok = j
            break
          end
        end
        #ok = findfirst(j -> !iszero(F[j, j]), (i + 1):n)
        if ok != 0 # ok !== nothing
          #j = ok + i # findfirst gives the index
          j = ok
          T[i,i] = 0
          T[j,j] = 0
          T[i,j] = 1
          T[j,i] = 1
        else
          ok = 0
          for j in (i + 1):n
            if !iszero(F[i, j])
              ok = j
              break
            end
          end
          if ok == 0
            error("Matrix is not of full rank")
          end
          #ok = findfirst(j -> !iszero(F[i, j]), (i + 1):n)
          #if ok === nothing
          #  error("Matrix is not of full rank")
          #end
          j = ok
          T[i, j] = 1 // (2 * F[j, i])
        end
        S = T * S
        F = T * F * transpose(_map(T, a))
      end
      T = identity_matrix(K, n)
      for j in (i + 1):n
        T[j, i] = divexact(-F[j, i], F[i, i])
      end
      F = T * F * transpose(_map(T, a))
      S = T * S
    end
    @assert isdiagonal(F)
  end
  return F, S
end

################################################################################
#
#  Hasse and Witt invariant
#
################################################################################

# Auxiliary function which works with a diagonal
function _hasse_invariant(D::Vector, p)
  h = 1
  n = length(D)
  for i in 1:n
    for j in (i + 1):n
      h = h * hilbert_symbol(D[i], D[j], p)
    end
  end
  return h
end

@doc Markdown.doc"""
    hasse_invariant(V::QuadSpace, p::Union{InfPlc, NfOrdIdl}) -> Int

Returns the Hasse invariant of the quadratic space `V` at `p`. This is equal
to the product of local Hilbert symbols $(a_i, a_j)_p$, $i < j$, where $V$ is
isometric to $\langle a_1,\dotsc,a_n\rangle$.
"""
function hasse_invariant(V::QuadSpace, p)
  return _hasse_invariant(diagonal(V), p)
end

function hasse_invariant(L::HermSpace, p)
  throw(error("The space must be quadratic"))
end

# This can be refactored to operate on the diagonal of a gram schmidt basis and
# the gram matrix.
# (Probably only on the diagonal of a gram schmidt basis)
function witt_invariant(L::QuadSpace, p::NfOrdIdl)
  h = hasse_invariant(L, p)
  F = gram_matrix(L)
  dett = det(F)
  K = base_ring(L)
  ncolsFmod8 = mod(ncols(F), 8)
  if ncolsFmod8 == 3 || ncolsFmod8 == 4
    c = -dett
  elseif ncolsFmod8 == 5 || ncolsFmod8 == 6
    c = K(-1)
  elseif ncolsFmod8 == 7 || ncolsFmod8 == 0
    c = dett
  else
    c = K(1)
  end
  return h * hilbert_symbol(K(-1), c, p)
end

function witt_invariant(L::QuadSpace, p::InfPlc)
  if iscomplex(p)
    return 1
  end

  h = hasse_invariant(L, p)
  F = gram_matrix(L)
  dett = det(F)
  K = base_ring(L)
  ncolsFmod8 = mod(ncols(F), 8)
  if ncolsFmod8 == 3 || ncolsFmod8 == 4
    c = -dett
  elseif ncolsFmod8 == 5 || ncolsFmod8 == 6
    c = K(-1)
  elseif ncolsFmod8 == 7 || ncolsFmod8 == 0
    c = dett
  else
    c = K(1)
  end
  @assert !iszero(c)
  if isnegative(c, p)
    return -h
  else
    return h
  end
end

@doc Markdown.doc"""
    witt_invariant(V::QuadSpace, p::Union{InfPlc, NfOrdIdl}) -> Int

Returns the Witt invariant of the quadratic space `V` at `p`.

See [Definition 3.2.1, Kir16].
"""
witt_invariant(V::QuadSpace, p)

function witt_invariant(L::HermSpace, p)
  throw(error("The space must be quadratic"))
end

################################################################################
#
#  Local equivalence
#
################################################################################

function isequivalent(L::QuadSpace, M::QuadSpace, p::NfOrdIdl)
  GL = gram_matrix(L)
  GM = gram_matrix(M)
  if GL == GM
    return true
  end

  return rank(GL) == rank(GM) && islocal_square(det(GL) * det(GM), p) && hasse_invariant(L, p) == hasse_invariant(M, p)
end

function isequivalent(L::QuadSpace, M::QuadSpace, p::InfPlc)
  if rank(L) != rank(M)
    return false
  end

  if iscomplex(p)
    return true
  end

  DL = diagonal(L)
  DM = diagonal(M)
  return count(x -> isnegative(x, p), DL) == count(x -> isnegative(x, p), DM)
end

# hermitian case

function isequivalent(L::HermSpace{AnticNumberField}, M::HermSpace{AnticNumberField}, p::fmpz)
  return _isequivalent(L, M, p)
end

function isequivalent(L::HermSpace, M::HermSpace, p::NfOrdIdl)
  return _isequivalent(L, M, p)
end

function _isequivalent(L::HermSpace, M::HermSpace, p)
  base_ring(L) != base_ring(M) && error("Both spaces must have the same base field")
  A = gram_matrix(L)
  B = gram_matrix(M)
  if A == B
    return true
  end

  if rank(L) != rank(M)
    return false
  end

  return islocal_norm(base_ring(L), det(L) * det(M), p)[1]
end

function isequivalent(L::HermSpace, M::HermSpace, P::InfPlc)
  if L == M
    return true
  end

  if iscomplex(P)
    return true
  end

  DL = diagonal(L)
  DM = diagonal(M)
  iL = count(d -> isnegative(d, P), DL)
  iM = count(d -> isnegative(d, P), DM)
  return iL == iM
end

@doc Markdown.doc"""
    isequivalent(L::AbsSpace, M::AbsSpace, p::Union{InfPlc, NfOrdIdl}) -> Bool

Returns whether `L` and `M` are equivalent over the completion at `p`.
"""
isequivalent(L::AbsSpace, M::AbsSpace, p)

################################################################################
#
#  Quadratic form with given invariants
#
################################################################################

function _quadratic_form_invariants(M; minimal = true)
  G, _ = _gram_schmidt(M, identity)
  D = diagonal(G)
  K = base_ring(M)
  O = maximal_order(K)
  sup = Dict{ideal_type(O), Bool}()
  for i in 1:length(D)
    f = factor(D[i] * O)
    for (P, e) in f
      if isodd(e)
        sup[P] = true
      end
    end
  end
  for (P, e) in prime_decomposition(O, 2)
    sup[P] = true
  end
  F = Dict{ideal_type(O), Int}()
  for P in keys(sup)
    e = _hasse_invariant(D, P)
    if e == -1 || !minimal
      F[P] = e
    end
  end
  I = [ (P, count(x -> isnegative(x, P), D)) for P in real_places(K) ];
  return prod(D), F, I
end

@doc Markdown.doc"""
    invariants(M::QuadSpace)
          -> FieldElem, Dict{NfOrdIdl, Int}, Vector{Tuple{InfPlc, Int}}

Returns a triple `(d, H, I)` of invariants of `M`, which determine the
equivalence class completely. The element `d` is the determinant of a Gram
matrix, `H` contains the non-trivial Hasse invariants and `I` contains for
each real place the negative index of inertia.

Note that `d` is determined only modulo squares.
"""
invariants(V::QuadSpace) = _quadratic_form_invariants(gram_matrix(V))

@doc Markdown.doc"""
    isequivalent(M::QuadSpace, L::QuadSpace) -> Bool

Tests if `M` and `L` are equivalent.
"""
function isequivalent(M::QuadSpace, L::QuadSpace)
  if gram_matrix(M) == gram_matrix(L)
    return true
  end
  d1, H1, I1 = invariants(M)
  d2, H2, I2 = invariants(L)
  return I1 == I2 && H1 == H2 && issquare(d1 * d2)[1]
end

function isequivalent(M::HermSpace, L::HermSpace)
  if gram_matrix(M) == gram_matrix(L)
    return true
  end

  if rank(M) != rank(L)
    return false
  end

  E = base_ring(M)
  # I could replace this with a islocal_norm at the ramified primes + primes
  # dividing right hand side
  return isnorm(E, det(M) * det(L))[1]
end

################################################################################
#
#  Definitness
#
################################################################################

# TODO: Add this functionality for Hermitian spaces

# Returns 0 if V is not definite
# Returns an element a != 0 such that a * canonical_basis of V has
# positive Gram matrix
function _isdefinite(V::AbsSpace)
  K = base_ring(V)
  R = maximal_order(K)
  if !istotally_real(K) || (ishermitian(V) && !istotally_complex(K))
    return zero(R)
  end
  D = diagonal(V)
  signs_to_consider = Tuple{InfPlc, Int}[]
  for v in real_places(K)
    S = Int[sign(d, v) for d in D]
    if length(unique(S)) != 1
      return zero(R)
    else
      push!(signs_to_consider, (v, S[1]))
    end
  end
  if length(signs_to_consider) == 1
    return R(signs_to_consider[1][2])
  else
    return element_with_signs(K, signs_to_consider)
  end
end

function ispositive_definite(V::AbsSpace)
  E = base_ring(V)
  if isquadratic(V)
    K = E
  else
    K = base_field(E)
  end
  R = maximal_order(K)
  if !istotally_real(K)
    return false
  end
  D = diagonal(V)
  for d in D
    if !istotally_positive(d)
      return false
    end
  end
  return true
end

function isnegative_definite(V::AbsSpace)
  E = base_ring(V)
  if isquadratic(V)
    K = E
  else
    K = base_field(E)
  end

  if !istotally_real(K)
    return false
  end
  D = diagonal(V)
  for d in D
    if !istotally_positive(-d)
      return false
    end
  end
  return true
end

function isdefinite(V::AbsSpace)
  return ispositive_definite(V) || isnegative_definite(V)
end

################################################################################
#
#  Isotropic
#
################################################################################

function isisotropic(V::QuadSpace, p)
  @assert base_ring(V) == nf(order(p))
  d = det(V)
  n = rank(V)
  K = base_ring(V)
  if d == 0
    return true
  elseif n <= 1
    return false
  elseif n == 2
    return islocal_square(-d, p)
  elseif n == 3
    return hasse_invariant(L, p) == hilbert_symbol(K(-1), K(-1), p)
  elseif n == 4
    return !islocal_square(d, p) || (hasse_invariant(L, p) == hilbert_symbol(K(-1), K(-1), p))
  else
    return true
  end
end

function isisotropic(V::HermSpace, q)
  if nf(order(q)) == base_ring(V)
    p = minimum(q)
  else
    p = q
  end
  @assert fixed_field(V) == nf(order(p))
  r = rank(V)
  if r >= 3
    return true
  elseif r == 0
    return false
  end
  d = det(V)
  if r == 1
    return d == 0
  end
  return islocal_norm(base_ring(V), -d, p)
end

function isisotropic(V::AbsSpace, p::InfPlc)
  n = rank(V)
  d = det(V)
  E = base_ring(L)
  if d == 0
    return true
  elseif n <= 1
    return false
  elseif iscomplex(p)
    return true
  else
    D = diagonal(V)
    return length(unique!([sign(evaluate(d, p)) for d in p])) == 2
  end
end

function _islocally_hyperbolic_hermitian_detclass(rk, d, E, K, p)
  if isodd(rk)
    return false
  end
  if d == 1
    if iseven(div(rk, 2))
      return true
    else
      return islocal_norm(E, K(-1), p)
    end
  else
    if iseven(div(rk, 2))
      return false
    else
      return !islocal_norm(E, K(-1), p)
    end
  end
end

function islocally_hyperbolic(V::HermSpace, p)
  rk = rank(V)
  if isodd(rk)
    return false
  end
  return islocal_norm(base_ring(V), det(V) * (-1)^(div(rk, 2)), p)
end

################################################################################
#
#  Helper functions (sort them later)
#
################################################################################

function image(f::NfRelToNfRelMor{T, T}, I::NfRelOrdFracIdl{T, S}) where {T, S}
  #S has to be an automorphism!!!!
  O = order(I)
  @assert ismaximal(O) # Otherwise the order might change
  K = nf(O)

  pb = pseudo_basis(I)

  z = sum(b * (f(a) * O) for (a, b) in pb)
  return z
end

# An element is locally a square if and only if the quadratic defect is 0, that is
# the valuation is inf.
# (see O'Meara, Introduction to quadratic forms, 3rd edition, p. 160)
function islocal_square(a, p)
  return quadratic_defect(a, p) isa PosInf
end

function _map(a::AbstractAlgebra.MatrixElem, f)
  z = similar(a)
  for i in 1:nrows(a)
    for j in 1:ncols(a)
      z[i, j] = f(a[i, j])
    end
  end
  return z
end

# I think I need a can_change_base_ring version

function element_with_signs(K, D::Dict{InfPlc, Int})
  return _element_with_signs(K, D)
end

function _element_with_signs(K, D)
  OK = maximal_order(K)
  G, mG = carlos_units(OK)
  r = real_places(K)
  z = id(G)
  for (v, s) in D
    for i in 1:length(r)
      if s == 1
        ss = 0
      else
        ss = 1
      end

      if v == r[i]
        z = z + ss * G[i]
      end
    end
  end
  return mG(z)
end

function element_with_signs(K, P::Vector{InfPlc}, S::Vector{Int})
  return _element_with_signs(K, zip(P, S))
end

function element_with_signs(K, A::Vector{Tuple{InfPlc, Int}})
  return _element_with_signs(K, A)
end

function prime_ideals_up_to(O::NfRelOrd, n::Union{Int, fmpz})
  p = 2
  z = ideal_type(O)[]
  while p <= n
    lp = prime_decomposition(base_ring(O), p)
    for q in lp
      if norm(q[1]) > n
        continue
      else
        lq = prime_decomposition(O, q[1])
        for Q in lq
          if absolute_norm(Q[1]) <= n
            push!(z, Q[1])
          end
        end
      end
    end
    p = next_prime(p)
  end
  return sort!(z, by = a -> absolute_norm(a))
end

################################################################################
#
#  Restriction of scalars
#
################################################################################

function restrict_scalars(V::AbsSpace, alpha = one(base_ring(V)))
  E = base_ring(V)
  n = rank(V)
  d = absolute_degree(E)
  B = absolute_basis(E)
  v = elem_type(E)[zero(E) for i in 1:n]
  w = elem_type(E)[zero(E) for i in 1:n]
  G = zero_matrix(FlintQQ, d * n, d * n)
  r = 1
  c = 1
  indices = Vector{Tuple{Int, Int}}(undef, n * d)
  for i in 1:n
    for bi in 1:length(B)
      v[i] = B[bi]
      c = 1
      for j in 1:n
        for bj in 1:length(B)
          w[j] = B[bj]
          G[r, c] = trace(alpha * inner_product(V, v, w), FlintQQ)
          w[j] = zero(E)
          c = c + 1
        end
      end
      v[i] = zero(E)
      indices[r] = (i, bi)
      r = r + 1
    end
  end

  VabstoV = function(v)
    @assert length(v) == d * n
    z = Vector{elem_type(E)}(undef, n)
    for i in 1:n
      z[i] = zero(E)
    end
    for k in 1:(d * n)
      (i, j) = indices[k]
      z[i] = z[i] + v[k] * B[j]
    end
    return z
  end

  VtoVabs = function(w)
    z = Vector{fmpq}(undef, d * n)
    k = 1
    for i in 1:n
      y = w[i]
      for j in 1:d
        z[k] = absolute_coeff(y, j - 1)
        k = k + 1
      end
    end
    return z
  end

  return quadratic_space(FlintQQ, G), VabstoV, VtoVabs
end

# Careful, starts at 0!
function absolute_coeff(z::nf_elem, i)
  return coeff(z, i)
end

function absolute_coeff(z::NumFieldElem, i)
  d = absolute_degree(base_field(parent(z)))
  rowindex = fld(i, d)
  colindex = (i % d)
  return absolute_coeff(coeff(z, rowindex), colindex)
end

function absolute_basis(K::NumField)
  k = base_field(K)
  kabs = absolute_basis(k)
  res = elem_type(K)[]
  for b in basis(K)
    for bb in kabs
      push!(res, bb * b)
    end
  end
  return res
end

function absolute_basis(K::NumField{fmpq})
  return basis(K)
end


text:
using FastIOBuffers: FastReadBuffer
using DataStructures: counter

""" A ROS Bag version X.Y """
abstract type Bag{X,Y} end

""" A chunked ROS bag V2.0+ """
struct ChunkedBag{X,Y,F} <: Bag{X,Y}
    file::F
    header::BagHeader
    indices::Vector{IndexData}
    chunkinfos::Vector{ChunkInfo}
    chunks::Vector{Chunk}
    connections::Vector{Connection}
    type_map::Dict{Int,DataType}
    topic_map::Dict{Int,String}
end
ChunkedBag(X, Y, file, header) = ChunkedBag{X,Y,typeof(file)}(file, header, [], [], [], [], Dict(), Dict())

Base.length(bag::ChunkedBag) = length(bag.chunks)

index!(bag::ChunkedBag, io::IO, rec::Record) = nothing
index!(bag::ChunkedBag, io::IO, rec::IndexData) = push!(bag.indices, rec)
index!(bag::ChunkedBag, io::IO, rec::ChunkInfo) = push!(bag.chunkinfos, rec)
index!(bag::ChunkedBag, io::IO, rec::Connection) = push!(bag.connections, rec)
index!(bag::ChunkedBag, io::IO, rec::Chunk) = begin
    push!(bag.chunks, rec)
    skip(io, rec.data_len)
end

function setup_maps!(bag::ChunkedBag, msg_types...)
    for conn in bag.connections
        # (a) map connection id -> topic
        bag.topic_map[conn.conn] = conn.topic

        # (b) map connection id -> message type
        type_str = conn.header[:type] # eg. "sensor_msgs/Imu"
        pkg, typ = Symbol.(split(type_str, "/"))
        # Do we have a module that implements this message type?
        # Look for eg. type_mod.sensor_msgs.Imu
        for type_mod in msg_types
            dtype = try
                package = getproperty(type_mod, pkg)
                getproperty(package, typ)
            catch
                continue
            end
            bag.type_map[conn.conn] = dtype # found it
            break
        end
    end
end

""" FileIO loader """
function load(f::File{format"ROSBAG"}, msg_types::Vararg{Module}; std_types=(RobotOSData.StdMsgs, RobotOSData.CommonMsgs))
    open(f) do s
        skipmagic(s)  # skip over the magic bytes
        version_str = readuntil(s.io, "\n")
        X,Y = parse.(Int, split(version_str, "."))
        if X == 2
            header = read(s.io, Record)::BagHeader
            bag = ChunkedBag(X, Y, f, header)
            load_bag!(bag, s.io, std_types..., msg_types...)
        else
            error("Bag version $X.$Y is not yet supported")
        end
    end
end

function load_bag!(bag::ChunkedBag, io::IO, msg_types...)
    # Skim through the first-level records (not recursing into chunks)
    # to build an inventory of the bag contents:
    while !eof(io)
        rec = read(io, Record)
        index!(bag, io, rec)
    end
    # Match connections to topics and messages to message types:
    setup_maps!(bag, msg_types...)
    if length(bag.chunks) != length(bag.chunkinfos)
        warn("Bad bag index: ChunkInfos do not match Chunks")
    end
    bag
end

function Base.show(io::IO, bag::ChunkedBag{X,Y}) where {X,Y}
    println(io, "ROS Bag v$X.$Y at $(bag.file.filename)")

    chunk_counter = counter(Symbol)
    chunksize = 0
    for chunk in bag.chunks
        chunksize += chunk.size
        comp = compression(chunk)
        push!(chunk_counter, comp == :none ? :uncompressed : comp)
    end
    unit = ""
    for outer unit in ("bytes", "kiB", "MiB", "GiB", "TiB")
        chunksize < 1024 && break
        chunksize /= 1024
    end
    print(io, join(["$c $k chunks" for (k,c) in chunk_counter], " / "))
    println(io, " with $(round(chunksize; digits=1)) $unit of message data")

    num_msgs = sum(length(index.msgs) for index in bag.indices)
    println(io, "$(num_msgs) messages in $(length(bag.connections)) topics:")
    rows = []
    for conn in bag.connections
        dtype = get(bag.type_map, conn.conn, Vector{UInt8})
        push!(rows, (" └─", conn.topic, ": ", conn.header[:type], " => ", dtype))
    end
    widths = [sum(length.(row[1:4])) for row in rows]
    pads = maximum(widths) .- widths
    for (row, pad) in zip(rows, pads)
        print(io, row[1])
        printstyled(io, row[2], color=:green)
        print(io, row[3])
        print(io, lpad("", pad))
        print(io, row[4])
        print(io, row[5])
        if row[6]==Vector{UInt8}
            printstyled(io, row[6], color=:red)
        else
            print(io, row[6])
        end
        println(io)
    end
end

""" Calls f(rec) for each record in chunk i. """
function replay(f, bag::ChunkedBag, i::Integer, filters=())
    open(bag.file) do s
        replay(f, bag, s.io, i, filters)
    end
end

""" Calls f(rec) for each record in the specified chunks. """
function replay(f, bag::ChunkedBag, indices::AbstractArray, filters=())
    open(bag.file) do s
        for i in indices
            replay(f, bag, s.io, i, filters)
        end
    end
end

function replay(f, bag::ChunkedBag, io::IO, chunk_idx::Integer, filters)
    infos = bag.chunkinfos[chunk_idx]
    if !select(filters, infos) return end
    seek(io, infos.chunk_pos)
    chunk = read(io, Record)::Chunk
    replay(f, bag, io, chunk, filters)
end

function replay(f, bag::ChunkedBag, io::IO, chunk::Chunk{:none}, filters)
    # Uncompressed chunk: we can simply parse the same io stream.
    replay_stream(f, bag, io, chunk.data_len, filters)
end

function replay(f, bag::ChunkedBag, io::IO, chunk::Chunk{:bz2}, filters)
    # Bzip2-compressed chunk: we need to expand the data first.
    compressed = read(io, chunk.data_len)
    expanded = Vector{UInt8}(undef, chunk.size)
    unbzip2!(compressed, expanded)
    chunk_io = FastReadBuffer(expanded)
    replay_stream(f, bag, chunk_io, chunk.size, filters)
end

function replay_stream(f, bag::ChunkedBag, io::IO, len::Integer, filters)
    # Parse and replay up to `len` bytes of records
    io_end = position(io) + len
    while position(io) < io_end
        rec = read(io, Record)
        replay(f, bag, io, rec, filters)
    end
end

function replay(f, bag::ChunkedBag, io::IO, msg::MessageData, filters)
    data_len = msg.data
    if select(filters, msg) # let's parse the message data then:
        dtype = get(bag.type_map, msg.conn, data_len)
        data = read(io, dtype) # defaults to reading data_len bytes
        f(MessageData(msg.conn, msg.time, data))
    else
        skip(io, data_len)
    end
end

function replay(f, bag::ChunkedBag, io::IO, rec::Record, filters)
    nothing # ignore all other record types
end

function read_messages(bag::ChunkedBag, chunk_indices, args...)
    messages = Dict{Int,Vector}()
    filters = make_filter.(Ref(bag), args)
    replay(bag, chunk_indices, filters) do msg
        out = get!(messages, msg.conn) do
            Vector{typeof(msg)}()
        end
        push!(out, msg)
    end
    for f in filters
        if isa(f, ConnectionFilter)
            return messages[f.conn]
        end
    end
    Dict(bag.topic_map[k] => v for (k,v) in pairs(messages))
end

select(filters::Tuple{}, x) = true # empty filters: include all
select(filters::Tuple, x) = all(f->select(f,x), filters)

struct ConnectionFilter
    conn::Int
end
select(filter::ConnectionFilter, msg::MessageData) = msg.conn == filter.conn
select(filter::ConnectionFilter, _) = true

select(filter::ROSTimespan, infos::ChunkInfo) =
    infos.start_time <= filter.end_time && infos.end_time >= filter.start_time
select(filter::ROSTimespan, msg::MessageData) = msg.time in filter
select(filter::ROSTimespan, _) = true

make_filter(bag, ts::ROSTimespan) = ts
make_filter(bag, s::String) = begin
    for (i,topic) in bag.topic_map
        if topic == s
            return ConnectionFilter(i)
        end
    end
    error("No such topic: $s")
end

Base.getindex(bag::ChunkedBag, i::Integer, filters...) = read_messages(bag, i, filters...)
Base.getindex(bag::ChunkedBag, i::AbstractArray, filters...) = read_messages(bag, i, filters...)
Base.getindex(bag::ChunkedBag, ::Colon, filters...) = read_messages(bag, eachindex(bag.chunkinfos), filters...)
Base.getindex(bag::ChunkedBag, filters...) = getindex(bag, :, filters...)

import TranscodingStreams
import TranscodingStreams: Memory, Error
import CodecBzip2: transcode, Bzip2Decompressor

""" Decode Bzip2 with minimal memory allocations """
function unbzip2!(input::Vector{UInt8}, output::Vector{UInt8})
    err = Error()
    codec = Bzip2Decompressor()
    TranscodingStreams.startproc(codec, :x, err)
    TranscodingStreams.haserror(err) && throw(err[])
    din, dout, status = TranscodingStreams.process(
        codec,
        Memory(pointer(input), length(input)),
        Memory(pointer(output), length(output)),
        err)
    TranscodingStreams.haserror(err) && throw(err[])
    status != :end && error("Bzip2: not the end")
    TranscodingStreams.finalize(codec)
    nothing
end

ext: jl


text:
using NamedArrays

n = NamedArray([1 2; 3 4], ([11,12], [13,14]))

s = sum(n, dims=1)

println( n ./ s )


text:
using Plots

Rotation(x, y, θ) = [cos(θ) -sin(θ);sin(θ) cos(θ)] * [x;y]

function KochCurve(; points, order, θ)
    for n in 1:order
        i = 1
        while i < length(points)
            insert!(points,i+1,(points[i+1] - points[i]) /3 + points[i])
            delta = (points[i+2] - points[i]) *2 /3 + points[i] - points[i+1]
            insert!(points, i+2, Rotation(delta[1], delta[2], θ)+points[i+1])
            insert!(points,i+3,(points[i+3] - points[i]) *2 /3 + points[i])
            i += 4
        end
    end
    return points
end

StartPoints = [[[0.0; 0.0],[10.0; 0.0]],
            [[0.0; 0.0],[10.0; 0.0],[5.0; 8.66],[0.0; 0.0]],
            [[0.0; 0.0],[10.0; 0.0],[5.0; 8.66],[0.0; 0.0]],
            [[0.0; 0.0],[10.0; 0.0],[5.0; 8.66],[0.0; 0.0]]
                ]
Adegs = [pi/3, pi/3, -pi/3, 2*pi/3,]

ylimlist = [(0,5), (0,10),(-3,9), (0,10)]
xlimlist = [(0,10), (0,10),(-1,11), (0,10)]

for j in 1:5
    for i in 1:4
        Parameters = Dict(
                    :points => StartPoints[i],
                    :order => j,
                    :θ => Adegs[i])
        points = KochCurve(;Parameters...)
        plot(hcat(points...)[1,:], hcat(points...)[2,:], legend = false, border=:none, fill =(0), ylims = ylimlist[i], xlims = xlimlist[i])
        savefig("C:\\Users\\Yaghoub\\Documents\\GitHub\\ComputationalPhysics-Fall2021\\ProblemSet1\\Figs\\Q1\\KochCurve$(i)O$(j).png")
    end
end


text:
module TrajectoryVisualization
import VegaLite

export visualize_trajectory, visualize_trajectory_batch

#=========================== Visualization Interface for Dynamics  Models ==========================#

function trajectory_data end
"Converts the state matrix x into a tabular format of `trajectory_data` that can be handed to
visualize_trajectory"
trajectory_data(control_system, x, player)

#===================================== generic implementations =====================================#

function visualize_trajectory(control_system, x; kwargs...)
    td = trajectory_data(control_system, x)
    visualize_trajectory(td; kwargs...)
end

"""
Visualizes a single trajectory `x` given as a matrix of states for a preset `control_system`.
"""
function visualize_trajectory(
    trajectory_data;
    canvas = VegaLite.@vlplot(),
    x_position_domain = extrema(s.px for s in trajectory_data) .+ (-0.01, 0.01),
    y_position_domain = extrema(s.py for s in trajectory_data) .+ (-0.01, 0.01),
    draw_line = true,
    legend = nothing,
)
    trajectory_visualizer =
        VegaLite.@vlplot(
            encoding = {
                x = {"px:q", scale = {domain = x_position_domain}, title = "Position x [m]"},
                y = {"py:q", scale = {domain = y_position_domain}, title = "Position y [m]"},
                order = "t:q",
                color = {"player:n", title = "Player", legend = legend},
            }
        ) +
        VegaLite.@vlplot(mark = {"point", shape = "circle", size = 25, clip = true, filled = true})

    if draw_line
        trajectory_visualizer += VegaLite.@vlplot(mark = {"line", clip = true})
    end

    canvas + (trajectory_data |> trajectory_visualizer)
end

"""
Generic visualization of a single tarjectory `trajectory_data`. Here, the `trajectory_data` is
provided in a tabular format and is self-contained. That is, the table-liek format must have columns
have positions `px` and `py`, time `t`, and player identifier `player`.
"""
function visualize_trajectory_batch(
    control_system,
    trajectory_batch;
    canvas = VegaLite.@vlplot(opacity = {value = 0.2}, width = 200, height = 200),
    kwargs...,
)
    mapreduce(+, trajectory_batch; init = canvas) do x
        visualize_trajectory(control_system, x; kwargs...)
    end
end

function visualize_trajectory_batch(
    trajectory_data_batch;
    canvas = VegaLite.@vlplot(opacity = {value = 0.2}, width = 200, height = 200),
    kwargs...,
)
    mapreduce(+, trajectory_data_batch; init = canvas) do trajectory_data
        visualize_trajectory(trajectory_data; kwargs...)
    end
end

end


text:
# All subtypes of AbstractAffectFeedForward must define the field affected_variables.
# TODO: make a unit test that checks for this.
abstract type AbstractAffectFeedForward end

abstract type AbstractCache end
abstract type FeedForwardChronology end

get_trigger(val::FeedForwardChronology) = val.trigger

abstract type AbstractOperationsProblem end
abstract type PowerSimulationsOperationsProblem <: AbstractOperationsProblem end


text:
function unit(::Type{T}, i::Integer) where T <: StaticVector
    T(ntuple(Val(length(T))) do j
        ifelse(i == j, 1, 0)
    end)
end

macro fixed_vector(name, parent)
    esc(quote

        struct $(name){S, T} <: $(parent){S, T}
            data::NTuple{S, T}

            function $(name){S, T}(x::NTuple{S,T}) where {S, T}
                new{S, T}(x)
            end

            function $(name){S, T}(x::NTuple{S,Any}) where {S, T}
                new{S, T}(StaticArrays.convert_ntuple(T, x))
            end
        end

        size_or(::Type{$(name)}, or) = or
        eltype_or(::Type{$(name)}, or) = or
        eltype_or(::Type{$(name){S, T} where S}, or) where {T} = T
        eltype_or(::Type{$(name){S, T} where T}, or) where {S} = or
        eltype_or(::Type{$(name){S, T}}, or) where {S, T} = T

        size_or(::Type{$(name){S, T} where S}, or) where {T} = or
        size_or(::Type{$(name){S, T} where T}, or) where {S} = Size{(S,)}()
        size_or(::Type{$(name){S, T}}, or) where {S, T} = (S,)

        # Array constructor
        function $(name){S}(x::AbstractVector{T}) where {S, T}
            @assert S <= length(x)
            $(name){S, T}(ntuple(i-> x[i], Val(S)))
        end

        function $(name){S, T1}(x::AbstractVector{T2}) where {S, T1, T2}
            @assert S <= length(x)
            $(name){S, T1}(ntuple(i-> T1(x[i]), Val(S)))
        end

        function $(name){S, T}(x) where {S, T}
            $(name){S, T}(ntuple(i-> T(x), Val(S)))
        end

        $(name){S}(x::T) where {S, T} = $(name){S, T}(ntuple(i-> x, Val(S)))
        $(name){1, T}(x::T) where T = $(name){1, T}((x,))
        $(name)(x::NTuple{S}) where {S} = $(name){S}(x)
        $(name)(x::T) where {S, T <: Tuple{Vararg{Any, S}}} = $(name){S, StaticArrays.promote_tuple_eltype(T)}(x)

        $(name){S}(x::T) where {S, T<:Tuple} = $(name){S, StaticArrays.promote_tuple_eltype(T)}(x)
        $(name){S, T}(x::StaticVector) where {S, T} = $(name){S, T}(Tuple(x))

        @generated function (::Type{$(name){S, T}})(x::$(name)) where {S, T}
            idx = [:(x[$i]) for i = 1:S]
            quote
                $($(name)){S, T}($(idx...))
            end
        end

        @generated function Base.convert(::Type{$(name){S, T}}, x::$(name)) where {S, T}
            idx = [:(x[$i]) for i = 1:S]
            quote
                $($(name)){S, T}($(idx...))
            end
        end

        @generated function (::Type{SV})(x::StaticVector) where SV <: $(name)
            len = size_or(SV, size(x))[1]
            if length(x) == len
                :(SV(Tuple(x)))
            elseif length(x) > len
                elems = [:(x[$i]) for i = 1:len]
                :(SV($(Expr(:tuple, elems...))))
            else
                error("Static Vector too short: $x, target type: $SV")
            end
        end

        Base.@pure StaticArrays.Size(::Type{$(name){S, Any}}) where {S} = Size(S)
        Base.@pure StaticArrays.Size(::Type{$(name){S, T}}) where {S,T} = Size(S)

        Base.@propagate_inbounds function Base.getindex(v::$(name){S, T}, i::Int) where {S, T}
            v.data[i]
        end

        Base.Tuple(v::$(name)) = v.data
        Base.convert(::Type{$(name){S, T}}, x::NTuple{S, T}) where {S, T} = $(name){S, T}(x)
        function Base.convert(::Type{$(name){S, T}}, x::Tuple) where {S, T}
            return $(name){S, T}(convert(NTuple{S, T}, x))
        end

        @generated function StaticArrays.similar_type(::Type{SV}, ::Type{T}, s::Size{S}) where {SV <: $(name), T, S}
            if length(S) === 1
                $(name){S[1], T}
            else
                StaticArrays.default_similar_type(T,s(),Val{length(S)})
            end
        end

        Base.:(*)(a::$name, b::$name) = a .* b
        Base.broadcasted(f, a::AbstractArray{T}, b::$name) where T <: $name = Base.broadcasted(f, a, (b,))
    end)
end

abstract type AbstractPoint{Dim, T} <: StaticVector{Dim, T} end
@fixed_vector Point AbstractPoint
@fixed_vector Vec StaticVector

const Mat = SMatrix
const VecTypes{N, T} = Union{StaticVector{N, T}, NTuple{N, T}}
const Vecf0{N} = Vec{N, Float32}
const Pointf0{N} = Point{N, Float32}
Base.isnan(p::StaticVector) = any(x-> isnan(x), p)

#Create constes like Mat4f0, Point2, Point2f0
for i=1:4
    for T=[:Point, :Vec]
        name = Symbol("$T$i")
        namef0 = Symbol("$T$(i)f0")
        @eval begin
            const $name = $T{$i}
            const $namef0 = $T{$i, Float32}
            export $name
            export $namef0
        end
    end
    name = Symbol("Mat$i")
    namef0 = Symbol("Mat$(i)f0")
    @eval begin
        const $name{T} = $Mat{$i,$i, T, $(i*i)}
        const $namef0 = $name{Float32}
        export $name
        export $namef0
    end
end

export Mat, Vec, Point, unit
export Vecf0, Pointf0


text:
using Sario
using Dates
using Test
using Documenter

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

@testset "loading" begin
    demfile = joinpath(@__DIR__, "elevation.dem")

    dem = Sario.load(demfile)
    rows, cols = (5, 6)

    @test dem == reshape(collect(Int16, 1:30), rows, cols)
    #  1   6  11  16  21  26
    #  2   7  12  17  22  27
    #  3   8  13  18  23  28
    #  4   9  14  19  24  29
    #  5  10  15  20  25  30

    @test Sario.load(demfile, (3, 4)) == dem[3, 4] == 18

    @test_throws DomainError Sario.load(demfile, (30, 40))

    @test Sario.load(demfile, (2:4, 2:4)) == dem[2:4, 2:4]
    @test Sario.load(demfile, (2:4, 2:4)) == Sario.load(demfile)[2:4, 2:4]

    dembig = joinpath(@__DIR__, "elevation_toobig.dem")
    @test_throws AssertionError Sario.load(dembig)

    expected_ints = [(Date(2015,1,1), Date(2016,1,1)),
                     (Date(2016,1,1), Date(2017,1,1))]

    @test Sario.parse_intlist_strings([("20150101", "20160101"),
                                       ("20160101", "20170101")]) == expected_ints
    @test Sario.parse_intlist_strings(["20150101_20160101.int", "20160101_20170101.int"]) == expected_ints

    # Unws
    unwfile = joinpath(@__DIR__, "testfile.unw")
    @test Sario.load(unwfile)[1:2, 1:3] ≈ [2.07206  1.57548  1.33183 ; 2.09438  2.00096  1.40108]

    intfile = joinpath(@__DIR__, "testfile.int")
    @test Sario.load(intfile)[1:2, 1:3] ≈ [-42739.2+77998.7im -386.675+82609.3im 23287.7+95590.0im;
                                           -53225.5+92191.6im -38782.8+84526.9im 12128.7+70777.0im];
    # TODO: add cc, other types
end

@testset "extras" begin
    # Take looks test
    a = [.1  0.01  2 ; 3  4  1 + 1im]

    @test all(Sario.take_looks(a, 2, 1) .≈ [1.55  2.005  1.5 + 0.5im])
    @test all(Sario.take_looks(a, 1, 2) .≈ [0.055; 3.5])

    out = zeros(eltype(a), 2, 1)
    @test all(Sario.take_looks!(out, a, 1, 2, 2, 1) .≈ [0.055; 3.5])

    b = reshape(collect(1:12), 3, 4)
    @test all(Sario.take_looks(b, 1, 2) .≈ [2   8 ; 4  10 ; 4  10 ])
end

@testset "DemRsc" begin
    demrsc = DemRsc(x_first=1.0, y_first=2.0, x_step=0.1, y_step=-0.2, file_length=6, width=5)
    @test demrsc.width == demrsc.cols == 5
    @test demrsc.file_length == demrsc.rows == 6
end


text:
using Test
using CausalDiscovery
using Random

aexpr = au"""(program
(= GRID_SIZE 16)

(object Particle (Cell 0 0 "blue"))

(: particles (List Particle))
(= particles
   (initnext (list)
             (updateObj (prev particles) (--> obj (Particle (uniformChoice (adjPositions (.. obj origin))))))))

(on clicked (= particles (addObj (prev particles) (Particle (Position (.. click x) (.. click y))))))
)"""
#
compiledMod = eval(compiletojulia(aexpr))

# time 0
state = compiledMod.init(nothing, nothing, nothing, nothing, nothing, MersenneTwister(0))
@test compiledMod.render(state.scene) == []

# time 1
state = compiledMod.next(state, nothing, nothing, nothing, nothing, nothing)
@test compiledMod.render(state.scene) == []

# time 2
state = compiledMod.next(state, compiledMod.Click(5,5), nothing, nothing, nothing, nothing)
@test compiledMod.render(state.scene) == [compiledMod.Cell(compiledMod.Position(5, 5), "blue", 0.8)]

# time 3
state = compiledMod.next(state, nothing, nothing, nothing, nothing, nothing)
@test compiledMod.render(state.scene) == [compiledMod.Cell(compiledMod.Position(4, 5), "blue", 0.8)]

# time 4
state = compiledMod.next(state, nothing, nothing, nothing, nothing, nothing)
@test compiledMod.render(state.scene) == [compiledMod.Cell(compiledMod.Position(3, 5), "blue", 0.8)]


text:
module LRHousing

using Chain: @chain
using Random
using StatsBase
using StatsBase.Statistics
using DataFrames
using CSV
using GLM
include("cor.jl")
include("shapirowilk.jl")
include("preset.jl")
include("visualize.jl")
include("anomalies.jl")
include("transformation.jl")

function read_data(path; missingstrings=["NA"], normalizenames=true)
	CSV.read(path, DataFrame, missingstrings=missingstrings, normalizenames=normalizenames)
end

"""
	fill_missing(X, stat)

Fill missing inside array with statistic.

## Example
```julia
using Statistics
fill_missing([1,2,3], mean)
```
"""
function fill_missing(X, stat)
	@chain begin
		X
		replace(_, missing => stat(filter(!ismissing, _)))
		convert.(nonmissingtype(eltype(_)), _)
	end
end

function fill_missing(df::AbstractDataFrame, numeric=detect_numeric(df))
	df = copy(df)
	for name in numeric
		df[!, name] = fill_missing(df[!, name], mean)
	end
	return df
end

function formula(columns, target; noconst=true)
	lhs = Term(Symbol(target))
	columns = filter(!isequal(target |> String), String.(columns))
	rhs = @. Term(Symbol(columns))
	if noconst
		FormulaTerm(lhs, (rhs...,))
	else
		FormulaTerm(lhs, (ConstantTerm(1), rhs...))
	end
end

function fit(data, columns, target; noconst=true)
	f = formula(columns, target; noconst=noconst)
	lm(f, data)
end

function rmse(pred, y; normalize=true)
	n = length(y)
	r = sqrt(mean(@. (pred - y)^2) / n)
	if normalize
		#= maxy = max(pred..., y...) =#
		#= miny = min(pred..., y...) =#
		#= r / (maxy - miny) =#
		r / iqr(y)
	else
		r
	end
end

function kfold_rmse(data, preset, k = 10; normalize=true, logging=false)
	n = size(data, 1)
	idx = collect(Iterators.partition(randperm(n), n÷(k-1)))
	tf = pp(data, preset)
	columns = target_columns(tf)
	f = formula(columns, preset.target)

	R = map(enumerate(idx)) do (k, I)
		train = derive(data[Not(I), :], tf)
		test = derive(data[I, :], tf)
		model = lm(f, train)
		pred = predict(model, test)
		r = rmse(pred, test[!, preset.target]; normalize=normalize)
		if logging
			@show (k, r)
		end
		r
	end
	mean(R)
end

function kfold_rmse(train_function, data, target, k=10; normalize=true, logging=false)
	n = size(data, 1)
	idx = collect(Iterators.partition(randperm(n), n÷(k-1)))
	RMSE = map(enumerate(idx)) do (k, I)
		train = data[Not(I), :]
		test = data[I, :]
		model = train_function(train)
		pred = predict(model, test)
		r = rmse(pred, test[!, target]; normalize=normalize)
		if logging
			@show (k, r)
		end
		r
	end
	mean(RMSE)
end

function rsquared(model, df::AbstractDataFrame)
	y = predict(model, df)
	rsquared(model, y)
end
function rsquared(model, y)
	my = mean(y)
	stot = sum(@. (y - my)^2)
	sres = sum(residuals(model).^2)
	1 - sres/stot
end

export lm

end # module


text:
### A Pluto.jl notebook ###
# v0.19.5

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ 82e0613e-b50e-11ec-0a7e-050ecf0c92ca
begin
	using Pkg
	using Revise
	Pkg.activate(dirname(Base.current_project(@__DIR__)))
	using PlutoUI, Portinari
end

# ╔═╡ dba679be-b646-48b2-ad36-e257fcc100f4
@bind sz Slider(1:100; show_value=true)

# ╔═╡ 6c52e239-e5bd-47c2-af70-19be6545acb3
x = collect(1:sz)

# ╔═╡ 4d754375-9fa1-43f6-880b-25c1b57e4ddf
y = rand(sz)

# ╔═╡ 78c862a6-a979-430b-8d46-1e1a94a394a9
line = Line(x, y, "base-line";
	attributes=D3Attr(style=(;fill="none", stroke="green", var"stroke-width"="3"))
)

# ╔═╡ 83870a8f-0f9b-4b0f-8d28-a8e1820e3b87
val = rand(sz) .* 500

# ╔═╡ 61ae25d3-e6c1-471b-b029-303f29c88970
shapes = Shape(x, y, val, "base-shapes";
	attributes=D3Attr(style=(;fill="rgba(255, 0, 0, 0.3)"))
)

# ╔═╡ 1149fcb4-909c-4250-bb69-39dff8a44c47
together = Context(
	(;domain=([extrema(x)...]), range=[0, 630]),
	(;domain=([extrema(y)...]), range=[0, 300]),
	[line, shapes],
	"together"
)

# ╔═╡ fa9c38b1-dbed-4fbd-a0a0-6ca80d38c855
together_rev = Context(
	(;domain=([extrema(x)...]), range=[0, 630]),
	(;domain=(reverse([extrema(y)...])), range=[0, 300]),
	[line, shapes],
	"together"
)

# ╔═╡ 6de965b1-62d8-491e-b632-386d8f7a0fca
merge = Context(
	(;domain=[0, 0], range=[50, 600]),
	(;domain=[0, 0], range=[50, 300]),
	[
		Scale(together, (0.0, 1.0), (0.0, 0.5), "together-scale"),
		Scale(together_rev, (0.0, 1.0), (0.5, 0.8), "together-rev-scale"),
	],
	"merge"
)

# ╔═╡ 6516c1ba-bbf3-466f-8619-4fa61838a5b0
Context(
	(;domain=([extrema(x)...]), range=[0, 630]),
	(;domain=([extrema(y)...]), range=[0, 300]),
	[
		line,
		Axis(Portinari.Right, "testaxis"; attributes=D3Attr(;
			style=(;transform="translateX(315px)")
		)),
		Axis(Portinari.Bottom, "testaxisy"; attributes=D3Attr(;
			style=(;transform="translateY(150px)")
		))
	],
	"test";
	drawAxis=false
)

# ╔═╡ Cell order:
# ╠═82e0613e-b50e-11ec-0a7e-050ecf0c92ca
# ╟─6c52e239-e5bd-47c2-af70-19be6545acb3
# ╟─4d754375-9fa1-43f6-880b-25c1b57e4ddf
# ╠═83870a8f-0f9b-4b0f-8d28-a8e1820e3b87
# ╠═78c862a6-a979-430b-8d46-1e1a94a394a9
# ╠═61ae25d3-e6c1-471b-b029-303f29c88970
# ╠═1149fcb4-909c-4250-bb69-39dff8a44c47
# ╠═fa9c38b1-dbed-4fbd-a0a0-6ca80d38c855
# ╠═6de965b1-62d8-491e-b632-386d8f7a0fca
# ╟─dba679be-b646-48b2-ad36-e257fcc100f4
# ╠═6516c1ba-bbf3-466f-8619-4fa61838a5b0


text:
function evalsingletsresonance(u::T, d::Vector{T},
    αs_singlets::Vector{T}, Ωs_singlets,
    βs_singlets, λ0::T, λ_multipliers::Vector{T}) where T <: Real

    u_rad = 2*π*u

    out = zeros(Complex{T}, length(αs_singlets))
    for i = 1:length(αs_singlets)
        τ = u_rad - d[i]

        λ = λ0*λ_multipliers[i]
        Ω = Ωs_singlets[i]
        out[i] = αs_singlets[i]*exp(im*βs_singlets[i])/(λ+im*(τ-Ω))
    end

    #return out # type stable.
    return tuple(out...) # type unstable.
end

# ::Vector{Vector{Complex{T}}} # if using type stable.
# type unstable.
function evalitpproxysysresonance(qs::Vector{Vector{Function}},
    u::T, d::Vector{T}, κs_λ::Vector{T},
    κs_β::Vector{Vector{T}}) where T <: Real

    @assert length(d) == length(qs)

    #out = collect( zeros(Complex{T}, length(qs[i])) for i = 1:length(qs) ) # type stable.

    N = sum( length(qs[i]) for i = 1:length(qs) )
    out = zeros(Complex{T}, N)
    n = 0

    u_rad = 2*π*u
    for i = 1:length(qs)
        r = u_rad - d[i]

        for k = 1:length(qs[i])
            #out[i][k] = qs[i][k](r, κs_λ[i], κs_β[i]) # type stable

            n += 1
            out[n] = qs[i][k](r, κs_λ[i])
        end
    end

    #return out # type stable.
    return tuple(out...) # type unstable.
end

function assembleqsevals(z1::NTuple{N, Vector{Complex{T}}}, qs::Vector{Vector{Function}}) where {N,T}
    #
    out_qs = Vector{Vector{Vector{Complex{T}}}}(undef, length(qs))

    n = 0
    for i = 1:length(qs)
        out_qs[i] = Vector{Vector{Complex{T}}}(undef, length(qs[i]))

        for k = 1:length(qs[i])

            n += 1
            out_qs[i][k] = z1[n]
        end
    end

    return out_qs
end

#
function visualizesubgroups(P, qs_evals, q_singlets_evals, f::Function)

    for i = 1:length(qs_evals)
        for k = 1:length(qs_evals[i])

            PyPlot.plot(P, f.(qs_evals[i][k]), label = "sub-group ($(i),$(k))")
        end
    end

    for k = 1:length(q_singlets_evals)

        PyPlot.plot(P, f.(q_singlets_evals[k]), label = "singlet $(k)")
    end
end

"""
save resonance groupings of a compound. f is real, imag, or abs.
"""
function savefiggroups(save_path::String,
    title_string::String,
    P,
    q_U::Vector{Complex{T}},
    qs_U,
    q_singlets_U,
    f::Function) where T

    plot_obj = Plots.plot( P,
        f.(q_U),
        title = title_string,
        label = "sum of sub-models",
        seriestype = :line,
        ticks = :native,
        xlims = (P[1],P[end]),
        hover = P,
        linewidth = 4,
        xflip = true,
        size = (1600, 900))

    #
    for i = 1:length(qs_evals)
        for k = 1:length(qs_evals[i])

            Plots.plot!(plot_obj, P, f.(qs_evals[i][k]), label = "sub-group ($(i),$(k))",
                seriestype = :line,
                linestyle = :dot,
                xflip = true,
                linewidth = 4)
        end
    end

    for k = 1:length(q_singlets_evals)

        Plots.plot!(plot_obj, P, f.(q_singlets_evals[k]), label = "singlet $(k)",
            seriestype = :line,
            linestyle = :dot,
            xflip = true,
            linewidth = 4)
    end

    Plots.savefig(plot_obj, save_path)
end


text:
module AWG5K
using ManualScraper: Parser

const SYM = :awg5k

# Directories to pull from inside ENA's help file.
const DIRS = ["/"]

# Files to pull from inside ENA's help file.
const FILTER = r"(?<=\/)[0-9a-z\-\.]+.htm$"

# This is necessary due to inconsistencies in the HTML files.
const INIT = [

    # Remove control chars and extended chars except CR and LF.
    (r"[^\x20-\x7e\x0a\x0d]", ""),

    # Remove attributes from tags
    (r"<([A-Za-z0-9]+) [^>]+>", s"<\1>"),

    # Remove extraneous tags
    ("", ""),
    ("", ""),
    ("
", ""),
    ("
", ""),
    ("
", ""),
    ("
", ""),

    # Contract whitespace and detritus
    (r":\s+", ":"),
    (r"\s+:", ":"),
    (r"\s+\|", "|"),
    (r"\|\s+", "|"),
    (r"=\s+","="),
    (r"\s+=","="),
    (" [1|2|3|4]", ""),
    (r"
SLISt:SUBSequence:TSTamp\?<\/p>\s+<\/Syntax>", "</Syntax>"),

    # convenience
    ("[,<StartIndex>[,<Size>]]", ""),
    ("[,", ",["),
    ("msus", "string"),
    ("subseq_name", "string"),
]

# Determine if the commands are for read/write, read, or write
const GETSET_REGEX = r""    # cannot be determined in idiomatic way

# How to parse the AWG's HTML files for commands and args.
# l1. Preceding line must end with <Syntax>
# l2. Start with whitespace and
: then follow with anything but whitespace
# l3. optionally some arguments
# l4. end with

# We will consume the lines as they are read to get multiple commands out
const CMD_REGEX = r"(?<=<Syntax>)(?:\s+
([^\s]+)(?: ([^<]+))?<\/p>)"

const CMD_PASS1 = [
    # Strip optional segments of command
    (r"\[SOURce\[([1n])\]\]", s":SOURce[\1]"),
    (r"\[:[:A-Za-z\|]+\]", s"")
]

const CMD_PASS2 = [

    # Save infixes
    (r"\[[a-z0-9\-\|]+\]", "[*]"),
    # but strip lowercase parts
    (r"[a-z]+", ""),

    # Replace infixes. Could be automated...
    (r"(MARK)\[\*\]", s"\1m"),
    (r"(SOUR)\[\*\]", s"\1ch"),
    (r"(DC)\[\*\]",   s"\1dc"),
    (r"(OUTP)\[\*\]", s"\1outp"),
    (r"(DOUT)\[\*\]", s"\1dout"),
    ("[:WIND[*]]",    ":WINDwind"),
    (r"(ELEM)\[\*\]", s"\1el"),
    (r"(WAV)\[\*\]",  s"\1wv"),

]

const ARG_REGEX = r"(?:<)?([A-Za-z0-9_\s]+)(?:>)?::=(?:<|\{)([A-Za-z0-9_\|]+)(?:>|\})"

const ARG_PASS = [
    (r"[\{\}]", ""),
    (r"&[lg]t;", ""),
    ("NR1", "v::Integer"),
    ("NRf", "v::Integer"),
    ("NR2", "v::Real"),
    ("NR3", "v::Real"),
    (r"\[?[Ss]tring\]?", "v::AbstractString"),
    (r"[Bb]oolean", "v::Bool"),
    (r"[Ll]ength", "v::Integer"),
]

const RETTYPE_REGEX = r""

const PARSER = Parser(INIT, GETSET_REGEX, CMD_REGEX, RETTYPE_REGEX,
    CMD_PASS1, CMD_PASS2, ARG_REGEX, ARG_PASS)

end


text:
module StaticJulia

# a ver
# Sugar defines a couple of utility functions we can use to work with typed ASTs
using LLVM, Sugar

using Sugar: LazyMethod, returntype, resolve_func, expr_type, getfunction,
    method_nargs, sugared

using LLVM: ref, Builder, Value, LLVMType, API

"""
Translates a julia type to a LLVM inbuild type
"""
julia_type_to_llvm{T}(t::T) = julia_type_to_llvm(T)

function julia_type_to_llvm(t::DataType)
    ref = ccall(:julia_type_to_llvm, LLVM.API.LLVMTypeRef, (Any, Ptr{Cint}), t, C_NULL)
    LLVM.LLVMType(ref)
end

# cannot run this example if we don't use Julia's llvm instance
LLVM.libllvm_exclusive && error("llvm is exclusive")

# just have these global for now!
const jlctx = LLVM.Context(cglobal(:jl_LLVMContext, Void))

immutable IRContext
    method::LazyMethod
    slotmap::Dict
    builder
    llvmfunc
    mod
    context
end

# Wasn't wrapped in LLVM, so I added it here for now!
# zero extend (widen)
zext!(builder::Builder, val::Value, dest_type::LLVMType, name::String="") =
    LLVM.Instruction(API.LLVMBuildZExt(ref(builder), ref(val), ref(dest_type), name))

zext!(builder::Builder, val::Value, dest_type::DataType, name::String="") =
    zext!(builder, val, julia_type_to_llvm(dest_type), name)

"""
Processes individual ast nodes and emits instructions for it
"""
function process_ast!(ctx, node)
    if Base.is_linenumber(node) || isa(node, Void)
        return # ignore
    elseif isa(node, Expr)
        head, args = node.head, node.args
        if head == :block
            for elem in args
                process_ast!(ctx, elem)
            end
            return
        elseif head == :call
            f = resolve_func(ctx.method, args[1])
            func_args = args[2:end]
            arg_types = map(x-> expr_type(ctx.method, x), func_args)

            # recursively process arguments
            llvm_args = map(x-> process_ast!(ctx, x), func_args)

            # Inbuilds:
            if (f == (+)) && length(func_args) == 2 && all(x-> x <: Integer, arg_types)
                return add!(ctx.builder, llvm_args...)
            elseif (f == (+)) && length(func_args) == 2 && all(x-> x <: AbstractFloat, arg_types)
                return fadd!(ctx.builder, func_args...)
            elseif f == (==) && length(func_args) == 2 && all(x-> x <: Integer, arg_types)
                tmp = icmp!(ctx.builder, LLVM.API.LLVMIntEQ, llvm_args...)

                return  zext!(ctx.builder, tmp, Bool)
            else
                # it's not an inbuild, so we need to compile it!
                llvm_fun = compile_llvm_func(f, arg_types, ctx.mod, jlctx)
                return call!(ctx.builder, llvm_fun, llvm_args) # emit call
            end
        elseif head == :unless
            # TODO implement
            error("Unless not implemented")
            #=
                [block1]
                <unless x> <goto 10>
                [block_unless]
                <label 10>
                [block3]
            =#
            condition = args[1]
            condition_llvm = process_ast!(ctx, condition)
            current_block = position(ctx.builder)
            ifblock_expr = args[2]
            elseblock_expr = length(args) == 3 ? args[3] : Expr(:block)
            ifblock = BasicBlock(ctx.llvmfunc, "then", ctx.context)
            continueblock = BasicBlock(ctx.llvmfunc, "continue", ctx.context)
            position!(ctx.builder, ifblock)
            process_ast!(ctx, ifblock_expr)
            position!(ctx.builder, current_block)
            return br!(ctx.builder, condition_llvm, ifblock, elseblock)
        elseif head == :(=)
            lhs, rhs = args
            @assert isa(lhs, Slot) # left hand should always be a slot
            llvm_value = process_ast!(ctx, rhs)
            ctx.slotmap[lhs] = llvm_value # update slot with the correct llvm value
            return llvm_value
        elseif head == :return
            ret_args = if length(args) == 1
                return ret!(ctx.builder, process_ast!(ctx, args[1]))
            elseif isempty(args)
                return #ingore empty return
            end
        end
    elseif isa(node, Slot)
        return ctx.slotmap[node] # get llvm value corresponding to slot
    end
    error("Node: $node not supported")
end

function compile_llvm_func(f::Function, param_types, mod, llvmctx = jlctx)
    # lazy method is a helper type from Sugar, which makes it easier
    # to get an AST and type information from a function
    method = LazyMethod((f, (param_types...)))
    compile_llvm_func(method, param_types, mod, llvmctx)
end

function compile_llvm_func(method::LazyMethod, param_types, mod, llvmctx = jlctx)

    ltypes = julia_type_to_llvm.(param_types)
    # needs to be a vector
    param_types_llvm = [ltypes...]
    println(method.signature)
    return_type = Core.Inference.return_type(method.signature...)
    return_type_llvm = julia_type_to_llvm(return_type)

    fun_type = LLVM.FunctionType(return_type_llvm, param_types_llvm)
    name = string(Symbol(getfunction(method)))
    llvmfun = LLVM.Function(mod, name, fun_type)

    # insert slot mappings for argument slots
    slotmap = Dict{Any, Any}(map(2:method_nargs(method)) do i
        SlotNumber(i) => parameters(llvmfun)[i - 1]
    end)
    # generate IR
    Builder(llvmctx) do builder
        ctx = IRContext(method, slotmap, builder, llvmfun, mod, llvmctx)
        entry = BasicBlock(llvmfun, "entry", llvmctx)
        position!(builder, entry)
        # get the body ast, in a decompiled form (we could also try
        # to directly work with Julia's lowered AST):
        body = sugared(method.signature...)
        process_ast!(ctx, body)
        println(mod)
    end
    verify(llvmfun)
    # optionally add attributes to function
    # push!(function_attributes(llvmfun), EnumAttribute("alwaysinline"))
    llvmfun
end

"""
Generates llvm ir, compiles it and returns it wrapped in a julia function
"""
function llvm_compile(f, param_types, mod)
    llvmfun = compile_llvm_func(f, param_types, mod)
    return_type = Core.Inference.return_type(f, Tuple{param_types...})
    argnames = [:($(Symbol(string("arg", i)))) for i in 1:length(param_types)]
    jlfun = quote
        function $(Symbol(name))($(argnames...))
            Base.llvmcall($(LLVM.ref(llvmfun)), $return_type, Tuple{$(param_types...)}, $(argnames...))
        end
    end
    eval(jlfun)
end

llvm_module(name::String) = LLVM.Module(name, jlctx)

export llvm_compile, llvm_module

end # module


text:
using CompScienceMeshes, BEAST
o, x, y, z = euclidianbasis(3)

sol = 5.0
#Δt, Nt = 0.08/sol,400
Δt, Nt = 0.12/sol,400

#D, Δx = 1.0, 0.25
D, Δx = 1.0, 0.35
Γ = meshsphere(D, Δx)
X = raviartthomas(Γ)

T = timebasisshiftedlagrange(Δt, Nt, 3)
U = timebasisdelta(Δt, Nt)

V = X ⊗ T
W = X ⊗ U

duration, delay, amplitude = 8.0/sol, 12.0/sol, 1.0
gaussian = creategaussian(duration, delay, duration)

direction, polarisation = z, x
E = planewave(polarisation, direction, derive(gaussian), sol)
T = MWSingleLayerTDIO(sol,-1/sol,-sol,2,0)

#Z = assemble(T,W,V);

@hilbertspace j
@hilbertspace j′
tdefie = @discretise T[j′,j] == -1E[j′]   j∈V  j′∈W
xefie = solve(tdefie)

# using PlotlyJS
# include(Pkg.dir("CompScienceMeshes","examples","plotlyjs_patches.jl"))
#
#
# Xefie, Δω, ω0 = fouriertransform(xefie, Δt, 0.0, 2)
# ω = collect(ω0 + (0:Nt-1)*Δω)
# _, i1 = findmin(abs.(ω-1.0*sol))
#
# ω1 = ω[i1]
# ue = Xefie[:,i1] / fouriertransform(gaussian)(ω1)
#
# fcre, geo = facecurrents(ue, X)
# t2 = patch(geo, real.(norm.(fcre)))
# PlotlyJS.plot(t2)


text:
using Base.Test
using NPZ, Compat

Debug = false

tmp = mktempdir()
if Debug
	println("temporary directory: $tmp")
end

TestArrays = Any[
	true,
	false,
	@compat(Int8(-42)),
	@compat(Int16(-42)),
	@compat(Int32(-42)),
	@compat(Int64(-42)),
	@compat(UInt8(42)),
	@compat(UInt16(42)),
	@compat(UInt32(42)),
	@compat(UInt64(42)),
	@compat(Float32(3.1415)),
	@compat(Float64(3.1415)),
	@compat(Complex64(1, 7)),
	@compat(Complex128(1, 7)),
	Bool[0, 1, 0, 1, 1, 0],
	Int8[-42, 0, 1, 2, 3, 4],
	Int16[-42, 0, 1, 2, 3, 4],
	Int32[-42, 0, 1, 2, 3, 4],
	Int64[-42, 0, 1, 2, 3, 4],
	UInt8[0, 1, 2, 3, 4],
	UInt16[0, 1, 2, 3, 4],
	UInt32[0, 1, 2, 3, 4],
	UInt64[0, 1, 2, 3, 4],
	Float64[-42, 0, 1, 2, 3.14, 4],
	Float32[-42, 0, 1, 2, 3.14, 4],
	Complex64[1+2im, 3, 4+5im, 6im, 7+8im],
	Complex128[1+2im, 3, 4+5im, 6im, 7+8im],
	[i-j for i in 1:3, j in 1:5],
	[i-j+k for i in 1:3, j in 1:4, k in 1:5],
]

# Write a NPZ file with all the test arrays and numbers,
# and read it back in.
old = Dict{UTF8String, Any}()
for (i, x) in enumerate(TestArrays)
	old["testvar_" * dec(i)] = x
end
filename = "$tmp/big.npz"
npzwrite(filename, old)
new = npzread(filename)
for (k, v) in old
	@test v == new[k]
end
@test old == new

# Write and then read NPY files for each array
for x in TestArrays
	filename = "$tmp/test.npy"
	npzwrite(filename, x)
	y = npzread(filename)
	@test typeof(x) == typeof(y)
	@test size(x) == size(y)
	@test x == y
end

if !Debug
	run(`rm -rf $tmp`)
end


text: include("units.jl")


text:
# Interface methods

function download_raster(T::Type{CHELSA{BioClim}}; layer::Integer=1)
    1 ≤ layer ≤ 19 || throw(ArgumentError("The layer must be between 1 and 19"))
    path = rasterpath(T, layer)
    url = rasterurl(T, layer)
    return _maybe_download(url, path)
end

rastername(::Type{CHELSA{BioClim}}, layer) = "CHELSA_bio10_$(lpad(layer, 2, "0")).tif"

rasterpath(::Type{CHELSA}) = joinpath(rasterpath(), "CHELSA")
rasterpath(::Type{CHELSA{BioClim}}) = joinpath(rasterpath(CHELSA), "BioClim")
rasterpath(T::Type{CHELSA{BioClim}}, layer) = joinpath(rasterpath(T), rastername(T, layer))

rasterurl(::Type{CHELSA}) = "ftp://envidatrepo.wsl.ch/uploads/chelsa/"
rasterurl(::Type{CHELSA{BioClim}}) = joinpath(rasterurl(CHELSA), "chelsa_V1/climatologies/bio/")
rasterurl(T::Type{CHELSA{BioClim}}, layer) = joinpath(rasterurl(T), rastername(T, layer))


text:
### A Pluto.jl notebook ###
# v0.15.1

using Markdown
using InteractiveUtils

# ╔═╡ f0a7e49a-e3e9-11eb-356b-3926060d8c13
md"""
# Chapter 8
#### Nicetak
 $(import Dates; Dates.format(Dates.today(), Dates.DateFormat("U d, Y")))
"""

# ╔═╡ 64d4f95a-6969-4b68-bdfa-e2fc308f2ec4
md"""
## 8.1
連結リストのコード8.6において, 連結リストの各ノードに格納された値を順に出力する関数`printList`の処理に要する計算量を評価してください.
"""

# ╔═╡ ca215632-c7c7-4ffc-a712-78868335644e
md"""
**答** 要素を一つずつ順番に出力するだけなので, $O(N)$
"""

# ╔═╡ 2ae58dfb-0086-44a2-a053-edeb6c8c5f54
md"""
## 8.2
サイズが$N$の連結リストにおいて, `get(i)`を`head`からスタートして$i$番目の要素を取得する関数とします. この時以下のコードの計算量を求めてください.

```
for (int i = 0; i < N; ++i){
	count << get(i) << endl;
}
```
"""

# ╔═╡ 9bb328a2-0fe2-48b9-8d9a-b71822bef6ae
md"""
**答** 連結リストにおいて$i$番目の要素を取得するためには, 最初の要素から順番に$i$回, 続く要素をたどる必要がある. したがって, `get(i)`は$O(N)$で実行され, それを$N$回繰り返すので$O(N^2)$. より正確に評価するなら$O\left(\frac{N(N + 1)}{2}\right)$と考えても良い.
"""

# ╔═╡ ead7ba60-3f49-47d5-bde1-cd4a73d38a8e
md"""
## 8.3
連結リストにおいて, サイズを$O(1)$で取得できるようにする方法を述べてください.
"""

# ╔═╡ 14b8b97b-39f9-417e-bd67-a625b8d05d4f
md"""
**答**
連結リストの実装において, 連結リストのサイズを常に保持するようにする. すなわち,

- 初期値は0
- 挿入(+1)と削除(-1)のたびにサイズを変更

となるように実装しておけば, サイズは$O(1)$で取得できる.
"""

# ╔═╡ df96aa07-7dfa-4d86-b11c-2471b9caa46a
md"""
## 8.4
単方向連結リストにおいて, 特定のノード$v$を削除する方法を述べてください. ただし, $O(N)$の計算量を要してもよいものとします.
"""

# ╔═╡ 73cdc535-703c-47df-a3dc-328ad3cd96a2
md"""
**答** ノード$v$の前の要素を取得すれば, 後のノードと繋ぎかえるだけである(自明に$O(1)$の操作).
ノード$v$の要素の取得は, 先頭のノードから順に次のノードが$v$であるかを判定していけばよい.
この繰り返す回数は$O(N)$であり, 全体の操作も$O(N)$となる.
"""

# ╔═╡ 0b8fe944-4a94-47f3-b0a3-42a35723cc9d
md"""
## 8.5
 $N$個の相異なる整数$a_0, a_1, \dots, a_{N-1}$と, $M$個の相異なる整数$b_0, b_1, \dots, b_{M-1}$が与えられます. $a$と$b$で共通する整数の個数を, 平均的に$O(N + M)$の計算量で求めるアルゴリズムを設計してください.

"""

# ╔═╡ 128f181f-6a33-43b7-b465-55ec0bdce845
md"""
**答**
1. ハッシュテーブルを$a_0, a_1, \dots, a_{N-1}$で構築する, $O(N)$
2. この中に, $b_0, b_1, \dots, b_{M-1}$ が存在するかを順番にチェクする. ハッシュテーブルの要素の検索は$O(1)$ であるから,この操作は $O(M)$

したがって, 全体の操作は$O(N + M)$.
"""

# ╔═╡ 34c3a14e-5e39-4465-b91a-eb403b2d9cda
function f₅(as::Vector{Int}, bs::Vector{Int})
	asₕ = Set(as)
	ans = 0
	for b in bs
		if b ∈ asₕ
			ans += 1
		end
	end

	return ans
end

# ╔═╡ 29fb9255-c233-46d0-a81f-c7b8c53d7ef0
f₅(collect(1:100), collect(1:3:100))

# ╔═╡ fba3be40-eec9-4489-9822-28c2bb2a3216
md"""
## 8.6
 $N$個の整数$a_0, a_1, \dots, a_{N-1}$と, $M$個の整数$b_0, b_1, \dots, b_{M-1}$が与えられます. $a_i = b_j$となるような添字$i, j$の組の個数を, 平均的に$O(N + M)$の計算量で求めるアルゴリズムを設計してください.
"""

# ╔═╡ 1a6f9771-cd95-4568-8a4f-3498d0cd55a1
md"""
**答**

1. 連想配列を用いて, $a_0, a_1, \dots, a_{N-1}$の中の各整数の個数を数える, $O(N)$
1. 各$b_i$に対応する上の連想配列の値を足し上げる, $O(M)$

よって, $O(N + M)$.
"""

# ╔═╡ 7e0eddba-4c9d-4c89-a3c5-abdd3170c232
function f₆(as::Vector{Int}, bs::Vector{Int})
	d = Dict{Int, Int}()
	for a in as
		if haskey(d, a)
			d[a] += 1
		else
			d[a] = 1
		end
	end

	ans = 0
	for b in bs
		if haskey(d, b)
			ans += d[b]
		end
	end

	return ans
end

# ╔═╡ 44358936-b446-420c-9db5-16df16b0de2b
f₆([1, 2, 3, 3], [2, 3])

# ╔═╡ 06d3fa80-a3f4-4c99-acf0-ee2d7ee030ee
md"""
## 8.7
 $N$個の整数$a_0, a_1, \dots, a_{N-1}$と, $N$個の整数$b_0, b_1, \dots, b_{N-1}$が与えられます. 2組の整数列からそれぞれ1個ずつ整数を選んで和を$K$とすることができるかどうかを, 平均的に
$O(N)$の計算量で判定するアルゴリズムを設計してください.　なお, 6.6節では, 類似の問題に対する, 二分探索法に基づいた$O(N \log N)$の計算量のアルゴリズムを示しました.
"""

# ╔═╡ d49be868-3522-452e-9615-9b096c713704
md"""
**答**
1. ハッシュテーブル $b_H$ を $b_0, b_1, \dots, b_{N-1}$ から作る ($O(N)$)
2.  $K - a_i$ が $b_H$の要素にあるかどうかを順番に確認する (一回あたり$O(1)$, 全体は$O(N)$)
3. ある $a_i$に対して2の条件が満たされた時, "Yes", そうでなければ "No" を出力

これは明らかに$O(N)$で動作する
"""

# ╔═╡ 1cd82896-c72b-4ad5-885d-a8a25382faa7
function f₇(as::Vector{Int}, bs::Vector{Int}, K::Int)
	bsₕ = Set(bs)

	ans = "No"
	for a in as
		if (K - a) ∈ bsₕ
			ans = "Yes"
			break
		end
	end

	return ans
end

# ╔═╡ c179229f-de71-4381-a0a6-562f9bab2f38
f₇([1, 2, 4], [8, 16, 32], 9)

# ╔═╡ Cell order:
# ╟─f0a7e49a-e3e9-11eb-356b-3926060d8c13
# ╟─64d4f95a-6969-4b68-bdfa-e2fc308f2ec4
# ╟─ca215632-c7c7-4ffc-a712-78868335644e
# ╟─2ae58dfb-0086-44a2-a053-edeb6c8c5f54
# ╟─9bb328a2-0fe2-48b9-8d9a-b71822bef6ae
# ╟─ead7ba60-3f49-47d5-bde1-cd4a73d38a8e
# ╟─14b8b97b-39f9-417e-bd67-a625b8d05d4f
# ╟─df96aa07-7dfa-4d86-b11c-2471b9caa46a
# ╟─73cdc535-703c-47df-a3dc-328ad3cd96a2
# ╟─0b8fe944-4a94-47f3-b0a3-42a35723cc9d
# ╟─128f181f-6a33-43b7-b465-55ec0bdce845
# ╠═34c3a14e-5e39-4465-b91a-eb403b2d9cda
# ╠═29fb9255-c233-46d0-a81f-c7b8c53d7ef0
# ╟─fba3be40-eec9-4489-9822-28c2bb2a3216
# ╟─1a6f9771-cd95-4568-8a4f-3498d0cd55a1
# ╠═7e0eddba-4c9d-4c89-a3c5-abdd3170c232
# ╠═44358936-b446-420c-9db5-16df16b0de2b
# ╟─06d3fa80-a3f4-4c99-acf0-ee2d7ee030ee
# ╟─d49be868-3522-452e-9615-9b096c713704
# ╠═1cd82896-c72b-4ad5-885d-a8a25382faa7
# ╠═c179229f-de71-4381-a0a6-562f9bab2f38


text:
# ------------------------------------------------------------------------------------------
# # Estructuras de datos
#
# Una vez que empecemos a trabajar con muchos datos a la vez, será conveniente guardar
# nuestros datos en estructuras como arreglos o diccionarios (más allá que sólo
# variables).

#
# Tipos de estructuras de datos que cubrimos
# 1. Tuplas
# 2. Diccionarios
# 3. Arreglos
#
#

# Como repaso, las tuplas y los arreglos ambos son secuencias ordenadas de elementos
# (entonces podemos accesarlos por medio de un índice).
# Los diccionarios y los arreglos son mutables.
#
# ¡Explicaremos más brevemente!
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Tuplas / Tuples
#
# Podemos crear una tupla encerrando una secuencia ordenada de elementos con `( )`.
#
# Sintaxis:

# ```julia
# (item1, item2, ...)```
# ------------------------------------------------------------------------------------------

mis_animales_favoritos = ("pingüino", "gato", "petauro_del_azúcar")

# ------------------------------------------------------------------------------------------
# Podemos accesar con el índice a esta tupla,
# ------------------------------------------------------------------------------------------

myfavoriteanimals[1]

# ------------------------------------------------------------------------------------------
# Pero como las tuplas con inmutables, no la podemos modificar
# ------------------------------------------------------------------------------------------

mis_animales_favoritos[1] = "nutria"

# ------------------------------------------------------------------------------------------
# ## Diccionarios
#
# Si tenemos conjuntos de datos relacionados entre sí, podemos guardar los datos en un
# diccionario. Un buen ejemplo es una lista de contactos, donde asociamos nombres a números
# de teléfono.
#
# Sintaxis:
# ```julia
# Dict(llave1 => valor1, llave2 => valor2, ....)```
# ------------------------------------------------------------------------------------------

miagenda = Dict("Jenny" => "867-5309", "Cazafantasmas" => "555-2368")

# ------------------------------------------------------------------------------------------
# En este ejemplo, cada nombre y número es un par de "llave" y "valor". Podemos tomar el
# númer de Jenny (un valor) usando la llave asociada.
# ------------------------------------------------------------------------------------------

miagenda["Jenny"]

# ------------------------------------------------------------------------------------------
# Podemos agregar otra entrada al diccionario de la manera siguiente
# ------------------------------------------------------------------------------------------

miagenda["Kramer"] = "555-FILK"

# ------------------------------------------------------------------------------------------
# Veamos como se ve nuestro diccionario ahora...
# ------------------------------------------------------------------------------------------

miagenda

# ------------------------------------------------------------------------------------------
# Para borrar a Kramer de nuestro diccionario - y simultáneamente tomar su número - usamos
# pop!
# ------------------------------------------------------------------------------------------

pop!(miagenda, "Kramer")

miagenda

# ------------------------------------------------------------------------------------------
# A diferencia de las tuplas y los arreglos, los diccionarios no están ordenados y no
# podemos accesarlos con un índice
# ------------------------------------------------------------------------------------------

miagenda[1]

# ------------------------------------------------------------------------------------------
# En el ejemplo anterior, `julia` piensa que estás tratando de accesar a un valor asociado a
# la llave `1`.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Arreglos
#
# A diferencia de las tuplas, los arreglos son mutables. A diferencia de los diccionarios,
# los arreglos contienen secuencias ordenadas de elementos.

# Podemos crear un arreglo encapsulando esta secuencia de elementos con `[ ]`.
#
# Sintaxis:

# ```julia
# [item1, item2, ...]```
#
#
# Por ejemplo, podemos usar un arreglo para recordar a mis amigos
# ------------------------------------------------------------------------------------------

myfriends = ["Ted", "Robyn", "Barney", "Lily", "Marshall"]

# ------------------------------------------------------------------------------------------
# O guardar una secuencia de números
# ------------------------------------------------------------------------------------------

fibonacci = [1, 1, 2, 3, 5, 8, 13]

mezcla = [1, 1, 2, 3, "Ted", "Robyn"]

# ------------------------------------------------------------------------------------------
# Una vez que tenenmos un arreglo, podemos tomar los datos individualmente dentro del
# arreglo accesándolos por su índice. Por ejemplo, si queremos al tercer amigo en myfriends,
# escribimos
# ------------------------------------------------------------------------------------------

myfriends[3]

# ------------------------------------------------------------------------------------------
# Podemos usar el índice para mutar la entrada del arreglo
# ------------------------------------------------------------------------------------------

myfriends[3] = "Baby Bop"

# ------------------------------------------------------------------------------------------
# También se puede editar con las fuciones de `push!` y `pop!`. `push!` agrega un elemento
# al final del arreglo y `pop!` quita al último elemento del arreglo.
#
# Podemos agregar otro número a la secuencia fibonacci
# ------------------------------------------------------------------------------------------

push!(fibonacci, 21)

# ------------------------------------------------------------------------------------------
# y quitarlo
# ------------------------------------------------------------------------------------------

pop!(fibonacci)

fibonacci

# ------------------------------------------------------------------------------------------
# Hasta ahora hemos dado ejemplos de arreglos de escalars unidimensionales, pero los
# arreglos pueden tener un número arbitrario de dimensiones y también pueden guardar otros
# arreglos.
#

# Por ejemplo, estos son arreglos de arreglos
# ------------------------------------------------------------------------------------------

favoritos = [["koobideh", "chocolate", "eggs"],["penguins", "cats", "sugargliders"]]

numbers = [[1, 2, 3], [4, 5], [6, 7, 8, 9]]

# ------------------------------------------------------------------------------------------
# Abajo hay arreglos de 2D y 3D poblados con valores aleatorios
# ------------------------------------------------------------------------------------------

rand(4, 3)

rand(4, 3, 2)

# ------------------------------------------------------------------------------------------
# ¡Cuidado copiando los arreglos!
# ------------------------------------------------------------------------------------------

fibonacci

somenumbers = fibonacci

somenumbers[1] = 404

fibonacci

# ------------------------------------------------------------------------------------------
# Editar `somenumbers` causa que `fibonacci` se edite también!
#
# En el ejemplo superior, en realidad no hicimos una copia de `fibonacci`. Sólo creamos una
# nueva manera de accesar las entradas del arreglo relacionado a `fibonacci`.
#
# Si queremos hacer una copia de un arreglo amarrado a `fibonacci`, usamos la función de
# `copy`.
# ------------------------------------------------------------------------------------------

# Primero restauramos a fibonnaci
fibonacci[1] = 1
fibonacci

somemorenumbers = copy(fibonacci)

somemorenumbers[1] = 404

fibonacci

# ------------------------------------------------------------------------------------------
# En el último ejemplo, no se editó a fibonacci. Entonces vemos que los arreglos amarrados a
# `somemorenumbers` y `fibonacci` son distintos.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### Exercises
#
# 3.1 Crea un arreglo, `arreglo`, que sea un arreglo de 1D de 2-elementos de 1-elemento de
# 1D, cada uno guardando el número 0.
# Accesa a `arreglo` para agregar un `1` a cada uno de los arreglos que contiene.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# 3.2 Trata de agregar "Emergencia" a `miagenda` con el valor `911`. Trata de agregar `911`
# como un entero y no como cadena. ¿Porqué no funciona?
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# 3.3 Crea un nuevo diccionario que se llame `agenda_flexible` que tenga el número de Jenny
# guardado como cadena y el de los Cazafantasmas como entero.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# 3.4 Agrega la llave de "Emergencia al valor (entero) `911` a `agenda_flexible`.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# 3.5 ¿Porqué podemos agregar un entero como valor a `agenda_flexible` pero no a
# `miagenda`? ¿Cómo pudimos haber inicializado `miagenda` para que aceptara enteros como
# valores?
# ------------------------------------------------------------------------------------------

max_stars_repo_path: .nbexports/introductory-tutorials/intro-to-julia-ES/ 3. Estructuras de datos.jl

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
module TestADT
using Test
using BitTwiddle

using BitTwiddle: ADTPlan, ConstructorPlan, emit, ctorid
import BitTwiddle: ctorid, argtypes

plan = ADTPlan(:MaybeInt, [
        ConstructorPlan(:Just, Tuple{Int}),
        ConstructorPlan(:Nil, Tuple{}),
        ])

println(@macroexpand @adt MaybeInt begin
    Just(::Int)
    Nil()
end)

@adt MaybeInt begin
    Just(::Int)
    Nil()
end

one = @inferred Just(1)
two = @inferred Just(2)
nil = @inferred Nil()
@test one isa MaybeInt
@test nil isa MaybeInt
@test one == one
@test one != two
@test nil == nil
@test one != nil
@test typeof(one) == typeof(nil)

out = @inferred destructure(Just, one)
@test out === (1,)
@test destructure(Just, two) === (2,)
@test () === destructure(Nil, nil)
@test_throws ArgumentError destructure(Nil, one)
@test_throws ArgumentError destructure(Just, nil)

@test occursin("Just", sprint(show, one))
@test occursin("Nil", sprint(show, nil))
for o in [nil, one, two]
    @test eval(parse(sprint(show, o))) === o
end

end#module


text:
using LinearAlgebra

################################################################################
# Evaluation functions for (unconstrained) ScSTOptiModel
# obj, grad!, prox!, objprox!
################################################################################
"Evaluate smooth objective"
function OptiMo.obj(p::ScSTOptiModel, x::Vector{Float64})
    chkfinite(x)
    precompMatJ!(p, x)
    fx = (p.data.x0' * p.eval.S[:, :, 1] * p.data.x0)[1]
    p.repo.nobjf += 1
    p.repo.objf = [p.repo.objf; fx]
    p.repo.delta = [p.repo.delta x]
    chkfinite(fx)
    return fx
end

"Evaluate gradient"
function OptiMo.grad!(p::ScSTOptiModel, x::Vector{Float64}, dfx::Vector{Float64})
    chkfinite(x)
    precompMatGradJ!(p, x)
    for i = 1:p.data.N
        dfx[i] = (p.eval.X[:, p.eval.tauIdx[i+1]]'*p.eval.C[:, :, i]*p.eval.X[:, p.eval.tauIdx[i+1]])[1]
    end
    p.repo.ngrad += 1
    chkfinite(dfx)
    return nothing
end

"Evaluate smooth objective and gradient"
function OptiMo.objgrad!(p::ScSTOptiModel, x::Vector{Float64}, dfx::Vector{Float64})
    chkfinite(x)
    precompMatGradJ!(p, x)
    fx = (p.data.x0'*p.eval.S[:, :, 1]*p.data.x0)[1]
    for i = 1:p.data.N
        dfx[i] = (p.eval.X[:, p.eval.tauIdx[i+1]]'*p.eval.C[:, :, i]*p.eval.X[:, p.eval.tauIdx[i+1] ])[1]
    end
    p.repo.nobjf += 1
    p.repo.ngrad += 1
    p.repo.objf = [p.repo.objf; fx]
    p.repo.delta = [p.repo.delta x]
    chkfinite(fx)
    chkfinite(dfx)
    return fx
end

"Constraints"
function OptiMo.cons!(p::ScSTOptiModel, x::Vector{Float64}, cx::Vector{Float64})
    chkfinite(x)
    tau, _ = gettau(p, x)
    p.data.constr(tau, cx)
    p.repo.ncons += 1
    chkfinite(cx)
    return nothing
end

"Transposed-Jacobian-vector multiplication"
function OptiMo.jtprod!(p::ScSTOptiModel, x::Vector{Float64}, v::Vector{Float64}, jtv::Vector{Float64})
    chkfinite(x)
    chkfinite(v)
    tau, _ = gettau(p, x)
    # x, jtv [N]
    # tau, jttv [N-1]
    # v [ncon]
    p.data.dconstr(tau, v, p.eval.jttv)
    jtv[p.data.N] = 0.0
    for i in p.data.N-1:-1:1
        jtv[i] = jtv[i+1] + p.eval.jttv[i]
    end
    p.repo.ncjtv += 1
    chkfinite(jtv)
    return nothing
end

"Project constraints"
function OptiMo.proj!(p::ScSTOptiModel, cx::Vector{Float64}, px::Vector{Float64})
    chkfinite(cx)
    p.data.pconstr(cx, px)
    p.repo.nproj += 1
    chkfinite(px)
    return nothing
end

"Prox operator"
function OptiMo.prox!(p::ScSTOptiModel, x::Vector{Float64}, a::Float64, z::Vector{Float64})
    chkfinite(x)
    chkfinite(a)
    proxl0simplex!(x, a * p.data.swc, p.data.dt, z)
    p.repo.nprox += 1
    chkfinite(z)
    return nothing
end

"Prox operator and objective"
function OptiMo.objprox!(p::ScSTOptiModel, x::Vector{Float64}, a::Float64, z::Vector{Float64})
    chkfinite(x)
    chkfinite(a)
    gz = proxl0simplex!(x, a * p.data.swc, p.data.dt, z)
    gz *= p.data.swc
    p.repo.nobjg += 1
    p.repo.nprox += 1
    chkfinite(z)
    chkfinite(gz)
    return gz
end

################################################################################
# Functions to precompute matrices
################################################################################
"Precompute matrices for cost function"
function precompMatJ!(p::ScSTOptiModel, x::Vector{Float64})
    chkfinite(x)
    if p.eval.delta_old_S != x
        copyto!(p.eval.delta_old_S, x)
        propagateDynamics!(p, x)
        computeS!(p.eval, p.data.E, p.data.nvec)
    end
    return nothing
end

"Precompute matrices for gradient"
function precompMatGradJ!(p::ScSTOptiModel, x::Vector{Float64})
    chkfinite(x)
    if p.eval.delta_old_C != x
        copyto!(p.eval.delta_old_C, x)
        precompMatJ!(p, x)
        computeC!(p.eval, p.data.Q, p.data.N)
    end
    return nothing
end

"Propagate dynamics"
function propagateDynamics!(p::ScSTOptiModel, x::Vector{Float64})
    chkfinite(x)
    p.repo.ndyna += 1
    # Get positive delta
    x = max.(x, 0.0)
    # Get switching times from delta
    tau, p.eval.tfd = delta2tau(x, p.data.t0)
    # Create grid from t0 to tfd
    #p.eval.tgrid[end] = p.eval.tfd
    #p.eval.tgrid = min.(p.eval.tgrid, p.eval.tgrid[end])
    p.eval.tgrid .= collect(range(p.data.t0, stop = p.eval.tfd, length=p.data.ngrid))
    chkfinite(p.eval.tgrid)
    # Fit switching times withing the grid
    p.eval.tvec, p.eval.tauIdx, p.eval.tgridIdx = mergeSortFindIndex(p.eval.tgrid, tau)
    p.eval.dvec = diff(p.eval.tvec)
    chkfinite(p.eval.dvec)
    # index for current input
    uIdx = 1
    # Compute Matrix Exponentials
    n = p.data.nx
    for i = 1:p.data.nvec-1
        # Verify which U input applies
        if uIdx < p.data.N && i >= p.eval.tauIdx[uIdx+1]
            uIdx += 1
        end
        # Linearize Dynamics
        chkfinite(p.eval.X[1:end-1, i])
        p.eval.A[:, :, i] = linearizeDyn(
            p.data.dynam,
            p.data.ddynam,
            p.eval.X[1:end-1, i],
            p.data.U[:, uIdx],
        )
        chkfinite(p.eval.A[:, :, i])
        # Compute Matrix Exponential
        tm = Matrix( [ - p.eval.A[:, :, i]' p.data.Q
                       zeros(n, n) p.eval.A[:, :, i] ] * p.eval.dvec[i] )
        chkfinite(tm)
        tm .= exp(tm)
        chkfinite(tm)
        # Assign \mathcal{E}_i
        p.eval.exm[:, :, i] .= tm[n+1:2*n, n+1:2*n]
        # Assign M_k
        p.eval.M[:, :, i] .= tm[n+1:2*n, n+1:2*n]' * tm[1:n, n+1:2*n]
        # Compute state at the next grid point
        p.eval.X[:, i+1] .= p.eval.exm[:, :, i] * p.eval.X[:, i]
    end
    return nothing
end

################################################################################
# Lower level functions
################################################################################
"Low level function for computing S matrices"
function computeS!(e::ScSTOeval, E::Matrix{Float64}, nvec::Int)
    e.S[:, :, nvec] .= E
    for i = nvec-1:-1:1
        e.S[:, :, i] .=
            e.M[:, :, i] .+ e.exm[:, :, i]' * e.S[:, :, i+1] * e.exm[:, :, i]
    end
end

"Low level function for computing C matrices"
function computeC!(e::ScSTOeval, Q::Matrix{Float64}, N::Int)
    for i = 1:N
        e.C[:, :, i] .=
            Q .+
            e.A[:, :, e.tauIdx[i+1]-1]' * e.S[:, :, e.tauIdx[i+1]] .+
            e.S[:, :, e.tauIdx[i+1]] * e.A[:, :, e.tauIdx[i+1]-1]
    end
end

"""
Linearize Dynamics
    A = linearizeDyn(dynam, ddynam, x, u)
"""
function linearizeDyn(dynam::Function, ddynam::Function, x::Vector{Float64}, u::Vector{Float64})
    chkfinite(x)
    chkfinite(u)
    f = dynam(x, u)
    df = ddynam(x, u)
    A = [df f - df * x; zeros(1, length(x) + 1)]
    chkfinite(A)
    return A
end

"""
Merge, sort, and find indices of grid and switching times
    tvec, tauIdx, tgridIdx = mergeSortFindIndex(tgrid, tau)
"""
function mergeSortFindIndex(tgrid::Vector{Float64}, tau::Vector{Float64})
    chkfinite(tgrid)
    chkfinite(tau)
    ngrid = length(tgrid)
    ntau = length(tau)
    tauIdx = Array{Int}(undef, ntau + 2)
    tauIdx[1] = 1
    tauIdx[end] = ntau + ngrid
    tgridIdx = Array{Int}(undef, ngrid)
    tgridIdx[1] = 1
    tgridIdx[end] = ntau + ngrid
    # Create merged and sorted time vector with grid and switching times
    ttmp = vcat(tgrid, tau)     # Concatenate grid and tau vector
    tidxtmp = sortperm(ttmp)    # Find permutation vector to sort ttemp
    tvec = ttmp[tidxtmp]        # Create full sorted tvec
    # Create index of the tau vector elements inside tvec
    for i = 1:ntau
        tauIdx[i+1] = findfirst(x -> x == ngrid + i, tidxtmp)
    end
    # Create index of the tgrid vector elements inside tvec
    for i = 1:ngrid
        tgridIdx[i] = findfirst(x -> x == i, tidxtmp)
    end
    chkfinite(tvec)
    chkfinite(tauIdx)
    chkfinite(tgridIdx)
    return tvec, tauIdx, tgridIdx
end

"""
    delta = tau2delta(tau, t0, tf)
converts from switching times to switching intervals.
"""
function tau2delta(tau::Vector{Float64}, t0::Float64, tf::Float64)
    chkfinite(tau)
    delta = diff([t0; tau; tf])
    chkfinite(delta)
    return delta
end

"""
    tau, tfd = delta2tau(delta, t0)
converts from switching intervals to switching times.
"""
function delta2tau(delta::Vector{Float64}, t0::Float64)
    chkfinite(delta)
    tau = Vector{Float64}(undef, length(delta) - 1)
    tau[1] = t0 + delta[1]
    for i = 2:length(tau)
        tau[i] = tau[i-1] + delta[i]
    end
    tfd = tau[end] + delta[end]
    chkfinite(tau)
    chkfinite(tfd)
    return tau, tfd
end

function chkfinite(x::R) where {R <: Real}
    if !isfinite(x)
        throw(ArgumentError("scalar contains Infs or NaNs"))
    end
    return true
end

function chkfinite(x::Vector{R}) where {R <: Real}
    for e in x
        if !isfinite(e)
            throw(ArgumentError("vector contains Infs or NaNs"))
        end
    end
    return true
end

function chkfinite(x::Matrix{R}) where {R <: Real}
    for e in x
        if !isfinite(e)
            throw(ArgumentError("matrix contains Infs or NaNs"))
        end
    end
    return true
end

function chkfinite(x::Array{R}) where {R <: Real}
    for e in x
        if !isfinite(e)
            throw(ArgumentError("array contains Infs or NaNs"))
        end
    end
    return true
end


text:
#計算の前提条件読込み
"""
    readprerequisite(ExcelFilePath::String,parameters::String,schedule::String;schedulcolumn::Int=6)

# 処理概要
スケジュール作成

# 引数
* `ExcelFilePath`:スケジュール作成に必要な情報が記載されたExcelファイルのフルパスを指定。
* `parameters`:パラメーターが記載されたシート名を指定
* `schedule`:ジョブ作成に必要な情報が記載されたシート名を指定
* `schedulcolumn`:スケジュールON,OFF開始列を指定（デフォルト：６）

# 結果（戻り値）
パラメータの読込み結果を出力
"""
function readprerequisite(ExcelFilePath::String,parameters::String,schedule::String;schedulcolumn::Int=6)
  parameters=DataFrames.DataFrame(XLSX.readtable(ExcelFilePath, parameters)...)
  scheduleplan=DataFrames.DataFrame(XLSX.readtable(ExcelFilePath, schedule)...)

  # 前提条件設定
  robotn=parameters[parameters[:,:parameter] .== "all_run_robot",:value][1] #ロボット
  run_unit_time=parameters[parameters[:,:parameter] .== "run_unit_time",:value][1] #実行単位時間、コマ割時間
  jobn=size(scheduleplan)[1] #ジョブ数
  timen=size(scheduleplan)[2]-schedulcolumn #時間コマ数

  blocktime_start=string(parameters[parameters[:,:parameter] .== "blocktime_start",:value][1])[1:5]  #スケジュール割当禁止開始時間
  blocktime_end=string(parameters[parameters[:,:parameter] .== "blocktime_end",:value][1])[1:5]  #スケジュール割当禁止終了時間
  blocktime_dow=parameters[parameters[:,:parameter] .== "blocktime_dow",:value][1]  #スケジュール割当禁止を設定するシート名

  return scheduleplan,robotn,run_unit_time,jobn,timen,schedule,blocktime_start,blocktime_end,blocktime_dow
end


text:
export iauStarpv
"""
Convert star catalog coordinates to position+velocity vector.

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

Status:  support function.

Given (Note 1):
   ra     double        right ascension (radians)
   dec    double        declination (radians)
   pmr    double        RA proper motion (radians/year)
   pmd    double        Dec proper motion (radians/year)
   px     double        parallax (arcseconds)
   rv     double        radial velocity (km/s, positive = receding)

Returned (Note 2):
   pv     double[2][3]  pv-vector (au, au/day)

Returned (function value):
         int           status:
                           0 = no warnings
                           1 = distance overridden (Note 6)
                           2 = excessive speed (Note 7)
                           4 = solution didn't converge (Note 8)
                        else = binary logical OR of the above

Notes:

   1. The star data accepted by this function are "observables" for an
   imaginary observer at the solar-system barycenter.  Proper motion
   and radial velocity are, strictly, in terms of barycentric
   coordinate time, TCB.  For most practical applications, it is
   permissible to neglect the distinction between TCB and ordinary
   "proper" time on Earth (TT/TAI).  The result will, as a rule, be
   limited by the intrinsic accuracy of the proper-motion and
   radial-velocity data;  moreover, the pv-vector is likely to be
   merely an intermediate result, so that a change of time unit
   would cancel out overall.

   In accordance with normal star-catalog conventions, the object's
   right ascension and declination are freed from the effects of
   secular aberration.  The frame, which is aligned to the catalog
   equator and equinox, is Lorentzian and centered on the SSB.

   2. The resulting position and velocity pv-vector is with respect to
   the same frame and, like the catalog coordinates, is freed from
   the effects of secular aberration.  Should the "coordinate
   direction", where the object was located at the catalog epoch, be
   required, it may be obtained by calculating the magnitude of the
   position vector pv[0][0-2] dividing by the speed of light in
   au/day to give the light-time, and then multiplying the space
   velocity pv[1][0-2] by this light-time and adding the result to
   pv[0][0-2].

   Summarizing, the pv-vector returned is for most stars almost
   identical to the result of applying the standard geometrical
   "space motion" transformation.  The differences, which are the
   subject of the Stumpff paper referenced below, are:

   (i) In stars with significant radial velocity and proper motion,
   the constantly changing light-time distorts the apparent proper
   motion.  Note that this is a classical, not a relativistic,
   effect.

   (ii) The transformation complies with special relativity.

   3. Care is needed with units.  The star coordinates are in radians
   and the proper motions in radians per Julian year, but the
   parallax is in arcseconds; the radial velocity is in km/s, but
   the pv-vector result is in au and au/day.

   4. The RA proper motion is in terms of coordinate angle, not true
   angle.  If the catalog uses arcseconds for both RA and Dec proper
   motions, the RA proper motion will need to be divided by cos(Dec)
   before use.

   5. Straight-line motion at constant speed, in the inertial frame,
   is assumed.

   6. An extremely small (or zero or negative) parallax is interpreted
   to mean that the object is on the "celestial sphere", the radius
   of which is an arbitrary (large) value (see the constant PXMIN).
   When the distance is overridden in this way, the status,
   initially zero, has 1 added to it.

   7. If the space velocity is a significant fraction of c (see the
   constant VMAX), it is arbitrarily set to zero.  When this action
   occurs, 2 is added to the status.

   8. The relativistic adjustment involves an iterative calculation.
   If the process fails to converge within a set number (IMAX) of
   iterations, 4 is added to the status.

   9. The inverse transformation is performed by the function
   iauPvstar.

Called:
   iauS2pv      spherical coordinates to pv-vector
   iauPm        modulus of p-vector
   iauZp        zero p-vector
   iauPn        decompose p-vector into modulus and direction
   iauPdp       scalar product of two p-vectors
   iauSxp       multiply p-vector by scalar
   iauPmp       p-vector minus p-vector
   iauPpp       p-vector plus p-vector

Reference:

   Stumpff, P., 1985, Astron.Astrophys. 144, 232-240.

This revision:  2017 March 16

SOFA release 2018-01-30

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

# int iauStarpv(double ra, double dec,
#               double pmr, double pmd, double px, double rv,
#               double pv[2][3])

function iauStarpv(ra::Real, dec::Real, pmr::Real, pmd::Real,
                   px::Real, rv::Real)

   pv = zeros(Float64, 3, 2)

   status = ccall((:iauStarpv, libsofa_c), Cint,
       (Cdouble, Cdouble, Cdouble, Cdouble,
       Cdouble, Cdouble, Ptr{Cdouble}),
       convert(Float64, ra), convert(Float64, dec),
       convert(Float64, pmr), convert(Float64, pmd),
       convert(Float64, px), convert(Float64, rv),
       pv)

   return status, SMatrix{2,3}(pv')
end