text:
################################################################################
#
#             EllCrv/QQ.jl : Minimal models and local information
#
# This file is part of Hecke.
#
# Copyright (c) 2015, 2016: Claus Fieker, Tommy Hofmann
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# * Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
# (C) 2016 Tommy Hofmann
# (C) 2016 Robin Ammon
# (C) 2016 Sofia Brenner
# (C) 2022 Jeroen Hanselman
#
################################################################################

export laska_kraus_connell, minimal_model, tates_algorithm_global, tates_algorithm_local, tidy_model,
       torsion_points, torsion_structure, torsion_bound, tamagawa_number, tamagawa_numbers,
       kodaira_symbol, kodaira_symbols, reduction_type, modp_reduction

################################################################################
#
#  Algorithm of Laska-Kraus-Connel
#
################################################################################

# algorithm of Laska-Kraus-Connell
@doc Markdown.doc"""
    laska_kraus_connell(E::EllCrv{fmpz}) -> Array{Nemo.fmpz}

Given an elliptic curve over $\mathbf Q$ with integral model, this returns an
isomorphic elliptic curve over $\mathbf Q$ with minimal discriminant.
"""
function laska_kraus_connell(E::EllCrv{fmpq})
  a1, a2, a3, a4, a6 = map(numerator,(a_invars(E)))

  b2, b4, b6, b8, c4, c6 = get_b_c_integral(E)

  delta = divexact(c4^3 - c6^2, 1728)

  u = fmpz(1)
  g = gcd(c6^2, delta)

  fac = factor(g)

  for (p, ord) in fac
    d = div(ord, 12)

    if p == 2
      a = divexact(c4, 2^(4*d))
      a = mod(a, 16)
      if a > 8
        a = a - 16
      end

      b = divexact(c6, 2^(6*d))
      b = mod(b, 32)
      if b > 16
        b = b - 32
      end

      if (mod(b, 4) != 3) && !((a == 0) && ((b == 0) || (b == 8)))
        d = d - 1
      end

    elseif p == 3
      ord1 = valuation(c6, 3)

      if (ord1 == 6*d + 2)
        d = d - 1
      end
    end
    u = u * p^d
  end

  c4 = divexact(c4, u^4)
  c6 = divexact(c6, u^6)

  b2 = mod(-c6, 12)
  if b2 > 6
      b2 = b2 - 12
  end

  b4 = divexact(b2^2 - c4, 24)
  b6 = divexact(-b2^3 + 36*b2*b4 - c6, 216)

  na1 = mod(b2, 2)
  na2 = divexact(b2 - na1, 4)
  na3 = mod(b6, 2)
  na4 = divexact(b4 - na1*na3, 2)
  na6 = divexact(b6 - na3, 4)

  return EllipticCurve([na1, na2, na3, na4, na6])::EllCrv{fmpq}
end

################################################################################
#
#  Tates algorithm
#
################################################################################

# Tate's algorithm over number fields, see Cremona, p. 66, Silverman p. 366
@doc Markdown.doc"""
    tates_algorithm_local(E::EllCrv{nf_elem}, pIdeal:: NfOrdIdl)
    -> EllipticCurve{nf_elem}, String, fmpz, fmpz, Bool

Returns a tuple $(\tilde E, K, m, f, c, s)$, where $\tilde E$ is a
minimal model for $E$ at the prime ideal $p$, $K$ is the Kodaira symbol,
$f$ is the conductor valuation at $p$, $c$ is the local Tamagawa number
at $p$ and s is false if and only if $E$ has non-split
multiplicative reduction.
"""
function tates_algorithm_local(E::EllCrv{nf_elem},pIdeal:: NfOrdIdl)

  #Assumption: Coefficients of E are in DVR with maximal ideal pIdeal of K.

  #Check if we have generators
  p = pIdeal.gen_one
  uniformizer = pIdeal.gen_two
  R = ring_of_integers(base_field(E))

  a1, a2, a3, a4, a6 = map(R,(a_invars(E)))

  b2, b4, b6, b8 = map(R,(b_invars(E)))
  c4, c6 = map(R,(c_invars(E)))

  delta = discriminant(E)
  delta = numerator(delta)

  n = valuation(delta, pIdeal)

  # test for type I0
  if n == 0
    return (E, "I0", FlintZZ(0), FlintZZ(1), true)::Tuple{EllCrv{nf_elem}, String, fmpz, fmpz, Bool}
  end

  # Maybe smods?
  # change coordinates so that p | a3, a4, a6
  if p == 2
    if mod(b2, pIdeal) == 0
      r = pth_root_mod(a4, pIdeal)
      t = mod(r*(r*(r + a2) + a4) + a6, pIdeal)
    else
      temp = invmod(a1, pIdeal)
      r = temp*a3
      t = temp*(a4 + r^2)
    end

  elseif p == 3
    if mod(b2, pIdeal) == 0
      r = pth_root_mod(-b6, pIdeal)
    else
      r = -invmod(b2, pIdeal)*b4
    end

    t = a1*r + a3
  else
    if mod(c4, pIdeal) == 0
      r = - invmod(FlintZZ(12), p)*b2
    else
      r = - invmod(FlintZZ(12)*c4, pIdeal)*(c6 + b2 * c4)
    end
      t = - invmod(FlintZZ(2), p)* (a1*r + a3)
      r = mod(R(r), pIdeal)
      t = mod(R(t), pIdeal)
  end

  trans = transform_rstu(E, [r, 0, t, 1])
  E = trans[1]

  a1, a2, a3, a4, a6 = map(R,(a_invars(E)))
  b2, b4, b6, b8 = map(R,(b_invars(E)))
  c4, c6 = map(R,(c_invars(E)))

  split = true
  # test for types In, II, III, IV
  if mod(c4, pIdeal) != 0
    if quadroots(1, a1, -a2, pIdeal)
      cp = FlintZZ(n)
    elseif mod(n, 2) == 0
      cp = FlintZZ(2)
      split = false
    else
      cp = FlintZZ(1)
      split = false
    end

    Kp = "I$(n)"
    fp = FlintZZ(1)

    return (E, Kp, fp, cp, split)::Tuple{EllCrv{nf_elem}, String, fmpz, fmpz, Bool}
  end

  if mod(a6, pIdeal^2) != 0
    Kp = "II"
    fp = FlintZZ(n)
    cp = FlintZZ(1)
    return (E, Kp, fp, cp, true)::Tuple{EllCrv{nf_elem}, String, fmpz, fmpz, Bool}
  end

  if mod(b8, pIdeal^3) != 0
    Kp = "III"
    fp = FlintZZ(n-1)
    cp = FlintZZ(2)
    return (E, Kp, fp, cp, true)::Tuple{EllCrv{nf_elem}, String, fmpz, fmpz, Bool}
  end

  if mod(b6, pIdeal^3) != 0
    if quadroots(1, divexact(a3, uniformizer), divexact(-a6, uniformizer^2), pIdeal)
      cp = FlintZZ(3)
    else
      cp = FlintZZ(1)
    end
    Kp = "IV"
    fp = FlintZZ(n - 2)
    return (E, Kp, fp, cp, true)::Tuple{EllCrv{nf_elem}, String, fmpz, fmpz, Bool}
  end

  # change coordinates so that p | a1, a2; p^2 | a3, a4; p^3 | a6
  # Taking square roots ok?
  if p == 2
    s = pth_root_mod(a2, pIdeal)
    t = uniformizer * pth_root(divexact(a6, uniformizer^2), pIdeal)
  elseif p ==3
    s = a1
    t = a3
  else
    s = -a1 * invmod(FlintZZ(2), p)
    t = -a3 * invmod(FlintZZ(2), p)
  end

  trans = transform_rstu(E, [0, s, t, 1])
  E = trans[1]

  a1, a2, a3, a4, a6 = map(R,(a_invars(E)))

  # set up auxililary cubic T^3 + bT^2 + cT + d
  b = divexact(a2, uniformizer)
  c = divexact(a4, uniformizer^2)
  d = divexact(a6, uniformizer^3)
  w = 27*d^2 - b^2*c^2 + 4*b^3*d - 18*b*c*d + 4*c^3
  x = 3*c - b^2
  # test for distinct roots: type I0*
  if mod(w, pIdeal) != 0
    Kp = "I0*"
    fp = FlintZZ(n - 4)
    cp = 1 + nrootscubic(b, c, d, pIdeal)
    return (E, Kp, fp, cp, true)::Tuple{EllCrv{nf_elem}, String, fmpz, fmpz, Bool}

  # test for double root: type Im*
  elseif mod(x, pIdeal) != 0
    # change coordinates so that the double root is T = 0
    if p == 2
      r = pth_root_mod(c, pIdeal)
    elseif p == 3
      r = c*inv_mod(b, pIdeal)
    else
      r = (b*c - 9*d) * invmod(FlintZZ(2)*x, pIdeal)
    end

    r = uniformizer * mod(r, pIdeal)

    trans = transform_rstu(E, [r, 0, 0, 1])
    E = trans[1]

    a1, a2, a3, a4, a6 = map(numerator,(a_invars(E)))

    # make a3, a4, a6 repeatedly more divisible by p
    m = 1
    mx = uniformizer^2
    my = uniformizer^2
    cp = FlintZZ(0)
    while cp == 0
      xa2 = R(divexact(a2, uniformizer))
      xa3 = R(divexact(a3, my))
      xa4 = R(divexact(a4, uniformizer*mx))
      xa6 = R(divexact(a6, mx*my))
      if mod(xa3^2 + 4*xa6, pIdeal) !=  0
        if quadroots(1, xa3, -xa6, pIdeal)
          cp = FlintZZ(4)
        else
          cp = FlintZZ(2)
        end

      else
        if p == 2
          t = my * pth_root_mod(xa6, pIdeal)
        else
          t = my * mod(R(-xa3*invmod(FlintZZ(2), p)), pIdeal)
        end

        trans = transform_rstu(E, [0, 0, t, 1])
        E = trans[1]

        a1, a2, a3, a4, a6 = map(numerator,(a_invars(E)))

        my = my*uniformizer
        m = m + 1
        xa2 = R(divexact(a2, uniformizer))
        xa3 = R(divexact(a3, my))
        xa4 = R(divexact(a4, uniformizer*mx))
        xa6 = R(divexact(a6, mx*my))

        if mod(xa4^2 - 4*xa2*xa6, pIdeal) != 0
          if quadroots(xa2, xa4, xa6, pIdeal)
            cp = FlintZZ(4)
          else
            cp = FlintZZ(2)
          end

        else
          if p == 2
            r = mx * pth_root_mod(xa6*inv_mod(xa2), pIdeal)
          else
            r = mx * mod(-xa4 * invmod(2*xa2, pIdeal), pIdeal)
          end

          trans = transform_rstu(E, [r, 0, 0, 1])
          E = trans[1]

          a1, a2, a3, a4, a6 = map(numerator,(a_invars(E)))

          mx = mx*uniformizer
          m = m + 1
        end
      end
    end

    fp = n - m - 4
    Kp = "I$(m)*"

    return (E, Kp, FlintZZ(fp), FlintZZ(cp), true)::Tuple{EllCrv{nf_elem}, String, fmpz, fmpz, Bool}

  else
    # Triple root case: types II*, III*, IV* or non-minimal
    # change coordinates so that the triple root is T = 0
    if p == 2
      r = b
    elseif p == 3
      r = pth_root_mod(-d, 3)
    else
      r = -b * invmod(FlintZZ(3), p)
    end

    r = p * mod(r, pIdeal)

    trans = transform_rstu(E, [r, 0, 0, 1])
    E = trans[1]

    a1, a2, a3, a4, a6 = map(numerator,(a_invars(E)))

    x3 = divexact(a3, uniformizer^2)
    x6 = divexact(a6, uniformizer^4)

    # Test for type IV*
    if mod(x3^2 + 4* x6, pIdeal) != 0
      if quadroots(1, x3, -x6, pIdeal)
        cp = FlintZZ(3)
      else
        cp = FlintZZ(1)
      end
      Kp = "IV*"
      fp = FlintZZ(n - 6)

      return (E, Kp, fp, FlintZZ(cp), true)::Tuple{EllCrv{nf_elem}, String, fmpz, fmpz, bool}
    else
      if p == 2
        t = -uniformizer^2 * pth_root_mod(x6, pIdeal)
      else
        t = uniformizer^2 * mod(x3 * invmod(FlintZZ(2), p), pIdeal)
      end

      trans = transform_rstu(E, [0, 0, t, 1])
      E = trans[1]

      a1, a2, a3, a4, a6 = map(numerator,(a_invars(E)))

      # Test for types III*, II*
      if mod(a4, uniformizer^4) != 0
        Kp = "III*"
        fp = FlintZZ(n - 7)
        cp = FlintZZ(2)

        return (E, Kp, fp, FlintZZ(cp), true)::Tuple{EllCrv{nf_elem}, String, fmpz, fmpz, Bool}
      elseif mod(a6, uniformizer^6) != 0
        Kp = "II*"
        fp = FlintZZ(n - 8)
        cp = FlintZZ(1)

        return (E, Kp, fp, FlintZZ(cp))::Tuple{EllCrv{nf_elem}, String,  fmpz, fmpz, Bool}
      else
        E = transform_rstu(E, [0, 0, 0, uniformizer])[1]
        return tates_algorithm_local(E, pIdeal)::Tuple{EllCrv{nf_elem}, String, fmpz, fmpz, Bool}
      end
    end
  end
end

@doc Markdown.doc"""
    tates_algorithm_local(E::EllCrv{fmpq}, p:: Int)
    -> EllipticCurve{fmpq}, String, fmpz, fmpz, Bool

Returns a tuple $(\tilde E, K, f, c, s)$, where $\tilde E$ is a
minimal model for $E$ at the prime ideal $p$, $K$ is the Kodaira symbol,
$f$ is the conductor valuation at $p$, $c$ is the local Tamagawa number
at $p$ and s is false if and only if $E$ has non-split
multiplicative reduction.
"""
function tates_algorithm_local(E::EllCrv{fmpq}, p)

  p = FlintZZ(p)

  a1, a2, a3, a4, a6 = map(numerator,(a_invars(E)))

  b2, b4, b6, b8, c4, c6 = get_b_c_integral(E)

  delta = discriminant(E)
  delta = numerator(delta)

  n = valuation(delta, p)

  # test for type I0
  if n == 0
    return (E, "I0", FlintZZ(0), FlintZZ(1), true)::Tuple{EllCrv{fmpq}, String, fmpz, fmpz, Bool}
  end

  # change coordinates so that p | a3, a4, a6
  if p == 2
    if mod(b2, p) == 0
      r = smod(a4, p)
      t = smod(r*(1 + a2 + a4) + a6, p)
    else
      r = smod(a3, p)
      t = smod(r + a4, p)
    end

  elseif p == 3
    if mod(b2, p) == 0
      r = smod(-b6, p)
    else
      r = smod(-b2*b4, p)
    end

    t = smod(a1*r + a3, p)
  else
    if mod(c4, p) == 0
      r = - invmod(FlintZZ(12), p)*b2
    else
      r = - invmod(FlintZZ(12)*c4, p)*(c6 + b2*c4)
    end
      t = - invmod(FlintZZ(2), p)* (a1*r + a3)
      r = smod(r, p)
      t = smod(t, p)
  end

  trans = transform_rstu(E, [r, 0, t, 1])
  E = trans[1]

  a1, a2, a3, a4, a6 = map(numerator,(a_invars(E)))

  b2, b4, b6, b8, c4, c6 = get_b_c_integral(E)

  split = true
  # test for types In, II, III, IV
  if mod(c4, p) != 0
    if quadroots(1, a1, -a2, p)
      cp = FlintZZ(n)
    elseif mod(n, 2) == 0
      cp = FlintZZ(2)
      split = false
    else
      cp = FlintZZ(1)
      split = false
    end

    Kp = "I$(n)"
    fp = FlintZZ(1)

    return (E, Kp, fp, cp, split)::Tuple{EllCrv{fmpq}, String, fmpz, fmpz, Bool}
  end

  if mod(a6, p^2) != 0
    Kp = "II"
    fp = FlintZZ(n)
    cp = FlintZZ(1)
    return (E, Kp, fp, cp, true)::Tuple{EllCrv{fmpq}, String, fmpz, fmpz, Bool}
  end

  if mod(b8, p^3) != 0
    Kp = "III"
    fp = FlintZZ(n-1)
    cp = FlintZZ(2)
    return (E, Kp, fp, cp, true)::Tuple{EllCrv{fmpq}, String, fmpz, fmpz, Bool}
  end

  if mod(b6, p^3) != 0
    if quadroots(1, divexact(a3, p), divexact(-a6, p^2), p)
      cp = FlintZZ(3)
    else
      cp = FlintZZ(1)
    end
    Kp = "IV"
    fp = n - 2
    return (E, Kp, FlintZZ(fp), cp, true)::Tuple{EllCrv{fmpq}, String, fmpz, fmpz, Bool}
  end

  # change coordinates so that p | a1, a2; p^2 | a3, a4; p^3 | a6
  if p == 2
    s = smod(a2, 2)
    t = 2 * smod(divexact(a6, 4), 2)
  else
    s = -a1 * invmod(FlintZZ(2), p)
    t = -a3 * invmod(FlintZZ(2), p)
  end

  trans = transform_rstu(E, fmpz[0, s, t, 1])
  E = trans[1]

  a1, a2, a3, a4, a6 = map(numerator,(a_invars(E)))

  b2, b4, b6, b8, c4, c6 = get_b_c_integral(E)

  # set up auxililary cubic T^3 + bT^2 + cT + d
  b = divexact(a2, p)
  c = divexact(a4, p^2)
  d = divexact(a6, p^3)
  w = 27*d^2 - b^2*c^2 + 4*b^3*d - 18*b*c*d + 4*c^3
  x = 3*c - b^2

  # test for distinct roots: type I0*
  if mod(w, p) != 0
    Kp = "I0*"
    fp = FlintZZ(n - 4)
    cp = 1 + nrootscubic(b, c, d, p)
    return (E, Kp, fp, FlintZZ(cp), true)::Tuple{EllCrv{fmpq}, String, fmpz, fmpz, Bool}

  # test for double root: type Im*
  elseif mod(x, p) != 0

    # change coordinates so that the double root is T = 0
    if p == 2
      r = c
    elseif p == 3
      r = b*c
    else
      r = (b*c - 9*d) * invmod(FlintZZ(2)*x, p)
    end

    r = p * smod(r, p)

    trans = transform_rstu(E, [r, 0, 0, 1])
    E = trans[1]

    a1, a2, a3, a4, a6 = map(numerator,(a_invars(E)))

    b2, b4, b6, b8, c4, c6 = get_b_c_integral(E)

    # make a3, a4, a6 repeatedly more divisible by p
    m = 1
    mx = p^2
    my = p^2
    cp = FlintZZ(0)

    while cp == 0
      xa2 = divexact(a2, p)
      xa3 = divexact(a3, my)
      xa4 = divexact(a4, p*mx)
      xa6 = divexact(a6, mx*my)

      if mod(xa3^2 + 4*xa6, p) !=  0
        if quadroots(1, xa3, -xa6, p)
          cp = FlintZZ(4)
        else
          cp = FlintZZ(2)
        end

      else
        if p == 2
          t = my * xa6
        else
          t = my * smod(-xa3*invmod(FlintZZ(2), p), p)
        end

        trans = transform_rstu(E, [0, 0, t, 1])
        E = trans[1]

        a1, a2, a3, a4, a6 = map(numerator,(a_invars(E)))

        b2, b4, b6, b8, c4, c6 = get_b_c_integral(E)

        my = my*p
        m = m + 1
        xa2 = divexact(a2, p)
        xa3 = divexact(a3, my)
        xa4 = divexact(a4, p*mx)
        xa6 = divexact(a6, mx*my)

        if mod(xa4^2 - 4*xa2*xa6, p) != 0
          if quadroots(xa2, xa4, xa6, p)
            cp = FlintZZ(4)
          else
            cp = FlintZZ(2)
          end

        else
          if p == 2
            r = mx * smod(xa6*xa2, 2)
          else
            r = mx * smod(-xa4 * invmod(2*xa2, p), p)
          end

          trans = transform_rstu(E, [r, 0, 0, 1])
          E = trans[1]

          a1, a2, a3, a4, a6 = map(numerator,(a_invars(E)))

          b2, b4, b6, b8, c4, c6 = get_b_c_integral(E)

          mx = mx*p
          m = m + 1
        end
      end
    end

    fp = n - m - 4
    Kp = "I$(m)*"

    return (E, Kp, FlintZZ(fp), FlintZZ(cp), true)::Tuple{EllCrv{fmpq}, String, fmpz, fmpz, Bool}

  else
    # Triple root case: types II*, III*, IV* or non-minimal
    # change coordinates so that the triple root is T = 0
    if p == 3
      rp = -d
    else
      rp = -b * invmod(FlintZZ(3), p)
    end

    r = p * smod(rp, p)

    trans = transform_rstu(E, [r, 0, 0, 1])
    E = trans[1]

    a1, a2, a3, a4, a6 = map(numerator,(a_invars(E)))

    b2, b4, b6, b8, c4, c6 = get_b_c_integral(E)

    x3 = divexact(a3, p^2)
    x6 = divexact(a6, p^4)

    # Test for type IV*
    if mod(x3^2 + 4* x6, p) != 0
      if quadroots(1, x3, -x6, p)
        cp = FlintZZ(3)
      else
        cp = FlintZZ(1)
      end
      Kp = "IV*"
      fp = FlintZZ(n - 6)

      return (E, Kp, fp, FlintZZ(cp), true)::Tuple{EllCrv{fmpq}, String, fmpz, fmpz, Bool}
    else
      if p == 2
        t = x6
      else
        t = x3 * invmod(FlintZZ(2), p)
      end

      t = -p^2 * smod(t, p)

      trans = transform_rstu(E, [0, 0, t, 1])
      E = trans[1]

      a1, a2, a3, a4, a6 = map(numerator,(a_invars(E)))

      b2, b4, b6, b8, c4, c6 = get_b_c_integral(E)

      # Test for types III*, II*
      if mod(a4, p^4) != 0
        Kp = "III*"
        fp = FlintZZ(n - 7)
        cp = FlintZZ(2)

        return (E, Kp, fp, FlintZZ(cp), true)::Tuple{EllCrv{fmpq}, String, fmpz, fmpz, Bool}
      elseif mod(a6, p^6) != 0
        Kp = "II*"
        fp = FlintZZ(n - 8)
        cp = FlintZZ(1)

        return (E, Kp, fp, FlintZZ(cp), true)::Tuple{EllCrv{fmpq}, String, fmpz, fmpz, Bool}
      else
        E = transform_rstu(E, [0, 0, 0, p])[1]
        return tates_algorithm_local(E, p)::Tuple{EllCrv{fmpq}, String, fmpz, fmpz, Bool}
      end
    end
  end
end

@doc Markdown.doc"""
    tates_algorithm_global(E::EllCrv{fmpq}) -> EllCrv{fmpq}

Return a global reduced minimal model for $E$ using Tate's algorithm.
"""
function tates_algorithm_global(E::EllCrv{fmpq})
  delta = abs(numerator(discriminant(E)))
  fac = factor(delta)

  p_list = [i[1] for i in fac]
  p_list = sort(p_list)

  output = []

  # apply tates algorithm successively for all primes dividing the discriminant
  for p in p_list
    E = tates_algorithm_local(E, p)[1]
  end

  # reduce coefficients (see tidy_model)
  E = tidy_model(E)

  return E::EllCrv{fmpq}
end

@doc Markdown.doc"""
    tamagawa number(E::EllCrv{fmpq}, p::Int) -> fmpz

Return the local Tamagawa number for E at p.
"""
function tamagawa_number(E::EllCrv{fmpq},p)
  return tates_algorithm_local(E,p)[4]
end

@doc Markdown.doc"""
    tamagawa numbers(E::EllCrv{fmpq}) -> Vector{(fmpz, fmpz)}

Return the sequence of Tamagawa numbers for $E$ at all the
bad primes $p$ of $E$.
"""
function tamagawa_numbers(E::EllCrv{fmpq})
  badp = bad_primes(E)
  return [tamagawa_number(E,p) for p in badp]
end

@doc Markdown.doc"""
    kodaira_symbol(E::EllCrv{fmpq}, p::Int) -> String

Return the reduction type of E at p using a Kodaira symbol.
"""
function kodaira_symbol(E::EllCrv{fmpq},p)
  return tates_algorithm_local(E,p)[2]
end

@doc Markdown.doc"""
    kodaira_symbols(E::EllCrv{fmpq}, p::Int) -> Vector{(fmpz, String)}

Return the reduction types of E at all bad primes as a sequence of
Kodaira symbols
"""
function kodaira_symbols(E::EllCrv{fmpq})
  badp = bad_primes(E)
  return [kodaira_symbol(E,p) for p in badp]
end

@doc Markdown.doc"""
    tamagawa number(E::EllCrv{fmpq}, p::NfOrdIdl) -> fmpz

Return the local Tamagawa number for E at p.
"""
function tamagawa_number(E::EllCrv{nf_elem},p::NfOrdIdl)
  return tates_algorithm_local(E,p)[4]
end

@doc Markdown.doc"""
    tamagawa numbers(E::EllCrv{fmpq}) -> Vector{(NfOrdIdl, fmpz)}

Return the sequence of Tamagawa numbers for $E$ at all the bad
prime ideals $p$ of $E$.
"""
function tamagawa_numbers(E::EllCrv{nf_elem})
  badp = bad_primes(E)
  return [tamagawa_number(E,p) for p in badp]
end

@doc Markdown.doc"""
    kodaira_symbol(E::EllCrv{nf_elem}, p::NfOrdIdl)
      -> String

Return the reduction type of E at the prime ideal p using
a Kodaira symbol.
"""
function kodaira_symbol(E::EllCrv{nf_elem},p::NfOrdIdl)
  return tates_algorithm_local(E,p)[2]
end

@doc Markdown.doc"""
    kodaira_symbols(E::EllCrv{nf_elem}, p::NfOrdIdl)
      -> Vector{(NfOrdIdl, String)}

Return the reduction types of E at all bad primes as a sequence of
Kodaira symbols.
"""
function kodaira_symbols(E::EllCrv{nf_elem})
  badp = bad_primes(E)
  return [kodaira_symbol(E,p) for p in badp]
end

@doc Markdown.doc"""
    reduction_type(E::EllCrv{fmpq}, p::fmpz) -> String

Return the reduction type of E at p. It can either be good, additive,
split multiplicative or nonsplit mutiplicative.
"""
function reduction_type(E::EllCrv{fmpq}, p)
  Ep, Kp, f, c, split = tates_algorithm_local(E, p)

  if Kp=="I0"
    return "Good"
  end

  if match(r"(I)([0-9]*)", Kp).match == Kp
    if split
      return "Split multiplicative"
    else
      return "Nonsplit multiplicative"
    end
  end

 return "Additive"

end

@doc Markdown.doc"""
    reduction_type(E::EllCrv{nf_elem}, p::NfOrdIdl) -> String

Return the reduction type of E at the prime ideal p.
It can either be good, additive, split multiplicative or
nonsplit mutiplicative.
"""
function reduction_type(E::EllCrv{nf_elem}, p::NfOrdIdl)
  Ep, Kp, f, c, split = tates_algorithm_local(E, p)

  if Kp=="I0"
    return "Good"
  end

  if match(r"(I)([0-9]*)", Kp).match == Kp
    if split
      return "Split multiplicative"
    else
      return "Nonsplit multiplicative"
    end
  end

 return "Additive"

end

################################################################################
#
#  Minimal model(s)
#
################################################################################

@doc Markdown.doc"""
    minimal_model(E::EllCrv{fmpq}) -> EllCrv{fmpq}

Returns the reduced global minimal model of $E$.
"""
function minimal_model(E::EllCrv{fmpq})
  return tates_algorithm_global(E)
end

@doc Markdown.doc"""
    minimal_model(E::EllCrv{fmpq}, p::Int) -> EllCrv{fmpq}

Returns a model of $E$, which is minimal at $p$. It is assumed that $p$
is prime.
"""
function minimal_model(E::EllCrv{fmpq}, p::Int)
  Ep = tates_algorithm_local(E, p)
  return Ep
end

@doc Markdown.doc"""
    minimal_model(E::EllCrv{fmpq}, p::NfOrdIdl) -> EllCrv{fmpq}

Returns a model of $E$, which is minimal at $p$. It is assumed that $p$
is a prime ideal.
"""
function minimal_model(E::EllCrv{fmpq}, p::NfOrdIdl)
  Ep = tates_algorithm_local(E, p)
  return Ep
end

@doc Markdown.doc"""
    tidy_model(E::EllCrv{fmpz}) -> EllCrv{fmpz}

Given an elliptic curve with minimal model, this functions returns an
isomorphic curve with reduced minimal model, that is, $a_1, a_3 \in \{0, 1\}$
and $a_2 \in \{-1,0,1\}$.
"""
function tidy_model(E::EllCrv{fmpq})

  a1, a2, a3, a4, a6 = map(numerator,(a_invars(E)))

  if mod(a1, 2) == 0
    s = -divexact(a1, 2)
  else
    s = -divexact(a1 - 1, 2)
  end

  if mod(-a2 + s*a1 + s^2, 3) == 0
    r = divexact(-a2 + s*a1 + s^2, 3)
  elseif mod(-a2 + s*a1 + s^2, 3) == 1
    r = divexact(-1 - a2 + s*a1 + s^2, 3)
  else
    r = divexact(1 - a2 + s*a1 + s^2, 3)
  end

  if mod(-a3 - r*a1, 2) == 0
    t = divexact(-a3 - r*a1, 2)
  else
    t = divexact(1 - a3 - r*a1, 2)
  end

  E = transform_rstu(E, [r, s, t, 1])[1]

  return E
end

################################################################################
#
#  Conductor
#
################################################################################

@doc Markdown.doc"""
    conductor(E::EllCrv{fmpq}) -> fmpz

Return the conductor of $E$ over QQ.
"""
function conductor(E::EllCrv{fmpq})
  badp = bad_primes(E)

  result = 1
  for p in badp
    result = result*(p^tates_algorithm_local(E,p)[3])
  end
  return result
end

@doc Markdown.doc"""
    conductor(E::EllCrv{nf_elem}) -> Nf)rdIdl

Return conductor of $E$ over a number field as an ideal.
"""
function conductor(E::EllCrv{nf_elem})
  badp = bad_primes(E)

  result = 1
  for p in badp
    result = result*(p^tates_algorithm_local(E,p)[3])
  end
  return result
end

#Magma returns the primes that divide the minimal discriminant
@doc Markdown.doc"""
    bad_primes(E::EllCrv{fmpq}) -> Vector{fmpz}

Return a list of the primes that divide the discriminant of $E$.
"""
function bad_primes(E::EllCrv{fmpq})

  d = ZZ(discriminant(E))
  L = factor(d)
  return [p for (p,e) in L]
end

@doc Markdown.doc"""
    bad_primes(E::EllCrv{fmpq}) -> Vector{NfOrdIdl}

Return a list of prime ideals that divide the discriminant of $E$.
"""
function bad_primes(E::EllCrv{nf_elem})
  R = ring_of_integers(base_field(E))
  d = R(discriminant(E))
  L = factor(d*R)
  return [p for (p,e) in L]
end

################################################################################
#
#  Reduction
#
################################################################################

#Magma also returns reduction map
@doc Markdown.doc"""
    bad_primes(E::EllCrv{nf_elem}, p::NfOrdIdl) -> EllCrv

Return the reduction of $E$ modulo the prime ideal p if p has good reduction
"""
function modp_reduction(E::EllCrv{nf_elem}, p::NfOrdIdl)
  if !isprime(p)
    throw(DomainError(p,"p is not a prime ideal"))
  end

  if p in bad_primes(E)
    throw(DomainError(p,"E has bad reduction at p"))
  end

  K, phi = ResidueField(order(p),p)

  a1, a2, a3, a4, a6 = map(phi,map(order(p), a_invars(E)))

  return EllipticCurve(K, [a1, a2, a3, a4, a6])

end

################################################################################
#
#  Integral invariants
#
################################################################################

# this needs to be rewritten
@doc Markdown.doc"""
    get_b_integral(E::EllCrv{fmpz}) -> Nemo.fmpz

Computes the invariants $b2$, $b4$, $b6$, $b8$ of an elliptic curve $E$ with integer coefficients.
"""
function get_b_integral(E)
  a1, a2, a3, a4, a6 = map(numerator,(a_invars(E)))

  b2 = a1^2 + 4*a2
  b4 = a1*a3 + 2*a4
  b6 = a3^2 + 4*a6
  b8 = a1^2*a6 - a1*a3*a4 + 4*a2*a6 + a2*a3^2 - a4^2

  return b2,b4,b6,b8
end

@doc Markdown.doc"""
    get_b_c_integral(E::EllCrv{fmpz}) -> Nemo.fmpz

Computes the invariants $b2$, $b4$, $b6$, $b8$, $c4$, $c6$ of an elliptic curve $E$ with integer coefficients.
"""
function get_b_c_integral(E)
    b2, b4, b6, b8 = get_b_integral(E)

    c4 = b2^2 - 24*b4
    c6 = -b2^3 + 36*b2*b4 - 216*b6

    return b2,b4,b6,b8,c4,c6
end


text:
using LightXML,TimeZones
smap = XMLDocument()

# create & attach a root node
smroot = create_root(smap, "urlset")

# set attributes
set_attribute(smroot,"xmlns","http://www.sitemaps.org/schemas/sitemap/0.9")
set_attribute(smroot,"xmlns:xsi","http://www.w3.org/2001/XMLSchema-instance")
set_attribute(smroot,"xsi:schemaLocation","http://www.sitemaps.org/schemas/sitemap/0.9 http://www.sitemaps.org/schemas/sitemap/0.9/sitemap.xsd")

# add html files
for fn in filter(s->occursin(".html",s),readdir(".."))
	url = new_child(smroot,"url")
	add_text(new_child(url,"loc"),"https://tobydriscoll.net/udmath305/"*fn)
	add_text(new_child(url,"lastmod"),"$(now(localzone()))")
	add_text(new_child(url,"priority"),"1.00")
end

save_file(smap, "../sitemap.xml")

ext: jl


text:
# This file was generated, do not modify it. # hide
e_tuned_iterated_pipe = evaluate(tuned_iterated_pipe, X, y,
                                 resampling=StratifiedCV(nfolds=6, rng=123),
                                 measures=[brier_loss, auc, accuracy])

max_stars_repo_name: alan-turing-institute/MLJTutorials

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
import Base: +, -, *, /, ^

+(x::RObject, y::RObject) = rcall(Symbol("+"), x, y)
-(x::RObject, y::RObject) = rcall(Symbol("-"), x, y)
*(x::RObject, y::RObject) = rcall(Symbol("*"), x, y)
/(x::RObject, y::RObject) = rcall(Symbol("/"), x, y)
^(x::RObject, y::RObject) = rcall(Symbol("^"), x, y)


text:
red_lighter         = "#FFEBEE"
red_light           = "#EF9A9A"
red                 = "#EF5350"
red_dark            = "#E53935"
red_darker          = "#C62828"

pink_lighter        = "#FCE4EC"
pink_light          = "#F48FB1"
pink                = "#EC407A"
pink_dark           = "#D81B60"
pink_darker         = "#AD1457"

purple_lighter      = "#F3E5F5"
purple_light        = "#CE93D8"
purple              = "#AB47BC"
purple_dark         = "#8E24AA"
purple_darker       = "#6A1B9A"

deep_purple_lighter = "#EDE7F6"
deep_purple_light   = "#B39DDB"
deep_purple         = "#7E57C2"
deep_purple_dark    = "#5E35B1"
deep_purple_darker  = "#4527A0"

indigo_lighter      = "#E8EAF6"
indigo_light        = "#9FA8DA"
indigo              = "#5C6BC0"
indigo_dark         = "#3949AB"
indigo_darker       = "#283593"

blue_lighter        = "#E3F2FD"
blue_light          = "#90CAF9"
blue                = "#42A5F5"
blue_dark           = "#1E88E5"
blue_darker         = "#1565C0"

light_blue_lighter  = "#E1F5FE"
light_blue_light    = "#81D4FA"
light_blue          = "#29B6F6"
light_blue_dark     = "#039BE5"
light_blue_darker   = "#0277BD"

cyan_lighter        = "#E0F7FA"
cyan_light          = "#80DEEA"
cyan                = "#26C6DA"
cyan_dark           = "#00ACC1"
cyan_darker         = "#00838F"

teal_lighter        = "#E0F2F1"
teal_light          = "#80CBC4"
teal                = "#26A69A"
teal_dark           = "#00897B"
teal_darker         = "#00695C"

green_lighter       = "#E8F5E9"
green_light         = "#A5D6A7"
green               = "#66BB6A"
green_dark          = "#43A047"
green_darker        = "#2E7D32"

light_green_lighter = "#F1F8E9"
light_green_light   = "#C5E1A5"
light_green         = "#9CCC65"
light_green_dark    = "#7CB342"
light_green_darker  = "#558B2F"

lime_lighter        = "#F9FBE7"
lime_light          = "#E6EE9C"
lime                = "#D4E157"
lime_dark           = "#C0CA33"
lime_darker         = "#9E9D24"

yellow_lighter      = "#FFFDE7"
yellow_light        = "#FFF59D"
yellow              = "#FFEE58"
yellow_dark         = "#FDD835"
yellow_darker       = "#F9A825"

amber_lighter       = "#FFF8E1"
amber_light         = "#FFE082"
amber               = "#FFCA28"
amber_dark          = "#FFB300"
amber_darker        = "#FF8F00"

orange_lighter      = "#FFF3E0"
orange_light        = "#FFCC80"
orange              = "#FFA726"
orange_dark         = "#FB8C00"
orange_darker       = "#EF6C00"

salmon_lighter      = "#FBE9E7"
salmon_light        = "#FFAB91"
salmon              = "#FF7043"
salmon_dark         = "#F4511E"
salmon_darker       = "#D84315"

brown_lighter       = "#EFEBE9"
brown_light         = "#BCAAA4"
brown               = "#8D6E63"
brown_dark          = "#6D4C41"
brown_darker        = "#4E342E"

grey_lighter        = "#FAFAFA"
grey_light          = "#EEEEEE"
grey                = "#BDBDBD"
grey_dark           = "#757575"
grey_darker         = "#424242"

blue_grey_lighter   = "#ECEFF1"
blue_grey_light     = "#B0BEC5"
blue_grey           = "#78909C"
blue_grey_dark      = "#546E7A"
blue_grey_darker    = "#37474F"

black               = "#000000"
white               = "#FFFFFF"


text:
@testset "combination: RGB" begin
    t1, t2 = RGBA{Float32}(1, 0.75, 0, 0.6), RGBA{Float64}(0, 0.5, 1, 0.4)
    c1, c2 = RGB{Float64}(t1), RGB{Float64}(t2)
    g1, g2 = AGray{Float32}(1, 0.6), AGray{Float64}(0, 0.4)

    for op in (CompositeClear,
               CompositeCopy,
               CompositeDestination,
               CompositeSourceOver,
               CompositeDestinationOver,
               CompositeSourceIn,
               CompositeDestinationIn,
               CompositeSourceOut,
               CompositeDestinationOut,
               CompositeSourceAtop,
               CompositeDestinationAtop,
               CompositeXor,
               CompositeLighter)
        @testset "op=$(keyword(op)), mode=$(keyword(mode))" for mode in (BlendNormal,
                                                                         BlendLighten,
                                                                         BlendLuminosity)
            b = blend(t1, t2, mode=mode, op=op)
            c = blend(c2, blend(c1, c2, mode=mode), opacity=alpha(t1))
            @test b === mode(t1, t2, op=op) === op(t1, t2, mode=mode)
            @test color(b) ≈ color(op(t1, RGBA(c, alpha(t2))))
            @test alpha(b) === alpha(blend(g1, g2, op=op))
        end
    end

end


text:
module TensorKitManifolds

export base, checkbase, isisometry, isunitary
export projecthermitian, projecthermitian!, projectantihermitian, projectantihermitian!
export projectcomplement, projectcomplement!, projectisometric, projectisometric!
export Grassmann, Stiefel, Unitary
export inner, retract, transport, transport!

using TensorKit

# Every submodule -- Grassmann, Stiefel, and Unitary -- implements their own methods for
# these. The signatures should be
# inner(W, Δ₁::Tangent, Δ₂::Tangent; metric)
# retract(W, Δ::Tangent, α::Real; alg)
# transport(Θ::Tangent, W, Δ::Tangent, α::Real, W′; alg)
# where the keyword arguments `alg` and `metric` should always be accepted, even if there is
# only one option for them and they are ignored. The `Tangent` is just a placeholder for the
# tangent type of each manifold. Similarly each submodule defines a `project!` function,
# which too should accept a keyword argument `metric`, even if it is ignored.
function inner end
function retract end
function transport end
function transport! end
function base end
function checkbase end
checkbase(x, y, z, args...) = checkbase(checkbase(x, y), z, args...)

# the machine epsilon for the elements of an object X, name inspired from eltype
eleps(X) = eps(real(eltype(X)))

function isisometry(W::AbstractTensorMap; tol = 10*eleps(W))
    WdW = W'*W
    s = zero(float(real(eltype(W))))
    for (c,b) in blocks(WdW)
        _subtractone!(b)
        s += dim(c)*length(b)
    end
    return norm(WdW) <= tol*sqrt(s)
end

isunitary(W::AbstractTensorMap; tol = 10*eleps(W)) =
    isisometry(W; tol = tol) && isisometry(W'; tol = tol)

function projecthermitian!(W::AbstractTensorMap)
    codomain(W) == domain(W) ||
        throw(DomainError("Tensor with distinct domain and codomain cannot be hermitian."))
    for (c,b) in blocks(W)
        _projecthermitian!(b)
    end
    return W
end
function projectantihermitian!(W::AbstractTensorMap)
    codomain(W) == domain(W) ||
        throw(DomainError("Tensor with distinct domain and codomain cannot be anithermitian."))
    for (c,b) in blocks(W)
        _projectantihermitian!(b)
    end
    return W
end

struct PolarNewton <: TensorKit.OrthogonalFactorizationAlgorithm
end
function projectisometric!(W::AbstractTensorMap; alg = Polar())
    if alg isa TensorKit.Polar || alg isa TensorKit.SDD
        foreach(blocks(W)) do (c,b)
            _polarsdd!(b)
        end
    elseif alg isa TensorKit.SVD
        foreach(blocks(W)) do (c,b)
            _polarsvd!(b)
        end
    elseif alg isa PolarNewton
        foreach(blocks(W)) do (c,b)
            _polarnewton!(b)
        end
    else
        throw(ArgumentError("unkown algorithm for projectisometric!: alg = $alg"))
    end
    return W
end

function projectcomplement!(X::AbstractTensorMap, W::AbstractTensorMap;
                            tol = 10*eleps(X))
    P = W'*X
    nP = norm(P)
    nX = norm(X)
    dP = dim(P)
    while nP > tol*max(dP, nX)
        X = mul!(X, W, P, -1, 1)
        P = W'*X
        nP = norm(P)
    end
    return X
end

projecthermitian(W::AbstractTensorMap) = projecthermitian!(copy(W))
projectantihermitian(W::AbstractTensorMap) = projectantihermitian!(copy(W))

function projectisometric(W::AbstractTensorMap;
                          alg::TensorKit.OrthogonalFactorizationAlgorithm = Polar())
    return projectisometric!(copy(W); alg=alg)
end
function projectcomplement(X::AbstractTensorMap, W::AbstractTensorMap,
                           tol = 10*eleps(X))
    return projectcomplement!(copy(X), W; tol=tol)
end

include("auxiliary.jl")
include("grassmann.jl")
include("stiefel.jl")
include("unitary.jl")

end


text:
Random.seed!(17)
N = 1000
net = CompleteBipartiteGraph(N, N, Network)
eprop!(net, "w", e -> rand())
res = Matching.run_bp(net, γ=0.1, maxiters=1000, verbose=false)
@test abs(res.energy - π^2/6) < 0.1
@test res.num_violations == 0 # no violated constraints

N = 1000
net = CompleteGraph(N, Network)
eprop!(net, "w", e -> rand())
res = Matching.run_bp(net, γ=0.1, maxiters=1000, verbose=false)
@test abs(res.energy - π^2/12) < 0.1
@test res.num_violations <= 5 # if graph is non-bipartite may fail


text:
using MyPlayground
using Test

@show Sys.CPU_THREADS
Sys.cpu_summary()

@testset begin
    @test true
    @test true
    @test true
    @test true
    @test true
end


text:
export ROCFFTError

struct ROCFFTError <: Exception
    code::rocfft_status
    msg::AbstractString
end
Base.show(io::IO, err::ROCFFTError) = print(io, "ROCFFTError(code $(err.code), $(err.msg))")

function ROCFFTError(code::rocfft_status)
    msg = status_message(code)
    return ROCFFTError(code, msg)
end

function status_message(status)
    if status == rocfft_status_success
        return "the operation completed successfully"
    elseif status == rocfft_status_failure
        return "an operation failed"
    elseif status == rocfft_status_invalid_arg_value
        return "invalid argument value"
    elseif status == rocfft_status_invalid_dimensions
        return "invalid array dimension"
    elseif status == rocfft_status_invalid_array_type
        return "invalid array type"
    elseif status == rocfft_status_invalid_strides
        return "invalid strides"
    elseif status == rocfft_status_invalid_distance
        return "invalid distance"
    elseif status == rocfft_status_invalid_offset
        return "invalid offset"
    else
        return "unknown error"
    end
end

macro check(fft_func)
    quote
        local err::rocfft_status
        err = $(esc(fft_func::Expr))
        if err != rocfft_status_success
            throw(ROCFFTError(err))
        end
        err
    end
end


text:
export weightedSample

function weightedSample(m::ConditionalModel, _data)
    return _weightedSample(getmoduletypencoding(m.model), m.model, m.args, _data)
end

export sourceWeightedSample

sourceWeightedSample(m::Model, data=NamedTuple()) = sourceWeightedSample(data)(m)

function sourceWeightedSample(_data)
    function(_m::Model)

        _datakeys = getntkeys(_data)
        proc(_m, st :: Assign)     = :($(st.x) = $(st.rhs))
        proc(_m, st :: Return)     = nothing
        proc(_m, st :: LineNumber) = nothing

        function proc(_m, st :: Sample)
            st.x ∈ _datakeys && return :(_ℓ += logdensity($(st.rhs), $(st.x)))
            return :($(st.x) = rand($(st.rhs)))
        end

        vals = map(x -> Expr(:(=), x,x),variables(_m))

        wrap(kernel) = @q begin
            _ℓ = 0.0
            $kernel

            return (_ℓ, $(Expr(:tuple, vals...)))
        end

        buildSource(_m, proc, wrap) |> MacroTools.flatten
    end
end

@gg function _weightedSample(M::Type{<:TypeLevel}, _m::Model, _args, _data)
    body = type2model(_m) |> sourceWeightedSample(_data) |> loadvals(_args, _data)
    @under_global from_type(_unwrap_type(M)) @q let M
        $body
    end
end


text:
for op in (:dirichlet,:neumann,:continuity,:ivp)
    @eval $op(d::PiecewiseSpace,k...) = InterlaceOperator($op(d.spaces,k...))
    @eval $op(d::UnionDomain,k...) = InterlaceOperator($op(d.domains,k...))
end

## diag provides a way to convert between DiagonalInterlaceOperator and bacn
# function blkdiagm{B<:Operator}(v::Vector{B})
#     if spacescompatible(map(domainspace,v)) && spacescompatible(map(rangespace,v))
#         # arrayspace
#         DiagonalInterlaceOperator(v)
#     elseif domainscompatible(map(domainspace,v)) && domainscompatible(map(rangespace,v))
#         # tuplespace
#         DiagonalInterlaceOperator(v)
#     else
#         DiagonalPiecewiseOperator(v)
#     end
# end

Base.blkdiag(A::DiagonalInterlaceOperator)=A.ops
Base.blkdiag(A::PlusOperator)=mapreduce(blkdiag,+,A.ops)
Base.blkdiag(A::TimesOperator)=mapreduce(blkdiag,.*,A.ops)

# TODO: general wrappers
for TYP in (:DerivativeWrapper,:ConversionWrapper)
    @eval Base.blkdiag{DT<:DiagonalInterlaceOperator}(A::($TYP{DT}))=A.op.ops
end

Base.blkdiag{FT<:PiecewiseSpace,OT<:DiagonalInterlaceOperator}(A::MultiplicationWrapper{FT,OT})=A.op.ops

Evaluation(S::TupleSpace,order::Number) =
    DiagonalInterlaceOperator(map(s->Evaluation(s,order),S),TupleSpace)

## Vector
# represents an operator applied to all spaces in an array space

immutable DiagonalArrayOperator{B<:Operator,T<:Number} <: Operator{T}
    op::B
    dimensions::Tuple{Vararg{Int}}
end

DiagonalArrayOperator{T}(op::Operator{T},dms::Tuple{Vararg{Int}})=DiagonalArrayOperator{typeof(op),T}(op,dms)
#DiagonalArrayOperator{T}(op::Operator{T},dms::Int)=DiagonalArrayOperator(op,(dms,))

function bandinds(D::DiagonalArrayOperator)
    bra,brb=bandinds(D.op)
    n=*(D.dimensions...)
    n*bra,n*brb
end

function getindex{B,T}(D::DiagonalArrayOperator{B,T},k::Integer,j::Integer)
    n=*(D.dimensions...)

    if mod(k,n) == mod(j,n)  # same block
        k=(k-1)÷n+1  # map k and j to block coordinates
        j=(j-1)÷n+1
        D.op[k,j]::T
    else
        zero(T)
    end
end

for op in (:domainspace,:rangespace)
    @eval $op(D::DiagonalArrayOperator)=ArraySpace($op(D.op),D.dimensions)
end

Derivative(AS::ArraySpace,k::Integer)=DerivativeWrapper(DiagonalArrayOperator(Derivative(AS.space,k),size(AS)),k)

function conversion_rule(AS::ArraySpace,BS::ArraySpace)
    if size(AS)==size(BS)
        ArraySpace(conversion_type(AS.space,BS.space),size(AS))
    else
        NoSpace()
    end
end

for OP in (:maxspace,:conversion_type)
    @eval function $OP(AS::ArraySpace,BS::ArraySpace)
        if size(AS)==size(BS)
            ArraySpace($OP(AS.space,BS.space),size(AS))
        else
            NoSpace()
        end
    end
end

function Conversion(AS::ArraySpace,BS::ArraySpace)
    @assert size(AS)==size(BS)
    ConversionWrapper(DiagonalArrayOperator(Conversion(AS.space,BS.space),size(AS)))
end

ToeplitzOperator{S,T,V,DD}(G::Fun{MatrixSpace{S,T,DD,1},V})=interlace(map(ToeplitzOperator,mat(G)))

## Sum Space

## Conversion

# TupleSpace maps down

function coefficients(v::AbstractVector,a::TupleSpace,b::TupleSpace)
    vs=vec(Fun(v,a))
    coefficients(detuple(map((f,s)->Fun(f,s),vs,b)))
end

function Conversion(a::TupleSpace,b::TupleSpace)
    m=findlast(s->dimension(s)==1,a.spaces)
    @assert all(s->dimension(s)==1,a.spaces[1:m-1]) &&
            a.spaces[1:m-1]==b.spaces[1:m-1]

    if m==0
        ConversionWrapper(DiagonalInterlaceOperator(map(Conversion,a.spaces,b.spaces),TupleSpace))
    elseif length(a)==m
        SpaceOperator(FiniteOperator(eye(m)),a,b)
    elseif length(a)==m+1
        ConversionWrapper(BlockOperator(eye(m),zeros(m,0),zeros(0,m),Conversion(a[end],b[end])))
    else
        ConversionWrapper(BlockOperator(eye(m),zeros(m,0),zeros(0,m),Conversion(TupleSpace(a[m+1:end]),TupleSpace(b[m+1:end]))))
    end
end

# Sum Space and PiecewiseSpace need to allow permutation of space orders
for TYP in (:SumSpace,:PiecewiseSpace)
    @eval function Conversion(S1::$TYP,S2::$TYP)
        if any(s->!isinf(dimension(s)),S1.spaces) || any(s->!isinf(dimension(s)),S2.spaces)
            error("Need to implement finite dimensional case")
        elseif sort([S1.spaces...])==sort([S2.spaces...])
            # swaps sumspace order
            ConversionWrapper(SpaceOperator(
            PermutationOperator(promote_type(eltype(domain(S1)),eltype(domain(S2))),S1.spaces,S2.spaces),
                          S1,S2))
        elseif all(map(hasconversion,S1.spaces,S2.spaces))
            # we can blocmk convert
            ConversionWrapper(DiagonalInterlaceOperator(map(Conversion,S1.spaces,S2.spaces),$TYP))
        elseif map(canonicalspace,S1.spaces)==map(canonicalspace,S2.spaces)
            error("Not implemented")
        elseif sort([map(canonicalspace,S1.spaces)...])==sort([map(canonicalspace,S2.spaces)...])
            # we can block convert after permuting
            P=PermutationOperator(promote_type(eltype(domain(S1)),eltype(domain(S2))),
                                  map(canonicalspace,S1.spaces),
                                  map(canonicalspace,S2.spaces))
            ds2=$TYP(S1.spaces[P.perm])
            ConversionWrapper(TimesOperator(Conversion(ds2,S2),SpaceOperator(P,S1,ds2)))
        elseif all(map(hasconversion,reverse(S1.spaces),S2.spaces))
            # special case that comes up, especially for two spaces
            rS1=SumSpace(reverse(S1.spaces))
            ConversionWrapper(TimesOperator(
                Conversion(rS1,S2),
                SpaceOperator(PermutationOperator(length(S1.spaces):-1:1),
                            S1,rS1)))
        elseif all(map(hasconversion,sort([map(canonicalspace,S1.spaces)...]),sort([map(canonicalspace,S2.spaces)...])))
            #TODO: general case
            @assert length(S1.spaces)==2
            ds2=$TYP(S1.spaces[[2,1]])
            ConversionWrapper(TimesOperator(Conversion(ds2,S2),Conversion(S1,ds2)))

        else
            # we don't know how to convert so go to default
            defaultConversion(S1,S2)
        end
    end
end

for (OPrule,OP) in ((:conversion_rule,:conversion_type),(:maxspace_rule,:maxspace),
                        (:union_rule,:union))
    for TYP in (:SumSpace,:PiecewiseSpace)
        @eval function $OPrule(S1::$TYP,S2::$TYP)
            cs1,cs2=map(canonicalspace,S1.spaces),map(canonicalspace,S2.spaces)
            if length(S1.spaces)!=length(S2.spaces)
                NoSpace()
            elseif canonicalspace(S1)==canonicalspace(S2)  # this sorts S1 and S2
                S1 ≤ S2?S1:S2  # choose smallest space by sorting
            elseif cs1==cs2
                # we can just map down
                # $TYP(map($OP,S1.spaces,S2.spaces))
                # this is commented out due to Issue #13261
                newspaces=[$OP(S1[k],S2[k]) for k=1:length(S1.spaces)]
                if any(b->b==NoSpace(),newspaces)
                    NoSpace()
                else
                    $TYP(newspaces)
                end
            elseif sort([cs1...])== sort([cs2...])
                # sort S1
                p=perm(cs1,cs2)
                $OP($TYP(S1.spaces[p]),S2)
            elseif length(S1)==length(S2)==2  &&
                    $OP(S1[1],S2[1])!=NoSpace() &&
                    $OP(S1[2],S2[2])!=NoSpace()
                #TODO: general length
                $TYP($OP(S1[1],S2[1]),
                     $OP(S1[2],S2[2]))
            elseif length(S1)==length(S2)==2  &&
                    $OP(S1[1],S2[2])!=NoSpace() &&
                    $OP(S1[2],S2[1])!=NoSpace()
                #TODO: general length
                $TYP($OP(S1[1],S2[2]),
                     $OP(S1[2],S2[1]))
            else
                NoSpace()
            end
        end
    end

    # TupleSpace doesn't allow reordering
    @eval function $OPrule(S1::TupleSpace,S2::TupleSpace)
        K=findfirst(s->!isa(s,ConstantSpace),S1)-1

        if length(S1)==length(S2)  &&
                all(s->isa(s,ConstantSpace),S2[1:K]) &&
                all(s->!isa(s,ConstantSpace),S2[K+1:end])
            newspaces=[$OP(S1[k],S2[k]) for k=K+1:length(S1)]

            if any(b->b==NoSpace(),newspaces)
                NoSpace()
            else
                TupleSpace(tuple(S1[1:K]...,newspaces...))
            end
        else
            NoSpace()
        end
    end

end

## Derivative

#TODO: do in @calculus_operator?

for TYP in (:PiecewiseSpace,:TupleSpace),(Op,OpWrap) in ((:Derivative,:DerivativeWrapper),
                                                         (:Integral,:IntegralWrapper))
    @eval $Op(S::$TYP,k::Integer)=$OpWrap(DiagonalInterlaceOperator(map(s->$Op(s,k),S.spaces),$TYP),k)
end

function Derivative(S::SumSpace,k::Integer)
    # we want to map before we decompose, as the map may introduce
    # mixed bases.
    if typeof(canonicaldomain(S))==typeof(domain(S))
        DerivativeWrapper(DiagonalInterlaceOperator(map(s->Derivative(s,k),S.spaces),SumSpace),k)
    else
        defaultDerivative(S,k)
    end
end

choosedomainspace(M::CalculusOperator{UnsetSpace},sp::SumSpace)=mapreduce(s->choosedomainspace(M,s),union,sp.spaces)

## Multiplcation for Array*Vector

function Multiplication{S,T,DD,S2,T2,DD2,dim}(f::Fun{MatrixSpace{S,T,DD,dim}},sp::VectorSpace{S2,T2,DD2,dim})
    @assert size(space(f),2)==length(sp)
    m=mat(f)
    MultiplicationWrapper(f,interlace(Operator{promote_type(eltype(f),eltype(sp))}[Multiplication(m[k,j],sp.space) for k=1:size(m,1),j=1:size(m,2)]))
end

## Multiply pieces

function Multiplication{PW<:PiecewiseSpace}(f::Fun{PW},sp::PiecewiseSpace)
    p=perm(domain(f).domains,domain(sp).domains)  # sort f
    vf=pieces(f)[p]

    MultiplicationWrapper(f,DiagonalInterlaceOperator(map(Multiplication,vf,sp.spaces),PiecewiseSpace))
end

Multiplication{SV1,SV2,T2,T1,D,d}(f::Fun{SumSpace{SV1,T1,D,d}},sp::SumSpace{SV2,T2,D,d}) =
    MultiplicationWrapper(f,mapreduce(g->Multiplication(g,sp),+,vec(f)))
Multiplication(f::Fun,sp::SumSpace) =
    MultiplicationWrapper(f,DiagonalInterlaceOperator(map(s->Multiplication(f,s),vec(sp)),SumSpace))

# we override coefficienttimes to split the multiplication down to components as union may combine spaes

coefficienttimes{S1<:SumSpace,S2<:SumSpace}(f::Fun{S1},g::Fun{S2}) = mapreduce(ff->ff*g,+,vec(f))
coefficienttimes{S1<:SumSpace}(f::Fun{S1},g::Fun) = mapreduce(ff->ff*g,+,vec(f))
coefficienttimes{S2<:SumSpace}(f::Fun,g::Fun{S2}) = mapreduce(gg->f*gg,+,vec(g))

coefficienttimes{S1<:PiecewiseSpace,S2<:PiecewiseSpace}(f::Fun{S1},g::Fun{S2})=depiece(map(coefficienttimes,pieces(f),pieces(g)))

## Definite Integral

# This makes sure that the defaults from a given Domain are respected for the UnionDomain.

DefiniteIntegral(d::UnionDomain) =
    ConcreteDefiniteIntegral(PiecewiseSpace(map(domainspace,map(DefiniteIntegral,d.domains))))
DefiniteLineIntegral(d::UnionDomain) =
    ConcreteDefiniteLineIntegral(PiecewiseSpace(map(domainspace,map(DefiniteLineIntegral,d.domains))))

####### This is a hack to get the Faraday Cage working.
function getindex{PWS<:PiecewiseSpace,T}(Σ::ConcreteDefiniteLineIntegral{PWS,T},k::Integer)
    n = length(domain(Σ))
    k ≤ n? one(T) : zero(T)
end
bandinds{PWS<:PiecewiseSpace,T}(Σ::ConcreteDefiniteLineIntegral{PWS,T}) =
    0,arclength(domain(Σ))-1
####### This is a hack to get the Faraday Cage working.

## TensorSpace of two PiecewiseSpaces

Base.getindex{PWS1<:PiecewiseSpace,PWS2<:PiecewiseSpace}(d::TensorSpace{Tuple{PWS1,PWS2}},i::Integer,j::Integer)=d[1][i]⊗d[2][j]
Base.getindex{PWS1<:PiecewiseSpace,PWS2<:PiecewiseSpace}(d::TensorSpace{Tuple{PWS1,PWS2}},i::Range,j::Range)=PiecewiseSpace(d[1][i])⊗PiecewiseSpace(d[2][j])

## ProductFun

##  Piecewise

function pieces{PS<:PiecewiseSpace}(U::ProductFun{PS})
    ps=space(U,1)
    sp2=space(U,2)
    m=length(ps)
    C=coefficients(U)
    [ProductFun(C[k:m:end,:],ps[k],sp2) for k=1:m]
end


text:
module LiquidsNEqProcesses

### Imports
using Reexport
using Parameters
using LinearAlgebra: I

@reexport using LiquidsDynamics

### Exports
export Interpolator, find_arrest, interpolate!

### Implementation
include("interpolation.jl")
include("arrest.jl")

end # module


text:
################################################################################
#  Copyright 2020, Tom Van Acker                                               #
################################################################################
# MultiStateSystems.jl                                                         #
# A Julia package to solve multi-state system models.                          #
# See http://github.com/timmyfaraday/MultiStateSystems.jl                      #
################################################################################

# load pkgs
using Test
using Unitful
using Measurements
using AdditionalUnits
using MultiStateSystems

# pkg const
const _MSM = Measurements
const _MSS = MultiStateSystems

@testset "MultiStateSystems" begin
    # form
    include("steady_state_process.jl")
    include("markov_process.jl")

    # io
    include("indices.jl")
    include("universal_generating_function.jl")

    # prob
    include("universal_generating_operator.jl")
end


text:
# ---
# title: 933. Number of Recent Calls
# id: problem933
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Queue
# link: <https://leetcode.com/problems/number-of-recent-calls/description/>
# hidden: true
# ---
#
# You have a `RecentCounter` class which counts the number of recent requests
# within a certain time frame.
#
# Implement the `RecentCounter` class:
#
#   * `RecentCounter()` Initializes the counter with zero recent requests.
#   * `int ping(int t)` Adds a new request at time `t`, where `t` represents some time in milliseconds, and returns the number of requests that has happened in the past `3000` milliseconds (including the new request). Specifically, return the number of requests that have happened in the inclusive range `[t - 3000, t]`.
#
# It is **guaranteed** that every call to `ping` uses a strictly larger value of
# `t` than the previous call.
#
#
#
# **Example 1:**
#
#
#
#     **Input**
#     ["RecentCounter", "ping", "ping", "ping", "ping"]
#     [[], [1], [100], [3001], [3002]]
#     **Output**
#     [null, 1, 2, 3, 3]
#
#     **Explanation**
#     RecentCounter recentCounter = new RecentCounter();
#     recentCounter.ping(1);     // requests = [ _1_ ], range is [-2999,1], return 1
#     recentCounter.ping(100);   // requests = [ _1_ , _100_ ], range is [-2900,100], return 2
#     recentCounter.ping(3001);  // requests = [ _1_ , _100_ , _3001_ ], range is [1,3001], return 3
#     recentCounter.ping(3002);  // requests = [1, _100_ , _3001_ , _3002_ ], range is [2,3002], return 3
#
#
#
#
# **Constraints:**
#
#   * `1 <= t <= 109`
#   * Each test case will call `ping` with **strictly increasing** values of `t`.
#   * At most `104` calls will be made to `ping`.
#
#
## @lc code=start
using LeetCode

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


text:
Pkg.checkout("POMDPs")

using POMDPs

POMDPs.add("MCTS") # for visualization


text:
# Day 1: Data Types

# Declared variables
i = 1
d = 2.0
s = "Hello "
# Declare second integer, double, and String variables.
# Read and save an integer, double, and String to your variables.
newI = parse(UInt64, readline())
newD = parse(Float64, readline())
newS = readline()

# Print the sum of both integer variables on a new line.
println(i + newI)
# Print the sum of the double variables on a new line.
println(d + newD)
# Concatenate and print the String variables on a new line
# The 's' variable above should be printed first.
println(string(s , newS))

max_stars_repo_path: Hackerrank/30 Days of Code/Julia/day 01.jl

max_stars_repo_name: Next-Gen-UI/Code-Dynamics

max_issues_repo_name: Next-Gen-UI/Code-Dynamics

max_forks_repo_name: Next-Gen-UI/Code-Dynamics


text:
## Translations for both the key and value
const KeywordDict = Dict{Symbol,Symbol}()

for s in (:complex, :basic, :modulus, :gaussian, :force, :deep, :none,
          :method, :fu, :combined, :default, :groebner, :measure,
          :double, :trig, :log, :denom, :numer, :check, :numerical,
          :minimal, :warn, :simplify, :rational, :manual, :implicit,
          :particular, :quick, :cubics, :quartics, :quintics)

    KeywordDict[Symbol(uppercasefirst(string(s)))] = s
end
KeywordDict[:Conditions] = :conds

translate_keyword(s) = get(KeywordDict,s,s)

@sjdoc Deep """
   Deep

an option for the functions `Together` and `Expand`.
"""

@sjdoc Conditions """
   Conditions

an option for `Integrate` and other integral transforms.
"""

@sjdoc Modulus """
    Modulus

an option for `Factor`.
"""

@sjdoc Gaussian """
    Gaussian

an option specifiying that solutions including Gaussian integers should be returned.
"""

# Convert Mxpr to sympy, pulling out Rule(a,b) to dict of keyword args.
# That is, we separate keyword args from positional argss
function sjtopy_kw(mx::T, kws) where T<:Mxpr
#    println("Looking for kw $kws")
    args = margs(mx)
    nargs = newargs()
    for i in 1:length(args)
        if is_Mxpr(args[i], :Rule)
            lhs = translate_keyword(args[i][1])
            rhs = translate_keyword(args[i][2])
            kws[lhs] = rhs
        else
            push!(nargs, sjtopy(args[i]))
        end
    end
    nargs
end

function sjtopy_kw(mx::T) where T<:Mxpr
    kws = Dict()  # type ? probably symbols
    nargs  = sjtopy_kw(mx, kws)
    return (nargs, kws)
end

## Called manually, not by mkapprules
## Not very efficient
function _separate_rules(args, kws; keylowercase=false)
    nargs = newargs()
    for i in 1:length(args)
        if isa(args[i], Mxpr{:Rule})
            length(args[i]) != 2 && error("Rule requires two arguments. " * length(args[i]) * " found.")
            k = translate_keyword(args[i][1])
            v = translate_keyword(args[i][2])
            if keylowercase
                if isa(k,Symbol)
                    k = Symbol(lowercase(string(k)))
                end
                if isa(v,Symbol)
                    v = Symbol(lowercase(string(v)))
                end
            end
            kws[k] = v
        else
            push!(nargs, args[i])
        end
    end
    nargs
end

## Called manually, not by mkapprules
function separate_rules1(args...; inkws...)
    kws = Dict()
    nargs  = _separate_rules(args, kws; inkws...)
    return (nargs, kws)
end

# Separate the Rule()'s from other arguments in an Mxpr expression
# Store keywords in a Dict so they can by passed as keword arguments.
# These do the same as above, but no conversion to sympy.
function separate_rules(mx::T, kws) where T<:Mxpr
    args = margs(mx)
    nargs = newargs()
    for i in 1:length(args)
        if is_Mxpr(args[i], :Rule)
            length(args[i]) != 2 && error("Rule requires two arguments. " * length(args[i]) * " found.")
            k = translate_keyword(args[i][1])
            v = translate_keyword(args[i][2])
            kws[k] = v
        else
            push!(nargs, args[i])
        end
    end
    nargs
end

function separate_rules(mx::T) where T<:Mxpr
    kws = Dict()  # type ? probably symbols
    nargs  = separate_rules(mx, kws)
    return (nargs, kws)
end

# Mma uses the expression Rule(a,b) to represent a keyword
# argument. It also uses Rule for many other things. This is awkward.
# Here, we separate keyword Rules from all other args including Rules
# that are not meant to be keywords.
# kws -- Dict of legal keywords with their default values.
# Only Rules with keywords in kws will be extracted
function separate_known_rules(mx::T, kws) where T<:Mxpr
    args = margs(mx)
    nargs = newargs()
    for i in 1:length(args)
        a = args[i]
        if is_Mxpr(a, :Rule)
            @checknargs a :Rule 2
            if haskey(kws,a[1])
                kws[a[1]] = a[2]
            else
                push!(nargs, a)
            end
        else
            push!(nargs, a)
        end
    end
    nargs
end


text:
module ORBIT

using ..ParticleSwarm
using ..Sequencing # Includes Clifford, Waveforms, and Qubit

# Fitness calls here are of the form
# "How well does <waveform> implement <pulse> on <qubit>?"

# The overarching "ORBIT" call will include a qubit and a gate, and define the
# fitness anonymously through the below.
function fitness(params::Vector{UInt16}, q::Qubit, p::Pulse)
  pulselength = div(length(params), 3)
  q[p] = ExactWaveform(params[1:pulselength], params[pulselength:2*pulselength],
                       params[2:pulselength:3*pulselength], true)
  pulseIndex = findfirst(SQPulse, p)

  # Every "1" in the readout detracts from the fitness score.
  1 - mean([sendSequenceAndReadout(q, benchmark1Qubit(
                                   pulseIndex, orbitSequenceLength))
        for _ in 1:orbitAveraging])
end

# Same but for FloatWaveforms, where we are optimizing vectors of Float64
function fitness(params::Vector{Float64}, q::Qubit, p::Pulse)
  q[p] = params
  pulseIndex = findfirst(SQPulse, p)
  1 - mean([sendSequenceAndReadout(q, benchmark1Qubit(
                                   pulseIndex, orbitSequenceLength))
        for _ in 1:orbitAveraging])
end

function optimize(q::Qubit, p::Pulse, numpoints::Integer, init::Waveform)
  if isa(init, ExactWaveform)
    if isa(q, QubitWithZ)
      varCount = 3*numpoints
      guess = [init.XYI; init.XYQ; init.Z]
    else
      varCount = 2*numpoints
      guess = [init.XYI; init.XYQ]
    end
  else
    varCount = numpoints
    guess = init.wavedata
  end

  if p[1] < 7 && isa(q.lineXYI[1], InstrumentControl.InsAWG5014C) # XY on AWG
end

end # End module ORBIT


text:
# -*- encoding: utf-8 -*-
#
# The MIT License (MIT)
#
# Copyright © 2021 Matteo Foglieni and Riccardo Gervasoni
#

@testset "test_Hit" begin
     cube = Cube()
     ray1 = Ray(Point(0.0, 0.0, 0.0), VEC_X)

     intersection1 = ray_intersection(cube, ray1)
     @test typeof(intersection1) == HitRecord
     @test HitRecord(
          Point(0.5, 0.0, 0.0),
          Normal(-1.0, 0.0, 0.0),
          Vec2d(0.5, 0.5),
          0.5,
          ray1
          ) ≈ intersection1

     ray2 = Ray( Point(0.0, 2.0, 0.0), -VEC_Y )
     intersection2 = ray_intersection(cube, ray2)
     @test typeof(intersection2) == HitRecord
     @test HitRecord(
          Point(0.0, 0.5, 0.0),
          Normal(0.0, 1.0, 0.0),
          Vec2d(0.5, 0.5),
          1.5,
          ray2
          ) ≈ intersection2

     ray3 = Ray( Point(0.0, 0.0, -1.5), VEC_Z )
     intersection3 = ray_intersection(cube, ray3)
     @test typeof(intersection3) == HitRecord
     @test HitRecord(
          Point(0.0, 0.0, -0.5),
          Normal(0.0, 0.0, -1.0),
          Vec2d(0.5, 0.5),
          1.0,
          ray3
          ) ≈ intersection3
end

@testset "test_noHit" begin
     cube = Cube()

     ray1 = Ray(Point(0, 0, 1), VEC_X)
     intersection1 = ray_intersection(cube, ray1)
     @test isnothing(intersection1)

     ray2 = Ray(Point(1.5, 1, 0), VEC_X+VEC_Y)
     intersection2 = ray_intersection(cube, ray2)
     @test isnothing(intersection2)

     ray3 = Ray(Point(0, 0, -1), -VEC_Z)
     intersection3 = ray_intersection(cube, ray3)
     @test isnothing(intersection3)
end

@testset "test_Transformation" begin
     cube1 = Cube( rotation_z(π/2.0) * translation(Vec(3.5, 0.0, 0.0)))

     ray1 = Ray(Point(0, 1, 0), VEC_Y)
     intersection1 = ray_intersection(cube1, ray1)
     @test typeof(intersection1) == HitRecord
     @test HitRecord(
          Point(0.0, 3.0, 0.0),
          Normal(0.0, -1.0, 0.0),
          Vec2d(0.5, 0.5),
          2.0,
          ray1
          ) ≈ intersection1

     cube2 = Cube(translation(Vec(0., 0., 1.5))*rotation_x(π/2.0))
     ray2 = Ray(Point(0.0, 0.0, 0.0), Vec(0., 0., 1.))
     intersection2 = ray_intersection(cube2, ray2)
     @test HitRecord(
          Point(0., 0., 1.),
          Normal(0.0, 0., -1.0),
          Vec2d(0.5, 0.5),
          1.0,
          ray2
          ) ≈ intersection2

     # Check if the plane failed to move by trying to hit the untransformed shape
     @test isnothing( ray_intersection(cube1, Ray( Point(0, 0, 0), -VEC_Z ) ) )
end

@testset "test_Normals" begin
     cube = Cube()

     P = Point(0, 0, 0.6)
     Q = Point(0.6, 0, 0)
     ray1 = Ray(P, Q-P)
     intersection1 = ray_intersection(cube, ray1)
     @test intersection1.normal ≈ Normal(0.0, 0.0, 1.0)

     ray2 = Ray(Point(-1.0, 0.0, 0.0), Vec(3.0, 0.25, 0.25))
     intersection2 = ray_intersection(cube, ray2)
     @test intersection2.normal ≈ Normal(-1.0, 0.0, 0.0)
end

@testset "test_UV_Coordinates" begin
     cube = Cube(translation(0.0, 0.5, 0.5))

     ray1 = Ray(Point(-1.0, 0.25, 0.5), VEC_X)
     @test ray_intersection(cube, ray1).surface_point ≈ Vec2d(0.25, 0.5)

     ray2 = Ray(Point(-1.0, 0.25, 0.75), VEC_X)
     @test ray_intersection(cube, ray2).surface_point ≈ Vec2d(0.25, 0.75)

     ray3 = Ray(Point(0.25, 0.75, 0.25), -VEC_X)
     @test ray_intersection(cube, ray3).surface_point ≈ Vec2d(0.75, 0.25)
end


text:
using Documenter, MeshAlgorithms

makedocs(;
    modules=[MeshAlgorithms],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/sjrodahl/MeshAlgorithms.jl/blob/{commit}{path}#L{line}",
    sitename="MeshAlgorithms.jl",
    authors="Sondre Rodahl",
#   assets=[],
)

#deploydocs(;
#    repo="github.com/sjrodahl/MeshAlgorithms.jl",
#)


text:
# We mark `isless()` as `@pure` because it's used in computing the correct
# (sorted) monomial order during type promotion, so its results must be
# assumed to be hyper-pure. Overloading this method in user code would be
# a bad idea.
@pure function isless(::Type{V1}, ::Type{V2}) where {V1 <: Variable, V2 <: Variable}
    name(V1) < name(V2)
end

# We should be able to get rid of this method as the fallback just
# redirects to `promote_rule(monomialtype(V1), monomialtype(V2))`
# but we should check first the implications with this `@pure` thing.
function promote_rule(::Type{V1}, ::Type{V2}) where {V1 <: Variable, V2 <: Variable}
    if V1 < V2
        Monomial{(V1(), V2()), 2}
    else
        Monomial{(V2(), V1()), 2}
    end
end

function _promote_monomial_noncommutative(::Type{Monomial{V1, N1}}, ::Type{Monomial{V2, N2}}) where {N1, V1, N2, V2}
    if V1 > V2
        _promote_monomial(Monomial{V2, N2}, Monomial{V1, N1})
    else
        varnames = shortest_common_supersequence(name.(V1), name.(V2))
        vars = Expr(:tuple, [Expr(:call, Expr(:curly, :Variable, Expr(:quote, n))) for n in varnames]...)
        quote
            Monomial{$vars, $(length(varnames))}
        end
    end
end

@pure merge(v1::Tuple{Vararg{Variable}}) = v1

@pure function merge(v1::Tuple{Vararg{Variable}}, v2::Variable)
    for i in 1:length(v1)
        if v2 === v1[i]
            return v1
        elseif v2 > v1[i]
            return tuple(v1[1:i-1]..., v2, v1[i:end]...)
        end
    end
    return tuple(v1..., v2)
end

@pure function merge(v1::Tuple{Vararg{Variable}}, v2::Variable, vars::Variable...)
    merge(merge(v1, v2), vars...)
end

@pure function promote_rule(::Type{<:Monomial{V1}}, ::Type{<:Monomial{V2}}) where {V1, V2}
    vars = merge(V1, V2...)
    Monomial{vars, length(vars)}
end

MP.promote_rule(::Type{Term{C,M1} where {C}}, M2::Type{<:MonomialLike}) where {M1} = (Term{C,promote_type(M1, M2)} where {C})
MP.promote_rule(M1::Type{<:MonomialLike}, ::Type{Term{C,M2} where {C}}) where {M2} = (Term{C,promote_type(M1, M2)} where {C})
MP.promote_rule(::Type{Term{C,M1} where {C}}, ::Type{Term{T,M2}}) where {T,M1,M2} = (Term{C,promote_type(M1, M2)} where {C})
MP.promote_rule(::Type{Term{T,M2}}, ::Type{Term{C,M1} where {C}}) where {T,M1,M2} = (Term{C,promote_type(M1, M2)} where {C})
MP.promote_rule(::Type{<:Polynomial}, ::Type{Term{C,M} where {C}}) where {M} = Polynomial
MP.promote_rule(::Type{Term{C,M} where {C}}, ::Type{<:Polynomial}) where {M} = Polynomial
MP.promote_rule(::Type{Polynomial}, ::Type{<:PolynomialLike}) = Polynomial
MP.promote_rule(::Type{<:PolynomialLike}, ::Type{Polynomial}) = Polynomial
MP.promote_rule_constant(::Type, ::Type{Polynomial}) = Any
MP.promote_rule_constant(::Type, ::Type{Term{C,M} where {C}}) where {M} = Any

Base.convert(::Type{Term{C,M} where {C}}, m::MonomialLike) where {M} = MP.term(convert(M, m))
Base.convert(::Type{Polynomial}, t::TermLike) = MP.polynomial(t)


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule GRASS_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("GRASS")
JLLWrappers.@generate_main_file("GRASS", UUID("58433f90-7b86-51a1-96c6-55447b2352ff"))
end  # module GRASS_jll


text:
using DataFrames, CSV, Images
using Base.Iterators: partition

data_cat = ["train", "valid"]

"""
  get_count_df(path_1, path_2, category)

Add a count of the number of images are patient column
and globally define a DataFrame.

Arguments:
1. `path_1`: Path to `labeled_studies` csv.
2. `path_2`: Path to the `image_paths` csv.
3. `category`: `train` or `valid`
"""

function get_count_df(path_1, path_2, category)
  df_1 = CSV.read(path_1, header = ["path", "label"])
  df_2 = CSV.read(path_2, header = ["path"])
  df_1[:count] = 0
  i = 1
  for (j, paths) in enumerate(df_1[:path])
    while(paths == df_2[:path][i][1:(findlast(isequal('/'), df_2[:path][i]))])
      df_1[:count][j] += 1
      i += 1
      if i == size(df_2)[1]
        break
      end
    end
  end
  if(category == "train")
    global df_train = df_1
  else
    global df_valid = df_1
  end
end

"""
  get_indices(study_type, category)

Get the starting and ending indices in the dataframe created
by calling `get_count_df` corresponding to the `study_type`

Arguments:
1. `study_type`: One of the 7 Study Types.
2. `category`: `train` or `valid`
"""

function get_indices(study_type, category)
  start_ind = 0
  end_ind = 0
  if category == "train"
    df = df_train
  else
    df = df_valid
  end
  for (i, path) in enumerate(df[:path])
    if split(path, "/")[3] == study_type
      if start_ind == 0
        start_ind = i
      end
    else
      if start_ind != 0 && end_ind == 0
        end_ind = i - 1
        break
      end
    end
  end
  if end_ind == 0
    end_ind = size(df)[1]
  end
  (start_ind, end_ind)
end

"""
  get_study_data(study_type)

Returns a dictionary with 2 DataFrames containing the image
paths for the `train` and `valid` set corresponding to the
`study_type`.

Arguments:
1. `study_type`: One of the 7 Study Types.
"""

function get_study_data(study_type)
  a1, b1 = get_indices(study_type, "train")
  a2, b2 = get_indices(study_type, "valid")
  Dict("train" => df_train[a1:b1, :], "valid" => df_valid[a2:b2, :])
end

im_mean = reshape([0.485, 0.456, 0.406], 1, 1, 3)
im_std = reshape([0.229, 0.224, 0.225], 1, 1, 3)

_normalize(img) = (img .- im_mean) ./ im_std

_denormalize(img) = img .* im_std .+ im_std

# Convert all the images to Grayscale since some are already grayscale and this might reduce train time
im2arr(img) = permutedims(reshape(float.(channelview(Gray.(imresize(img, (224, 224))))), 224, 224, 1), (2, 1, 3))

"""
  get_batched_images(study_type, batch_size; path_t1 = "",
                     path_t2 = "", path_v1 = "", path_v2 = "")

Load the images of the train and validation set and get it into proper
shape for running the classifier.

Arguments:
1. `study_type`: One of the 7 study types.
2. `batch_size`: Batch Size for training and validation.
3. `path_t1`: Optional path to the `train_labeled_studies.csv`.
4. `path_t2`: Optional path to the `train_image_paths.csv`.
5. `path_v1`: Optional path to the `valid_labeled_studies.csv`.
6. `path_v2`: Optional path to the `valid_image_paths.csv`.
"""

function get_batched_images(study_type, batch_size; path_t1 = "",
                            path_t2 = "", path_v1 = "", path_v2 = "")
  try # Check if df_train is loaded
    df_train[1]
  catch
    warn("df_train was not loaded")
    get_count_df(path_t1, path_t2, "train")
  end
  try # Check if df_valid is loaded
    df_valid[1]
  catch
    warn("df_valid was not loaded")
    get_count_df(path_v1, path_v2, "valid")
  end
  dict = get_study_data(study_type)
  images = Dict("train_imgs" => [], "train_labs" => [],
                "valid_imgs" => [], "valid_labs" => [])
  c = 0
  for cate in data_cat
    for (i, path) in enumerate(dict[cate][:path])
      for j in 1:dict[cate][:count][i]
        img = im2arr(load(dict[cate][:path][i] * "image$(j).png"))
        push!(images[cate * "_imgs"], 255 * img)
        push!(images[cate * "_labs"], dict[cate][:label][i])
        c += 1
        if(c%1000 == 0)
          info("$c images loaded")
          gc() # Without manual garbage collection cache resources might be exhausted
        end
      end
    end
  end
  gc() # Clear cache once images have been loaded
  # Randomly shuffle the train images
  perm_train_inds = randperm(length(images["train_imgs"]))
  images["train_imgs"] = images["train_imgs"][perm_train_inds]
  images["train_labs"] = images["train_labs"][perm_train_inds]
  Dict("train" => [(cat(4, images["train_imgs"][i]...), reshape(Float32.(images["train_labs"][i]), 1, :))
        for i in partition(1:length(images["train_imgs"]), batch_size)],
       "valid" => [(cat(4, images["valid_imgs"][i]...), reshape(Float32.(images["valid_labs"][i]), 1, :))
        for i in partition(1:length(images["valid_imgs"]), batch_size)])
end


text:
# =========== temp

println(Base.Docs.doc(scatter))
println(Main.Docs.doc(Makie.convert_arguments, Tuple{Any, AbstractMatrix}))
println(Base.Docs.doc(Makie.convert_arguments, Tuple{Any, AbstractVector, AbstractVector, AbstractMatrix}))
println(Base.Docs.doc(Makie.convert_arguments, Tuple{Any, AbstractVector, AbstractVector, Function}))

const atomics = (
    heatmap,
    image,
    lines,
    linesegments,
    mesh,
    meshscatter,
    scatter,
    surface,
    text,
    Makie.volume
)

# =============================================
# Trying out heatmap plot from all_samples.jl
path = joinpath(@__DIR__, "..", "docs", "test")
scene = Scene(resolution = (500, 500))
scene = heatmap!(scene, rand(32, 32))
save("heatmap.png", scene)

path = joinpath(@__DIR__, "..", "docs", "test")
# save("heatmap.png", scene)
img = Makie.scene2image(scene)
save(path, img)

# =============================================
# Print code for when database search returns multiple results

# First find those damn entries!
entries = example_database(lines)
entries = findall(x -> x.title == database_key, database)
entries = findall(x -> x.tags == "scatter", database)
entries = findall(x -> occursin("scatter", collect(x.tags)), database)

tgs = database[2].tags

occursin("scatter", tgs)

findall(database) do entry
    # find tags
    tags_found = any(tag -> string(tag) in entry.tags, "scatter")
end

len_entries = length(entries)
println("len_entries = ", len_entries)

for i = 1:len_entries
    # println(x.title)
    # println(x.source)
    println("i = ", i)
    sprint() do io
        print_code(
            io, database, entries[1],
            scope_start = "",
            scope_end = "",
            indent = "",
            resolution = (entry)-> "resolution = (500, 500)",
            outputfile = (entry, ending)-> Pkg.dir("Makie", "docs", "media", string(entry.unique_name, ending))
        )
    end
end

sprint(print_code(stdout, database, entries[1],
scope_start = "",
scope_end = "",
indent = "",
resolution = (entry)-> "resolution = (500, 500)",
outputfile = (entry, ending)-> Pkg.dir("Makie", "docs", "media", string(entry.unique_name, ending))
))

# autogenerate example library including plots test
# TODO: not working right now since a lot of plots are broken
# TODO: need to figure out if there is a way to "close" a Makie plot window
pathroot = joinpath(@__DIR__, "..", "docs", "src", "library-test")
imgpath = joinpath(pathroot, "plots")
mdpath = joinpath(pathroot, "library-test.md")
open(mdpath, "w") do io
    println(io, "# Examples library test -- autogenerated")
    isempty(database) && error("database is empty!")
    counter = 1
    for idx in 1:length(database)
        println(io, "## index = $idx")
        # println(io, "## `$(Makie.to_string(func))`")
        try
            println(io, "Example $counter, \"$(database[idx].title)\"")
            _print_source(io, idx; style = "julia")
            println(io, "`plot goes here`\n")
            # TODO: add code to generate + embed plots
            counter += 1
        catch
            println("ERROR: Didn't work with $tag at index $idx\n")
        end
        println(io, "\n")
    end
end


text:
# Epochs more-or-less need to define offsets relative to one another (as a Duration), and that's it!

"""
    abstract Clock

An `Clock` defines the instant of time to which a clock measures its time
relative to, i.e. the instant at time = 0. A `Time` object is simply a
combination of a `Duration` since this instant and a `Clock`.

If the relative difference between two clocks is known, then the two clocks can
be subtracted to get this difference as a duration. If the ticking rates of the
two clocks are different, the result should be in the units of the first clock
(i.e. `c1` in `c1 - c2`). This duration will be used automatically when
converting `Time` from one clock to another.
"""
abstract type Clock end

########################
## TAI time standards ##
########################

"""
     abstract TAIClock <: Clock

International Atomic Time (TAI) is used to define mulitple other standard
measurements of time, such as UTC and GPS time. Time standards defined against
TAI can be easily transformed between one another by knowing the relative
difference between their epochs (where the clock time = 0).

The relative difference between standard TAI time (at a given date) and a
custom-defined `TAIClock` class can be given by defining the `epoch()` of the
clock (as a `Base.Dates.Date`) and by defining the subtraction operation between
the clock and `TAIDate`. This will allow times to be converted between your
custom clock and any other TAI-derived clock.
"""
abstract type TAIClock <: Clock end

"""
    TAIDate <: TAIClock
    TAIDate(date)
    TAIDate(year, month, day)

An International Atomic Time (TAI) clock with epoch (time = 0) set at the
beginning of the given `date`.
"""
struct TAIDate <: TAIClock
    date::Date
end
TAIDate(year, month, day) = TAIDate(Date(year, month, day))

"""
    epoch(clock)

Returns the date corresponding to the epoch (starting time, t = 0) of the clock.
E.g. `epoch(TAIDate(Dates.Date(1980, 1, 6)))` is 1980-01-06.
"""
epoch(tai_date::TAIDate) = tai_date.date

function Base.:-(c1::TAIDate, c2::TAIDate)
    return (epoch(c1) - epoch(c2)).value * days
end

# By default, the TAI-derived times will go through TAIDate as an intermediatary
function Base.:-(c1::C1, c2::C2) where {C1 <: TAIClock, C2 <: TAIClock}
    return (c1 - TAIDate(epoch(c1))) - (c2 - TAIDate(epoch(c1)))
end

########################
## GPS time standards ##
########################

"""
    GPSEpoch <: TAIClock
    GPSEpoch()

The GPS time standard is defined relative to the start of 6th January, 1980.
Equivalent to `GPSWeek(0)`.
"""
struct GPSEpoch <: TAIClock; end

#Base.:-(::GPSEpoch, ::GPSEpoch) = Duration{1, ZeroUnit}(zerounit)

Base.show(io::IO, ::GPSEpoch) = print(io, "GPS epoch")

function Base.:-(c1::GPSEpoch, c2::GPSEpoch)
    return Duration{1}(zerounit)
end

function Base.:-(c1::GPSEpoch, c2::TAIDate)
    return (epoch(c1) - epoch(c2)).value * days - gps_offset
end

function Base.:-(c1::TAIDate, c2::GPSEpoch)
    return (epoch(c1) - epoch(c2)).value * days + gps_offset
end

Base.@pure epoch(::GPSEpoch) = Date(1980, 1, 6)

"""
    GPSWeek <: TAIClock
    GPSWeek(week)

The GPS epoch began on Sunday 6th January, 1980, and the GPS broadcast signal
resets to zero every Sunday at midnight. For instance, week 1908 corresponds to
the week beginning 2016-07-31 (a Sunday).
"""
struct GPSWeek <: TAIClock
    week::Int
end

function epoch(gps_week::GPSWeek)
    epoch(GPSEpoch()) + Week(gps_week.week)
end

Base.show(io::IO, c::GPSWeek) = print(io, "GPS week $(c.week) ($(epoch(c)))")

"""
    gpsweek(date)

Get the GPS week for the given date, i.e. how many Sundays have passed between
6th January 1980 and `date` (including `date`).
"""
gpsweek(date::Date) = fld(Duration(date - epoch(GPSEpoch())), weeks)

"""
    GPSWeek(date)
    GPSWeek(day, year, month)

Create a `GPSWeek` clock corresponding to the given `date`.
"""
GPSWeek(date::Date) = GPSWeek(gpsweek(date))
GPSWeek(year, month, day) = GPSWeek(gpsweek(Date(year, month, day)))

const gps_offset = 19seconds

function Base.:-(c1::GPSWeek, c2::GPSWeek)
    return (c1.week - c2.week)*weeks
end

function Base.:-(c1::GPSWeek, c2::TAIDate)
    return (epoch(c1) - epoch(c2)).value * days - gps_offset
end

function Base.:-(c1::TAIDate, c2::GPSWeek)
    return (epoch(c1) - epoch(c2)).value * days + gps_offset
end

########################
## UTC time standards ##
########################

"""
    UTCDate <: TAIClock
    UTCDate(date)
    UTCDate(year, month, day)

A clock which starts (defines it epoch) at the beggining of the specificied
Coordinated Universal Time (UTC) day.
"""
struct UTCDate <: TAIClock
    date::Date
end

UTCDate(year, month, day) = UTCDate(Date(year, month, day))

epoch(utc_date::UTCDate) = utc_date.date

Base.show(io::IO, utc_date::UTCDate) = (print(io, utc_date.date); print(io, " (UTC)"))

include("../gen/leap_seconds.jl")

function lookup_leapsecs(date::UTCDate)
    d = epoch(date)
    if d > leap_secs_expiry_date
        throw(ErrorException("Unknown number of leap seconds for date $d"))
    end
    (searchsortedlast(utc_offsets, d) + leap_index_offset)seconds
end

function Base.:-(c1::UTCDate, c2::TAIDate)
    return (epoch(c1) - epoch(c2)).value * days - lookup_leapsecs(c1)
end

function Base.:-(c1::TAIDate, c2::UTCDate)
    return (epoch(c1) - epoch(c2)).value * days + lookup_leapsecs(c2)
end

function Base.:-(c1::UTCDate, c2::UTCDate)
    return (epoch(c1) - epoch(c2)).value * days - (lookup_leapsecs(c1) - lookup_leapsecs(c2))
end

#=
immutable OffsetEpoch{E <: Epoch} <: Epoch
    offset::Duration
    epoch::E
end

Base.:-(e1::OffsetEpoch, e2::Epoch)       = e1.duration + (e1.epoch - e2)
Base.:-(e1::Epoch, e2::OffsetEpoch)       = e1.duration + (e1 - e2.epoch)
Base.:-(e1::OffsetEpoch, e2::OffsetEpoch) = e1.duration - e2.duration + (e1.epoch - e2.epoch)

function Base.show(io::IO, e::OffsetEpoch)
    print(io, e.duration)
    print(io, " after ")
    print(io, e.epoch)
end
=#


text:
####################################
# The AbstractBlockArray interface #
####################################

"""
    abstract AbstractBlockArray{T, N} <: AbstractArray{T, N}

The abstract type that represents a blocked array. Types that implement
the `AbstractBlockArray` interface should subtype from this type.

** Typealiases **

* `AbstractBlockMatrix{T}` -> `AbstractBlockArray{T, 2}`

* `AbstractBlockVector{T}` -> `AbstractBlockArray{T, 1}`

* `AbstractBlockVecOrMat{T}` -> `Union{AbstractBlockMatrix{T}, AbstractBlockVector{T}}`
"""
abstract type AbstractBlockArray{T, N} <: LayoutArray{T, N} end
const AbstractBlockMatrix{T} = AbstractBlockArray{T, 2}
const AbstractBlockVector{T} = AbstractBlockArray{T, 1}
const AbstractBlockVecOrMat{T} = Union{AbstractBlockMatrix{T}, AbstractBlockVector{T}}

block2string(b, s) = string(join(map(string,b), '×'), "-blocked ", Base.dims2string(s))
_block_summary(a) = string(block2string(blocksize(a), size(a)), " ", typeof(a))
Base.summary(a::AbstractBlockArray) = _block_summary(a)
_show_typeof(io, a) = show(io, typeof(a))
function _block_summary(io, a)
    print(io, block2string(blocksize(a), size(a)))
    print(io, ' ')
    _show_typeof(io, a)
end
Base.summary(io::IO, a::AbstractBlockArray) = _block_summary(io, a)

# avoid to_shape which complicates axes
Base.similar(a::AbstractBlockArray{T}) where {T}                             = similar(a, T)
Base.similar(a::AbstractBlockArray, ::Type{T}) where {T}                     = similar(a, T, axes(a))
Base.similar(a::AbstractBlockArray{T}, dims::Tuple) where {T}                = similar(a, T, dims)

Base.IndexStyle(::Type{<:AbstractBlockArray}) = IndexCartesian()

# need to overload axes to return BlockAxis
@inline size(block_array::AbstractBlockArray) = map(length, axes(block_array))
@inline axes(block_array::AbstractBlockArray) = throw(error("axes for ", typeof(block_array), " is not implemented"))

"""
    getblock(A, inds...)

Returns the block at blockindex `inds...`. An alternative syntax is `A[Block(inds...)].
Throws a `BlockBoundsError` if this block is out of bounds.

```jldoctest; setup = quote using BlockArrays end
julia> v = Array(reshape(1:6, (2, 3)))
2×3 Array{Int64,2}:
 1  3  5
 2  4  6

julia> A = BlockArray(v, [1,1], [2,1])
2×2-blocked 2×3 BlockArray{Int64,2}:
 1  3  │  5
 ──────┼───
 2  4  │  6

julia> getblock(A, 2, 1)
1×2 Array{Int64,2}:
 2  4

julia> A[Block(1, 2)]
1×1 Array{Int64,2}:
 5
```
"""
function getblock(A::AbstractBlockArray{T,N}, ::Vararg{Integer, N}) where {T,N}
    throw(error("getblock for ", typeof(A), " is not implemented"))
end

"""
    getblock!(X, A, inds...)

Stores the block at blockindex `inds` in `X` and returns it. Throws a `BlockBoundsError` if the
attempted assigned block is out of bounds.

```jldoctest; setup = quote using BlockArrays end
julia> A = PseudoBlockArray(ones(2, 3), [1, 1], [2, 1])
2×2-blocked 2×3 PseudoBlockArray{Float64,2}:
 1.0  1.0  │  1.0
 ──────────┼─────
 1.0  1.0  │  1.0

julia> x = zeros(1, 2);

julia> getblock!(x, A, 2, 1);

julia> x
1×2 Array{Float64,2}:
 1.0  1.0
```
"""
getblock!(X, A::AbstractBlockArray{T,N}, ::Vararg{Integer, N}) where {T,N} = throw(error("getblock! for ", typeof(A), " is not implemented"))

@inline getblock!(X, A::AbstractBlockArray{T,N}, block::Block{N}) where {T,N}             = getblock!(X, A, block.n...)
@inline getblock!(X, A::AbstractBlockVector, block::Block{1})                       = getblock!(X, A, block.n[1])
@inline getblock!(X, A::AbstractBlockArray{T, N}, block::Vararg{Block{1}, N}) where {T,N} = getblock!(X, A, (Block(block).n)...)

"""
    setblock!(A, v, inds...)

Stores the block `v` in the block at block index `inds` in `A`. An alternative syntax is `A[Block(inds...)] = v`.
Throws a `BlockBoundsError` if this block is out of bounds.

```jldoctest; setup = quote using BlockArrays end
julia> A = PseudoBlockArray(zeros(2, 3), [1, 1], [2, 1]);

julia> setblock!(A, [1 2], 1, 1);

julia> A[Block(2, 1)] = [3 4];

julia> A
2×2-blocked 2×3 PseudoBlockArray{Float64,2}:
 1.0  2.0  │  0.0
 ──────────┼─────
 3.0  4.0  │  0.0
```
"""
setblock!(A::AbstractBlockArray{T,N}, v, ::Vararg{Integer, N}) where {T,N} = throw(error("setblock! for ", typeof(A), " is not implemented"))

@inline setblock!(A::AbstractBlockArray{T, N}, v, block::Block{N}) where {T,N}      = setblock!(A, v, block.n...)
@inline setblock!(A::AbstractBlockVector, v, block::Block{1})                       = setblock!(A, v, block.n[1])
@inline setblock!(A::AbstractBlockArray{T, N}, v, block::Vararg{Block{1}, N}) where {T,N} = setblock!(A, v, (Block(block).n)...)

"""
    BlockBoundsError([A], [inds...])

Thrown when a block indexing operation into a block array, `A`, tried to access an out-of-bounds block, `inds`.
"""
struct BlockBoundsError <: Exception
    a::Any
    i::Any
    BlockBoundsError() = new()
    BlockBoundsError(a::AbstractArray) = new(a)
    BlockBoundsError(a::AbstractArray, @nospecialize(i)) = new(a,i)
end

BlockBoundsError(a::AbstractArray, I::Block) = BlockBoundsError(a, I.n)

function Base.showerror(io::IO, ex::BlockBoundsError)
    print(io, "BlockBoundsError")
    if isdefined(ex, :a)
        print(io, ": attempt to access ")
            print(io, summary(ex.a))
        if isdefined(ex, :i)
            print(io, " at block index [")
            join(io, ex.i, ',')
            print(io, ']')
        end
    end
end

"""
    blockcheckbounds(A, inds...)

Throw a `BlockBoundsError` if the specified block indexes are not in bounds for the given block array.
Subtypes of `AbstractBlockArray` should
specialize this method if they need to provide custom block bounds checking behaviors.

```jldoctest; setup = quote using BlockArrays end
julia> A = BlockArray(rand(2,3), [1,1], [2,1]);

julia> blockcheckbounds(A, 3, 2)
ERROR: BlockBoundsError: attempt to access 2×2-blocked 2×3 BlockArray{Float64,2,Array{Array{Float64,2},2},Tuple{BlockedUnitRange{Array{Int64,1}},BlockedUnitRange{Array{Int64,1}}}} at block index [3,2]
[...]
```
"""
@inline function blockcheckbounds(A::AbstractArray{T, N}, i::Vararg{Integer, N}) where {T,N}
    if blockcheckbounds(Bool, A, i...)
        return
    else
        throw(BlockBoundsError(A, i))
    end
end

@inline function blockcheckbounds(::Type{Bool}, A::AbstractArray{T, N}, i::Vararg{Integer, N}) where {T,N}
    n = blockaxes(A)
    k = 0
    for idx in 1:N # using enumerate here will allocate
        k += 1
        @inbounds _i = i[idx]
        Block(_i) in n[k] || return false
    end
    return true
end

blockcheckbounds(A::AbstractArray{T, N}, i::Block{N}) where {T,N} = blockcheckbounds(A, i.n...)

# Convert to @generated...
@propagate_inbounds Base.getindex( block_arr::AbstractBlockArray{T, N}, block::Block{N}) where {T,N}       =  getblock(block_arr, block.n...)
@propagate_inbounds Base.setindex!(block_arr::AbstractBlockArray{T, N}, v, block::Block{N}) where {T,N}    =  setblock!(block_arr, v, block.n...)
@propagate_inbounds Base.getindex( block_arr::AbstractBlockVector, block::Block{1})                  =  getblock(block_arr, block.n[1])
@propagate_inbounds Base.setindex!(block_arr::AbstractBlockVector, v, block::Block{1})               =  setblock!(block_arr, v, block.n[1])
@inline Base.getindex(block_arr::AbstractBlockArray{T,N}, block::Vararg{Block{1}, N}) where {T,N}     =  getblock(block_arr, (Block(block).n)...)
@inline Base.setindex!(block_arr::AbstractBlockArray{T,N}, v, block::Vararg{Block{1}, N}) where {T,N} =  setblock!(block_arr, v, (Block(block).n)...)

"""
    eachblock(A::AbstractBlockArray)

Create a generator that iterates over each block of an `AbstractBlockArray`
returning views.

```jldoctest; setup = quote using BlockArrays end
julia> v = Array(reshape(1:6, (2, 3)))
2×3 Array{Int64,2}:
 1  3  5
 2  4  6

julia> A = BlockArray(v, [1,1], [2,1])
2×2-blocked 2×3 BlockArray{Int64,2}:
 1  3  │  5
 ──────┼───
 2  4  │  6

julia> Matrix.(collect(eachblock(A)))
2×2 Array{Array{Int64,2},2}:
 [1 3]  [5]
 [2 4]  [6]

julia> sum.(eachblock(A))
2×2 Array{Int64,2}:
 4  5
 6  6
```
"""
function eachblock(A::AbstractBlockArray)
    # blockinds = CartesianIndices(blocksize(A))
    blockinds = CartesianIndices(axes.(blocklasts.(axes(A)),1))
    (view(A, Block(Tuple(I))) for I in blockinds)
end


text:
"""

transformation between locale identifiers in differnet forms
1. form used by POSIX:  <lang>_<region>[.<charset>][@<cextension>]
2. form used by BCP47:  <lang>[_<script>][_<region>][-<variant>][-lextensions]
"""
module LocaleIdTranslations

using ResourceBundles
using Unicode
import ResourceBundles: SEP

export cloc_to_loc, loc_to_cloc, all_locales

"""

    cloc_to_loc(posix::String) -> unicode-locale-id-string

Posix string has the general form `<lang>[_<country>][.<charset>][@<extension>]`.
We transform this to the following string:
`<lang[_<script>][_country>][_<variant>][-x-posix-<extension>]`.
The `@extension` may be transformed to a <script>, a variant, or a private extension.
The charset is ignored. The extension is optional in input and output.
"""
function cloc_to_loc(cloc::String)
    cloc == "" && return cloc
    if cloc == "C" || uppercase(cloc) == "POSIX"
        return "C"
    end
    cloc = lowercase(cloc)
    m = match(REG_POSIX, cloc)
    if m == nothing
        return cloc
    end
    mc = m.captures
    lang = mc[1]
    reg = nostring(mc[3])
    charset = nostring(mc[5])
    ext = nostring(mc[7])
    script = ""
    var = ""
    if !isempty(ext)
        langext = string(lang, mc[6])
        if haskey(EXTENSION_TO_SCRIPT, ext)
            script = EXTENSION_TO_SCRIPT[ext]
            ext = ""
        elseif haskey(LANGEXT_TO_VAR, langext)
            var = LANGEXT_TO_VAR[langext]
            ext = ""
        elseif haskey(LANGEXT_TO_NEW, langext)
            lang = LANGEXT_TO_NEW[key]
            ext = ""
        elseif ext in EXTENSIONS_IGNORE
            ext = ""
        end
    end
    io = IOBuffer()
    write(io, lang)
    isempty(script) || write(io, SEP, script)
    isempty(reg) || write(io, SEP, reg)
    isempty(var) || write(io, SEP, var)
    isempty(ext) || write(io, "-x-posix-", ext)
    String(take!(io))
end

nostring(s::AbstractString) = s
nostring(::Nothing) = ""

const REG_POSIX = r"(^[[:alpha:]]+)([_-]([[:alpha:]]+))?(\.([[:alnum:]_-]+))?(@([[:alnum:]]+))?$"

const EXTENSION_TO_SCRIPT = Dict(
        "cyrillic" => "Cyrl",
        "latin" => "Latn",
        "devanagari" => "Deva")

const LANGEXT_TO_VAR = Dict(
        "es@valencia" => "valencia"
)

const LANGEXT_TO_NEW = Dict(
        "aa@saaho" => "ssy",
        "gez@abegede" => "aa")

const EXTENSIONS_IGNORE = ["euro"]

const S0 = Symbol("")

"""
    loc_to_cloc(loc::LocaleId) -> POSIX string

Translate from LocaleId to String in POSIX Format.
Comparable with cloc_to_loc.
"""
function loc_to_cloc(loc::LocaleId)
    s = string(loc)
    ( s == "C" || uppercase(s) == "POSIX" ) && return "C"
    s == "" && return s
    lang = loc.language
    reg = loc.region
    script = loc.script
    # only interested in extensions of the form "-x-posix-<ext>"
    extd = loc.extensions === nothing ? Symbol[] : get(loc.extensions, 'x', Symbol[])
    ext = length(extd) == 2 && extd[1] == :posix ? extd[2] : S0

    pext = get(LANG_TO_EXT, string(lang), "")

    allo = all_locales() # in the operating system storage (locale -a)
    if reg != S0 # all posix locales have a region
        if script != S0
            key = string(script)
            if haskey(SCRIPT_TO_EXT, key)
                pext = SCRIPT_TO_EXT[key]
            end
        end
        if ext != S0 && isempty(pext)
            pext = string(ext)
        end
        su = locale_name(lang, reg, pext)
        if findfirst(x -> x == su, allo) != 0
            return su
        end
        su = locale_name(loc.language, reg)
        if findfirst(x -> x == su, allo) != 0
            return su
        else
            s = loc.language
        end
    else
        s = loc.language
    end

    s = string(s)
    su = locale_name(s, uppercase(s)) # check the case, when lang and reg are the same
    if findfirst( x ->x == su, allo) != 0
        return su
    end

    reg = get(PROBABLE_SUBTAGS, s, "")
    if reg != ""
        su = locale_name(s, reg)
        if findfirst( x ->x == su, allo) != 0
            return su
        end
    end

    ix = findfirst(x -> x == s || startswith(x, s * '_'), allo)
    return ix != 0 ? allo[ix] : "C"
end

locale_name(lang, reg) = string(lang, '_', reg, UTF8)
locale_name(lang, reg, ext) = isempty(ext) ? locale_name(lang, reg) : string(lang, '_', reg, UTF8, '@', ext)

const PROBABLE_SUBTAGS = Dict("en" => "US", "zh" => "CN", "sv" => "SE")
const SCRIPT_TO_EXT = Dict("Deva" => "devanagari", "Latn" => "latin", "Cyrl" => "cyrillic")
const LANG_TO_EXT = Dict("ssy" => "saaho")

allloc = String[]

function all_locales()
    global allloc
    if isempty(allloc)
        list = String[]
        push!(list, "C", "POSIX")
        for loc in eachline(`locale -a`)
            loca, ext = splitext(loc)
            simple = ext == "" && '_' ∉ loca
            standard = loca == "C" || loca == "POSIX"
            if ( simple || ext == UTF8 ) && !standard
                push!(list, loc)
            end
        end
        allloc = list
    end
    allloc
end

const UTF8 = ".utf8"

end # module


text: flattenandsort = sort ∘ collect ∘ Iterators.flatten

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas