text:
# From here on, this file contains specializations of DG methods on the
# curved 3D meshes `StructuredMesh{3}, P4estMesh{3}` to the compressible
# Euler equations.
#
# The specialized methods contain relatively verbose and ugly code in the sense
# that we do not use the common "pointwise physics" interface of Trixi.jl.
# However, this is currently (November 2021) necessary to get a significant
# speed-up by using SIMD instructions via LoopVectorization.jl.
#
# TODO: SIMD. We could get even more speed-up if we did not need to permute
#             array dimensions, i.e., if we changed the basic memory layout.
#
# We do not wrap this code in `@muladd begin ... end` block. Optimizations like
# this are handled automatically by LoopVectorization.jl.

# We specialize on `PtrArray` since these will be returned by `Trixi.wrap_array`
# if LoopVectorization.jl can handle the array types. This ensures that `@turbo`
# works efficiently here.
@inline function split_form_kernel!(_du::PtrArray, u_cons::PtrArray,
                                    element, mesh::Union{StructuredMesh{3}, P4estMesh{3}},
                                    nonconservative_terms::Val{false},
                                    equations::CompressibleEulerEquations3D,
                                    volume_flux::typeof(flux_shima_etal_turbo),
                                    dg::DGSEM, cache, alpha)
  @unpack derivative_split = dg.basis
  @unpack contravariant_vectors = cache.elements

  # Create a temporary array that will be used to store the RHS with permuted
  # indices `[i, j, k, v]` to allow using SIMD instructions.
  # `StrideArray`s with purely static dimensions do not allocate on the heap.
  du = StrideArray{eltype(u_cons)}(undef,
    (ntuple(_ -> StaticInt(nnodes(dg)), ndims(mesh))...,
     StaticInt(nvariables(equations))))

  # Convert conserved to primitive variables on the given `element`.
  u_prim = StrideArray{eltype(u_cons)}(undef,
    (ntuple(_ -> StaticInt(nnodes(dg)), ndims(mesh))...,
     StaticInt(nvariables(equations))))

  @turbo for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
    rho    = u_cons[1, i, j, k, element]
    rho_v1 = u_cons[2, i, j, k, element]
    rho_v2 = u_cons[3, i, j, k, element]
    rho_v3 = u_cons[4, i, j, k, element]
    rho_e  = u_cons[5, i, j, k, element]

    v1 = rho_v1 / rho
    v2 = rho_v2 / rho
    v3 = rho_v3 / rho
    p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3))

    u_prim[i, j, k, 1] = rho
    u_prim[i, j, k, 2] = v1
    u_prim[i, j, k, 3] = v2
    u_prim[i, j, k, 4] = v3
    u_prim[i, j, k, 5] = p
  end

  # x direction
  # At first, we create new temporary arrays with permuted memory layout to
  # allow using SIMD instructions along the first dimension (which is contiguous
  # in memory).
  du_permuted = StrideArray{eltype(u_cons)}(undef,
    (StaticInt(nnodes(dg)^2), StaticInt(nnodes(dg)),
     StaticInt(nvariables(equations))))

  u_prim_permuted = StrideArray{eltype(u_cons)}(undef,
    (StaticInt(nnodes(dg)^2), StaticInt(nnodes(dg)),
     StaticInt(nvariables(equations))))

  @turbo for v in eachvariable(equations),
             k in eachnode(dg),
             j in eachnode(dg),
             i in eachnode(dg)
    jk = j + nnodes(dg) * (k- 1)
    u_prim_permuted[jk, i, v] = u_prim[i, j, k, v]
  end
  fill!(du_permuted, zero(eltype(du_permuted)))

  # We must also permute the contravariant vectors.
  contravariant_vectors_x = StrideArray{eltype(contravariant_vectors)}(undef,
    (StaticInt(nnodes(dg)^2), StaticInt(nnodes(dg)),
     StaticInt(ndims(mesh))))

  @turbo for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
    jk = j + nnodes(dg) * (k- 1)
    contravariant_vectors_x[jk, i, 1] = contravariant_vectors[1, 1, i, j, k, element]
    contravariant_vectors_x[jk, i, 2] = contravariant_vectors[2, 1, i, j, k, element]
    contravariant_vectors_x[jk, i, 3] = contravariant_vectors[3, 1, i, j, k, element]
  end

  # Next, we basically inline the volume flux. To allow SIMD vectorization and
  # still use the symmetry of the volume flux and the derivative matrix, we
  # loop over the triangular part in an outer loop and use a plain inner loop.
  for i in eachnode(dg), ii in (i+1):nnodes(dg)
    @turbo for jk in Base.OneTo(nnodes(dg)^2)
      rho_ll = u_prim_permuted[jk, i, 1]
      v1_ll  = u_prim_permuted[jk, i, 2]
      v2_ll  = u_prim_permuted[jk, i, 3]
      v3_ll  = u_prim_permuted[jk, i, 4]
      p_ll   = u_prim_permuted[jk, i, 5]

      rho_rr = u_prim_permuted[jk, ii, 1]
      v1_rr  = u_prim_permuted[jk, ii, 2]
      v2_rr  = u_prim_permuted[jk, ii, 3]
      v3_rr  = u_prim_permuted[jk, ii, 4]
      p_rr   = u_prim_permuted[jk, ii, 5]

      normal_direction_1 = 0.5 * (
        contravariant_vectors_x[jk, i, 1] + contravariant_vectors_x[jk, ii, 1])
      normal_direction_2 = 0.5 * (
        contravariant_vectors_x[jk, i, 2] + contravariant_vectors_x[jk, ii, 2])
      normal_direction_3 = 0.5 * (
        contravariant_vectors_x[jk, i, 3] + contravariant_vectors_x[jk, ii, 3])

      v_dot_n_ll = v1_ll * normal_direction_1 + v2_ll * normal_direction_2 + v3_ll * normal_direction_3
      v_dot_n_rr = v1_rr * normal_direction_1 + v2_rr * normal_direction_2 + v3_rr * normal_direction_3

      # Compute required mean values
      rho_avg = 0.5 * (rho_ll + rho_rr)
      v1_avg  = 0.5 * ( v1_ll +  v1_rr)
      v2_avg  = 0.5 * ( v2_ll +  v2_rr)
      v3_avg  = 0.5 * ( v3_ll +  v3_rr)
      v_dot_n_avg = 0.5 * (v_dot_n_ll + v_dot_n_rr)
      p_avg   = 0.5 * (  p_ll +   p_rr)
      velocity_square_avg = 0.5 * (v1_ll * v1_rr + v2_ll * v2_rr + v3_ll * v3_rr)

      # Calculate fluxes depending on normal_direction
      f1 = rho_avg * v_dot_n_avg
      f2 = f1 * v1_avg + p_avg * normal_direction_1
      f3 = f1 * v2_avg + p_avg * normal_direction_2
      f4 = f1 * v3_avg + p_avg * normal_direction_3
      f5 = ( f1 * velocity_square_avg + p_avg * v_dot_n_avg * equations.inv_gamma_minus_one
            + 0.5 * (p_ll * v_dot_n_rr + p_rr * v_dot_n_ll) )

      # Add scaled fluxes to RHS
      factor_i = alpha * derivative_split[i, ii]
      du_permuted[jk, i, 1] += factor_i * f1
      du_permuted[jk, i, 2] += factor_i * f2
      du_permuted[jk, i, 3] += factor_i * f3
      du_permuted[jk, i, 4] += factor_i * f4
      du_permuted[jk, i, 5] += factor_i * f5

      factor_ii = alpha * derivative_split[ii, i]
      du_permuted[jk, ii, 1] += factor_ii * f1
      du_permuted[jk, ii, 2] += factor_ii * f2
      du_permuted[jk, ii, 3] += factor_ii * f3
      du_permuted[jk, ii, 4] += factor_ii * f4
      du_permuted[jk, ii, 5] += factor_ii * f5
    end
  end

  @turbo for v in eachvariable(equations),
             k in eachnode(dg),
             j in eachnode(dg),
             i in eachnode(dg)
    jk = j + nnodes(dg) * (k- 1)
    du[i, j, k, v] = du_permuted[jk, i, v]
  end

  # y direction
  # We must also permute the contravariant vectors.
  contravariant_vectors_y = StrideArray{eltype(contravariant_vectors)}(undef,
    (StaticInt(nnodes(dg)), StaticInt(nnodes(dg)), StaticInt(nnodes(dg)),
     StaticInt(ndims(mesh))))

  @turbo for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
    contravariant_vectors_y[i, j, k, 1] = contravariant_vectors[1, 2, i, j, k, element]
    contravariant_vectors_y[i, j, k, 2] = contravariant_vectors[2, 2, i, j, k, element]
    contravariant_vectors_y[i, j, k, 3] = contravariant_vectors[3, 2, i, j, k, element]
  end

  # A possible permutation of array dimensions with improved opportunities for
  # SIMD vectorization appeared to be slower than the direct version used here
  # in preliminary numerical experiments on an AVX2 system.
  for j in eachnode(dg), jj in (j+1):nnodes(dg)
    @turbo for k in eachnode(dg), i in eachnode(dg)
      rho_ll = u_prim[i, j, k, 1]
      v1_ll  = u_prim[i, j, k, 2]
      v2_ll  = u_prim[i, j, k, 3]
      v3_ll  = u_prim[i, j, k, 4]
      p_ll   = u_prim[i, j, k, 5]

      rho_rr = u_prim[i, jj, k, 1]
      v1_rr  = u_prim[i, jj, k, 2]
      v2_rr  = u_prim[i, jj, k, 3]
      v3_rr  = u_prim[i, jj, k, 4]
      p_rr   = u_prim[i, jj, k, 5]

      normal_direction_1 = 0.5 * (
        contravariant_vectors_y[i, j, k, 1] + contravariant_vectors_y[i, jj, k, 1])
      normal_direction_2 = 0.5 * (
        contravariant_vectors_y[i, j, k, 2] + contravariant_vectors_y[i, jj, k, 2])
      normal_direction_3 = 0.5 * (
        contravariant_vectors_y[i, j, k, 3] + contravariant_vectors_y[i, jj, k, 3])

      v_dot_n_ll = v1_ll * normal_direction_1 + v2_ll * normal_direction_2 + v3_ll * normal_direction_3
      v_dot_n_rr = v1_rr * normal_direction_1 + v2_rr * normal_direction_2 + v3_rr * normal_direction_3

      # Compute required mean values
      rho_avg = 0.5 * (rho_ll + rho_rr)
      v1_avg  = 0.5 * ( v1_ll +  v1_rr)
      v2_avg  = 0.5 * ( v2_ll +  v2_rr)
      v3_avg  = 0.5 * ( v3_ll +  v3_rr)
      v_dot_n_avg = 0.5 * (v_dot_n_ll + v_dot_n_rr)
      p_avg   = 0.5 * (  p_ll +   p_rr)
      velocity_square_avg = 0.5 * (v1_ll * v1_rr + v2_ll * v2_rr + v3_ll * v3_rr)

      # Calculate fluxes depending on normal_direction
      f1 = rho_avg * v_dot_n_avg
      f2 = f1 * v1_avg + p_avg * normal_direction_1
      f3 = f1 * v2_avg + p_avg * normal_direction_2
      f4 = f1 * v3_avg + p_avg * normal_direction_3
      f5 = ( f1 * velocity_square_avg + p_avg * v_dot_n_avg * equations.inv_gamma_minus_one
            + 0.5 * (p_ll * v_dot_n_rr + p_rr * v_dot_n_ll) )

      # Add scaled fluxes to RHS
      factor_j = alpha * derivative_split[j, jj]
      du[i, j, k, 1] += factor_j * f1
      du[i, j, k, 2] += factor_j * f2
      du[i, j, k, 3] += factor_j * f3
      du[i, j, k, 4] += factor_j * f4
      du[i, j, k, 5] += factor_j * f5

      factor_jj = alpha * derivative_split[jj, j]
      du[i, jj, k, 1] += factor_jj * f1
      du[i, jj, k, 2] += factor_jj * f2
      du[i, jj, k, 3] += factor_jj * f3
      du[i, jj, k, 4] += factor_jj * f4
      du[i, jj, k, 5] += factor_jj * f5
    end
  end

  # z direction
  # The memory layout is already optimal for SIMD vectorization in this loop.
  # We just squeeze the first two dimensions to make the code slightly faster.
  GC.@preserve u_prim begin
    u_prim_reshaped = PtrArray(pointer(u_prim),
      (StaticInt(nnodes(dg)^2), StaticInt(nnodes(dg)),
       StaticInt(nvariables(equations))))

    du_reshaped = PtrArray(pointer(du),
      (StaticInt(nnodes(dg)^2), StaticInt(nnodes(dg)),
       StaticInt(nvariables(equations))))

    # We must also permute the contravariant vectors.
    contravariant_vectors_z = StrideArray{eltype(contravariant_vectors)}(undef,
      (StaticInt(nnodes(dg)^2), StaticInt(nnodes(dg)),
       StaticInt(ndims(mesh))))

    @turbo for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
      ij = i + nnodes(dg) * (j- 1)
      contravariant_vectors_z[ij, k, 1] = contravariant_vectors[1, 3, i, j, k, element]
      contravariant_vectors_z[ij, k, 2] = contravariant_vectors[2, 3, i, j, k, element]
      contravariant_vectors_z[ij, k, 3] = contravariant_vectors[3, 3, i, j, k, element]
    end

    for k in eachnode(dg), kk in (k+1):nnodes(dg)
      @turbo for ij in Base.OneTo(nnodes(dg)^2)
        rho_ll = u_prim_reshaped[ij, k, 1]
        v1_ll  = u_prim_reshaped[ij, k, 2]
        v2_ll  = u_prim_reshaped[ij, k, 3]
        v3_ll  = u_prim_reshaped[ij, k, 4]
        p_ll   = u_prim_reshaped[ij, k, 5]

        rho_rr = u_prim_reshaped[ij, kk, 1]
        v1_rr  = u_prim_reshaped[ij, kk, 2]
        v2_rr  = u_prim_reshaped[ij, kk, 3]
        v3_rr  = u_prim_reshaped[ij, kk, 4]
        p_rr   = u_prim_reshaped[ij, kk, 5]

        normal_direction_1 = 0.5 * (
          contravariant_vectors_z[ij, k, 1] + contravariant_vectors_z[ij, kk, 1])
        normal_direction_2 = 0.5 * (
          contravariant_vectors_z[ij, k, 2] + contravariant_vectors_z[ij, kk, 2])
        normal_direction_3 = 0.5 * (
          contravariant_vectors_z[ij, k, 3] + contravariant_vectors_z[ij, kk, 3])

        v_dot_n_ll = v1_ll * normal_direction_1 + v2_ll * normal_direction_2 + v3_ll * normal_direction_3
        v_dot_n_rr = v1_rr * normal_direction_1 + v2_rr * normal_direction_2 + v3_rr * normal_direction_3

        # Compute required mean values
        rho_avg = 0.5 * (rho_ll + rho_rr)
        v1_avg  = 0.5 * ( v1_ll +  v1_rr)
        v2_avg  = 0.5 * ( v2_ll +  v2_rr)
        v3_avg  = 0.5 * ( v3_ll +  v3_rr)
        v_dot_n_avg = 0.5 * (v_dot_n_ll + v_dot_n_rr)
        p_avg   = 0.5 * (  p_ll +   p_rr)
        velocity_square_avg = 0.5 * (v1_ll * v1_rr + v2_ll * v2_rr + v3_ll * v3_rr)

        # Calculate fluxes depending on normal_direction
        f1 = rho_avg * v_dot_n_avg
        f2 = f1 * v1_avg + p_avg * normal_direction_1
        f3 = f1 * v2_avg + p_avg * normal_direction_2
        f4 = f1 * v3_avg + p_avg * normal_direction_3
        f5 = ( f1 * velocity_square_avg + p_avg * v_dot_n_avg * equations.inv_gamma_minus_one
              + 0.5 * (p_ll * v_dot_n_rr + p_rr * v_dot_n_ll) )

        # Add scaled fluxes to RHS
        factor_k = alpha * derivative_split[k, kk]
        du_reshaped[ij, k, 1] += factor_k * f1
        du_reshaped[ij, k, 2] += factor_k * f2
        du_reshaped[ij, k, 3] += factor_k * f3
        du_reshaped[ij, k, 4] += factor_k * f4
        du_reshaped[ij, k, 5] += factor_k * f5

        factor_kk = alpha * derivative_split[kk, k]
        du_reshaped[ij, kk, 1] += factor_kk * f1
        du_reshaped[ij, kk, 2] += factor_kk * f2
        du_reshaped[ij, kk, 3] += factor_kk * f3
        du_reshaped[ij, kk, 4] += factor_kk * f4
        du_reshaped[ij, kk, 5] += factor_kk * f5
      end
    end
  end # GC.@preserve u_prim begin

  # Finally, we add the temporary RHS computed here to the global RHS in the
  # given `element`.
  @turbo for v in eachvariable(equations),
             k in eachnode(dg),
             j in eachnode(dg),
             i in eachnode(dg)
    _du[v, i, j, k, element] += du[i, j, k, v]
  end
end

@inline function split_form_kernel!(_du::PtrArray, u_cons::PtrArray,
                                    element, mesh::Union{StructuredMesh{3}, P4estMesh{3}},
                                    nonconservative_terms::Val{false},
                                    equations::CompressibleEulerEquations3D,
                                    volume_flux::typeof(flux_ranocha_turbo),
                                    dg::DGSEM, cache, alpha)
  @unpack derivative_split = dg.basis
  @unpack contravariant_vectors = cache.elements

  # Create a temporary array that will be used to store the RHS with permuted
  # indices `[i, j, k, v]` to allow using SIMD instructions.
  # `StrideArray`s with purely static dimensions do not allocate on the heap.
  du = StrideArray{eltype(u_cons)}(undef,
    (ntuple(_ -> StaticInt(nnodes(dg)), ndims(mesh))...,
     StaticInt(nvariables(equations))))

  # Convert conserved to primitive variables on the given `element`. In addition
  # to the usual primitive variables, we also compute logarithms of the density
  # and pressure to increase the performance of the required logarithmic mean
  # values.
  u_prim = StrideArray{eltype(u_cons)}(undef,
    (ntuple(_ -> StaticInt(nnodes(dg)), ndims(mesh))...,
     StaticInt(nvariables(equations) + 2))) # We also compute "+ 2" logs

  @turbo for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
    rho    = u_cons[1, i, j, k, element]
    rho_v1 = u_cons[2, i, j, k, element]
    rho_v2 = u_cons[3, i, j, k, element]
    rho_v3 = u_cons[4, i, j, k, element]
    rho_e  = u_cons[5, i, j, k, element]

    v1 = rho_v1 / rho
    v2 = rho_v2 / rho
    v3 = rho_v3 / rho
    p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3))

    u_prim[i, j, k, 1] = rho
    u_prim[i, j, k, 2] = v1
    u_prim[i, j, k, 3] = v2
    u_prim[i, j, k, 4] = v3
    u_prim[i, j, k, 5] = p
    u_prim[i, j, k, 6] = log(rho)
    u_prim[i, j, k, 7] = log(p)
  end

  # x direction
  # At first, we create new temporary arrays with permuted memory layout to
  # allow using SIMD instructions along the first dimension (which is contiguous
  # in memory).
  du_permuted = StrideArray{eltype(u_cons)}(undef,
    (StaticInt(nnodes(dg)^2), StaticInt(nnodes(dg)),
     StaticInt(nvariables(equations))))

  u_prim_permuted = StrideArray{eltype(u_cons)}(undef,
    (StaticInt(nnodes(dg)^2), StaticInt(nnodes(dg)),
     StaticInt(nvariables(equations) + 2)))

  @turbo for v in indices(u_prim, 4), # v in eachvariable(equations) misses +2 logs
             k in eachnode(dg),
             j in eachnode(dg),
             i in eachnode(dg)
    jk = j + nnodes(dg) * (k- 1)
    u_prim_permuted[jk, i, v] = u_prim[i, j, k, v]
  end
  fill!(du_permuted, zero(eltype(du_permuted)))

  # We must also permute the contravariant vectors.
  contravariant_vectors_x = StrideArray{eltype(contravariant_vectors)}(undef,
    (StaticInt(nnodes(dg)^2), StaticInt(nnodes(dg)),
     StaticInt(ndims(mesh))))

  @turbo for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
    jk = j + nnodes(dg) * (k- 1)
    contravariant_vectors_x[jk, i, 1] = contravariant_vectors[1, 1, i, j, k, element]
    contravariant_vectors_x[jk, i, 2] = contravariant_vectors[2, 1, i, j, k, element]
    contravariant_vectors_x[jk, i, 3] = contravariant_vectors[3, 1, i, j, k, element]
  end

  # Next, we basically inline the volume flux. To allow SIMD vectorization and
  # still use the symmetry of the volume flux and the derivative matrix, we
  # loop over the triangular part in an outer loop and use a plain inner loop.
  for i in eachnode(dg), ii in (i+1):nnodes(dg)
    @turbo for jk in Base.OneTo(nnodes(dg)^2)
      rho_ll     = u_prim_permuted[jk, i, 1]
      v1_ll      = u_prim_permuted[jk, i, 2]
      v2_ll      = u_prim_permuted[jk, i, 3]
      v3_ll      = u_prim_permuted[jk, i, 4]
      p_ll       = u_prim_permuted[jk, i, 5]
      log_rho_ll = u_prim_permuted[jk, i, 6]
      log_p_ll   = u_prim_permuted[jk, i, 7]

      rho_rr     = u_prim_permuted[jk, ii, 1]
      v1_rr      = u_prim_permuted[jk, ii, 2]
      v2_rr      = u_prim_permuted[jk, ii, 3]
      v3_rr      = u_prim_permuted[jk, ii, 4]
      p_rr       = u_prim_permuted[jk, ii, 5]
      log_rho_rr = u_prim_permuted[jk, ii, 6]
      log_p_rr   = u_prim_permuted[jk, ii, 7]

      normal_direction_1 = 0.5 * (
        contravariant_vectors_x[jk, i, 1] + contravariant_vectors_x[jk, ii, 1])
      normal_direction_2 = 0.5 * (
        contravariant_vectors_x[jk, i, 2] + contravariant_vectors_x[jk, ii, 2])
      normal_direction_3 = 0.5 * (
        contravariant_vectors_x[jk, i, 3] + contravariant_vectors_x[jk, ii, 3])

      v_dot_n_ll = v1_ll * normal_direction_1 + v2_ll * normal_direction_2 + v3_ll * normal_direction_3
      v_dot_n_rr = v1_rr * normal_direction_1 + v2_rr * normal_direction_2 + v3_rr * normal_direction_3

      # Compute required mean values
      # We inline the logarithmic mean to allow LoopVectorization.jl to optimize
      # it efficiently. This is equivalent to
      #   rho_mean = ln_mean(rho_ll, rho_rr)
      x1 = rho_ll
      log_x1 = log_rho_ll
      y1 = rho_rr
      log_y1 = log_rho_rr
      x1_plus_y1 = x1 + y1
      y1_minus_x1 = y1 - x1
      z1 = y1_minus_x1^2 / x1_plus_y1^2
      special_path1 = x1_plus_y1 / (2 + z1*(2/3 + z1*(2/5 + 2/7*z1)))
      regular_path1 = y1_minus_x1 / (log_y1 - log_x1)
      rho_mean = ifelse(z1 < 1.0e-4, special_path1, regular_path1)

      # Algebraically equivalent to `inv_ln_mean(rho_ll / p_ll, rho_rr / p_rr)`
      # in exact arithmetic since
      #     log((ϱₗ/pₗ) / (ϱᵣ/pᵣ)) / (ϱₗ/pₗ - ϱᵣ/pᵣ)
      #   = pₗ pᵣ log((ϱₗ pᵣ) / (ϱᵣ pₗ)) / (ϱₗ pᵣ - ϱᵣ pₗ)
      # inv_rho_p_mean = p_ll * p_rr * inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll)
      x2 = rho_ll * p_rr
      log_x2 = log_rho_ll + log_p_rr
      y2 = rho_rr * p_ll
      log_y2 = log_rho_rr + log_p_ll
      x2_plus_y2 = x2 + y2
      y2_minus_x2 = y2 - x2
      z2 = y2_minus_x2^2 / x2_plus_y2^2
      special_path2 = (2 + z2*(2/3 + z2*(2/5 + 2/7*z2))) / x2_plus_y2
      regular_path2 = (log_y2 - log_x2) / y2_minus_x2
      inv_rho_p_mean = p_ll * p_rr * ifelse(z2 < 1.0e-4, special_path2, regular_path2)

      v1_avg = 0.5 * (v1_ll + v1_rr)
      v2_avg = 0.5 * (v2_ll + v2_rr)
      v3_avg = 0.5 * (v3_ll + v3_rr)
      p_avg  = 0.5 * (p_ll + p_rr)
      velocity_square_avg = 0.5 * (v1_ll*v1_rr + v2_ll*v2_rr + v3_ll*v3_rr)

      # Calculate fluxes depending on normal_direction
      f1 = rho_mean * 0.5 * (v_dot_n_ll + v_dot_n_rr)
      f2 = f1 * v1_avg + p_avg * normal_direction_1
      f3 = f1 * v2_avg + p_avg * normal_direction_2
      f4 = f1 * v3_avg + p_avg * normal_direction_3
      f5 = ( f1 * ( velocity_square_avg + inv_rho_p_mean * equations.inv_gamma_minus_one )
          + 0.5 * (p_ll * v_dot_n_rr + p_rr * v_dot_n_ll) )

      # Add scaled fluxes to RHS
      factor_i = alpha * derivative_split[i, ii]
      du_permuted[jk, i, 1] += factor_i * f1
      du_permuted[jk, i, 2] += factor_i * f2
      du_permuted[jk, i, 3] += factor_i * f3
      du_permuted[jk, i, 4] += factor_i * f4
      du_permuted[jk, i, 5] += factor_i * f5

      factor_ii = alpha * derivative_split[ii, i]
      du_permuted[jk, ii, 1] += factor_ii * f1
      du_permuted[jk, ii, 2] += factor_ii * f2
      du_permuted[jk, ii, 3] += factor_ii * f3
      du_permuted[jk, ii, 4] += factor_ii * f4
      du_permuted[jk, ii, 5] += factor_ii * f5
    end
  end

  @turbo for v in eachvariable(equations),
             k in eachnode(dg),
             j in eachnode(dg),
             i in eachnode(dg)
    jk = j + nnodes(dg) * (k- 1)
    du[i, j, k, v] = du_permuted[jk, i, v]
  end

  # y direction
  # We must also permute the contravariant vectors.
  contravariant_vectors_y = StrideArray{eltype(contravariant_vectors)}(undef,
    (StaticInt(nnodes(dg)), StaticInt(nnodes(dg)), StaticInt(nnodes(dg)),
     StaticInt(ndims(mesh))))

  @turbo for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
    contravariant_vectors_y[i, j, k, 1] = contravariant_vectors[1, 2, i, j, k, element]
    contravariant_vectors_y[i, j, k, 2] = contravariant_vectors[2, 2, i, j, k, element]
    contravariant_vectors_y[i, j, k, 3] = contravariant_vectors[3, 2, i, j, k, element]
  end

  # A possible permutation of array dimensions with improved opportunities for
  # SIMD vectorization appeared to be slower than the direct version used here
  # in preliminary numerical experiments on an AVX2 system.
  for j in eachnode(dg), jj in (j+1):nnodes(dg)
    @turbo for k in eachnode(dg), i in eachnode(dg)
      rho_ll     = u_prim[i, j, k, 1]
      v1_ll      = u_prim[i, j, k, 2]
      v2_ll      = u_prim[i, j, k, 3]
      v3_ll      = u_prim[i, j, k, 4]
      p_ll       = u_prim[i, j, k, 5]
      log_rho_ll = u_prim[i, j, k, 6]
      log_p_ll   = u_prim[i, j, k, 7]

      rho_rr     = u_prim[i, jj, k, 1]
      v1_rr      = u_prim[i, jj, k, 2]
      v2_rr      = u_prim[i, jj, k, 3]
      v3_rr      = u_prim[i, jj, k, 4]
      p_rr       = u_prim[i, jj, k, 5]
      log_rho_rr = u_prim[i, jj, k, 6]
      log_p_rr   = u_prim[i, jj, k, 7]

      normal_direction_1 = 0.5 * (
        contravariant_vectors_y[i, j, k, 1] + contravariant_vectors_y[i, jj, k, 1])
      normal_direction_2 = 0.5 * (
        contravariant_vectors_y[i, j, k, 2] + contravariant_vectors_y[i, jj, k, 2])
      normal_direction_3 = 0.5 * (
        contravariant_vectors_y[i, j, k, 3] + contravariant_vectors_y[i, jj, k, 3])

      v_dot_n_ll = v1_ll * normal_direction_1 + v2_ll * normal_direction_2 + v3_ll * normal_direction_3
      v_dot_n_rr = v1_rr * normal_direction_1 + v2_rr * normal_direction_2 + v3_rr * normal_direction_3

      # Compute required mean values
      # We inline the logarithmic mean to allow LoopVectorization.jl to optimize
      # it efficiently. This is equivalent to
      #   rho_mean = ln_mean(rho_ll, rho_rr)
      x1 = rho_ll
      log_x1 = log_rho_ll
      y1 = rho_rr
      log_y1 = log_rho_rr
      x1_plus_y1 = x1 + y1
      y1_minus_x1 = y1 - x1
      z1 = y1_minus_x1^2 / x1_plus_y1^2
      special_path1 = x1_plus_y1 / (2 + z1*(2/3 + z1*(2/5 + 2/7*z1)))
      regular_path1 = y1_minus_x1 / (log_y1 - log_x1)
      rho_mean = ifelse(z1 < 1.0e-4, special_path1, regular_path1)

      # Algebraically equivalent to `inv_ln_mean(rho_ll / p_ll, rho_rr / p_rr)`
      # in exact arithmetic since
      #     log((ϱₗ/pₗ) / (ϱᵣ/pᵣ)) / (ϱₗ/pₗ - ϱᵣ/pᵣ)
      #   = pₗ pᵣ log((ϱₗ pᵣ) / (ϱᵣ pₗ)) / (ϱₗ pᵣ - ϱᵣ pₗ)
      # inv_rho_p_mean = p_ll * p_rr * inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll)
      x2 = rho_ll * p_rr
      log_x2 = log_rho_ll + log_p_rr
      y2 = rho_rr * p_ll
      log_y2 = log_rho_rr + log_p_ll
      x2_plus_y2 = x2 + y2
      y2_minus_x2 = y2 - x2
      z2 = y2_minus_x2^2 / x2_plus_y2^2
      special_path2 = (2 + z2*(2/3 + z2*(2/5 + 2/7*z2))) / x2_plus_y2
      regular_path2 = (log_y2 - log_x2) / y2_minus_x2
      inv_rho_p_mean = p_ll * p_rr * ifelse(z2 < 1.0e-4, special_path2, regular_path2)

      v1_avg = 0.5 * (v1_ll + v1_rr)
      v2_avg = 0.5 * (v2_ll + v2_rr)
      v3_avg = 0.5 * (v3_ll + v3_rr)
      p_avg  = 0.5 * (p_ll + p_rr)
      velocity_square_avg = 0.5 * (v1_ll*v1_rr + v2_ll*v2_rr + v3_ll*v3_rr)

      # Calculate fluxes depending on normal_direction
      f1 = rho_mean * 0.5 * (v_dot_n_ll + v_dot_n_rr)
      f2 = f1 * v1_avg + p_avg * normal_direction_1
      f3 = f1 * v2_avg + p_avg * normal_direction_2
      f4 = f1 * v3_avg + p_avg * normal_direction_3
      f5 = ( f1 * ( velocity_square_avg + inv_rho_p_mean * equations.inv_gamma_minus_one )
          + 0.5 * (p_ll * v_dot_n_rr + p_rr * v_dot_n_ll) )

      # Add scaled fluxes to RHS
      factor_j = alpha * derivative_split[j, jj]
      du[i, j, k, 1] += factor_j * f1
      du[i, j, k, 2] += factor_j * f2
      du[i, j, k, 3] += factor_j * f3
      du[i, j, k, 4] += factor_j * f4
      du[i, j, k, 5] += factor_j * f5

      factor_jj = alpha * derivative_split[jj, j]
      du[i, jj, k, 1] += factor_jj * f1
      du[i, jj, k, 2] += factor_jj * f2
      du[i, jj, k, 3] += factor_jj * f3
      du[i, jj, k, 4] += factor_jj * f4
      du[i, jj, k, 5] += factor_jj * f5
    end
  end

  # z direction
  # The memory layout is already optimal for SIMD vectorization in this loop.
  # We just squeeze the first two dimensions to make the code slightly faster.
  GC.@preserve u_prim begin
    u_prim_reshaped = PtrArray(pointer(u_prim),
      (StaticInt(nnodes(dg)^2), StaticInt(nnodes(dg)),
       StaticInt(nvariables(equations) + 2)))

    du_reshaped = PtrArray(pointer(du),
      (StaticInt(nnodes(dg)^2), StaticInt(nnodes(dg)),
       StaticInt(nvariables(equations))))

    # We must also permute the contravariant vectors.
    contravariant_vectors_z = StrideArray{eltype(contravariant_vectors)}(undef,
      (StaticInt(nnodes(dg)^2), StaticInt(nnodes(dg)),
       StaticInt(ndims(mesh))))

    @turbo for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
      ij = i + nnodes(dg) * (j- 1)
      contravariant_vectors_z[ij, k, 1] = contravariant_vectors[1, 3, i, j, k, element]
      contravariant_vectors_z[ij, k, 2] = contravariant_vectors[2, 3, i, j, k, element]
      contravariant_vectors_z[ij, k, 3] = contravariant_vectors[3, 3, i, j, k, element]
    end

    for k in eachnode(dg), kk in (k+1):nnodes(dg)
      @turbo for ij in Base.OneTo(nnodes(dg)^2)
        rho_ll     = u_prim_reshaped[ij, k, 1]
        v1_ll      = u_prim_reshaped[ij, k, 2]
        v2_ll      = u_prim_reshaped[ij, k, 3]
        v3_ll      = u_prim_reshaped[ij, k, 4]
        p_ll       = u_prim_reshaped[ij, k, 5]
        log_rho_ll = u_prim_reshaped[ij, k, 6]
        log_p_ll   = u_prim_reshaped[ij, k, 7]

        rho_rr     = u_prim_reshaped[ij, kk, 1]
        v1_rr      = u_prim_reshaped[ij, kk, 2]
        v2_rr      = u_prim_reshaped[ij, kk, 3]
        v3_rr      = u_prim_reshaped[ij, kk, 4]
        p_rr       = u_prim_reshaped[ij, kk, 5]
        log_rho_rr = u_prim_reshaped[ij, kk, 6]
        log_p_rr   = u_prim_reshaped[ij, kk, 7]

        normal_direction_1 = 0.5 * (
          contravariant_vectors_z[ij, k, 1] + contravariant_vectors_z[ij, kk, 1])
        normal_direction_2 = 0.5 * (
          contravariant_vectors_z[ij, k, 2] + contravariant_vectors_z[ij, kk, 2])
        normal_direction_3 = 0.5 * (
          contravariant_vectors_z[ij, k, 3] + contravariant_vectors_z[ij, kk, 3])

        v_dot_n_ll = v1_ll * normal_direction_1 + v2_ll * normal_direction_2 + v3_ll * normal_direction_3
        v_dot_n_rr = v1_rr * normal_direction_1 + v2_rr * normal_direction_2 + v3_rr * normal_direction_3

        # Compute required mean values
        # We inline the logarithmic mean to allow LoopVectorization.jl to optimize
        # it efficiently. This is equivalent to
        #   rho_mean = ln_mean(rho_ll, rho_rr)
        x1 = rho_ll
        log_x1 = log_rho_ll
        y1 = rho_rr
        log_y1 = log_rho_rr
        x1_plus_y1 = x1 + y1
        y1_minus_x1 = y1 - x1
        z1 = y1_minus_x1^2 / x1_plus_y1^2
        special_path1 = x1_plus_y1 / (2 + z1*(2/3 + z1*(2/5 + 2/7*z1)))
        regular_path1 = y1_minus_x1 / (log_y1 - log_x1)
        rho_mean = ifelse(z1 < 1.0e-4, special_path1, regular_path1)

        # Algebraically equivalent to `inv_ln_mean(rho_ll / p_ll, rho_rr / p_rr)`
        # in exact arithmetic since
        #     log((ϱₗ/pₗ) / (ϱᵣ/pᵣ)) / (ϱₗ/pₗ - ϱᵣ/pᵣ)
        #   = pₗ pᵣ log((ϱₗ pᵣ) / (ϱᵣ pₗ)) / (ϱₗ pᵣ - ϱᵣ pₗ)
        # inv_rho_p_mean = p_ll * p_rr * inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll)
        x2 = rho_ll * p_rr
        log_x2 = log_rho_ll + log_p_rr
        y2 = rho_rr * p_ll
        log_y2 = log_rho_rr + log_p_ll
        x2_plus_y2 = x2 + y2
        y2_minus_x2 = y2 - x2
        z2 = y2_minus_x2^2 / x2_plus_y2^2
        special_path2 = (2 + z2*(2/3 + z2*(2/5 + 2/7*z2))) / x2_plus_y2
        regular_path2 = (log_y2 - log_x2) / y2_minus_x2
        inv_rho_p_mean = p_ll * p_rr * ifelse(z2 < 1.0e-4, special_path2, regular_path2)

        v1_avg = 0.5 * (v1_ll + v1_rr)
        v2_avg = 0.5 * (v2_ll + v2_rr)
        v3_avg = 0.5 * (v3_ll + v3_rr)
        p_avg  = 0.5 * (p_ll + p_rr)
        velocity_square_avg = 0.5 * (v1_ll*v1_rr + v2_ll*v2_rr + v3_ll*v3_rr)

        # Calculate fluxes depending on normal_direction
        f1 = rho_mean * 0.5 * (v_dot_n_ll + v_dot_n_rr)
        f2 = f1 * v1_avg + p_avg * normal_direction_1
        f3 = f1 * v2_avg + p_avg * normal_direction_2
        f4 = f1 * v3_avg + p_avg * normal_direction_3
        f5 = ( f1 * ( velocity_square_avg + inv_rho_p_mean * equations.inv_gamma_minus_one )
            + 0.5 * (p_ll * v_dot_n_rr + p_rr * v_dot_n_ll) )

        # Add scaled fluxes to RHS
        factor_k = alpha * derivative_split[k, kk]
        du_reshaped[ij, k, 1] += factor_k * f1
        du_reshaped[ij, k, 2] += factor_k * f2
        du_reshaped[ij, k, 3] += factor_k * f3
        du_reshaped[ij, k, 4] += factor_k * f4
        du_reshaped[ij, k, 5] += factor_k * f5

        factor_kk = alpha * derivative_split[kk, k]
        du_reshaped[ij, kk, 1] += factor_kk * f1
        du_reshaped[ij, kk, 2] += factor_kk * f2
        du_reshaped[ij, kk, 3] += factor_kk * f3
        du_reshaped[ij, kk, 4] += factor_kk * f4
        du_reshaped[ij, kk, 5] += factor_kk * f5
      end
    end
  end # GC.@preserve u_prim begin

  # Finally, we add the temporary RHS computed here to the global RHS in the
  # given `element`.
  @turbo for v in eachvariable(equations),
             k in eachnode(dg),
             j in eachnode(dg),
             i in eachnode(dg)
    _du[v, i, j, k, element] += du[i, j, k, v]
  end
end


text:
"""
MyAdaGrad implementation.
See: http://jmlr.org/papers/v12/duchi11a.html
"""
type MyAdaGrad
  alpha::Float64
  states::ObjectIdDict
end

MyAdaGrad(alpha::Float64) = MyAdaGrad(alpha, ObjectIdDict())

@compat function (opt::MyAdaGrad){T}(value::Array{T}, grad::Array{T})
  state = get!(opt.states, value, nothing)
  if state == nothing
    sqgrad = zeros(T, length(value))
    opt.states[value] = sqgrad
  else
    sqgrad = state::Array{T}
  end
  @inbounds @simd for i = 1:length(grad)
    sqgrad[i] += grad[i] * grad[i]
    value[i] -= T(opt.alpha) * grad[i] / sqrt(sqgrad[i] + T(1e-6))
  end
  fill!(grad, T(0.0))
end

ext: jl

max_stars_repo_name: hshindo/JukaiNLP

max_issues_repo_name: hshindo/JukaiNLP

max_forks_repo_name: hshindo/JukaiNLP


text:
using SpecialFunctions, Documenter

makedocs(modules=[SpecialFunctions],
         sitename="SpecialFunctions.jl",
         authors="Jeff Bezanson, Stefan Karpinski, Viral B. Shah, et al.",
         pages=["Home" => "index.md",
                "Overview" => "functions_overview.md",
                "List" => "functions_list.md"])

deploydocs(repo="github.com/JuliaMath/SpecialFunctions.jl.git")


text:
function observation_function(sns_pose::Array, obj_pose::Array)
  diff_x = obj_pose[1] - sns_pose[1]
  diff_y = obj_pose[2] - sns_pose[2]
  phi = atan(diff_y, diff_x) - sns_pose[3]
  while phi >= pi
    phi -= 2 * pi
  end
  while phi < -pi
    phi += 2 * pi
  end
  return [hypot(diff_x, diff_y), phi]
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
## Run Metropolis-Hastings using a Bayesian probit model with a Normal prior
## N(0, priorVar*I) on the vasoconstriction data

using Distributions, GeometricMCMC

include("vaso.jl")
include("probitNormalPrior.jl")

# Create Model instance
model =
  Model(nPars, data, logPrior, logLikelihood, gradLogPosterior, randPrior);

# Create instance mhOpts of Metropolis-Hastings options
mhOpts = MhOpts(55000, 5000, 0.1);

# Run Metropolis-Hastings simulation
mcmc, z = mh(model, mhOpts);

# Compute ZV-MH mean estimators based on linear polynomial
linearZvMcmc, linearCoef = linearZv(mcmc, z);

# Compute ZV-MH mean estimators based on quadratic polynomial
quadraticZvMcmc, quadraticCoef = quadraticZv(mcmc, z);


text:
# # What is Parsing?

# When you read texts, you understand what the words mean.
# Your smart brain analyzes the series of letter symbols.
# Like you move your muscles without being aware of every detail, intuitively,
# your brain understands the meaning of text after assigning syntax to each word.
#
# This post explains how a computer can read a text (parsing) with
# a simple example, reading a name in a letter, a list of students, or citizens.
# The post also explains how to program a [`CombinedParser`](@ref) in julia using [`@syntax`](@ref).
#
# Let's get to it:
# These examples represent how a name should be parsed into a Julia `NamedTuple` representation.
# You can write a program to parse such names and addresses
using CombinedParsers
using CombinedParsers.Regexp
# ## Names
# can be written in many ways.
# The example uses parser building blocks:
import CombinedParsers.Regexp: word, words
import CombinedParsers.Regexp: whitespace_horizontal, whitespace_maybe, newline
# A name can be written in terms of these building blocks,
@syntax name =
    Either(
        Sequence(
            :lastname => word,
            re" *, *", # comma separating preceding lastname
            :name     => words),
        Sequence(
            v -> (lastname = v[3], name= v[1]),
            words, *(" "), word)
    )

# `CombinedParsers.@syntax` defines a string macro `name""` for using the syntax parser in your Julia program:
# See also [`@syntax`](@ref)
name"Person, Called"

# `CombinedParsers.@syntax` also defines a syntax parser function `name()` to use in your Julia program:
(re"Best, *" * name)("Best, yours truly")

# ## Addresses
# On a letter, you see the
@syntax for street_address in texts
    examples = "Am Hang 19" => (street="Am Hang", no=19)
    Sequence(:street => !!Repeat(AnyChar()),
             " ",
             :no =>Numeric(Int))
end;

# In programming, regular expressions are the standard language for matching sequences of letters, a `String`.
# The street is a sequence of alphabetical symbols or space (not numbers or `AnyChar()`)
# can be written concisely as `re"Regex"`:
@syntax for street_address in texts ## todo: override in Either
    Sequence(:street => !!re"[[:alpha:] ]+", whitespace_horizontal, :no =>Numeric(Int))
end;

street_address"Allee 47"

# Defined `@syntax`'s can be combined.
whitenewline = Repeat1(Either(whitespace_horizontal,newline));
@syntax for address in texts
    Sequence(
        Optional("Adresse:",whitenewline),
        :door => street_address, whitenewline,
        :zip =>Numeric(Int), whitespace_horizontal,
        :city => !re"[[:alpha:] ]+",
        whitespace_maybe
    )
end

# A `@syntax` can be applied as a julia String macro
address"""Allee 47
80000 Augsburg
"""

# ## Person's address data
# For person entries in texts a `CombinedParser` can be intuitively written as
@syntax for person_adresses in texts
    Sequence(Optional("Name:",whitenewline),
             :person => name, Repeat1(whitenewline),
             :adresses => join(address, Repeat1(whitenewline)));
end;

person_adresses"""
Name: Gottfried Mutbürger

Adresse:
Am Hang 19
86653 Glauberg

Adresse:
Allee 47
80650 Zweistadt
"""

# This example demonstrated the terms "syntax" and "parsing".
# Here, a syntax defines the ways in which information is written.
# Parsing means the computer reads the information back from the way it was written:
# The `CombinedParser.jl` definitions create a program that parses a formal `NamedTuple` representation from texts written in that way.


text:
# Selective includes.
include("./HessianRep/HessianDense.jl")
include("./HessianRep/HessianSparse.jl")
include("./HessianRep/HessianSparseNoTril.jl")
include("./HessianRep/HessianOp.jl")

include("./Types/Types.jl")
include("./Utilities/Utilities.jl")

include("./SolveModel/SolveModelKARC.jl")
include("./SolveModel/SolveModelKTR.jl")
include("./SolveModel/SolveModelARCDiagAbs.jl")
include("./SolveModel/SolveModelARCDiag.jl")
include("./SolveModel/SolveModelARCDiag_HO.jl")
include("./SolveModel/SolveModelARCDiag_HO_vs_Nwt.jl")
include("./SolveModel/SolveModelTRDiag.jl")
include("./SolveModel/SolveModelTRDiagAbs.jl")
include("./SolveModel/SolveModelTRDiag-HO-lambda-vs-Nwt.jl")
include("./SolveModel/SolveDiag.jl")
include("./SolveModel/SolveDiagTR.jl")
include("./SolveModel/SolveModelST_TR.jl")
include("./SolveModel/SolveModelST_ARC.jl")
include("./SolveModel/cgARC.jl")
#include("./SolveModel/krylov_aux.jl")

include("./Utilities/ldlt_symm.jl")
include("./PreProcess/TParamsKARC.jl")
include("./PreProcess/TParamsKTR.jl")
include("./PreProcess/TParamsST.jl")
include("./PreProcess/PreProcessLDLt.jl")
include("./PreProcess/PreProcessSpectral.jl")
# include("./PreProcess/PreProcessMA97.jl")
# include("./PreProcess/PreProcessMA57.jl")
include("./PreProcess/PreProcessKARC.jl")
include("./PreProcess/PreProcessKTR.jl")
include("./PreProcess/PreProcessST_TR.jl")

################################################################################
# Functions needed for high order corrections
################################################################################
# include("./SolveModel/SolveModelTRDiag-HO.jl")
# include("./SolveModel/SolveModelTRDiag-HO-lambda.jl")
include("./SolveModel/high-order-correction/shamanskii.jl")
include("./SolveModel/high-order-correction/shamanskii-lambda.jl")
# include("./SolveModel/high-order-correction/shamanskii-ma57.jl")
# include("./SolveModel/high-order-correction/shamanskii-ma57-bfgs.jl")
include("./SolveModel/high-order-correction/shamanskii-bfgs.jl")
include("./SolveModel/high-order-correction/chebyshev.jl")
include("./SolveModel/high-order-correction/halley.jl")
include("./SolveModel/high-order-correction/super-halley.jl")

################################################################################
## Model, temporary, shoudln't be here
################################################################################
include("autodiff_high_order_model.jl")


text:
"""
```julia
legend_pos_from_angle(theta, xmin, xcenter, xmax, ymin, ycenter, ymax, inout)
```

Return `(x,y)` at an angle `theta` degrees from
`(xcenter,ycenter)` on a rectangle defined by (`xmin`, `xmax`, `ymin`, `ymax`).
"""
function legend_pos_from_angle(theta, xmin, xcenter, xmax, ymin, ycenter, ymax)
    (s, c) = sincosd(theta)
    x = c < 0 ? (xmin - xcenter) / c : (xmax - xcenter) / c
    y = s < 0 ? (ymin - ycenter) / s : (ymax - ycenter) / s
    A = min(x, y)
    return (xcenter + A * c, ycenter + A * s)
end

"""
Split continuous range `[-1,1]` evenly into an integer `[1,2,3]`
"""
function legend_anchor_index(x)
    x < -1 // 3 && return 1
    x < 1 // 3 && return 2
    return 3
end

"""
Turn legend argument into a (theta, :inner) or (theta, :outer) tuple.
For backends where legend position is given in normal coordinates (0,0) -- (1,1),
so :topleft exactly corresponds to (45, :inner) etc.

If `leg` is a (::Real,::Real) tuple, keep it as is.
"""
legend_angle(leg::Real) = (leg, :inner)
legend_angle(leg::Tuple{S,T}) where {S<:Real,T<:Real} = leg
legend_angle(leg::Tuple{S,Symbol}) where {S<:Real} = leg
legend_angle(leg::Symbol) = get(
    (
        topleft          = (135, :inner),
        top              = (90, :inner),
        topright         = (45, :inner),
        left             = (180, :inner),
        right            = (0, :inner),
        bottomleft       = (225, :inner),
        bottom           = (270, :inner),
        bottomright      = (315, :inner),
        outertopleft     = (135, :outer),
        outertop         = (90, :outer),
        outertopright    = (45, :outer),
        outerleft        = (180, :outer),
        outerright       = (0, :outer),
        outerbottomleft  = (225, :outer),
        outerbottom      = (270, :outer),
        outerbottomright = (315, :outer),
    ),
    leg,
    (45, :inner),
)


text:
module Selections

import Base: !, ~, -, |, &
using Base: Callable
using Base.Broadcast: broadcasted, materialize, materialize!
using Tables: columns, columntable, materializer, rowtable, schema

export select #, select_colnames, select_renames, rename, subset, transform
export col, cols, not
export if_keys, if_values, if_pairs, if_eltype, if_matches, colrange
export all_cols, other_cols, else_cols
export key_map, key_prefix, key_replace, key_suffix
export bycol, bycol!, byrow, byrow!, bytab, bytab!

include("utils.jl")
include("types/function_composition.jl")

include("types/abstract.jl")
include("types/context_selection.jl")
include("types/renaming.jl")
include("types/intermediate_containers.jl")
include("types/chaining.jl")
include("types/selection.jl")
include("types/transformation.jl")

include("resolving/input_conversions.jl")
include("resolving/chaining.jl")
include("resolving/resolving.jl")
include("resolving/select.jl")

end # module


text:
struct PSRate{P} <: PopulationState{P}
  n::Int64 # population size
  neurontype::P
  state_now::Vector{Float64}
  alloc_du::Vector{Float64}
  alloc_r::Vector{Float64}
  input::Vector{Float64}
end
function PSRate(p::NeuronType,n::Integer)
  PSRate(n,p,[zeros(Float64,n) for _ in (1:4) ]... )
end

iofunction(x,ps::PSRate) = iofunction(x,ps.neurontype)
ioinv(x,ps::PSRate) =ioinv(x,ps.neurontype)

function local_update!(t::Real,dt::Float64,ps::PSRate)
  copy!(ps.alloc_du,ps.state_now)  # u
  ps.alloc_du .-= ps.input  # (u-I)
  lmul!(-dt/ps.neurontype.τ,ps.alloc_du) # du =  dt/τ (-u+I)
  ps.state_now .+= ps.alloc_du # u_t+1 = u_t + du
  return nothing
end

# threshold-linear input-output function
struct NTReLU <: NeuronType
  τ::Float64 # time constant
  α::Float64 # gain modulation
end
@inline function iofunction(x::Float64,p::NTReLU)
  if x<=0.0
    return 0.0
  else
    return x*p.α
  end
end
@inline function ioinv(y::Float64,p::NTReLU)
  y<=0.0 && return 0.0
  return y/p.α
end

# threshold-quadratic input-output function
struct NTQuadratic <: NeuronType
  τ::Float64 # time constant
  α::Float64 # gain modulation
end
@inline function iofunction(x::Float64,p::NTQuadratic)
  if x<=0.0
    return 0.0
  else
    return x*x*p.α
  end
end
@inline function ioinv(y::Float64,p::NTQuadratic)
  y<=0.0 && return 0.0
  return sqrt(y/p.α)
end

struct NTStep <: NeuronType
  θ::Float64 # step boundary
end
@inline function iofunction(x::Float64,p::NTStep)
  return (x<p.θ ? zero(Float64) : one(Float64))
end

"""
    send_signal(conn::ConnectionStateRate)

Computes the input to postsynaptic population, given the current state of presynaptic population.
For a rate model, it applies the iofunction to the neuron potentials, gets the rate values
then multiplies rates by weights, adding the result to the input of the postsynaptic population.
"""
function forward_signal!(t::Real,dt::Real,pspost::PSRate,conn::BaseConnection,pspre::PSRate)
  # convert a state to rates r = iofun(u)
  r = pspre.alloc_r
  broadcast!(x->iofunction(x,pspre.neurontype),r,pspre.state_now)
  # multiply by weights, add to input of b .  input_b += W * r
  mul!(pspost.input,conn.weights,r,1,1)
  return nothing
end


text:
"""
Called with an integer `n > 0` to genrate a `d x n` matrix of contexts.
"""
abstract type AbstractContextSampler end

"""
Driver to manage how the bandit policy interacts with its environment.
"""
abstract type AbstractDriver end

"""
A metric that will be called as `metric(X, a, r)` after each batch outputted by a driver.
"""
abstract type AbstractMetric end

"""
A callable policy that outputs actions via `policy(X)` and is updated via [`update!`](@ref)
"""
abstract type AbstractPolicy end

"""
A callable object to ouput rewards via `sampler(X, a)`.
"""
abstract type AbstractRewardSampler end

"""
    BanditDataset(d::Int64)

Stores the trajectory of a bandit simulation.
"""
mutable struct BanditDataset
    X::Matrix{Float64}
    a::Vector{Int64}
    r::Matrix{Float64}

    function BanditDataset(d::Int64)
        if d <= 0
            throw(ArgumentError("d must be positive"))
        end
        X = Matrix{Float64}(undef, d, 0)
        a = Int64[]
        r = Matrix{Float64}(undef, 1, 0)
        return new(X, a, r)
    end
end

"""
    rm_data(data::BanditDataset, a::Int64)

Return the data associated with arm `a`.
"""
function arm_data(data::BanditDataset, a::Int64)
    if a <= 0
        throw(ArgumentError("a must be positive"))
    end
    idx = data.a .== a
    return data.X[:, idx], data.r[:, idx]
end

"""
    append_data!(data::BanditDataset, X::AbstractMatrix, a::AbstractVector{<:Int},
              r::AbstractMatrix)

Add a batch of data to the dataset.
"""
function append_data!(
    data::BanditDataset, X::AbstractMatrix, a::AbstractVector{<:Int}, r::AbstractMatrix
)
    check_regression_data(X, r)
    if size(X, 1) != size(data.X, 1)
        throw(DimensionMismatch("X does not match the dimension of the dataset"))
    end
    data.X = hcat(data.X, X)
    data.r = hcat(data.r, r)
    return data.a = vcat(data.a, a)
end


text:
"""
    @init_parallel_kernel(package, numbertype)

Initialize the package ParallelKernel, giving access to its main functionality. Creates a module `Data` in the module where `@init_parallel_kernel` is called from. The module `Data` contains the types `Data.Number`, `Data.Array` and `Data.DeviceArray` (type `?Data` *after* calling `@init_parallel_kernel` to see the full description of the module).

# Arguments
- `package::Module`: the package used for parallelization (CUDA or Threads).
- `numbertype::DataType`: the type of numbers used by @zeros, @ones and @rand and in all array types of module `Data` (e.g. Float32 or Float64). It is contained in `Data.Number` after @init_parallel_stencil.

See also: [`Data`](@ref)
"""
macro init_parallel_kernel(package, numbertype)
    check_already_initialized()
    numbertype_val = eval_arg(__module__, numbertype)
    checkargs_init(package, numbertype_val)
    esc(init_parallel_kernel(__module__, package, numbertype_val))
end
macro init_parallel_kernel(args...) @ArgumentError("wrong number of arguments."); end

function init_parallel_kernel(caller::Module, package::Symbol, numbertype::DataType; datadoc_call=:())
    if package == PKG_CUDA
        @assert CUDA.functional(true)
        data_module = Data_cuda(numbertype)
        import_cmd  = :(import CUDA)
    elseif package == PKG_THREADS
        data_module = Data_threads(numbertype)
        import_cmd  = :()
    end
    if !isdefined(caller, :Data) || (@eval(caller, isa(Data, Module)) &&  length(symbols(caller, :Data)) == 1)  # Only if the module Data does not exist in the caller or is empty, create it.
        if (datadoc_call==:()) datadoc_call = :(@doc ParallelStencil.ParallelKernel.DATA_DOC Data) end
        @eval(caller, $data_module)
        @eval(caller, $datadoc_call)
    elseif isdefined(caller, :Data) && isdefined(caller.Data, :DeviceArray)
        @warn "Module Data from previous module initialization found in caller module ($caller); module Data not created. If you are working interactively in the REPL, then you can ignore this message."
    else
        @warn "Module Data cannot be created in caller module ($caller) as there is already a user defined symbol (module/variable...) with this name. ParallelStencil is still usable but without the features of the Data module."
    end
    @eval(caller, $import_cmd)
    set_package(package)
    set_numbertype(numbertype)
    set_initialized(true)
    return nothing
end

macro is_initialized() is_initialized() end
macro get_package() get_package() end
macro get_numbertype() get_numbertype() end
let
    global is_initialized, set_initialized, set_package, get_package, set_numbertype, get_numbertype, check_initialized, check_already_initialized
    _is_initialized::Bool       = false
    package::Symbol             = PKG_NONE
    numbertype::DataType        = NUMBERTYPE_NONE
    set_initialized(flag::Bool) = (_is_initialized = flag)
    is_initialized()            = _is_initialized
    set_package(pkg::Symbol)    = (package = pkg)
    get_package()               = package
    set_numbertype(T::DataType) = (numbertype = T)
    get_numbertype()            = numbertype
    check_initialized()         = if !is_initialized() @NotInitializedError("no macro or function of the module can be called before @init_parallel_kernel.") end
    check_already_initialized() = if is_initialized() @IncoherentCallError("ParallelKernel has already been initialized.") end
end

checkargs_init(package, numbertype) = (check_package(package); check_numbertype(numbertype);)


text:
@testset "Fitting" begin
  wl = geostatsimage("WalkerLake")
  TI = reshape(wl[:Z], size(domain(wl)))[1:20,1:20]
  d = georef((z=TI,))
  γwalker = EmpiricalVariogram(d, :z, maxlag=15.)

  # variogram types to fit
  Vs = (GaussianVariogram, SphericalVariogram,
        ExponentialVariogram, MaternVariogram)

  # all fits lead to same sill
  γs = [fit(V, γwalker) for V in Vs]
  for γ in γs
    @test isapprox(sill(γ), 0.054, atol=1e-3)
  end

  # best fit is a Gaussian variogram
  γbest = fit(Variogram, γwalker)
  @test γbest isa GaussianVariogram
  @test isapprox(sill(γbest), 0.054, atol=1e-3)

  # make sure convenient methods work
  γ₁ = fit(GaussianVariogram, γwalker, h -> 1/h)
  γ₂ = fit(Variogram, γwalker, h -> 1/h)
  @test sill(γ₁) > 0
  @test sill(γ₂) > 0

  if visualtests
    plts = map(γs) do γ
      plot(γwalker, legend=false)
      plot!(γ, 0., 15.)
    end
    @test_reference "data/fitting.png" plot(plts...)
  end
end


text:
using ModelingToolkit
using Base.Test

@Var a=1.0 b
a1 = Variable(:a,1.0)
@test a1 == a

@Var begin
    a = 1.0
    b
end

@IVar t
@DVar x(t)
@DVar y(t)=sin(1)+exp(1)
@DVar z(t)
x1 = DependentVariable(:x,dependents = [t])
y1 = DependentVariable(:y, sin(1) + exp(1),dependents = [t])
z1 = DependentVariable(:z,dependents = [t])
@test x1 == x
@test y1 == y
@test z1 == z
@IVar begin
    t
    s = cos(2.5)
end
t1 = IndependentVariable(:t)
s1 = IndependentVariable(:s, cos(2.5))
@test t1 == t
@test s1 == s
@Deriv D''~t
D1 = Differential(t, 2)
@test D1 == D
@Const c=0 v=2
c1 = Constant(0)
v1 = Constant(2)
@test c1 == c
@test v1 == v


text:
function count_terminal_grabbed_boxes!(cbelt::ConveyorBelt, server::ArmAllocationServer)

    for (arm_name, arm_props) in server.agent_prop_set
        if arm_props.obj_grabbed != ""
            cbelt.deposited_boxes += 1
        end
    end
end

### Conveyor Belt functions ###

"""
Simulate the next time-step for the conveyor belt. New boxes are added by the mirror generator
in the top-level script.
"""
function propagate_belt!(cbelt::ConveyorBelt, server::ArmAllocationServer)

    diagnostic_info = ""

    to_delete = Vector{String}(undef,0)

    # First push forward boxes according to route
    for (box_name, route) in cbelt.box_routes

        # Remove first entry of route
        popfirst!(route)

        # If at end of belt, remove
        if length(route) == 0 || route[1].coordinates[1] > cbelt.belt_end
            @info box_name," was not picked up!"
            diagnostic_info = string(box_name," lost!")
            push!(to_delete, box_name)
            cbelt.lost_boxes += 1

            # Delete from last arm bookkeeping
            delete!(server.agent_task_windows, (server.agent_ordering[end], box_name))
            delete!(server.agent_prop_set[server.agent_ordering[end]].workspace_objs, box_name)

            if haskey(server.agent_task_allocation, server.agent_ordering[end]) &&
                server.agent_task_allocation[server.agent_ordering[end]].name == box_name
                delete!(server.agent_task_allocation, server.agent_ordering[end])
            end
        end
    end

    # Now delete boxes to delete
    for box_to_del in to_delete
        remove_box!(cbelt, box_to_del)
    end

    # Update time
    cbelt.current_time += 1.0

    # IMP: New boxes added in top-level script

    return diagnostic_info
end

function remove_box!(cbelt::ConveyorBelt, box_name::String)
    delete!(cbelt.box_routes, box_name)
    cbelt.num_active_boxes -= 1
end

"""
Prior to the next allocation, update the arm-object task windows.
"""
function update_time_windows!(cbelt::ConveyorBelt, server::ArmAllocationServer)

    @assert cbelt.current_time == server.current_time

    # For root arm, all objs before start are also to be considered
    for (box_name, route) in cbelt.box_routes

        first_route_tp = route[1]
        if first_route_tp.coordinates[1] <= server.agent_set[server.agent_ordering[1]].reach[1]
            push!(server.agent_prop_set[server.agent_ordering[1]].workspace_objs, box_name)
        end
    end

    for (arm_name, arm) in server.agent_set

        agent_prop_set = server.agent_prop_set[arm_name]
        interaction_events = agent_prop_set.interaction_events

        # Get rid of stale interaction_events
        # AND update the start time for ongoing interaction events
        interaction_events = Vector{InteractionEvent}(undef, 0)
        for ie in agent_prop_set.interaction_events
            if ie.timestamps[FINISH] >= cbelt.current_time+1.0
                if ie.timestamps[START] < cbelt.current_time+1.0
                    ie.timestamps[START] = cbelt.current_time+1.0
                end
                push!(interaction_events, ie)
            end
        end

        # Iterate over boxes and insert in dict if valid
        # TODO - Change to handle update to existing routes
        for (box_name, route) in cbelt.box_routes

            # First check if boxes in or out of workspace already
            first_route_tp = route[1]
            @assert first_route_tp.timestamp >= server.current_time+1.0

            if first_route_tp.coordinates[1] >= arm.reach[2]
                # Out of workspace
                if box_name in agent_prop_set.workspace_objs
                    pop!(agent_prop_set.workspace_objs, box_name)
                    delete!(server.agent_task_windows, (arm_name, box_name))
                end
            else
                # Insert in workspace_objs IF the current position is within the workspace
                if first_route_tp.coordinates[1] >= arm.reach[1] && first_route_tp.coordinates[1] < arm.reach[2]
                    push!(agent_prop_set.workspace_objs, box_name)
                end

                # Regardless, update the interaction events
                if haskey(server.agent_task_windows, (arm_name, box_name)) == false

                    # now compute time windows
                    min_time = Inf
                    max_time = -Inf

                    for tpoint in route

                        tm = tpoint.timestamp
                        coords = tpoint.coordinates

                        if coords[1] >= arm.reach[1] && coords[1] <= arm.reach[2]
                            # Within window
                            if tm < min_time
                                min_time = tm
                            end
                            if tm > max_time
                                max_time = tm
                            end
                        end
                    end

                    # Assign windows and update temporal event list
                    if min_time != Inf && max_time != -Inf

                        server.agent_task_windows[(arm_name, box_name)] = [min_time, max_time+arm.downtime/2.0, max_time + arm.downtime]
                        ie_timestamps = Dict(START=>min_time, FINISH=>max_time+arm.downtime/2.0, SUCCESS=>max_time + arm.downtime)
                        push!(interaction_events, InteractionEvent(arm_name, box_name, ie_timestamps))

                    end
                end
            end
        end

        # Re-sort the temporal event list
        sort!(interaction_events, by=x->x.timestamps[START])
        agent_prop_set.interaction_events = interaction_events

    end
end

"""
Given the starting x-coordinate of the box and belt properties (speed etc), map out the
trajectory of the box over the planning horizon.
"""
function generate_route_from_init_point(cbelt::ConveyorBelt, box_start_time::Float64, box_start_point::Float64)

    box_route = Vector{TimePoint}(undef, 0)

    while box_start_point < cbelt.belt_end

        push!(box_route, TimePoint(box_start_time, [box_start_point, 0.5]))

        box_start_point += cbelt.curr_speed
        box_start_time += 1.0
    end

    return box_route
end

function add_box_to_belt!(cbelt::ConveyorBelt, box_start_point::Float64)

    cbelt.num_active_boxes += 1
    cbelt.num_total_boxes += 1
    new_box_name = string("box",cbelt.num_total_boxes)

    box_route = generate_route_from_init_point(cbelt, cbelt.current_time+1.0, box_start_point)
    cbelt.box_routes[new_box_name] = box_route

end

"""
Get box x-coordinate that is closest to the reference time.
"""
function get_box_point(cbelt::ConveyorBelt, box_name::String, ref_time::Float64)

    box_route = cbelt.box_routes[box_name]

    for box_tp in box_route
        if abs(box_tp.timestamp - ref_time) < 0.5
            return box_tp.coordinates
        end
    end

    # Shouldn't reach this point
    return [Inf, Inf]
end

"""
Sets up the mirrored feasible task sequence generator by reflecting the
simulator about the Y-axis.
"""
function setup_mirror_gen!(cbelt::ConveyorBelt)

    # Copy over everything while mirroring the x-axes of the positions
    n_arms = length(cbelt.mirror_gen.arm_set)

    mirror_arm_set = Vector{Arm}(undef, n_arms)
    cbelt.mirror_gen.arm_stages = Vector{GENSTAGE}(undef, n_arms)

    for (i, arm) in enumerate(cbelt.mirror_gen.arm_set)

        new_arm = Arm(base_pos = [-arm.base_pos[1], arm.base_pos[2]],
                      bin_pos = [-arm.bin_pos[1], arm.bin_pos[2]],
                      reach = [-arm.reach[1], -arm.reach[2]],
                      orig_pos = [-arm.orig_pos[1], arm.orig_pos[2]],
                      grasp_model = arm.grasp_model,
                      downtime = arm.downtime
                      )

        mirror_arm_set[i] = new_arm
        cbelt.mirror_gen.arm_stages[i] = PREDROP
    end

    cbelt.mirror_gen.arm_set = mirror_arm_set
end

"""
Simulate one more timestep of the mirrored feasible sequence generator.
"""
function mirror_gen_step!(cbelt::ConveyorBelt, rng::RNG=Random.GLOBAL_RNG) where {RNG <: AbstractRNG}

    mirror_gen = cbelt.mirror_gen

    # First propagate boxes forward with speed and add to belt the new boxes
    keys_to_delete = Set{Int64}()

    for (i, pos) in mirror_gen.box_pos

        new_pos = pos + cbelt.curr_speed
        if new_pos > 0.0
            push!(keys_to_delete, i)
            add_box_to_belt!(cbelt, -pos)
        else
            mirror_gen.box_pos[i] = new_pos
        end
    end

    for kd in keys_to_delete
        delete!(mirror_gen.box_pos, kd)
    end

    # Now propagate arms and generate boxes if PREDROP
    for (i, arm) in enumerate(mirror_gen.arm_set)

        stage = mirror_gen.arm_stages[i]

        if stage == PREDROP::GENSTAGE

            # Drop box with NEW BOX PROB
            if rand(rng) <= cbelt.new_box_prob

                # Generate a new box position at random
                dist = Distributions.Uniform(arm.reach[2], arm.reach[1])
                new_pos = rand(rng, dist)

                mirror_gen.gen_box_id += 1
                mirror_gen.box_pos[mirror_gen.gen_box_id] = new_pos

                new_stage = POSTDROP::GENSTAGE

            else
                new_stage = stage
            end

        elseif stage == POSTDROP::GENSTAGE
            new_stage = ATBIN::GENSTAGE
        else
            new_stage = PREDROP
        end

        mirror_gen.arm_stages[i] = new_stage
    end
    @debug mirror_gen.arm_stages
    @debug mirror_gen.box_pos
end

"""
A simplified model of the trajectory planning and execution problem.
The arm can move from any place in its workspace to another in unit time.
"""
function propagate_arms!(server::ArmAllocationServer, cbelt::ConveyorBelt, rng::RNG=Random.GLOBAL_RNG) where {RNG <: AbstractRNG}

    diagnostic_info = ""

    # Update time
    server.current_time += 1.0

    @assert isequal(server.current_time,cbelt.current_time) "Server-Belt time mismatch!"

    # iterate through arms
    for (arm_name, arm) in server.agent_set

        arm_props = server.agent_prop_set[arm_name]

        if haskey(server.agent_task_allocation, arm_name)

            arm_assgn = server.agent_task_allocation[arm_name]

            # If box not in box routes, it is a stale entry that has already been picked up or lost - delete
            if ~(haskey(cbelt.box_routes, arm_assgn.name))
                delete!(server.agent_task_allocation, arm_name)
            else
                # Attempt a grab IF the noted time of the attempt is now
                if abs(server.current_time - arm_assgn.time) < 1.0

                    # Take arm to object position at reference point
                    arm_props.endeff_pos = get_box_point(cbelt, arm_assgn.name, arm_assgn.time)

                    @debug "$(arm_name) about to attempt $(arm_assgn.name)"
                    if rand(rng) < arm.grasp_model.bernoulli_prob  # Successful grab

                        @info arm_name, "has grabbed ",arm_assgn.name,"!"
                        arm_props.obj_grabbed = arm_assgn.name

                        # remove from arm-obj mutual assignment, workspace_objs, and windows
                        delete!(server.agent_task_allocation, arm_name)
                        delete!(server.agent_task_windows, (arm_name, arm_assgn.name))
                        delete!(arm_props.workspace_objs, arm_assgn.name)

                        # Remove box from belt
                        remove_box!(cbelt, arm_assgn.name)
                    else    # Unsuccessful grab

                        @info arm_name, "has missed ",arm_assgn.name,"!"
                        diagnostic_info = string(arm_name," misses ",arm_assgn.name,";",diagnostic_info)
                    end
                else
                    # Just go to rest pos
                    arm_props.endeff_pos = arm.orig_pos
                end # end abs()
            end # cbelt has box key
        else
            # First check if going to deposit object
            if arm_props.obj_grabbed != ""

                arm_props.endeff_pos = arm.bin_pos
                @info arm_name, " is depositing", arm_props.obj_grabbed
                arm_props.obj_grabbed = ""

                cbelt.deposited_boxes += 1

            else

                @info arm_name, "is at origin!"
                arm_props.endeff_pos = arm.orig_pos

            end
        end
    end

    return diagnostic_info
end

"""
Sets up the conveyor belt and arm allocation server with three hardcoded arms.
The other problem parameters are taken as named arguments and have defaults.
"""
function setup_conveyor_scenario(;grasp_probs::Vector{Float64}=[0.75, 0.75, 0.75],
                                 belt_speed::Float64=0.07,
                                 new_box_prob::Float64=0.75)

    # Set up the three arms
    arm1 = Arm(base_pos=[0.2, 0.75], bin_pos=[0.4, 0.75], reach=[0.06, 0.36],
               orig_pos=[0.2, 0.65], grasp_model=DiscreteGPM(grasp_probs[1]), downtime=2.0)
    arm1_props = ArmProperties(endeff_pos = arm1.orig_pos)

    arm2 = Arm(base_pos=[0.5, 0.25], bin_pos=[0.3, 0.25], reach=[0.36, 0.66],
               orig_pos=[0.5, 0.35], grasp_model=DiscreteGPM(grasp_probs[2]), downtime=2.0)
    arm2_props = ArmProperties(endeff_pos = arm2.orig_pos)

    arm3 = Arm(base_pos=[0.8, 0.75], bin_pos=[0.6, 0.75], reach=[0.66, 0.96],
               orig_pos=[0.8, 0.65], grasp_model=DiscreteGPM(grasp_probs[3]), downtime=2.0)
    arm3_props = ArmProperties(endeff_pos = arm3.orig_pos)

    arm_set = Dict("arm1"=>arm1, "arm2"=>arm2, "arm3"=>arm3)
    agent_prop_set = Dict("arm1"=>arm1_props, "arm2"=>arm2_props, "arm3"=>arm3_props)
    arm_ordering = ["arm1", "arm2", "arm3"]

    # Set up belt
    cbelt = ConveyorBelt(curr_speed=belt_speed, belt_start=0., belt_end=1., belt_width=0.2,
                      new_box_prob=new_box_prob, mirror_gen = MirrorGenerator(arm_set = [arm1, arm2, arm3]))
    cbelt.current_time = 0.0
    setup_mirror_gen!(cbelt)

    # Set up assignment server
    server = ArmAllocationServer(current_time=0.0, agent_set=arm_set,
                              agent_prop_set=agent_prop_set, agent_ordering = ["arm1", "arm2", "arm3"])

    return (server, cbelt)
end # function

function parse_conveyor_commandline()

    s = ArgParseSettings()

    @add_arg_table s begin
        "--trials", "-n"
            arg_type = Int64
            default = 100
        "--timesteps", "-t"
            arg_type = Int64
            default = 500
        "--belt_speed", "-s"
            arg_type = Float64
            default = 0.07
        "--new_box_prob", "-p"
            arg_type = Float64
            default = 0.75
        "--grasp_probs", "-g"
            arg_type = String
            default = "0.75,0.75,0.75"
        "--qpol_name", "-q"
            arg_type = String
            default = ""
        "baseline"
            help = "Which method to evaluate (dispatch/hungarian/mcts/qlearning/scoba)"
            required = true
        "out_file_name"
            required = true
    end

    return parse_args(s)
end


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

# Universal Generating Operator
"""
    solve!(ntw::MultiStateSystems.AbstractNetwork)

This function determines the universal generating function related to the output
of all users `usr` of a network `ntw`.
"""
function solve!(ntw::AbstractNetwork)
    for nn in ntws(ntw) solve_network!(nn) end
    for usr in ntw.usr if haskey(usr,:ind)
        for ind in usr[:ind]
            if ind == :GRA usr[:GRA] = GRA(usr) end
    end end end
    if get_info(ntw, :dependent_sources)
        s_ugf = ntw.props[:source_ugf]
        s_prb, s_val = [pr[end] for pr in s_ugf.prb], s_ugf.val
        for usr in ntw.usr
            n_ugf = usr[:ugf]
            n_prb, n_val = [pr[end] for pr in n_ugf.prb], n_ugf.val
            usr[:mat] = s_prb * n_prb'
            val, prb = reduce(kron(n_val,ustrip.(s_val)), kron(n_prb,s_prb))

            msr = Expr(:kw, get_msr(ntw)[1], n_val)
            usr[:std] = eval(:(solvedSTD(prob = $(n_prb), $msr)))
            usr[:ugf] = UGF(get_msr(ntw)[1],val,prb)
    end end
    for usr in ntw.usr if haskey(usr,:ind)
        for ind in usr[:ind]
            if ind == :EENS usr[:EENS] = EENS(usr) end
    end end end
end

## Probability Function
probability_function(pr::Vector) =
    vec([prod(np) for np in Iterators.product(pr...)])

## Structure Function
function user_structure_function(ntw::AbstractNetwork, u_node::Int)
    npaths, cpaths = paths(ntw, u_node)
    cpaths != [[]] || return nothing
    while length(npaths[1]) ≠ 0
        has_duplicate_paths(npaths) ? parallel_reduction!(npaths, cpaths) : ~ ;
        has_unique_sequence(npaths) ? series_reduction!(npaths, cpaths) : ~ ;
        while has_bridge_sequence(npaths) bridge_reduction!(npaths, cpaths) end
    end
    return cpaths[1][1]
end
function set_structure_function!(ntw::Network)
    for u_node in usr_nodes(ntw)
        expr = user_structure_function(ntw, u_node)
        for nu in ntw.ulib[u_node] ntw.usr[nu][:str] = :($expr) end
    end
end

### Flow
#### Parallel
par(args) = :(+($(args...)))
has_duplicate_paths(paths::Vector{<:Vector}) =
    !(length(paths) == length(unique(paths)))
parallel_combination(paths::Vector{<:Vector}) =
    [par_seg(unique([paths[np][ni] for np in 1:length(paths)]))
                                   for ni in 1:length(paths[1])]
par_seg(segment::Vector) = length(segment) > 1 ? par(segment) : segment[1] ;
function parallel_reduction!(npaths::Vector{<:Vector}, cpaths::Vector{<:Vector})
    for unique_npath in unique(npaths)
        idx = broadcast(==,Ref(unique_npath),npaths)
        if sum(idx) > 1
            # clean-up of the component path
            new_cpath = parallel_combination(cpaths[idx])
            deleteat!(cpaths,idx)
            push!(cpaths,new_cpath)
            # clean-up of the nodal path
            deleteat!(npaths,idx)
            push!(npaths,unique_npath)
    end end
end

#### Series
ser(args) = :(min($(args...)))
has_unique_sequence(npaths::Vector{<:Vector}) = !isempty(unique_sequence(npaths))
check_sequence(seq::Vector, npath::Vector) =
    length(seq) <= length(npath) && !has_subvector(npath, seq) &&
    has_subvector(npath, seq[2:end-1])
unique_elements(id::Int, idx::Vector, paths::Vector{<:Vector}) =
    setdiff(paths[id],(paths[ni] for ni in 1:length(paths) if !in(ni,idx))...)
function unique_sequence(npaths::Vector{<:Vector})
    length(npaths) > 1 || return [1]

    all_seqs, nogo_seqs = [], []
    for npath in npaths
        for ni in 1:length(npath)-3+1
            seq = npath[ni:ni+3-1]
            push!(all_seqs, seq)
            for xpath in npaths
                if check_sequence(seq, xpath) push!(nogo_seqs, seq) end
            end
        end
    end
    seqs = setdiff(all_seqs, nogo_seqs)
    !isempty(seqs) || return []
    return [np for (np,npath) in enumerate(npaths)
               if  has_subvector(npath,combine_seqs(seqs))]
end
function combine_seqs(seqs)
    while true
        new_seqs = []
        for seq in seqs for alt_seq in seqs
            if seq[2:end] == alt_seq[1:end-1]
                push!(new_seqs,vcat(seq,alt_seq[end]))
            end end
        end
        !isempty(new_seqs) || break
        seqs = new_seqs
    end
    return seqs[1]
end
function series_reduction!(npaths::Vector{<:Vector}, cpaths::Vector{<:Vector})
    path_idx = unique_sequence(npaths)
    for np in path_idx
        # clean-up nodal paths
        idx = findall(in(unique_elements(np, path_idx, npaths)), npaths[np])
        deleteat!(npaths[np], idx)
        # clean-up of the component paths
        idx = findall(in(unique_elements(np, path_idx, cpaths)), cpaths[np])
        if length(idx) > 1
            cpaths[np][idx[1]] = ser([cpaths[np][id] for id in idx])
        end
        deleteat!(cpaths[np], idx[2:end])
    end
end

#### bridge
delta(arga,argb) = :(-($arga,$argb))
function bridge(paths)
    left_path = ser(setdiff(paths[1],paths[4]))
    right_path = ser(setdiff(paths[4],paths[1]))

    left_bottom = setdiff(paths[1],paths[2])[1]
    right_bottom = setdiff(paths[4],paths[3])[1]

    temp = setdiff(paths[1],paths[4],[left_bottom])
    left_top = length(temp) == 1 ? temp[1] : ser(temp) ;
    temp = setdiff(paths[4],paths[1],[right_bottom])
    right_top = length(temp) == 1 ? temp[1] : ser(temp) ;

    bridge_comp = setdiff(paths[2],union(paths[1],paths[4]))[1]

    delta_left = delta(left_top,left_bottom)
    delta_right = delta(right_top,right_bottom)

    bridge_path = ser([:(abs($delta_left)), :(abs($delta_right)), bridge_comp])
    cond = :(1 * (_UF.ustrip($delta_left) * _UF.ustrip($delta_right) < 0))

    idx_pre = findfirst(x -> x == temp[1], paths[4])
    idx_end = findlast(x -> x == right_bottom, paths[4])

    return vcat(paths[4][1:idx_pre-1],
                par([left_path, right_path, :(*($bridge_path,$cond))]),
                par([left_bottom,right_bottom]),
                paths[4][idx_end+1:end])
end
swapped_paths(npaths) = [(ni,nj) for ni in 1:length(npaths) for nj in ni:length(npaths)
                                 if is_swapped_path(npaths[ni],npaths[nj])]
has_bridge_sequence(npaths::Vector{<:Vector}) = !isempty(bridge_sequence(npaths))
function is_swapped_path(xpath::Vector, ypath::Vector)
    length(xpath) == length(ypath) || return false
    idx = findall(x -> !iszero(x), xpath .- ypath)
    return  sum(xpath .- ypath) == 0 &&
            sum(.!iszero.(xpath .- ypath)) == 2 &&
            diff(idx)[1] in [1,-1] &&
            sort(xpath[idx]) == sort(ypath[idx])
end
function bridge_sequence(npaths::Vector{<:Vector})
    !isempty(swapped_paths(npaths)) || return []
    (ni,nj) = swapped_paths(npaths)[1]
    long_left_path, long_right_path = npaths[ni], npaths[nj]
    id = findfirst(x -> !iszero(x), npaths[ni] .- npaths[nj])
    short_left_path = long_right_path[1:end .!= id]
    short_right_path = long_left_path[1:end .!= id]
    short_left_path[1:end .!= id] == short_right_path[1:end .!= id] || return []
    selected_paths = [short_left_path, long_left_path, long_right_path, short_right_path]
    idx = [findfirst(x -> x == 1,broadcast(==,Ref(path),npaths)) for path in selected_paths]
    return idx, short_left_path[1:end .!= id]
end
function bridge_reduction!(npaths::Vector{<:Vector}, cpaths::Vector{<:Vector})
    idx, npath = bridge_sequence(npaths)
    # clean-up nodal paths
    push!(npaths, npath)
    deleteat!(npaths, sort(idx))
    # clean-up component paths
    push!(cpaths, bridge(cpaths[idx]))
    deleteat!(cpaths, sort(idx))
end

function solve_network!(ntw::AbstractNetwork)
    set_msr!(ntw)
    set_ugf!(ntw)
    set_structure_function!(ntw)

    msr = get_msr(ntw)[1]
    pr, vl, idx_itr = get_prb(ntw), get_val(ntw), get_idx_itr(ntw)
    Prb = probability_function(pr)

    skip = Int[]
    for (nu,usr) in enumerate(ntw.usr) if !in(nu,skip)
        if get_info(usr,:eval_dep)
            Val = zeros(Number, length(Prb), length(usr[:eval_dep_ids]))

            for (nc,nu) in enumerate(usr[:eval_dep_ids])
                expr = ntw.usr[nu][:str]
                exrp = quote function structure_function(idx,val) $expr end end
                eval(exrp)

                for (ni,id) in enumerate(idx_itr)
                    Val[ni,nc] = Base.invokelatest(structure_function,id,vl)
                end
            end

            rVal, rPrb = reduce(Val, Prb)

            for (nc,nu) in enumerate(usr[:eval_dep_ids])
                ntw.usr[nu][:ugf] = UGF(msr, rVal[:,nc], rPrb, rdc=false)
                ntw.usr[nu][:std] = STD(ntw.usr[nu][:ugf])
            end

            push!(skip, usr[:eval_dep_ids]...)
        else
            Val = zeros(Number, length(Prb))

            expr = usr[:str]
            exrp = quote function structure_function(idx,val) $expr end end
            eval(exrp)

            for (ni,id) in enumerate(idx_itr)
                Val[ni] = Base.invokelatest(structure_function,id,vl)
            end

            usr[:ugf] = UGF(msr, Val, Prb)                                      # Currently, only one measure is supported
            usr[:std] = STD(usr[:ugf])
    end end end
    set_info!(ntw, :solved, true)
end

# ## Structure Function
# function user_structure_function(ntw::AbstractNetwork, u_node::Int)
#     npaths, cpaths = paths(ntw, u_node)
#     cpaths != [[]] || return nothing
#     while length(npaths[1]) ≠ 0
#         has_duplicate_paths(npaths) ? vertical_reduction!(npaths,cpaths) :
#                                       horizontal_reduction!(npaths,cpaths) ;
#     end
#     return cpaths[1][1]
# end
# function set_structure_function!(ntw::Network)
#     for u_node in usr_nodes(ntw)
#         expr = user_structure_function(ntw, u_node)
#         for nu in ntw.ulib[u_node] ntw.usr[nu][:str] = :($expr) end
#     end
# end

# ### Flow
# par(args) = :(+($(args...)))
# ser(args) = :(min($(args...)))

# unique_elements(id,paths::Array) =
#     setdiff(paths[id],(paths[ni] for ni in CartesianIndices(paths) if ni≠id)...)
# has_unique_elements(id,paths::Array) = !isempty(unique_elements(id,paths))
# has_duplicate_paths(paths::Array) = !(length(paths)==length(unique(paths)))
# vertical_combination(paths::Array) =
#     [par_seg(unique([paths[np][ni] for np in CartesianIndices(paths)]))
#                                    for ni in CartesianIndices(paths[1])]
# par_seg(segment::Array) = length(segment) > 1 ? par(segment) : segment[1] ;

# function vertical_reduction!(npaths::Vector, cpaths::Vector)
#     for unique_npath in unique(npaths)
#         idx = broadcast(==,Ref(unique_npath),npaths)
#         if sum(idx) > 1
#             # clean-up of the component path
#             new_cpath = vertical_combination(cpaths[idx])
#             deleteat!(cpaths,idx)
#             push!(cpaths,new_cpath)
#             # clean-up of the nodal path
#             deleteat!(npaths,idx)
#             push!(npaths,unique_npath)
#     end end
# end
# function horizontal_reduction!(npaths::Vector, cpaths::Vector)
#     for ni in CartesianIndices(cpaths)
#         if has_unique_elements(ni,npaths)
#             # clean-up of the nodal path
#             idx = findall(in(unique_elements(ni,npaths)),npaths[ni])
#             deleteat!(npaths[ni],idx)
#             # clean-up of the component paths TODO enumerate over seperate sequences larger than one
#             idx = findall(in(unique_elements(ni,cpaths)),cpaths[ni])
#             if length(idx) > 1
#                 cpaths[ni][idx[1]] = ser([cpaths[ni][id] for id in idx])
#             end
#             deleteat!(cpaths[ni],idx[2:end])
#     end end
# end
# function unique_sequence(paths::Vector{<:Vector})
#     length(paths) > 1 || return [1]
#     nogo_seqs, all_seqs = [], []
#     for path in paths
#         for ni in 1:length(path)-3+1
#             seq = path[ni:ni+3-1]
#             push!(all_seqs,seq)
#             for npath in paths
#                 if length(npath) >= length(seq) && has_subvector(npath, seq[2:end-1]) && !has_subvector(npath, seq)
#                     push!(nogo_seqs,seq)
#                 end
#             end
#         end
#     end
#     seqs = setdiff(all_seqs,nogo_seqs)
#     res = [[np for (np,path) in enumerate(paths) if has_subvector(path,seq)] for seq in seqs]
#     sort!(res, by=length)
#     return first(res)
# end


text:
using FIRRTL
using Test

include("types.jl")


text:
using ImageMagick,FileIO,JLD,FixedPointNumbers
const imgdir = "/mnt/ai/work/askin/150723_exp3_GCAMP_Images"

# Save all images as a list of images
function img2jld()
    global img = Any[]
    for i=1:7002
        i%100==0 && print(".")
        push!(img, load("$imgdir/$i.tiff"))
    end
    println(); info("Saving...")
    JLD.save("img.jld","img",img)
end

# Save all images as a 512x512x7002 array of UFixed16 values
function img3jld(img)
    global img3=Array(UFixed16,512,512,7002)
    for i=1:7002; copy!(img3,1+(i-1)*512*512,img[i].data,1,512*512); end
    JLD.save("img3.jld","img3",img3)
end

# Find centers as local maxima of pixel values: doesn't work very well
function centers{T}(x::Matrix{T}; r=30)
    c = Any[]
    for i=1:size(x,1)
        for j=1:size(x,2)
            localmax(x,i,j,r) && push!(c, (i,j,x[i,j]))
        end
    end
    return sort(c, by=(x->x[3]), rev=true)
end

function findmax2(a)
    v,i = findmax(a)
    v,ind2sub(a,i)
end

function centers32{T}(x::Matrix{T}; r=32)
    c = Any[]
    for j=1:r:size(x,2)
        jj=min(j+r-1,size(x,2))
        for i=1:r:size(x,1)
            ii=min(i+r-1,size(x,1))
            push!(c, ((i+ii)/2,(j+jj)/2,mean(sub(x,i:ii,j:jj))))
        end
    end
    return sort(c, by=(x->x[3]), rev=true)
end

function localmax{T}(x::Matrix{T}, i, j, r)
    for ii=max(1,i-r):min(size(x,1),i+r)
        for jj=max(1,j-r):min(size(x,2),j+r)
            if x[i,j] < x[ii,jj]
                return false
            end
        end
    end
    return true
end

# Find the image with the max pixel for every 20 image block
function findmax20()
    zz = Int[]
    for i=0:20:size(img3,3)
        d = img3[:,:,i+1:min(i+20,7002)]
        (v,n) = findmax(d)
        (x,y,z) = ind2sub(d,n)
        push!(zz,z)
    end
    zz
end

# Find variation in activation looking at the value of the max pixel for every 20 image block
function findmaxpix()
    zz = Float32[]
    for i=0:20:size(img3,3)
        d = img3[:,:,i+1:min(i+20,7002)]
        push!(zz,maximum(d))
    end
    zz
end

# Find variation based on distance around the max pixel in every block
function pixeldistance(; r=50)
    dr = Float32[]; dv = Float32[]
    for i=0:20:size(img3,3)
        d = img3[:,:,i+1:min(i+20,7002)]
        (v,n) = findmax(d)
        (v > 0.5) && continue
        (x,y,z) = ind2sub(d,n)
        for xx = max(1,x-r):min(512,x+r)
            for yy = max(1,y-r):min(512,y+r)
                vv = d[xx,yy,z]
                push!(dr, sqrt((x-xx)^2 + (y-yy)^2))
                push!(dv, Float32(vv)-Float32(v))
            end
        end
    end
    (dr, dv)
end

function brighten(img)
    c = copy(img)
    for i=1:length(c.data)
        c.data[i] = min(1, 5*c.data[i])
        # c.data[i] = max(0, 1-4*c.data[i])
    end
    return c
end

max_stars_repo_name: JuliaTagBot/NeuronTracking

max_issues_repo_name: JuliaTagBot/NeuronTracking

max_forks_repo_name: JuliaTagBot/NeuronTracking


text:
# visualize the final states of vocab charlm model after the roots of words

using TSne, Statistics

include("charlm.jl")

# load model
if !@isdefined(model)
    include(Knet.dir("src/models21/deepmap.jl"))
    include(Knet.dir("src/models21/fileio.jl"))
    model = charlm_model()
    charlm_loss(model, tst) # necessary to initialize weights, fix this
    setweights!(model, "charlm.jld2")
end
#@show charlm_loss(model,tst)

# load unique roots from data: 8864 of them
if !@isdefined(rootpos)
    posdict = Dict{String,Int}()
    rootpos = Dict{String,Set{String}}()
    for line in eachline("trmor2018.filtered")
        local sur, lex, root, pos
        sur,lex = split(line,'\t')
        root, pos = split(lex,['+','^'])[1:2]
        if !haskey(posdict, pos)
            posdict[pos] = 1+length(posdict)
        end
        push!(get!(rootpos,root,Set{String}()), pos)
    end
    roots = collect(keys(rootpos))
end

# construct a dataset with batchsize = 1
if !@isdefined(rootdata)
    rootdata = map(roots) do root
        chars = split(root,"")
        ints = map(chars) do char
            if char == "þ"; char = "ş"; end # a bug in dataset
            if char == "ý"; char = "ı"; end
            if char == "ð"; char = "ğ"; end
            surface_vocab.w2i[char]
        end
        permutedims(ints)       # batch needs (B,T) shape
    end
end

# try clustering the final hidden state
# model(charlm_input(tst[1])) |> "37×128×8 CUDA.CuArray{Float32, 3, CUDA.Mem.DeviceBuffer}"
# model[1:end-1](charlm_input(tst[1])) |> "350×128×8 CUDA.CuArray{Float32, 3, CUDA.Mem.DeviceBuffer}"

hidden_states(model, batch) = model[1:end-1](charlm_input(batch)) |> Array{Float64} # Warning: Float32 makes tsne much slower!
final_hidden_states(model, batch) = hidden_states(model, batch)[:,:,end]
final_hidden_states(model, data::Vector{Matrix{Int64}}) = reduce(hcat, (x->final_hidden_states(model,x)).(data))
#@show summary(final_hidden_states(model, tst))
#surface_labels(batch) = reduce(*, surface_vocab.i2w[batch]; dims=2)
#surface_labels(data::Vector{Matrix{Int64}}) = reduce(vcat, surface_labels.(data))

if !@isdefined rootfhs
    rootfhs = final_hidden_states(model, rootdata) |> permutedims # (8864, 350)
end

if !@isdefined root_tsne
    root_tsne = tsne(rootfhs)
end

using Plots, ColorSchemes
plotly()
labels = map(roots) do r; posdict[first(rootpos[r])]; end
poscolor = get(ColorSchemes.tab20, labels, :extrema)
scatter(root_tsne[:,1], root_tsne[:,2]; hover=roots, color=poscolor)

# From: https://github.com/nassarhuda/JuliaTutorials/blob/master/TSNE/TSNE.ipynb
# function my_tsne_plot(Y,d,labels; marker=5)
#     colorsvec = get(ColorSchemes.tab10, labels, :extrema)
#     plot(legend=false,axis=false,grid=false)
#     tsneplot = scatter!(Y[:,1], Y[:,2],
#                         marker=(marker,marker,:auto,stroke(0)),alpha=0.8,
#                         background_color=:black,
#                         color=colorsvec)
#     for i = 0:9
#         mx,my = mean(Y[findall(labels.==i),:],dims=1)
#         s = Shape(mx.+[-d,d,d,-d],my.+[-d,-d,d,d])
#         plot!(s,color=ColorSchemes.tab10[i+1],linewidth=0,alpha=0.9)
#         annotate!(mx,my,text(string(i),color=:white))
#     end
#     return tsneplot
# end

# function tsne_mnist(N=3000; plot=false)
#     # using MLDatasets, ColorSchemes, Plots
#     alldata, labels = MNIST.traindata(Float64)
#     data = reshape(permutedims(alldata[:, :, 1:N], (3, 1, 2)),
#                    N, size(alldata, 1)*size(alldata, 2))
#     labels = labels[1:N]
#     rescale(A; dims=1) = (A .- mean(A, dims=dims)) ./ max.(std(A, dims=dims), eps())
#     # Normalize the data, this should be done if there are large scale differences in the dataset
#     X = rescale(data, dims=1)
#     @show summary(X)
#     Y1 = tsne(X, 2, 50, 1000, 20.0)
#     plot ? my_tsne_plot(Y1, 5, labels) : Y1
# end

# Timing: tsne(tstfhs[1:1000,:], 2, 50, 1000, 20.0);

# mnist tsne:        charlm tsne: (why is it taking so much longer with Float32?)
# N    secs  KL      secs  KL
# 1000 3     0.7470  3     0.3175
# 2000 14    1.0122  13    0.3681
# 3000 31    1.1529  25    0.4125
# 4000 48    1.2079

nothing

max_stars_repo_name: denizyuret/research

max_issues_repo_name: denizyuret/research

max_forks_repo_name: denizyuret/research


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

using Markdown
using InteractiveUtils

# ╔═╡ 7357f256-ead0-401c-9cff-1aa099b67aa0
begin
    import Pkg
    Pkg.activate(mktempdir())
    Pkg.add([
        Pkg.PackageSpec(name="Plots", version="1"),
        Pkg.PackageSpec(name="PlutoUI", version="0.7"),
	    Pkg.PackageSpec(name="DifferentialEquations", version="6"),
	    Pkg.PackageSpec(name="Parameters", version="0.12"),
        Pkg.PackageSpec(name="LabelledArrays", version="1"),
        Pkg.PackageSpec(name="Setfield", version="0.7"),
    ])
    using Plots, PlutoUI, LinearAlgebra, DifferentialEquations, Parameters, LabelledArrays, Setfield

	Plots.gr(fmt=:png, lw=2)

	PlutoUI.TableOfContents()
end

# ╔═╡ 0e455e50-bdcb-11eb-08b9-8d2d26788d14
md"""
# Fig 1.07 Collins toggle switch

For Figures 1.7, 7.13, 7.14, 7.15
"""

# ╔═╡ 30cdc452-ff17-4b08-939a-98c7dbc0c674
# Initial conditions
u0 = LVector(s1=0.075, s2=2.5)

# ╔═╡ 520f36c3-e164-459e-8869-99a8421774f6
# Simulation time
tend = 50.0

# ╔═╡ f47e84f9-7ff5-4c32-996a-47c600b1a71a
# Parameters
p = (a1 = 3.0, a2 = 2.5, β = 4.0, γ = 4.0, i1 = 0.0, i2 = 0.0)

# ╔═╡ 463ef6e7-73db-4b55-b97f-71e90c2f14f4
begin
	# Events and callbacks
    on_i1!(integrator) = integrator.p = @set p.i1 = 10.0
    off_i1!(integrator) = integrator.p = @set p.i1 = 0.0
    on_i2!(integrator) = integrator.p = @set p.i2 = 10.0
    off_i2!(integrator) = integrator.p = @set p.i2 = 0.0

    cbs = CallbackSet(PresetTimeCallback([10.0], on_i2!),
                  PresetTimeCallback([20.0], off_i2!),
                  PresetTimeCallback([30.0], on_i1!),
                  PresetTimeCallback([40.0], off_i1!))
end

# ╔═╡ d29832bb-74f9-443a-af7f-a12e7ae652c3
md"""
# Appendix
Running environment and some auxillary functions.

**You can use [this helper tool](https://fonsp.com/article-test-3/pkghelper.html) to generate these commands!**
"""

# ╔═╡ 796ff12c-b74d-4fcf-b228-4e3c6fed2a89
begin
	# Convenience functions
    hill(x, k) = x / (x + k)
    hill(x, k, n) = hill(x^n, k^n)
    exprel(x) = ifelse(x≈zero(x), one(x), x / expm1(x))
end

# ╔═╡ 34178e95-cff5-400e-a5a6-bdbcf8326214
# The model
	function collins!(du, u, p, t)
		@unpack a1, a2, β, γ, i1, i2 = p
		@unpack s1, s2 = u
		du.s1 = a1 * hill(1 + i2, s2, β) - s1
		du.s2 = a2 * hill(1 + i1, s1, γ) - s2
    	return du
	end

# ╔═╡ 759580e6-0910-43a2-b81a-cdaa3bf34441
sol = solve(ODEProblem(collins!, u0, tend, p), callback=cbs)

# ╔═╡ b922fac1-e43c-45eb-a6db-0a9b18f541aa
plot(sol, legend=:right, xlabel="Time", ylabel="Concentration", title="Figure 1.7")

# ╔═╡ Cell order:
# ╠═0e455e50-bdcb-11eb-08b9-8d2d26788d14
# ╠═34178e95-cff5-400e-a5a6-bdbcf8326214
# ╠═463ef6e7-73db-4b55-b97f-71e90c2f14f4
# ╠═30cdc452-ff17-4b08-939a-98c7dbc0c674
# ╠═520f36c3-e164-459e-8869-99a8421774f6
# ╠═f47e84f9-7ff5-4c32-996a-47c600b1a71a
# ╠═759580e6-0910-43a2-b81a-cdaa3bf34441
# ╠═b922fac1-e43c-45eb-a6db-0a9b18f541aa
# ╠═d29832bb-74f9-443a-af7f-a12e7ae652c3
# ╠═796ff12c-b74d-4fcf-b228-4e3c6fed2a89
# ╠═7357f256-ead0-401c-9cff-1aa099b67aa0

max_stars_repo_name: stevengogogo/pluto-notebooks

max_issues_repo_name: stevengogogo/pluto-notebooks

max_forks_repo_name: stevengogogo/pluto-notebooks


text:
@testset "Topology: topology_to_graphs" begin
    @testset "Directed" begin
        nodes_matrix = permutedims([[1, "A"] [2, "B"] [3, "C"] [4, "D"] [5, "E"] [6, "F"]], [2, 1])
        edges_matrix = [[1, 2, 1.0] [2, 3, 2.0] [2, 6, 2.0] [3, 4, 3.0] [3, 6, 3.0] [5, 6, 5.0]]'
        traffic_matrix = [[1, 3, .5] [1, 6, .5] [4, 1, .5] [4, 5, .5] [5, 2, .5] [6, 4, .5]]'
        t = Topology(nodes_matrix, edges_matrix, traffic_matrix)

        g, k = topology_to_graphs(t)

        @test nv(g) == 6
        @test ne(g) == 12 # Directed graph
        @test nv(k) == 6
        @test ne(k) == 6

        for i in 1:6
            @test get_prop(g, i, :name) == string('A' + i - 1)
            @test get_prop(k, i, :name) == string('A' + i - 1)
        end
        for (s, t) in [(1, 2), (2, 3), (2, 6), (3, 4), (3, 6), (5, 6)]
            @test has_edge(g, s, t)
            @test get_prop(g, s, t, :capacity) == s

            @test has_edge(g, t, s)
            @test get_prop(g, t, s, :capacity) == s
        end
        for (s, t) in [(1, 3), (1, 6), (4, 1), (4, 5), (5, 2), (6, 4)]
            @test has_edge(k, s, t)
            @test ! has_edge(k, t, s)

            @test get_prop(k, s, t, :demand) == 0.5
        end
    end

    @testset "Undirected" begin
        nodes_matrix = permutedims([[1, "A"] [2, "B"] [3, "C"] [4, "D"] [5, "E"] [6, "F"]], [2, 1])
        edges_matrix = [[1, 2, 1.0] [2, 3, 2.0] [2, 6, 2.0] [3, 4, 3.0] [3, 6, 3.0] [5, 6, 5.0]]'
        traffic_matrix = [[1, 3, .5] [1, 6, .5] [4, 1, .5] [4, 5, .5] [5, 2, .5] [6, 4, .5]]'
        t = Topology(nodes_matrix, edges_matrix, traffic_matrix)

        g, k = topology_to_graphs(t, make_network_undirected=false)

        @test nv(g) == 6
        @test ne(g) == 6 # Undirected graph
        @test nv(k) == 6
        @test ne(k) == 6

        for i in 1:6
            @test get_prop(g, i, :name) == string('A' + i - 1)
            @test get_prop(k, i, :name) == string('A' + i - 1)
        end
        for (s, t) in [(1, 2), (2, 3), (2, 6), (3, 4), (3, 6), (5, 6)]
            @test has_edge(g, s, t)
            @test get_prop(g, s, t, :capacity) == s

            @test ! has_edge(g, t, s)
        end
        for (s, t) in [(1, 3), (1, 6), (4, 1), (4, 5), (5, 2), (6, 4)]
            @test has_edge(k, s, t)
            @test ! has_edge(k, t, s)

            @test get_prop(k, s, t, :demand) == 0.5
        end
    end
end


text:
using DiffLinearAlgebra, Compat.Test

import FDM: assert_approx_equal, central_fdm
import DiffLinearAlgebra: AA, AM, AVM, AS, ASVM, Arg1, Arg2

@testset "DiffLinearAlgebra" begin
    include("test_util.jl")
    include("util.jl")
    include("generic.jl")
    include("blas.jl")
    include("diagonal.jl")
    include("triangular.jl")
    include("uniformscaling.jl")
    include("factorization/cholesky.jl")
    include("test_imports.jl")
end


text:
@testset "LSIF" begin
  for (i, (pair, rtol)) in enumerate([(pair₁, 2e-1), (pair₂, 4e-1)])
    d_nu, d_de = pair
    Random.seed!(123)
    x_nu, x_de = rand(d_nu, 1000), rand(d_de, 500)

    # estimated density ratio
    r̂ = densratio(x_nu, x_de, LSIF(σ=1.0, b=100))

    # simplex constraints
    @test abs(mean(r̂) - 1) ≤ 2e-2
    @test all(r̂ .≤ Inf)

    r = pdf.(d_nu, x_de) ./ pdf.(d_de, x_de)
    @test r ≈ r̂ rtol=rtol

    if visualtests
      gr(size=(800,800))
      @plottest plot_d_nu(pair,x_de,r̂) joinpath(datadir,"LSIF-$i.png") !istravis
    end
  end
end


text:
module Raven

println("Raven: Small polaron mobility by hopping")

include("types.jl")
include("rates.jl")
include("ToFeT_io.jl")

end # module