From 744ce59840ad00d35c6f66aa9e4198e76893cfe3 Mon Sep 17 00:00:00 2001
From: chaoming <adaduo@outlook.com>
Date: Tue, 9 Jan 2024 11:18:50 +0800
Subject: [PATCH 1/4] [math] add ``ein_rearrange``, ``ein_reduce``, and
 ``ein_repeat`` inspired by `einops` pckage

---
 brainpy/_src/math/einops.py                   |  728 ++++++++
 brainpy/_src/math/einops_parsing.py           |  153 ++
 brainpy/_src/math/interoperability.py         |   10 +-
 brainpy/_src/math/others.py                   |   27 +-
 brainpy/_src/math/tests/test_einops.py        |  332 ++++
 .../_src/math/tests/test_einops_parsing.py    |  111 ++
 brainpy/math/__init__.py                      |    1 +
 brainpy/math/einops.py                        |    6 +
 brainpy/math/interoperability.py              |    1 +
 docs/tutorial_math/einops_in_brainpy.ipynb    | 1509 +++++++++++++++++
 docs/tutorial_math/index.rst                  |    1 +
 docs/tutorial_math/test_images.npy            |  Bin 0 -> 1327232 bytes
 12 files changed, 2876 insertions(+), 3 deletions(-)
 create mode 100644 brainpy/_src/math/einops.py
 create mode 100644 brainpy/_src/math/einops_parsing.py
 create mode 100644 brainpy/_src/math/tests/test_einops.py
 create mode 100644 brainpy/_src/math/tests/test_einops_parsing.py
 create mode 100644 brainpy/math/einops.py
 create mode 100644 docs/tutorial_math/einops_in_brainpy.ipynb
 create mode 100644 docs/tutorial_math/test_images.npy

diff --git a/brainpy/_src/math/einops.py b/brainpy/_src/math/einops.py
new file mode 100644
index 000000000..d42026974
--- /dev/null
+++ b/brainpy/_src/math/einops.py
@@ -0,0 +1,728 @@
+import functools
+import itertools
+from collections import OrderedDict
+from typing import Set, Tuple, List, Dict, Union, Callable, Optional, cast
+
+import jax
+import numpy as np
+
+from . import compat_numpy as bnp
+from . import others as bnp2
+from .einops_parsing import ParsedExpression, _ellipsis, AnonymousAxis, EinopsError
+from .ndarray import Array
+
+__all__ = [
+  'ein_reduce', 'ein_rearrange', 'ein_repeat', 'ein_shape',
+]
+
+Tensor = Union[Array, jax.Array]
+ReductionCallable = Callable[[Tensor, Tuple[int, ...]], Tensor]
+Reduction = Union[str, ReductionCallable]
+
+_reductions = ("min", "max", "sum", "mean", "prod", "any", "all")
+
+# magic integers are required to stay within
+# traceable subset of language
+_unknown_axis_length = -999999
+_expected_axis_length = -99999
+
+
+def _product(sequence: List[int]) -> int:
+  """minimalistic product that works both with numbers and symbols. Supports empty lists"""
+  result = 1
+  for element in sequence:
+    result *= element
+  return result
+
+
+def _reduce_axes(tensor, reduction_type: Reduction, reduced_axes: List[int]):
+  if callable(reduction_type):
+    # custom callable
+    return reduction_type(tensor, tuple(reduced_axes))
+  else:
+    # one of built-in operations
+    assert reduction_type in _reductions
+    if reduction_type == "mean":
+      if not bnp2.is_float_type(tensor):
+        raise NotImplementedError("reduce_mean is not available for non-floating tensors")
+    return __reduce(tensor, reduction_type, tuple(reduced_axes))
+
+
+def __reduce(x: Union[Array, jax.Array], operation: str, reduced_axes):
+  if operation == "min":
+    return x.min(axis=reduced_axes)
+  elif operation == "max":
+    return x.max(axis=reduced_axes)
+  elif operation == "sum":
+    return x.sum(axis=reduced_axes)
+  elif operation == "mean":
+    return x.mean(axis=reduced_axes)
+  elif operation == "prod":
+    return x.prod(axis=reduced_axes)
+  elif operation == "any":
+    return x.any(axis=reduced_axes)
+  elif operation == "all":
+    return x.all(axis=reduced_axes)
+  else:
+    raise NotImplementedError("Unknown reduction ", operation)
+
+
+def _optimize_transformation(init_shapes, reduced_axes, axes_reordering, final_shapes):
+  # 'collapses' neighboring axes if those participate in the result pattern in the same order
+  # TODO add support for added_axes
+  assert len(axes_reordering) + len(reduced_axes) == len(init_shapes)
+  # joining consecutive axes that will be reduced
+  # possibly we can skip this if all backends can optimize this (not sure)
+  reduced_axes = tuple(sorted(reduced_axes))
+  for i in range(len(reduced_axes) - 1)[::-1]:
+    if reduced_axes[i] + 1 == reduced_axes[i + 1]:
+      removed_axis = reduced_axes[i + 1]
+      removed_length = init_shapes[removed_axis]
+      init_shapes = init_shapes[:removed_axis] + init_shapes[removed_axis + 1:]
+      init_shapes[removed_axis - 1] *= removed_length
+      reduced_axes = reduced_axes[: i + 1] + tuple(axis - 1 for axis in reduced_axes[i + 2:])
+
+  # removing axes that are moved together during reshape
+  def build_mapping():
+    init_to_final = {}
+    for axis in range(len(init_shapes)):
+      if axis in reduced_axes:
+        init_to_final[axis] = None
+      else:
+        after_reduction = sum(x is not None for x in init_to_final.values())
+        init_to_final[axis] = list(axes_reordering).index(after_reduction)
+    return init_to_final
+
+  init_axis_to_final_axis = build_mapping()
+
+  for init_axis in range(len(init_shapes) - 1)[::-1]:
+    if init_axis_to_final_axis[init_axis] is None:
+      continue
+    if init_axis_to_final_axis[init_axis + 1] is None:
+      continue
+    if init_axis_to_final_axis[init_axis] + 1 == init_axis_to_final_axis[init_axis + 1]:
+      removed_axis = init_axis + 1
+      removed_length = init_shapes[removed_axis]
+      removed_axis_after_reduction = sum(x not in reduced_axes for x in range(removed_axis))
+
+      reduced_axes = tuple(axis if axis < removed_axis else axis - 1 for axis in reduced_axes)
+      init_shapes = init_shapes[:removed_axis] + init_shapes[removed_axis + 1:]
+      init_shapes[removed_axis - 1] *= removed_length
+      old_reordering = axes_reordering
+      axes_reordering = []
+      for axis in old_reordering:
+        if axis == removed_axis_after_reduction:
+          pass
+        elif axis < removed_axis_after_reduction:
+          axes_reordering.append(axis)
+        else:
+          axes_reordering.append(axis - 1)
+      init_axis_to_final_axis = build_mapping()
+
+  return init_shapes, reduced_axes, axes_reordering, final_shapes
+
+
+CookedRecipe = Tuple[Optional[List[int]], Optional[List[int]], List[int], Dict[int, int], Optional[List[int]], int]
+
+# Actual type is tuple[tuple[str, int], ...]
+# However torch.jit.script does not "understand" the correct type,
+# and torch_specific will use list version.
+HashableAxesLengths = Tuple[Tuple[str, int], ...]
+FakeHashableAxesLengths = List[Tuple[str, int]]
+
+
+class TransformRecipe:
+  """
+  Recipe describes actual computation pathway.
+  Recipe can be applied to a tensor or variable.
+  """
+
+  # structure is non-mutable. In future, this can be non-mutable dataclass (python 3.7+)
+  # update: pytorch 2.0 torch.jit.script seems to have problems with dataclasses unless they were explicitly provided
+
+  def __init__(
+      self,
+      # list of sizes (or just sizes) for elementary axes as they appear in left expression.
+      # this is what (after computing unknown parts) will be a shape after first transposition.
+      # This does not include any ellipsis dimensions.
+      elementary_axes_lengths: List[int],
+      # if additional axes are provided, they should be set in prev array
+      # This shows mapping from name to position
+      axis_name2elementary_axis: Dict[str, int],
+      # each dimension in input can help to reconstruct length of one elementary axis
+      # or verify one of dimensions. Each element points to element of elementary_axes_lengths.
+      input_composition_known_unknown: List[Tuple[List[int], List[int]]],
+      # permutation applied to elementary axes, if ellipsis is absent
+      axes_permutation: List[int],
+      # permutation puts reduced axes in the end, we only need to know the first position.
+      first_reduced_axis: int,
+      # at which positions which of elementary axes should appear. Axis position -> axis index.
+      added_axes: Dict[int, int],
+      # ids of axes as they appear in result, again pointers to elementary_axes_lengths,
+      # only used to infer result dimensions
+      output_composite_axes: List[List[int]],
+  ):
+    self.elementary_axes_lengths: List[int] = elementary_axes_lengths
+    self.axis_name2elementary_axis: Dict[str, int] = axis_name2elementary_axis
+    self.input_composition_known_unknown: List[Tuple[List[int], List[int]]] = input_composition_known_unknown
+    self.axes_permutation: List[int] = axes_permutation
+
+    self.first_reduced_axis: int = first_reduced_axis
+    self.added_axes: Dict[int, int] = added_axes
+    self.output_composite_axes: List[List[int]] = output_composite_axes
+
+
+def _reconstruct_from_shape_uncached(
+    self: TransformRecipe, shape: List[int], axes_dims: FakeHashableAxesLengths
+) -> CookedRecipe:
+  """
+  Reconstruct all actual parameters using shape.
+  Shape is a tuple that may contain integers, shape symbols (tf, theano) and UnknownSize (tf, previously mxnet)
+  known axes can be integers or symbols, but not Nones.
+  """
+  # magic number
+  need_init_reshape = False
+
+  # last axis is allocated for collapsed ellipsis
+  axes_lengths: List[int] = list(self.elementary_axes_lengths)
+  for axis, dim in axes_dims:
+    axes_lengths[self.axis_name2elementary_axis[axis]] = dim
+
+  for input_axis, (known_axes, unknown_axes) in enumerate(self.input_composition_known_unknown):
+    length = shape[input_axis]
+    if len(known_axes) == 0 and len(unknown_axes) == 1:
+      # shortcut for the most common case
+      axes_lengths[unknown_axes[0]] = length
+      continue
+
+    known_product = 1
+    for axis in known_axes:
+      known_product *= axes_lengths[axis]
+
+    if len(unknown_axes) == 0:
+      if isinstance(length, int) and isinstance(known_product, int) and length != known_product:
+        raise EinopsError(f"Shape mismatch, {length} != {known_product}")
+    else:
+      # assert len(unknown_axes) == 1, 'this is enforced when recipe is created, so commented out'
+      if isinstance(length, int) and isinstance(known_product, int) and length % known_product != 0:
+        raise EinopsError(f"Shape mismatch, can't divide axis of length {length} in chunks of {known_product}")
+
+      unknown_axis = unknown_axes[0]
+      inferred_length: int = length // known_product
+      axes_lengths[unknown_axis] = inferred_length
+
+    if len(known_axes) + len(unknown_axes) != 1:
+      need_init_reshape = True
+
+  # at this point all axes_lengths are computed (either have values or variables, but not Nones)
+
+  # elementary axes are ordered as they appear in input, then all added axes
+  init_shapes: Optional[List[int]] = axes_lengths[: len(self.axes_permutation)] if need_init_reshape else None
+
+  need_final_reshape = False
+  final_shapes: List[int] = []
+  for grouping in self.output_composite_axes:
+    lengths = [axes_lengths[elementary_axis] for elementary_axis in grouping]
+    final_shapes.append(_product(lengths))
+    if len(lengths) != 1:
+      need_final_reshape = True
+
+  added_axes: Dict[int, int] = {
+    pos: axes_lengths[pos_in_elementary] for pos, pos_in_elementary in self.added_axes.items()
+  }
+
+  # this list can be empty
+  reduced_axes = list(range(self.first_reduced_axis, len(self.axes_permutation)))
+
+  n_axes_after_adding_axes = len(added_axes) + len(self.axes_permutation)
+
+  axes_reordering: Optional[List[int]] = self.axes_permutation
+  if self.axes_permutation == list(range(len(self.axes_permutation))):
+    axes_reordering = None
+
+  _final_shapes = final_shapes if need_final_reshape else None
+  return init_shapes, axes_reordering, reduced_axes, added_axes, _final_shapes, n_axes_after_adding_axes
+
+
+_reconstruct_from_shape = functools.lru_cache(1024)(_reconstruct_from_shape_uncached)
+
+
+def _apply_recipe(
+    recipe: TransformRecipe, tensor: Tensor, reduction_type: Reduction, axes_lengths: HashableAxesLengths
+) -> Tensor:
+  # this method implements actual work for all backends for 3 operations
+  try:
+    init_shapes, axes_reordering, reduced_axes, added_axes, final_shapes, n_axes_w_added = (
+      _reconstruct_from_shape(recipe, bnp.shape(tensor), axes_lengths))
+  except TypeError:
+    # shape or one of passed axes lengths is not hashable (i.e. they are symbols)
+    _result = _reconstruct_from_shape_uncached(recipe, bnp.shape(tensor), axes_lengths)
+    (init_shapes, axes_reordering, reduced_axes, added_axes, final_shapes, n_axes_w_added) = _result
+  if init_shapes is not None:
+    tensor = bnp.reshape(bnp.as_jax(tensor), init_shapes)
+  if axes_reordering is not None:
+    tensor = bnp.transpose(bnp.as_jax(tensor), axes_reordering)
+  if len(reduced_axes) > 0:
+    tensor = _reduce_axes(bnp.as_jax(tensor), reduction_type=reduction_type, reduced_axes=reduced_axes)
+  if len(added_axes) > 0:
+    tensor = bnp2.add_axes(tensor, n_axes=n_axes_w_added, pos2len=added_axes)
+  if final_shapes is not None:
+    tensor = bnp.reshape(bnp.as_jax(tensor), final_shapes)
+  return tensor
+
+
+def _apply_recipe_array_api(
+    xp, recipe: TransformRecipe, tensor: Tensor, reduction_type: Reduction, axes_lengths: HashableAxesLengths
+) -> Tensor:
+  # completely-inline implementation
+  init_shapes, axes_reordering, reduced_axes, added_axes, final_shapes, n_axes_w_added = _reconstruct_from_shape(
+    recipe, tensor.shape, axes_lengths
+  )
+  if init_shapes is not None:
+    tensor = xp.reshape(tensor, init_shapes)
+  if axes_reordering is not None:
+    tensor = xp.permute_dims(tensor, axes_reordering)
+  if len(reduced_axes) > 0:
+    if callable(reduction_type):
+      # custom callable
+      tensor = reduction_type(tensor, tuple(reduced_axes))
+    else:
+      # one of built-in operations
+      assert reduction_type in _reductions
+      tensor = getattr(xp, reduction_type)(tensor, axis=tuple(reduced_axes))
+  if len(added_axes) > 0:
+    # we use broadcasting
+    for axis_position, axis_length in added_axes.items():
+      tensor = xp.expand_dims(tensor, axis=axis_position)
+
+    final_shape = list(tensor.shape)
+    for axis_position, axis_length in added_axes.items():
+      final_shape[axis_position] = axis_length
+
+    tensor = xp.broadcast_to(tensor, final_shape)
+  if final_shapes is not None:
+    tensor = xp.reshape(tensor, final_shapes)
+  return tensor
+
+
+@functools.lru_cache(256)
+def _prepare_transformation_recipe(
+    pattern: str,
+    operation: Reduction,
+    axes_names: Tuple[str, ...],
+    ndim: int,
+) -> TransformRecipe:
+  """Perform initial parsing of pattern and provided supplementary info
+  axes_lengths is a tuple of tuples (axis_name, axis_length)
+  """
+  left_str, rght_str = pattern.split("->")
+  left = ParsedExpression(left_str)
+  rght = ParsedExpression(rght_str)
+
+  # checking that axes are in agreement - new axes appear only in repeat, while disappear only in reduction
+  if not left.has_ellipsis and rght.has_ellipsis:
+    raise EinopsError("Ellipsis found in right side, but not left side of a pattern {}".format(pattern))
+  if left.has_ellipsis and left.has_ellipsis_parenthesized:
+    raise EinopsError("Ellipsis inside parenthesis in the left side is not allowed: {}".format(pattern))
+  if operation == "rearrange":
+    if left.has_non_unitary_anonymous_axes or rght.has_non_unitary_anonymous_axes:
+      raise EinopsError("Non-unitary anonymous axes are not supported in rearrange (exception is length 1)")
+    difference = set.symmetric_difference(left.identifiers, rght.identifiers)
+    if len(difference) > 0:
+      raise EinopsError("Identifiers only on one side of expression (should be on both): {}".format(difference))
+  elif operation == "repeat":
+    difference = set.difference(left.identifiers, rght.identifiers)
+    if len(difference) > 0:
+      raise EinopsError("Unexpected identifiers on the left side of repeat: {}".format(difference))
+    axes_without_size = set.difference(
+      {ax for ax in rght.identifiers if not isinstance(ax, AnonymousAxis)},
+      {*left.identifiers, *axes_names},
+    )
+    if len(axes_without_size) > 0:
+      raise EinopsError("Specify sizes for new axes in repeat: {}".format(axes_without_size))
+  elif operation in _reductions or callable(operation):
+    difference = set.difference(rght.identifiers, left.identifiers)
+    if len(difference) > 0:
+      raise EinopsError("Unexpected identifiers on the right side of reduce {}: {}".format(operation, difference))
+  else:
+    raise EinopsError("Unknown reduction {}. Expect one of {}.".format(operation, _reductions))
+
+  if left.has_ellipsis:
+    n_other_dims = len(left.composition) - 1
+    if ndim < n_other_dims:
+      raise EinopsError(f"Wrong shape: expected >={n_other_dims} dims. Received {ndim}-dim tensor.")
+    ellipsis_ndim = ndim - n_other_dims
+    ell_axes = [_ellipsis + str(i) for i in range(ellipsis_ndim)]
+    left_composition = []
+    for composite_axis in left.composition:
+      if composite_axis == _ellipsis:
+        for axis in ell_axes:
+          left_composition.append([axis])
+      else:
+        left_composition.append(composite_axis)
+
+    rght_composition = []
+    for composite_axis in rght.composition:
+      if composite_axis == _ellipsis:
+        for axis in ell_axes:
+          rght_composition.append([axis])
+      else:
+        group = []
+        for axis in composite_axis:
+          if axis == _ellipsis:
+            group.extend(ell_axes)
+          else:
+            group.append(axis)
+        rght_composition.append(group)
+
+    left.identifiers.update(ell_axes)
+    left.identifiers.remove(_ellipsis)
+    if rght.has_ellipsis:
+      rght.identifiers.update(ell_axes)
+      rght.identifiers.remove(_ellipsis)
+  else:
+    if ndim != len(left.composition):
+      raise EinopsError(f"Wrong shape: expected {len(left.composition)} dims. Received {ndim}-dim tensor.")
+    left_composition = left.composition
+    rght_composition = rght.composition
+
+  # parsing all dimensions to find out lengths
+  axis_name2known_length: Dict[Union[str, AnonymousAxis], int] = OrderedDict()
+  for composite_axis in left_composition:
+    for axis_name in composite_axis:
+      if isinstance(axis_name, AnonymousAxis):
+        axis_name2known_length[axis_name] = axis_name.value
+      else:
+        axis_name2known_length[axis_name] = _unknown_axis_length
+
+  # axis_ids_after_first_reshape = range(len(axis_name2known_length)) at this point
+
+  repeat_axes_names = []
+  for axis_name in rght.identifiers:
+    if axis_name not in axis_name2known_length:
+      if isinstance(axis_name, AnonymousAxis):
+        axis_name2known_length[axis_name] = axis_name.value
+      else:
+        axis_name2known_length[axis_name] = _unknown_axis_length
+      repeat_axes_names.append(axis_name)
+
+  axis_name2position = {name: position for position, name in enumerate(axis_name2known_length)}
+
+  # axes provided as kwargs
+  for elementary_axis in axes_names:
+    if not ParsedExpression.check_axis_name(elementary_axis):
+      raise EinopsError("Invalid name for an axis", elementary_axis)
+    if elementary_axis not in axis_name2known_length:
+      raise EinopsError("Axis {} is not used in transform".format(elementary_axis))
+    axis_name2known_length[elementary_axis] = _expected_axis_length
+
+  input_axes_known_unknown = []
+  # some shapes are inferred later - all information is prepared for faster inference
+  for i, composite_axis in enumerate(left_composition):
+    known: Set[str] = {axis for axis in composite_axis if axis_name2known_length[axis] != _unknown_axis_length}
+    unknown: Set[str] = {axis for axis in composite_axis if axis_name2known_length[axis] == _unknown_axis_length}
+    if len(unknown) > 1:
+      raise EinopsError("Could not infer sizes for {}".format(unknown))
+    assert len(unknown) + len(known) == len(composite_axis)
+    input_axes_known_unknown.append(
+      ([axis_name2position[axis] for axis in known], [axis_name2position[axis] for axis in unknown])
+    )
+
+  axis_position_after_reduction: Dict[str, int] = {}
+  for axis_name in itertools.chain(*left_composition):
+    if axis_name in rght.identifiers:
+      axis_position_after_reduction[axis_name] = len(axis_position_after_reduction)
+
+  result_axes_grouping: List[List[int]] = [
+    [axis_name2position[axis] for axis in composite_axis] for i, composite_axis in enumerate(rght_composition)
+  ]
+
+  ordered_axis_left = list(itertools.chain(*left_composition))
+  ordered_axis_rght = list(itertools.chain(*rght_composition))
+  reduced_axes = [axis for axis in ordered_axis_left if axis not in rght.identifiers]
+  order_after_transposition = [axis for axis in ordered_axis_rght if axis in left.identifiers] + reduced_axes
+  axes_permutation = [ordered_axis_left.index(axis) for axis in order_after_transposition]
+  added_axes = {
+    i: axis_name2position[axis_name]
+    for i, axis_name in enumerate(ordered_axis_rght)
+    if axis_name not in left.identifiers
+  }
+
+  first_reduced_axis = len(order_after_transposition) - len(reduced_axes)
+
+  return TransformRecipe(
+    elementary_axes_lengths=list(axis_name2known_length.values()),
+    axis_name2elementary_axis={axis: axis_name2position[axis] for axis in axes_names},
+    input_composition_known_unknown=input_axes_known_unknown,
+    axes_permutation=axes_permutation,
+    first_reduced_axis=first_reduced_axis,
+    added_axes=added_axes,
+    output_composite_axes=result_axes_grouping,
+  )
+
+
+def _prepare_recipes_for_all_dims(
+    pattern: str, operation: Reduction, axes_names: Tuple[str, ...]
+) -> Dict[int, TransformRecipe]:
+  """
+  Internal function, used in layers.
+  Layer makes all recipe creation when it is initialized, thus to keep recipes simple we pre-compute for all dims
+  """
+  left_str, rght_str = pattern.split("->")
+  left = ParsedExpression(left_str)
+  dims = [len(left.composition)]
+  if left.has_ellipsis:
+    dims = [len(left.composition) - 1 + ellipsis_dims for ellipsis_dims in range(8)]
+  return {ndim: _prepare_transformation_recipe(pattern, operation, axes_names, ndim=ndim) for ndim in dims}
+
+
+def ein_reduce(tensor: Union[Tensor, List[Tensor]], pattern: str, reduction: Reduction, **axes_lengths: int) -> Tensor:
+  """
+  ``ein_reduce`` provides combination of reordering and reduction using reader-friendly notation.
+
+  Examples for reduce operation:
+
+  ```python
+  >>> x = np.random.randn(100, 32, 64)
+
+  # perform max-reduction on the first axis
+  >>> y = ein_reduce(x, 't b c -> b c', 'max')
+
+  # same as previous, but with clearer axes meaning
+  >>> y = ein_reduce(x, 'time batch channel -> batch channel', 'max')
+
+  >>> x = np.random.randn(10, 20, 30, 40)
+
+  # 2d max-pooling with kernel size = 2 * 2 for image processing
+  >>> y1 = ein_reduce(x, 'b c (h1 h2) (w1 w2) -> b c h1 w1', 'max', h2=2, w2=2)
+
+  # if one wants to go back to the original height and width, depth-to-space trick can be applied
+  >>> y2 = ein_rearrange(y1, 'b (c h2 w2) h1 w1 -> b c (h1 h2) (w1 w2)', h2=2, w2=2)
+  >>> assert ein_shape(x, 'b _ h w') == ein_shape(y2, 'b _ h w')
+
+  # Adaptive 2d max-pooling to 3 * 4 grid
+  >>> ein_reduce(x, 'b c (h1 h2) (w1 w2) -> b c h1 w1', 'max', h1=3, w1=4).shape
+  (10, 20, 3, 4)
+
+  # Global average pooling
+  >>> ein_reduce(x, 'b c h w -> b c', 'mean').shape
+  (10, 20)
+
+  # Subtracting mean over batch for each channel
+  >>> y = x - ein_reduce(x, 'b c h w -> () c () ()', 'mean')
+
+  # Subtracting per-image mean for each channel
+  >>> y = x - ein_reduce(x, 'b c h w -> b c () ()', 'mean')
+
+  ```
+
+  Parameters:
+      tensor: tensor: tensor of any supported library (e.g. numpy.ndarray, tensorflow, pytorch).
+          list of tensors is also accepted, those should be of the same type and shape
+      pattern: string, reduction pattern
+      reduction: one of available reductions ('min', 'max', 'sum', 'mean', 'prod'), case-sensitive
+          alternatively, a callable f(tensor, reduced_axes) -> tensor can be provided.
+          This allows using various reductions, examples: np.max, tf.reduce_logsumexp, torch.var, etc.
+      axes_lengths: any additional specifications for dimensions
+
+  Returns:
+      tensor of the same type as input
+  """
+  try:
+    hashable_axes_lengths = tuple(axes_lengths.items())
+    shape = bnp.shape(tensor)
+    recipe = _prepare_transformation_recipe(pattern, reduction, axes_names=tuple(axes_lengths), ndim=len(shape))
+    return _apply_recipe(recipe,
+                         cast(Tensor, tensor),
+                         reduction_type=reduction,
+                         axes_lengths=hashable_axes_lengths)
+  except EinopsError as e:
+    message = ' Error while processing {}-reduction pattern "{}".'.format(reduction, pattern)
+    if not isinstance(tensor, list):
+      message += "\n Input tensor shape: {}. ".format(shape)
+    else:
+      message += "\n Input is list. "
+    message += "Additional info: {}.".format(axes_lengths)
+    raise EinopsError(message + "\n {}".format(e))
+
+
+def ein_rearrange(tensor: Union[Tensor, List[Tensor]], pattern: str, **axes_lengths) -> Tensor:
+  """
+  ``ein_rearrange`` is a reader-friendly smart element reordering for multidimensional tensors.
+  This operation includes functionality of transpose (axes permutation), reshape (view), squeeze, unsqueeze,
+  stack, concatenate and other operations.
+
+  Examples for rearrange operation:
+
+  ```python
+  # suppose we have a set of 32 images in "h w c" format (height-width-channel)
+  >>> images = [np.random.randn(30, 40, 3) for _ in range(32)]
+
+  # stack along first (batch) axis, output is a single array
+  >>> ein_rearrange(images, 'b h w c -> b h w c').shape
+  (32, 30, 40, 3)
+
+  # concatenate images along height (vertical axis), 960 = 32 * 30
+  >>> ein_rearrange(images, 'b h w c -> (b h) w c').shape
+  (960, 40, 3)
+
+  # concatenated images along horizontal axis, 1280 = 32 * 40
+  >>> ein_rearrange(images, 'b h w c -> h (b w) c').shape
+  (30, 1280, 3)
+
+  # reordered axes to "b c h w" format for deep learning
+  >>> ein_rearrange(images, 'b h w c -> b c h w').shape
+  (32, 3, 30, 40)
+
+  # flattened each image into a vector, 3600 = 30 * 40 * 3
+  >>> ein_rearrange(images, 'b h w c -> b (c h w)').shape
+  (32, 3600)
+
+  # split each image into 4 smaller (top-left, top-right, bottom-left, bottom-right), 128 = 32 * 2 * 2
+  >>> ein_rearrange(images, 'b (h1 h) (w1 w) c -> (b h1 w1) h w c', h1=2, w1=2).shape
+  (128, 15, 20, 3)
+
+  # space-to-depth operation
+  >>> ein_rearrange(images, 'b (h h1) (w w1) c -> b h w (c h1 w1)', h1=2, w1=2).shape
+  (32, 15, 20, 12)
+
+  ```
+
+  When composing axes, C-order enumeration used (consecutive elements have different last axis)
+  Find more examples in einops tutorial.
+
+  Parameters:
+      tensor: tensor of any supported library (e.g. numpy.ndarray, tensorflow, pytorch).
+              list of tensors is also accepted, those should be of the same type and shape
+      pattern: string, rearrangement pattern
+      axes_lengths: any additional specifications for dimensions
+
+  Returns:
+      tensor of the same type as input. If possible, a view to the original tensor is returned.
+
+  """
+  return ein_reduce(tensor, pattern, reduction="rearrange", **axes_lengths)
+
+
+def ein_repeat(tensor: Union[Tensor, List[Tensor]], pattern: str, **axes_lengths) -> Tensor:
+  """
+  ``ein_repeat`` allows reordering elements and repeating them in arbitrary combinations.
+  This operation includes functionality of repeat, tile, broadcast functions.
+
+  Examples for repeat operation:
+
+  ```python
+  # a grayscale image (of shape height x width)
+  >>> image = np.random.randn(30, 40)
+
+  # change it to RGB format by repeating in each channel
+  >>> ein_repeat(image, 'h w -> h w c', c=3).shape
+  (30, 40, 3)
+
+  # repeat image 2 times along height (vertical axis)
+  >>> ein_repeat(image, 'h w -> (repeat h) w', repeat=2).shape
+  (60, 40)
+
+  # repeat image 2 time along height and 3 times along width
+  >>> ein_repeat(image, 'h w -> (h2 h) (w3 w)', h2=2, w3=3).shape
+  (60, 120)
+
+  # convert each pixel to a small square 2x2. Upsample image by 2x
+  >>> ein_repeat(image, 'h w -> (h h2) (w w2)', h2=2, w2=2).shape
+  (60, 80)
+
+  # pixelate image first by downsampling by 2x, then upsampling
+  >>> downsampled = ein_reduce(image, '(h h2) (w w2) -> h w', 'mean', h2=2, w2=2)
+  >>> ein_repeat(downsampled, 'h w -> (h h2) (w w2)', h2=2, w2=2).shape
+  (30, 40)
+
+  ```
+
+  When composing axes, C-order enumeration used (consecutive elements have different last axis)
+  Find more examples in einops tutorial.
+
+  Parameters:
+      tensor: tensor of any supported library (e.g. numpy.ndarray, tensorflow, pytorch).
+          list of tensors is also accepted, those should be of the same type and shape
+      pattern: string, rearrangement pattern
+      axes_lengths: any additional specifications for dimensions
+
+  Returns:
+      Tensor of the same type as input. If possible, a view to the original tensor is returned.
+
+  """
+  return ein_reduce(tensor, pattern, reduction="repeat", **axes_lengths)
+
+
+def ein_shape(x, pattern: str) -> dict:
+  """
+  Parse a tensor shape to dictionary mapping axes names to their lengths.
+
+  ```python
+  # Use underscore to skip the dimension in parsing.
+  >>> x = np.zeros([2, 3, 5, 7])
+  >>> ein_shape(x, 'batch _ h w')
+  {'batch': 2, 'h': 5, 'w': 7}
+
+  # `parse_shape` output can be used to specify axes_lengths for other operations:
+  >>> y = np.zeros([700])
+  >>> ein_rearrange(y, '(b c h w) -> b c h w', **ein_shape(x, 'b _ h w')).shape
+  (2, 10, 5, 7)
+
+  ```
+
+  For symbolic frameworks may return symbols, not integers.
+
+  Parameters:
+      x: tensor of any supported framework
+      pattern: str, space separated names for axes, underscore means skip axis
+
+  Returns:
+      dict, maps axes names to their lengths
+  """
+  exp = ParsedExpression(pattern, allow_underscore=True)
+  shape = bnp.shape(x)
+  if exp.has_composed_axes():
+    raise RuntimeError(f"Can't parse shape with composite axes: {pattern} {shape}")
+  if len(shape) != len(exp.composition):
+    if exp.has_ellipsis:
+      if len(shape) < len(exp.composition) - 1:
+        raise RuntimeError(f"Can't parse shape with this number of dimensions: {pattern} {shape}")
+    else:
+      raise RuntimeError(f"Can't parse shape with different number of dimensions: {pattern} {shape}")
+  if exp.has_ellipsis:
+    ellipsis_idx = exp.composition.index(_ellipsis)
+    composition = (
+        exp.composition[:ellipsis_idx]
+        + ["_"] * (len(shape) - len(exp.composition) + 1)
+        + exp.composition[ellipsis_idx + 1:]
+    )
+  else:
+    composition = exp.composition
+  result = {}
+  for (axis_name,), axis_length in zip(composition, shape):  # type: ignore
+    if axis_name != "_":
+      result[axis_name] = axis_length
+  return result
+
+
+# _enumerate_directions is not exposed in the public API
+def _enumerate_directions(x):
+  """
+  For an n-dimensional tensor, returns tensors to enumerate each axis.
+  ```python
+  x = np.zeros([2, 3, 4]) # or any other tensor
+  i, j, k = _enumerate_directions(x)
+  result = i + 2*j + 3*k
+  ```
+
+  `result[i, j, k] = i + 2j + 3k`, and also has the same shape as result
+  Works very similarly to numpy.ogrid (open indexing grid)
+  """
+  shape = bnp.shape(x)
+  result = []
+  for axis_id, axis_length in enumerate(shape):
+    shape = [1] * len(shape)
+    shape[axis_id] = axis_length
+    result.append(bnp.reshape(bnp.arange(0, axis_length), shape))
+  return result
diff --git a/brainpy/_src/math/einops_parsing.py b/brainpy/_src/math/einops_parsing.py
new file mode 100644
index 000000000..6ce055bdb
--- /dev/null
+++ b/brainpy/_src/math/einops_parsing.py
@@ -0,0 +1,153 @@
+import keyword
+import warnings
+from typing import List, Optional, Set, Tuple, Union
+
+_ellipsis: str = '…'  # NB, this is a single unicode symbol. String is used as it is not a list, but can be iterated
+
+
+class EinopsError(Exception):
+  pass
+
+
+class AnonymousAxis(object):
+  """Important thing: all instances of this class are not equal to each other """
+
+  def __init__(self, value: str):
+    self.value = int(value)
+    if self.value <= 1:
+      if self.value == 1:
+        raise EinopsError('No need to create anonymous axis of length 1. Report this as an issue')
+      else:
+        raise EinopsError('Anonymous axis should have positive length, not {}'.format(self.value))
+
+  def __repr__(self):
+    return "{}-axis".format(str(self.value))
+
+
+class ParsedExpression:
+  """
+  non-mutable structure that contains information about one side of expression (e.g. 'b c (h w)')
+  and keeps some information important for downstream
+  """
+
+  def __init__(self, expression: str, *, allow_underscore: bool = False,
+               allow_duplicates: bool = False):
+    self.has_ellipsis: bool = False
+    self.has_ellipsis_parenthesized: Optional[bool] = None
+    self.identifiers: Set[str] = set()
+    # that's axes like 2, 3, 4 or 5. Axes with size 1 are exceptional and replaced with empty composition
+    self.has_non_unitary_anonymous_axes: bool = False
+    # composition keeps structure of composite axes, see how different corner cases are handled in tests
+    self.composition: List[Union[List[str], str]] = []
+    if '.' in expression:
+      if '...' not in expression:
+        raise EinopsError('Expression may contain dots only inside ellipsis (...)')
+      if str.count(expression, '...') != 1 or str.count(expression, '.') != 3:
+        raise EinopsError(
+          'Expression may contain dots only inside ellipsis (...); only one ellipsis for tensor ')
+      expression = expression.replace('...', _ellipsis)
+      self.has_ellipsis = True
+
+    bracket_group: Optional[List[str]] = None
+
+    def add_axis_name(x):
+      if x in self.identifiers:
+        if not (allow_underscore and x == "_") and not allow_duplicates:
+          raise EinopsError('Indexing expression contains duplicate dimension "{}"'.format(x))
+      if x == _ellipsis:
+        self.identifiers.add(_ellipsis)
+        if bracket_group is None:
+          self.composition.append(_ellipsis)
+          self.has_ellipsis_parenthesized = False
+        else:
+          bracket_group.append(_ellipsis)
+          self.has_ellipsis_parenthesized = True
+      else:
+        is_number = str.isdecimal(x)
+        if is_number and int(x) == 1:
+          # handling the case of anonymous axis of length 1
+          if bracket_group is None:
+            self.composition.append([])
+          else:
+            pass  # no need to think about 1s inside parenthesis
+          return
+        is_axis_name, reason = self.check_axis_name_return_reason(x, allow_underscore=allow_underscore)
+        if not (is_number or is_axis_name):
+          raise EinopsError('Invalid axis identifier: {}\n{}'.format(x, reason))
+        if is_number:
+          x = AnonymousAxis(x)
+        self.identifiers.add(x)
+        if is_number:
+          self.has_non_unitary_anonymous_axes = True
+        if bracket_group is None:
+          self.composition.append([x])
+        else:
+          bracket_group.append(x)
+
+    current_identifier = None
+    for char in expression:
+      if char in '() ':
+        if current_identifier is not None:
+          add_axis_name(current_identifier)
+        current_identifier = None
+        if char == '(':
+          if bracket_group is not None:
+            raise EinopsError("Axis composition is one-level (brackets inside brackets not allowed)")
+          bracket_group = []
+        elif char == ')':
+          if bracket_group is None:
+            raise EinopsError('Brackets are not balanced')
+          self.composition.append(bracket_group)
+          bracket_group = None
+      elif str.isalnum(char) or char in ['_', _ellipsis]:
+        if current_identifier is None:
+          current_identifier = char
+        else:
+          current_identifier += char
+      else:
+        raise EinopsError("Unknown character '{}'".format(char))
+
+    if bracket_group is not None:
+      raise EinopsError('Imbalanced parentheses in expression: "{}"'.format(expression))
+    if current_identifier is not None:
+      add_axis_name(current_identifier)
+
+  def flat_axes_order(self) -> List:
+    result = []
+    for composed_axis in self.composition:
+      assert isinstance(composed_axis, list), 'does not work with ellipsis'
+      for axis in composed_axis:
+        result.append(axis)
+    return result
+
+  def has_composed_axes(self) -> bool:
+    # this will ignore 1 inside brackets
+    for axes in self.composition:
+      if isinstance(axes, list) and len(axes) > 1:
+        return True
+    return False
+
+  @staticmethod
+  def check_axis_name_return_reason(name: str, allow_underscore: bool = False) -> Tuple[bool, str]:
+    if not str.isidentifier(name):
+      return False, 'not a valid python identifier'
+    elif name[0] == '_' or name[-1] == '_':
+      if name == '_' and allow_underscore:
+        return True, ''
+      return False, 'axis name should should not start or end with underscore'
+    else:
+      if keyword.iskeyword(name):
+        warnings.warn("It is discouraged to use axes names that are keywords: {}".format(name), RuntimeWarning)
+      if name in ['axis']:
+        warnings.warn("It is discouraged to use 'axis' as an axis name "
+                      "and will raise an error in future", FutureWarning)
+      return True, ''
+
+  @staticmethod
+  def check_axis_name(name: str) -> bool:
+    """
+    Valid axes names are python identifiers except keywords,
+    and additionally should not start or end with underscore
+    """
+    is_valid, _reason = ParsedExpression.check_axis_name_return_reason(name)
+    return is_valid
diff --git a/brainpy/_src/math/interoperability.py b/brainpy/_src/math/interoperability.py
index 22fe25caf..948538371 100644
--- a/brainpy/_src/math/interoperability.py
+++ b/brainpy/_src/math/interoperability.py
@@ -7,7 +7,10 @@
 
 
 __all__ = [
-  'as_device_array', 'as_jax', 'as_ndarray', 'as_numpy', 'as_variable', 'is_bp_array'
+  'as_device_array', 'as_jax', 'as_ndarray', 'as_numpy', 'as_variable',
+  'from_numpy',
+
+  'is_bp_array'
 ]
 
 
@@ -99,3 +102,8 @@ def as_variable(tensor, dtype=None):
   """
   from .object_transform.variables import Variable
   return Variable(tensor, dtype=dtype)
+
+
+def from_numpy(arr, dtype=None):
+  return as_ndarray(arr, dtype=dtype)
+
diff --git a/brainpy/_src/math/others.py b/brainpy/_src/math/others.py
index f3cf4f516..94aeebb16 100644
--- a/brainpy/_src/math/others.py
+++ b/brainpy/_src/math/others.py
@@ -1,22 +1,27 @@
 # -*- coding: utf-8 -*-
 
 
-from typing import Optional
+from typing import Optional, Union
 
+import jax
 import jax.numpy as jnp
 from jax.tree_util import tree_map
 
 from brainpy import check, tools
 from .compat_numpy import fill_diagonal
 from .environment import get_dt, get_int
-from .ndarray import Array
 from .interoperability import as_jax
+from .ndarray import Array
 
 __all__ = [
   'shared_args_over_time',
   'remove_diag',
   'clip_by_norm',
   'exprel',
+  'is_float_type',
+  # 'reduce',
+  'add_axis',
+  'add_axes',
 ]
 
 
@@ -119,3 +124,21 @@ def exprel(x, threshold: float = None):
     else:
       threshold = 1e-5
   return _exprel(x, threshold)
+
+
+def is_float_type(x: Union[Array, jax.Array]):
+  return x.dtype in ("float16", "float32", "float64", "float128", "bfloat16")
+
+
+def add_axis(x: Union[Array, jax.Array], new_position: int):
+  x = as_jax(x)
+  return jnp.expand_dims(x, new_position)
+
+
+def add_axes(x: Union[Array, jax.Array], n_axes, pos2len):
+  x = as_jax(x)
+  repeats = [1] * n_axes
+  for axis_position, axis_length in pos2len.items():
+    x = add_axis(x, axis_position)
+    repeats[axis_position] = axis_length
+  return jnp.tile(x, repeats)
diff --git a/brainpy/_src/math/tests/test_einops.py b/brainpy/_src/math/tests/test_einops.py
new file mode 100644
index 000000000..e6738009e
--- /dev/null
+++ b/brainpy/_src/math/tests/test_einops.py
@@ -0,0 +1,332 @@
+import numpy
+import pytest
+
+import brainpy.math as bm
+from brainpy._src.math.einops import ein_rearrange, ein_reduce, ein_repeat, _enumerate_directions
+from brainpy._src.math.einops_parsing import EinopsError
+
+REDUCTIONS = ("min", "max", "sum", "mean", "prod")
+
+identity_patterns = [
+  "...->...",
+  "a b c d e-> a b c d e",
+  "a b c d e ...-> ... a b c d e",
+  "a b c d e ...-> a ... b c d e",
+  "... a b c d e -> ... a b c d e",
+  "a ... e-> a ... e",
+  "a ... -> a ... ",
+  "a ... c d e -> a (...) c d e",
+]
+
+equivalent_rearrange_patterns = [
+  ("a b c d e -> (a b) c d e", "a b ... -> (a b) ... "),
+  ("a b c d e -> a b (c d) e", "... c d e -> ... (c d) e"),
+  ("a b c d e -> a b c d e", "... -> ... "),
+  ("a b c d e -> (a b c d e)", "... ->  (...)"),
+  ("a b c d e -> b (c d e) a", "a b ... -> b (...) a"),
+  ("a b c d e -> b (a c d) e", "a b ... e -> b (a ...) e"),
+]
+
+equivalent_reduction_patterns = [
+  ("a b c d e -> ", " ... ->  "),
+  ("a b c d e -> (e a)", "a ... e -> (e a)"),
+  ("a b c d e -> d (a e)", " a b c d e ... -> d (a e) "),
+  ("a b c d e -> (a b)", " ... c d e  -> (...) "),
+]
+
+
+def test_collapsed_ellipsis_errors_out():
+  x = numpy.zeros([1, 1, 1, 1, 1])
+  ein_rearrange(x, "a b c d ... ->  a b c ... d")
+  with pytest.raises(EinopsError):
+    ein_rearrange(x, "a b c d (...) ->  a b c ... d")
+
+  ein_rearrange(x, "... ->  (...)")
+  with pytest.raises(EinopsError):
+    ein_rearrange(x, "(...) -> (...)")
+
+
+def test_ellipsis_ops_numpy():
+  x = numpy.arange(2 * 3 * 4 * 5 * 6).reshape([2, 3, 4, 5, 6])
+  for pattern in identity_patterns:
+    assert numpy.array_equal(x, ein_rearrange(x, pattern)), pattern
+
+  for pattern1, pattern2 in equivalent_rearrange_patterns:
+    assert numpy.array_equal(ein_rearrange(x, pattern1), ein_rearrange(x, pattern2))
+
+  for reduction in ["min", "max", "sum"]:
+    for pattern1, pattern2 in equivalent_reduction_patterns:
+      assert numpy.array_equal(ein_reduce(x, pattern1, reduction=reduction),
+                               ein_reduce(x, pattern2, reduction=reduction))
+
+  # now just check coincidence with numpy
+  all_rearrange_patterns = [*identity_patterns]
+  for pattern_pairs in equivalent_rearrange_patterns:
+    all_rearrange_patterns.extend(pattern_pairs)
+
+
+def test_rearrange_consistency_numpy():
+  shape = [1, 2, 3, 5, 7, 11]
+  x = numpy.arange(numpy.prod(shape)).reshape(shape)
+  for pattern in [
+    "a b c d e f -> a b c d e f",
+    "b a c d e f -> a b d e f c",
+    "a b c d e f -> f e d c b a",
+    "a b c d e f -> (f e) d (c b a)",
+    "a b c d e f -> (f e d c b a)",
+  ]:
+    result = ein_rearrange(x, pattern)
+    assert len(numpy.setdiff1d(x, result)) == 0
+    assert result.dtype == x.dtype
+
+  result = ein_rearrange(x, "a b c d e f -> a (b) (c d e) f")
+  assert numpy.array_equal(x.flatten(), result.flatten())
+
+  result = ein_rearrange(x, "a aa aa1 a1a1 aaaa a11 -> a aa aa1 a1a1 aaaa a11")
+  assert numpy.array_equal(x, result)
+
+  result1 = ein_rearrange(x, "a b c d e f -> f e d c b a")
+  result2 = ein_rearrange(x, "f e d c b a -> a b c d e f")
+  assert numpy.array_equal(result1, result2)
+
+  result = ein_rearrange(ein_rearrange(x, "a b c d e f -> (f d) c (e b) a"), "(f d) c (e b) a -> a b c d e f", b=2, d=5)
+  assert numpy.array_equal(x, result)
+
+  sizes = dict(zip("abcdef", shape))
+  temp = ein_rearrange(x, "a b c d e f -> (f d) c (e b) a", **sizes)
+  result = ein_rearrange(temp, "(f d) c (e b) a -> a b c d e f", **sizes)
+  assert numpy.array_equal(x, result)
+
+  x2 = numpy.arange(2 * 3 * 4).reshape([2, 3, 4])
+  result = ein_rearrange(x2, "a b c -> b c a")
+  assert x2[1, 2, 3] == result[2, 3, 1]
+  assert x2[0, 1, 2] == result[1, 2, 0]
+
+
+def test_rearrange_permutations_numpy():
+  # tests random permutation of axes against two independent numpy ways
+  for n_axes in range(1, 10):
+    input = numpy.arange(2 ** n_axes).reshape([2] * n_axes)
+    permutation = numpy.random.permutation(n_axes)
+    left_expression = " ".join("i" + str(axis) for axis in range(n_axes))
+    right_expression = " ".join("i" + str(axis) for axis in permutation)
+    expression = left_expression + " -> " + right_expression
+    result = ein_rearrange(input, expression)
+
+    for pick in numpy.random.randint(0, 2, [10, n_axes]):
+      assert input[tuple(pick)] == result[tuple(pick[permutation])]
+
+  for n_axes in range(1, 10):
+    input = numpy.arange(2 ** n_axes).reshape([2] * n_axes)
+    permutation = numpy.random.permutation(n_axes)
+    left_expression = " ".join("i" + str(axis) for axis in range(n_axes)[::-1])
+    right_expression = " ".join("i" + str(axis) for axis in permutation[::-1])
+    expression = left_expression + " -> " + right_expression
+    result = ein_rearrange(input, expression)
+    assert result.shape == input.shape
+    expected_result = numpy.zeros_like(input)
+    for original_axis, result_axis in enumerate(permutation):
+      expected_result |= ((input >> original_axis) & 1) << result_axis
+
+    assert numpy.array_equal(result, expected_result)
+
+
+def test_reduction_imperatives():
+  for reduction in REDUCTIONS:
+    # slight redundancy for simpler order - numpy version is evaluated multiple times
+    input = numpy.arange(2 * 3 * 4 * 5 * 6, dtype="int64").reshape([2, 3, 4, 5, 6])
+    if reduction in ["mean", "prod"]:
+      input = input / input.astype("float64").mean()
+    test_cases = [
+      ["a b c d e -> ", {}, getattr(input, reduction)()],
+      ["a ... -> ", {}, getattr(input, reduction)()],
+      ["(a1 a2) ... (e1 e2) -> ", dict(a1=1, e2=2), getattr(input, reduction)()],
+      [
+        "a b c d e -> (e c) a",
+        {},
+        getattr(input, reduction)(axis=(1, 3)).transpose(2, 1, 0).reshape([-1, 2]),
+      ],
+      [
+        "a ... c d e -> (e c) a",
+        {},
+        getattr(input, reduction)(axis=(1, 3)).transpose(2, 1, 0).reshape([-1, 2]),
+      ],
+      [
+        "a b c d e ... -> (e c) a",
+        {},
+        getattr(input, reduction)(axis=(1, 3)).transpose(2, 1, 0).reshape([-1, 2]),
+      ],
+      ["a b c d e -> (e c a)", {}, getattr(input, reduction)(axis=(1, 3)).transpose(2, 1, 0).reshape([-1])],
+      ["(a a2) ... -> (a2 a) ...", dict(a2=1), input],
+    ]
+    for pattern, axes_lengths, expected_result in test_cases:
+      result = ein_reduce(bm.from_numpy(input.copy()), pattern, reduction=reduction, **axes_lengths)
+      result = bm.as_numpy(result)
+      print(reduction, pattern, expected_result, result)
+      assert numpy.allclose(result, expected_result), f"Failed at {pattern}"
+
+
+def test_enumerating_directions():
+  for shape in [[], [1], [1, 1, 1], [2, 3, 5, 7]]:
+    x = numpy.arange(numpy.prod(shape)).reshape(shape)
+    axes1 = _enumerate_directions(x)
+    axes2 = _enumerate_directions(bm.from_numpy(x))
+    assert len(axes1) == len(axes2) == len(shape)
+    for ax1, ax2 in zip(axes1, axes2):
+      ax2 = bm.as_numpy(ax2)
+      assert ax1.shape == ax2.shape
+      assert numpy.allclose(ax1, ax2)
+
+
+def test_concatenations_and_stacking():
+  for n_arrays in [1, 2, 5]:
+    shapes = [[], [1], [1, 1], [2, 3, 5, 7], [1] * 6]
+    for shape in shapes:
+      arrays1 = [numpy.arange(i, i + numpy.prod(shape)).reshape(shape) for i in range(n_arrays)]
+      arrays2 = [bm.from_numpy(array) for array in arrays1]
+      result0 = numpy.asarray(arrays1)
+      result1 = ein_rearrange(arrays1, "...->...")
+      result2 = ein_rearrange(arrays2, "...->...")
+      assert numpy.array_equal(result0, result1)
+      assert numpy.array_equal(result1, bm.as_numpy(result2))
+
+      result1 = ein_rearrange(arrays1, "b ... -> ... b")
+      result2 = ein_rearrange(arrays2, "b ... -> ... b")
+      assert numpy.array_equal(result1, bm.as_numpy(result2))
+
+
+def test_gradients_imperatives():
+  # lazy - just checking reductions
+  for reduction in REDUCTIONS:
+    if reduction in ("any", "all"):
+      continue  # non-differentiable ops
+    x = numpy.arange(1, 1 + 2 * 3 * 4).reshape([2, 3, 4]).astype("float32")
+    y0 = bm.from_numpy(x)
+    if not hasattr(y0, "grad"):
+      continue
+
+    y1 = ein_reduce(y0, "a b c -> c a", reduction=reduction)
+    y2 = ein_reduce(y1, "c a -> a c", reduction=reduction)
+    y3 = ein_reduce(y2, "a (c1 c2) -> a", reduction=reduction, c1=2)
+    y4 = ein_reduce(y3, "... -> ", reduction=reduction)
+
+    y4.backward()
+    grad = bm.as_numpy(y0.grad)
+
+
+def test_tiling_imperatives():
+  input = numpy.arange(2 * 3 * 5, dtype="int64").reshape([2, 1, 3, 1, 5])
+  test_cases = [
+    (1, 1, 1, 1, 1),
+    (1, 2, 1, 3, 1),
+    (3, 1, 1, 4, 1),
+  ]
+  for repeats in test_cases:
+    expected = numpy.tile(input, repeats)
+    converted = bm.from_numpy(input)
+    repeated = bm.tile(converted, repeats)
+    result = bm.as_numpy(repeated)
+    assert numpy.array_equal(result, expected)
+
+
+repeat_test_cases = [
+  # all assume that input has shape [2, 3, 5]
+  ("a b c -> c a b", dict()),
+  ("a b c -> (c copy a b)", dict(copy=2, a=2, b=3, c=5)),
+  ("a b c -> (a copy) b c ", dict(copy=1)),
+  ("a b c -> (c a) (copy1 b copy2)", dict(a=2, copy1=1, copy2=2)),
+  ("a ...  -> a ... copy", dict(copy=4)),
+  ("... c -> ... (copy1 c copy2)", dict(copy1=1, copy2=2)),
+  ("...  -> ... ", dict()),
+  (" ...  -> copy1 ... copy2 ", dict(copy1=2, copy2=3)),
+  ("a b c  -> copy1 a copy2 b c () ", dict(copy1=2, copy2=1)),
+]
+
+
+def check_reversion(x, repeat_pattern, **sizes):
+  """Checks repeat pattern by running reduction"""
+  left, right = repeat_pattern.split("->")
+  reduce_pattern = right + "->" + left
+  repeated = ein_repeat(x, repeat_pattern, **sizes)
+  reduced_min = ein_reduce(repeated, reduce_pattern, reduction="min", **sizes)
+  reduced_max = ein_reduce(repeated, reduce_pattern, reduction="max", **sizes)
+  assert numpy.array_equal(x, reduced_min)
+  assert numpy.array_equal(x, reduced_max)
+
+
+def test_repeat_numpy():
+  # check repeat vs reduce. Repeat works ok if reverse reduction with min and max work well
+  x = numpy.arange(2 * 3 * 5).reshape([2, 3, 5])
+  x1 = ein_repeat(x, "a b c -> copy a b c ", copy=1)
+  assert numpy.array_equal(x[None], x1)
+  for pattern, axis_dimensions in repeat_test_cases:
+    check_reversion(x, pattern, **axis_dimensions)
+
+
+test_cases_repeat_anonymous = [
+  # all assume that input has shape [1, 2, 4, 6]
+  ("a b c d -> c a d b", dict()),
+  ("a b c d -> (c 2 d a b)", dict(a=1, c=4, d=6)),
+  ("1 b c d -> (d copy 1) 3 b c ", dict(copy=3)),
+  ("1 ...  -> 3 ... ", dict()),
+  ("() ... d -> 1 (copy1 d copy2) ... ", dict(copy1=2, copy2=3)),
+  ("1 b c d -> (1 1) (1 b) 2 c 3 d (1 1)", dict()),
+]
+
+
+def test_anonymous_axes():
+  x = numpy.arange(1 * 2 * 4 * 6).reshape([1, 2, 4, 6])
+  for pattern, axis_dimensions in test_cases_repeat_anonymous:
+    check_reversion(x, pattern, **axis_dimensions)
+
+
+def test_list_inputs():
+  x = numpy.arange(2 * 3 * 4 * 5 * 6).reshape([2, 3, 4, 5, 6])
+
+  assert numpy.array_equal(
+    ein_rearrange(list(x), "... -> (...)"),
+    ein_rearrange(x, "... -> (...)"),
+  )
+  assert numpy.array_equal(
+    ein_reduce(list(x), "a ... e -> (...)", "min"),
+    ein_reduce(x, "a ... e -> (...)", "min"),
+  )
+  assert numpy.array_equal(
+    ein_repeat(list(x), "...  -> b (...)", b=3),
+    ein_repeat(x, "...  -> b (...)", b=3),
+  )
+
+
+def bit_count(x):
+  return sum((x >> i) & 1 for i in range(20))
+
+
+def test_reduction_imperatives_booleans():
+  """Checks that any/all reduction works in all frameworks"""
+  x_np = numpy.asarray([(bit_count(x) % 2) == 0 for x in range(2 ** 6)]).reshape([2] * 6)
+
+  for axis in range(6):
+    expected_result_any = numpy.any(x_np, axis=axis, keepdims=True)
+    expected_result_all = numpy.all(x_np, axis=axis, keepdims=True)
+    assert not numpy.array_equal(expected_result_any, expected_result_all)
+
+    axes = list("abcdef")
+    axes_in = list(axes)
+    axes_out = list(axes)
+    axes_out[axis] = "1"
+    pattern = (" ".join(axes_in)) + " -> " + (" ".join(axes_out))
+
+    res_any = ein_reduce(bm.from_numpy(x_np), pattern, reduction="any")
+    res_all = ein_reduce(bm.from_numpy(x_np), pattern, reduction="all")
+
+    assert numpy.array_equal(expected_result_any, bm.as_numpy(res_any))
+    assert numpy.array_equal(expected_result_all, bm.as_numpy(res_all))
+
+  # expected result: any/all
+  expected_result_any = numpy.any(x_np, axis=(0, 1), keepdims=True)
+  expected_result_all = numpy.all(x_np, axis=(0, 1), keepdims=True)
+  pattern = "a b ... -> 1 1 ..."
+  res_any = ein_reduce(bm.from_numpy(x_np), pattern, reduction="any")
+  res_all = ein_reduce(bm.from_numpy(x_np), pattern, reduction="all")
+  assert numpy.array_equal(expected_result_any, bm.as_numpy(res_any))
+  assert numpy.array_equal(expected_result_all, bm.as_numpy(res_all))
diff --git a/brainpy/_src/math/tests/test_einops_parsing.py b/brainpy/_src/math/tests/test_einops_parsing.py
new file mode 100644
index 000000000..069c7bbac
--- /dev/null
+++ b/brainpy/_src/math/tests/test_einops_parsing.py
@@ -0,0 +1,111 @@
+import pytest
+
+from brainpy._src.math.einops_parsing import EinopsError, ParsedExpression, AnonymousAxis, _ellipsis
+
+
+class AnonymousAxisPlaceholder:
+  def __init__(self, value: int):
+    self.value = value
+    assert isinstance(self.value, int)
+
+  def __eq__(self, other):
+    return isinstance(other, AnonymousAxis) and self.value == other.value
+
+
+def test_anonymous_axes():
+  a, b = AnonymousAxis('2'), AnonymousAxis('2')
+  assert a != b
+  c, d = AnonymousAxisPlaceholder(2), AnonymousAxisPlaceholder(3)
+  assert a == c and b == c
+  assert a != d and b != d
+  assert [a, 2, b] == [c, 2, c]
+
+
+def test_elementary_axis_name():
+  for name in ['a', 'b', 'h', 'dx', 'h1', 'zz', 'i9123', 'somelongname',
+               'Alex', 'camelCase', 'u_n_d_e_r_score', 'unreasonablyLongAxisName']:
+    assert ParsedExpression.check_axis_name(name)
+
+  for name in ['', '2b', '12', '_startWithUnderscore', 'endWithUnderscore_', '_', '...', _ellipsis]:
+    assert not ParsedExpression.check_axis_name(name)
+
+
+def test_invalid_expressions():
+  # double ellipsis should raise an error
+  ParsedExpression('... a b c d')
+  with pytest.raises(EinopsError):
+    ParsedExpression('... a b c d ...')
+  with pytest.raises(EinopsError):
+    ParsedExpression('... a b c (d ...)')
+  with pytest.raises(EinopsError):
+    ParsedExpression('(... a) b c (d ...)')
+
+  # double/missing/enclosed parenthesis
+  ParsedExpression('(a) b c (d ...)')
+  with pytest.raises(EinopsError):
+    ParsedExpression('(a)) b c (d ...)')
+  with pytest.raises(EinopsError):
+    ParsedExpression('(a b c (d ...)')
+  with pytest.raises(EinopsError):
+    ParsedExpression('(a) (()) b c (d ...)')
+  with pytest.raises(EinopsError):
+    ParsedExpression('(a) ((b c) (d ...))')
+
+  # invalid identifiers
+  ParsedExpression('camelCase under_scored cApiTaLs ß ...')
+  with pytest.raises(EinopsError):
+    ParsedExpression('1a')
+  with pytest.raises(EinopsError):
+    ParsedExpression('_pre')
+  with pytest.raises(EinopsError):
+    ParsedExpression('...pre')
+  with pytest.raises(EinopsError):
+    ParsedExpression('pre...')
+
+
+def test_parse_expression():
+  parsed = ParsedExpression('a1  b1   c1    d1')
+  assert parsed.identifiers == {'a1', 'b1', 'c1', 'd1'}
+  assert parsed.composition == [['a1'], ['b1'], ['c1'], ['d1']]
+  assert not parsed.has_non_unitary_anonymous_axes
+  assert not parsed.has_ellipsis
+
+  parsed = ParsedExpression('() () () ()')
+  assert parsed.identifiers == set()
+  assert parsed.composition == [[], [], [], []]
+  assert not parsed.has_non_unitary_anonymous_axes
+  assert not parsed.has_ellipsis
+
+  parsed = ParsedExpression('1 1 1 ()')
+  assert parsed.identifiers == set()
+  assert parsed.composition == [[], [], [], []]
+  assert not parsed.has_non_unitary_anonymous_axes
+  assert not parsed.has_ellipsis
+
+  aap = AnonymousAxisPlaceholder
+
+  parsed = ParsedExpression('5 (3 4)')
+  assert len(parsed.identifiers) == 3 and {i.value for i in parsed.identifiers} == {3, 4, 5}
+  assert parsed.composition == [[aap(5)], [aap(3), aap(4)]]
+  assert parsed.has_non_unitary_anonymous_axes
+  assert not parsed.has_ellipsis
+
+  parsed = ParsedExpression('5 1 (1 4) 1')
+  assert len(parsed.identifiers) == 2 and {i.value for i in parsed.identifiers} == {4, 5}
+  assert parsed.composition == [[aap(5)], [], [aap(4)], []]
+
+  parsed = ParsedExpression('name1 ... a1 12 (name2 14)')
+  assert len(parsed.identifiers) == 6
+  assert parsed.identifiers.difference({'name1', _ellipsis, 'a1', 'name2'}).__len__() == 2
+  assert parsed.composition == [['name1'], _ellipsis, ['a1'], [aap(12)], ['name2', aap(14)]]
+  assert parsed.has_non_unitary_anonymous_axes
+  assert parsed.has_ellipsis
+  assert not parsed.has_ellipsis_parenthesized
+
+  parsed = ParsedExpression('(name1 ... a1 12) name2 14')
+  assert len(parsed.identifiers) == 6
+  assert parsed.identifiers.difference({'name1', _ellipsis, 'a1', 'name2'}).__len__() == 2
+  assert parsed.composition == [['name1', _ellipsis, 'a1', aap(12)], ['name2'], [aap(14)]]
+  assert parsed.has_non_unitary_anonymous_axes
+  assert parsed.has_ellipsis
+  assert parsed.has_ellipsis_parenthesized
diff --git a/brainpy/math/__init__.py b/brainpy/math/__init__.py
index cf7a766b4..02f671345 100644
--- a/brainpy/math/__init__.py
+++ b/brainpy/math/__init__.py
@@ -8,6 +8,7 @@
 from .compat_numpy import *
 from .compat_tensorflow import *
 from .compat_pytorch import *
+from .einops import *
 
 # functions
 from .activations import *
diff --git a/brainpy/math/einops.py b/brainpy/math/einops.py
new file mode 100644
index 000000000..5dcb4ce67
--- /dev/null
+++ b/brainpy/math/einops.py
@@ -0,0 +1,6 @@
+from brainpy._src.math.einops import (
+  ein_repeat as ein_repeat,
+  ein_shape as ein_shape,
+  ein_reduce as ein_reduce,
+  ein_rearrange as ein_rearrange,
+)
diff --git a/brainpy/math/interoperability.py b/brainpy/math/interoperability.py
index f6356bca7..6956f9ba2 100644
--- a/brainpy/math/interoperability.py
+++ b/brainpy/math/interoperability.py
@@ -6,6 +6,7 @@
   as_ndarray as as_ndarray,
   as_numpy as as_numpy,
   as_variable as as_variable,
+  from_numpy as from_numpy,
   is_bp_array as is_bp_array,
 )
 
diff --git a/docs/tutorial_math/einops_in_brainpy.ipynb b/docs/tutorial_math/einops_in_brainpy.ipynb
new file mode 100644
index 000000000..2489d6bae
--- /dev/null
+++ b/docs/tutorial_math/einops_in_brainpy.ipynb
@@ -0,0 +1,1509 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Array operations with ``ein_rearrange``, ``ein_reduce``, and ``ein_repeat``\n",
+    "\n",
+    "We don't write \n",
+    "```python\n",
+    "y = x.transpose(0, 2, 3, 1)\n",
+    "```\n",
+    "We write comprehensible code\n",
+    "```python\n",
+    "y = bm.ein_rearrange(x, 'b c h w -> b h w c')\n",
+    "```\n",
+    "\n",
+    "\n",
+    "## What's in this tutorial?\n",
+    "\n",
+    "- fundamentals: reordering, composition and decomposition of axes\n",
+    "- operations: `ein_rearrange`, `ein_reduce`, `ein_repeat`\n",
+    "- how much you can do with a single operation!\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Preparations"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2024-01-09T03:16:51.896023200Z",
+     "start_time": "2024-01-09T03:16:49.966551200Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "# Examples are given for numpy. This code also setups ipython/jupyter\n",
+    "# so that numpy arrays in the output are displayed as images\n",
+    "import numpy\n",
+    "\n",
+    "import brainpy.math as bm"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Load a batch of images to play with"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "source": [
+    "Please download [the data](./test_images.npy)."
+   ],
+   "metadata": {
+    "collapsed": false
+   }
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2024-01-09T03:16:51.903282300Z",
+     "start_time": "2024-01-09T03:16:51.898250400Z"
+    }
+   },
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "(6, 96, 96, 3) float64\n"
+     ]
+    }
+   ],
+   "source": [
+    "ims = numpy.load('./test_images.npy', allow_pickle=False)\n",
+    "# There are 6 images of shape 96x96 with 3 color channels packed into tensor\n",
+    "print(ims.shape, ims.dtype)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {
+    "pycharm": {
+     "name": "#%%\n"
+    },
+    "ExecuteTime": {
+     "end_time": "2024-01-09T03:16:51.910514400Z",
+     "start_time": "2024-01-09T03:16:51.905419300Z"
+    }
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": "(96, 96, 3)"
+     },
+     "execution_count": 3,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# display the first image (whole 4d tensor can't be rendered)\n",
+    "ims[0].shape"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {
+    "pycharm": {
+     "name": "#%%\n"
+    },
+    "ExecuteTime": {
+     "end_time": "2024-01-09T03:16:51.916049400Z",
+     "start_time": "2024-01-09T03:16:51.912295Z"
+    }
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": "(96, 96, 3)"
+     },
+     "execution_count": 4,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# second image in a batch\n",
+    "ims[1].shape"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2024-01-09T03:16:51.987415500Z",
+     "start_time": "2024-01-09T03:16:51.917288700Z"
+    }
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": "(96, 96, 3)"
+     },
+     "execution_count": 5,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# rearrange, as its name suggests, rearranges elements\n",
+    "# below we swapped height and width.\n",
+    "# In other words, transposed first two axes (dimensions)\n",
+    "bm.ein_rearrange(ims[0], 'h w c -> w h c').shape"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Composition of axes\n",
+    "transposition is very common and useful, but let's move to other capabilities provided by einops"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2024-01-09T03:16:52.001062900Z",
+     "start_time": "2024-01-09T03:16:51.984159900Z"
+    }
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": "(576, 96, 3)"
+     },
+     "execution_count": 6,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# einops allows seamlessly composing batch and height to a new height dimension\n",
+    "# We just rendered all images by collapsing to 3d tensor!\n",
+    "bm.ein_rearrange(ims, 'b h w c -> (b h) w c').shape"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2024-01-09T03:16:52.043645400Z",
+     "start_time": "2024-01-09T03:16:52.002184500Z"
+    }
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": "(96, 576, 3)"
+     },
+     "execution_count": 7,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# or compose a new dimension of batch and width\n",
+    "bm.ein_rearrange(ims, 'b h w c -> h (b w) c').shape"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "metadata": {
+    "pycharm": {
+     "name": "#%%\n"
+    },
+    "ExecuteTime": {
+     "end_time": "2024-01-09T03:16:52.044717500Z",
+     "start_time": "2024-01-09T03:16:52.032578100Z"
+    }
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": "(96, 576, 3)"
+     },
+     "execution_count": 8,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# resulting dimensions are computed very simply\n",
+    "# length of newly composed axis is a product of components\n",
+    "# [6, 96, 96, 3] -> [96, (6 * 96), 3]\n",
+    "bm.ein_rearrange(ims, 'b h w c -> h (b w) c').shape"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "metadata": {
+    "pycharm": {
+     "name": "#%%\n"
+    },
+    "ExecuteTime": {
+     "end_time": "2024-01-09T03:16:52.059635400Z",
+     "start_time": "2024-01-09T03:16:52.039293900Z"
+    }
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": "(165888,)"
+     },
+     "execution_count": 9,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# we can compose more than two axes. \n",
+    "# let's flatten 4d array into 1d, resulting array has as many elements as the original\n",
+    "bm.ein_rearrange(ims, 'b h w c -> (b h w c)').shape"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "pycharm": {
+     "name": "#%% md\n"
+    }
+   },
+   "source": [
+    "## Decomposition of axis"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 10,
+   "metadata": {
+    "pycharm": {
+     "name": "#%%\n"
+    },
+    "ExecuteTime": {
+     "end_time": "2024-01-09T03:16:52.104413Z",
+     "start_time": "2024-01-09T03:16:52.056324200Z"
+    }
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": "(2, 3, 96, 96, 3)"
+     },
+     "execution_count": 10,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# decomposition is the inverse process - represent an axis as a combination of new axes\n",
+    "# several decompositions possible, so b1=2 is to decompose 6 to b1=2 and b2=3\n",
+    "bm.ein_rearrange(ims, '(b1 b2) h w c -> b1 b2 h w c ', b1=2).shape"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 11,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2024-01-09T03:16:52.136340300Z",
+     "start_time": "2024-01-09T03:16:52.073847300Z"
+    }
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": "(192, 288, 3)"
+     },
+     "execution_count": 11,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# finally, combine composition and decomposition:\n",
+    "bm.ein_rearrange(ims, '(b1 b2) h w c -> (b1 h) (b2 w) c ', b1=2).shape"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 12,
+   "metadata": {
+    "pycharm": {
+     "name": "#%%\n"
+    },
+    "ExecuteTime": {
+     "end_time": "2024-01-09T03:16:52.165079200Z",
+     "start_time": "2024-01-09T03:16:52.106539200Z"
+    }
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": "(288, 192, 3)"
+     },
+     "execution_count": 12,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# slightly different composition: b1 is merged with width, b2 with height\n",
+    "# ... so letters are ordered by w then by h\n",
+    "bm.ein_rearrange(ims, '(b1 b2) h w c -> (b2 h) (b1 w) c ', b1=2).shape"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 13,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2024-01-09T03:16:52.199903Z",
+     "start_time": "2024-01-09T03:16:52.144629900Z"
+    }
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": "(192, 288, 3)"
+     },
+     "execution_count": 13,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# move part of width dimension to height. \n",
+    "# we should call this width-to-height as image width shrunk by 2 and height doubled. \n",
+    "# but all pixels are the same!\n",
+    "# Can you write reverse operation (height-to-width)?\n",
+    "bm.ein_rearrange(ims, 'b h (w w2) c -> (h w2) (b w) c', w2=2).shape"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Order of axes matters"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 14,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2024-01-09T03:16:52.200972800Z",
+     "start_time": "2024-01-09T03:16:52.190142300Z"
+    }
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": "(96, 576, 3)"
+     },
+     "execution_count": 14,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# compare with the next example\n",
+    "bm.ein_rearrange(ims, 'b h w c -> h (b w) c').shape"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 15,
+   "metadata": {
+    "pycharm": {
+     "name": "#%%\n"
+    },
+    "ExecuteTime": {
+     "end_time": "2024-01-09T03:16:52.250337300Z",
+     "start_time": "2024-01-09T03:16:52.196592800Z"
+    }
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": "(96, 576, 3)"
+     },
+     "execution_count": 15,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# order of axes in composition is different\n",
+    "# rule is just as for digits in the number: leftmost digit is the most significant, \n",
+    "# while neighboring numbers differ in the rightmost axis.\n",
+    "\n",
+    "# you can also think of this as lexicographic sort\n",
+    "bm.ein_rearrange(ims, 'b h w c -> h (w b) c').shape"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 16,
+   "metadata": {
+    "pycharm": {
+     "name": "#%%\n"
+    },
+    "ExecuteTime": {
+     "end_time": "2024-01-09T03:16:52.277698500Z",
+     "start_time": "2024-01-09T03:16:52.228269800Z"
+    }
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": "(96, 576, 3)"
+     },
+     "execution_count": 16,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# what if b1 and b2 are reordered before composing to width?\n",
+    "bm.ein_rearrange(ims, '(b1 b2) h w c -> h (b1 b2 w) c ', b1=2).shape "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "outputs": [
+    {
+     "data": {
+      "text/plain": "(96, 576, 3)"
+     },
+     "execution_count": 17,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "bm.ein_rearrange(ims, '(b1 b2) h w c -> h (b2 b1 w) c ', b1=2).shape "
+   ],
+   "metadata": {
+    "collapsed": false,
+    "ExecuteTime": {
+     "end_time": "2024-01-09T03:16:52.314368100Z",
+     "start_time": "2024-01-09T03:16:52.262594800Z"
+    }
+   },
+   "execution_count": 17
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "pycharm": {
+     "name": "#%% md\n"
+    }
+   },
+   "source": [
+    "## Meet einops.reduce\n",
+    "\n",
+    "In einops-land you don't need to guess what happened\n",
+    "```python\n",
+    "x.mean(-1)\n",
+    "```\n",
+    "Because you write what the operation does\n",
+    "```python\n",
+    "bm.ein_reduce(x, 'b h w c -> b h w', 'mean')\n",
+    "```\n",
+    "\n",
+    "if axis is not present in the output — you guessed it — axis was reduced."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 18,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2024-01-09T03:16:52.354728900Z",
+     "start_time": "2024-01-09T03:16:52.298014600Z"
+    }
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": "(96, 96, 3)"
+     },
+     "execution_count": 18,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# average over batch\n",
+    "bm.ein_reduce(ims, 'b h w c -> h w c', 'mean').shape"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 19,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2024-01-09T03:16:52.355832600Z",
+     "start_time": "2024-01-09T03:16:52.340237700Z"
+    }
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": "(96, 96, 3)"
+     },
+     "execution_count": 19,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# the previous is identical to familiar:\n",
+    "ims.mean(axis=0).shape\n",
+    "# but is so much more readable"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 20,
+   "metadata": {
+    "pycharm": {
+     "name": "#%%\n"
+    },
+    "ExecuteTime": {
+     "end_time": "2024-01-09T03:16:52.408044400Z",
+     "start_time": "2024-01-09T03:16:52.345070800Z"
+    }
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": "(96, 96)"
+     },
+     "execution_count": 20,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# Example of reducing of several axes \n",
+    "# besides mean, there are also min, max, sum, prod\n",
+    "bm.ein_reduce(ims, 'b h w c -> h w', 'min').shape"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 21,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2024-01-09T03:16:52.438192700Z",
+     "start_time": "2024-01-09T03:16:52.365121Z"
+    }
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": "(48, 288, 3)"
+     },
+     "execution_count": 21,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# this is mean-pooling with 2x2 kernel\n",
+    "# image is split into 2x2 patches, each patch is averaged\n",
+    "bm.ein_reduce(ims, 'b (h h2) (w w2) c -> h (b w) c', 'mean', h2=2, w2=2).shape"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 22,
+   "metadata": {
+    "pycharm": {
+     "name": "#%%\n"
+    },
+    "ExecuteTime": {
+     "end_time": "2024-01-09T03:16:52.466068200Z",
+     "start_time": "2024-01-09T03:16:52.429666600Z"
+    }
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": "(48, 288, 3)"
+     },
+     "execution_count": 22,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# max-pooling is similar\n",
+    "# result is not as smooth as for mean-pooling\n",
+    "bm.ein_reduce(ims, 'b (h h2) (w w2) c -> h (b w) c', 'max', h2=2, w2=2).shape"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 23,
+   "metadata": {
+    "pycharm": {
+     "name": "#%%\n"
+    },
+    "ExecuteTime": {
+     "end_time": "2024-01-09T03:16:52.508614800Z",
+     "start_time": "2024-01-09T03:16:52.453429200Z"
+    }
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": "(288, 192)"
+     },
+     "execution_count": 23,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# yet another example. Can you compute result shape?\n",
+    "bm.ein_reduce(ims, '(b1 b2) h w c -> (b2 h) (b1 w)', 'mean', b1=2).shape"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "pycharm": {
+     "name": "#%%\n"
+    }
+   },
+   "source": [
+    "## Stack and concatenate"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 24,
+   "metadata": {
+    "pycharm": {
+     "name": "#%%\n"
+    },
+    "ExecuteTime": {
+     "end_time": "2024-01-09T03:16:52.509704200Z",
+     "start_time": "2024-01-09T03:16:52.486964100Z"
+    }
+   },
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "<class 'list'> with 6 tensors of shape (96, 96, 3)\n"
+     ]
+    },
+    {
+     "data": {
+      "text/plain": "[(96, 96, 3), (96, 96, 3), (96, 96, 3), (96, 96, 3), (96, 96, 3), (96, 96, 3)]"
+     },
+     "execution_count": 24,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# rearrange can also take care of lists of arrays with the same shape\n",
+    "x = list(ims)\n",
+    "print(type(x), 'with', len(x), 'tensors of shape', x[0].shape)\n",
+    "# that's how we can stack inputs\n",
+    "# \"list axis\" becomes first (\"b\" in this case), and we left it there\n",
+    "res = bm.ein_rearrange(x, 'b h w c -> b h w c')\n",
+    "\n",
+    "[r.shape for r in res]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 25,
+   "metadata": {
+    "pycharm": {
+     "name": "#%%\n"
+    },
+    "ExecuteTime": {
+     "end_time": "2024-01-09T03:16:52.524732200Z",
+     "start_time": "2024-01-09T03:16:52.495686100Z"
+    }
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": "(96, 96, 3, 6)"
+     },
+     "execution_count": 25,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# but new axis can appear in the other place:\n",
+    "bm.ein_rearrange(x, 'b h w c -> h w c b').shape"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 26,
+   "metadata": {
+    "collapsed": false,
+    "pycharm": {
+     "name": "#%%\n"
+    },
+    "ExecuteTime": {
+     "end_time": "2024-01-09T03:16:52.528015200Z",
+     "start_time": "2024-01-09T03:16:52.511870500Z"
+    }
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": "False"
+     },
+     "execution_count": 26,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# that's equivalent to numpy stacking, but written more explicitly\n",
+    "numpy.array_equal(bm.ein_rearrange(x, 'b h w c -> h w c b'), numpy.stack(x, axis=3))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 27,
+   "metadata": {
+    "pycharm": {
+     "name": "#%%\n"
+    },
+    "ExecuteTime": {
+     "end_time": "2024-01-09T03:16:52.586497800Z",
+     "start_time": "2024-01-09T03:16:52.517938100Z"
+    }
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": "(96, 576, 3)"
+     },
+     "execution_count": 27,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# ... or we can concatenate along axes\n",
+    "bm.ein_rearrange(x, 'b h w c -> h (b w) c').shape"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 28,
+   "metadata": {
+    "pycharm": {
+     "name": "#%%\n"
+    },
+    "ExecuteTime": {
+     "end_time": "2024-01-09T03:16:52.589607600Z",
+     "start_time": "2024-01-09T03:16:52.524732200Z"
+    }
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": "False"
+     },
+     "execution_count": 28,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# which is equivalent to concatenation\n",
+    "numpy.array_equal(bm.ein_rearrange(x, 'b h w c -> h (b w) c'), numpy.concatenate(x, axis=1))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Addition or removal of axes\n",
+    "\n",
+    "You can write 1 to create a new axis of length 1. Similarly you can remove such axis.\n",
+    "\n",
+    "There is also a synonym `()` that you can use. That's a composition of zero axes and it also has a unit length."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 29,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2024-01-09T03:16:52.601830300Z",
+     "start_time": "2024-01-09T03:16:52.531696500Z"
+    }
+   },
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "(6, 1, 96, 96, 1, 3)\n",
+      "(6, 96, 96, 3)\n"
+     ]
+    }
+   ],
+   "source": [
+    "x = bm.ein_rearrange(ims, 'b h w c -> b 1 h w 1 c') # functionality of numpy.expand_dims\n",
+    "print(x.shape)\n",
+    "print(bm.ein_rearrange(x, 'b 1 h w 1 c -> b h w c').shape) # functionality of numpy.squeeze"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 30,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2024-01-09T03:16:52.652283400Z",
+     "start_time": "2024-01-09T03:16:52.562649Z"
+    }
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": "(96, 576, 3)"
+     },
+     "execution_count": 30,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# compute max in each image individually, then show a difference \n",
+    "x = bm.ein_reduce(ims, 'b h w c -> b () () c', 'max') - ims\n",
+    "bm.ein_rearrange(x, 'b h w c -> h (b w) c').shape"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Repeating elements\n",
+    "\n",
+    "Third operation we introduce is `repeat`"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 31,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2024-01-09T03:16:52.708988500Z",
+     "start_time": "2024-01-09T03:16:52.634965400Z"
+    }
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": "(96, 5, 96, 3)"
+     },
+     "execution_count": 31,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# repeat along a new axis. New axis can be placed anywhere\n",
+    "bm.ein_repeat(ims[0], 'h w c -> h new_axis w c', new_axis=5).shape"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 32,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2024-01-09T03:16:52.714789300Z",
+     "start_time": "2024-01-09T03:16:52.710069Z"
+    }
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": "(96, 5, 96, 3)"
+     },
+     "execution_count": 32,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# shortcut\n",
+    "bm.ein_repeat(ims[0], 'h w c -> h 5 w c').shape"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 33,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2024-01-09T03:16:52.757633Z",
+     "start_time": "2024-01-09T03:16:52.714789300Z"
+    }
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": "(96, 288, 3)"
+     },
+     "execution_count": 33,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# repeat along w (existing axis)\n",
+    "bm.ein_repeat(ims[0], 'h w c -> h (repeat w) c', repeat=3).shape"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 34,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2024-01-09T03:16:52.853440Z",
+     "start_time": "2024-01-09T03:16:52.757633Z"
+    }
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": "(192, 192, 3)"
+     },
+     "execution_count": 34,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# repeat along two existing axes\n",
+    "bm.ein_repeat(ims[0], 'h w c -> (2 h) (2 w) c').shape"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 35,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2024-01-09T03:16:52.935098900Z",
+     "start_time": "2024-01-09T03:16:52.853440Z"
+    }
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": "(96, 288, 3)"
+     },
+     "execution_count": 35,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# order of axes matters as usual - you can repeat each element (pixel) 3 times \n",
+    "# by changing order in parenthesis\n",
+    "bm.ein_repeat(ims[0], 'h w c -> h (w repeat) c', repeat=3).shape"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Note: `repeat` operation covers functionality identical to `numpy.repeat`, `numpy.tile` and actually more than that."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "pycharm": {
+     "name": "#%% md\n"
+    }
+   },
+   "source": [
+    "## Reduce ⇆ repeat\n",
+    "\n",
+    "reduce and repeat are like opposite of each other: first one reduces amount of elements, second one increases.\n",
+    "\n",
+    "In the following example each image is repeated first, then we reduce over new axis to get back original tensor. Notice that operation patterns are \"reverse\" of each other"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 36,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2024-01-09T03:16:53.086847800Z",
+     "start_time": "2024-01-09T03:16:52.936595200Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "repeated = bm.ein_repeat(ims, 'b h w c -> b h new_axis w c', new_axis=2)\n",
+    "reduced = bm.ein_reduce(repeated, 'b h new_axis w c -> b h w c', 'min')\n",
+    "\n",
+    "\n",
+    "assert bm.allclose(ims, reduced)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Fancy examples in random order\n",
+    "\n",
+    "(a.k.a. mad designer gallery)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 37,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2024-01-09T03:16:53.124865300Z",
+     "start_time": "2024-01-09T03:16:53.089018Z"
+    }
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": "(192, 288, 3)"
+     },
+     "execution_count": 37,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# interweaving pixels of different pictures\n",
+    "# all letters are observable\n",
+    "bm.ein_rearrange(ims, '(b1 b2) h w c -> (h b1) (w b2) c ', b1=2).shape"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 38,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2024-01-09T03:16:53.139588200Z",
+     "start_time": "2024-01-09T03:16:53.123858300Z"
+    }
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": "(192, 288, 3)"
+     },
+     "execution_count": 38,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# interweaving along vertical for couples of images\n",
+    "bm.ein_rearrange(ims, '(b1 b2) h w c -> (h b1) (b2 w) c', b1=2).shape"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 39,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2024-01-09T03:16:53.186247700Z",
+     "start_time": "2024-01-09T03:16:53.140592800Z"
+    }
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": "(96, 288, 3)"
+     },
+     "execution_count": 39,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# interweaving lines for couples of images\n",
+    "# exercise: achieve the same result without einops in your favourite framework\n",
+    "bm.ein_reduce(ims, '(b1 b2) h w c -> h (b2 w) c', 'max', b1=2).shape"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 40,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2024-01-09T03:16:53.232730900Z",
+     "start_time": "2024-01-09T03:16:53.178674500Z"
+    }
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": "(144, 288)"
+     },
+     "execution_count": 40,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# color can be also composed into dimension\n",
+    "# ... while image is downsampled\n",
+    "bm.ein_reduce(ims, 'b (h 2) (w 2) c -> (c h) (b w)', 'mean').shape"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 41,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2024-01-09T03:16:53.302503900Z",
+     "start_time": "2024-01-09T03:16:53.236495100Z"
+    }
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": "(24, 192)"
+     },
+     "execution_count": 41,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# disproportionate resize\n",
+    "bm.ein_reduce(ims, 'b (h 4) (w 3) c -> (h) (b w)', 'mean').shape"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 42,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2024-01-09T03:16:53.365480400Z",
+     "start_time": "2024-01-09T03:16:53.303630100Z"
+    }
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": "(48, 576)"
+     },
+     "execution_count": 42,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# spilt each image in two halves, compute mean of the two\n",
+    "bm.ein_reduce(ims, 'b (h1 h2) w c -> h2 (b w)', 'mean', h1=2).shape"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 43,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2024-01-09T03:16:53.413333100Z",
+     "start_time": "2024-01-09T03:16:53.364414400Z"
+    }
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": "(96, 576, 3)"
+     },
+     "execution_count": 43,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# split in small patches and transpose each patch\n",
+    "bm.ein_rearrange(ims, 'b (h1 h2) (w1 w2) c -> (h1 w2) (b w1 h2) c', h2=8, w2=8).shape"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 44,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2024-01-09T03:16:53.499062100Z",
+     "start_time": "2024-01-09T03:16:53.407925200Z"
+    }
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": "(96, 576, 3)"
+     },
+     "execution_count": 44,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# stop me someone!\n",
+    "bm.ein_rearrange(ims, 'b (h1 h2 h3) (w1 w2 w3) c -> (h1 w2 h3) (b w1 h2 w3) c', h2=2, w2=2, w3=2, h3=2).shape"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 45,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2024-01-09T03:16:53.546329400Z",
+     "start_time": "2024-01-09T03:16:53.459186600Z"
+    }
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": "(192, 288, 3)"
+     },
+     "execution_count": 45,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "bm.ein_rearrange(ims, '(b1 b2) (h1 h2) (w1 w2) c -> (h1 b1 h2) (w1 b2 w2) c', h1=3, w1=3, b2=3).shape"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 46,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2024-01-09T03:16:53.587041200Z",
+     "start_time": "2024-01-09T03:16:53.505732100Z"
+    }
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": "(96, 576, 3)"
+     },
+     "execution_count": 46,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# patterns can be arbitrarily complicated\n",
+    "bm.ein_reduce(ims, '(b1 b2) (h1 h2 h3) (w1 w2 w3) c -> (h1 w1 h3) (b1 w2 h2 w3 b2) c', 'mean', \n",
+    "       h2=2, w1=2, w3=2, h3=2, b2=2).shape"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 47,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2024-01-09T03:16:53.608899300Z",
+     "start_time": "2024-01-09T03:16:53.556416400Z"
+    }
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": "(96, 576, 3)"
+     },
+     "execution_count": 47,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# subtract background in each image individually and normalize\n",
+    "# pay attention to () - this is composition of 0 axis, a dummy axis with 1 element.\n",
+    "im2 = bm.ein_reduce(ims, 'b h w c -> b () () c', 'max') - ims\n",
+    "im2 /= bm.ein_reduce(im2, 'b h w c -> b () () c', 'max')\n",
+    "bm.ein_rearrange(im2, 'b h w c -> h (b w) c').shape"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 48,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2024-01-09T03:16:53.742684900Z",
+     "start_time": "2024-01-09T03:16:53.578494900Z"
+    }
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": "(96, 576, 3)"
+     },
+     "execution_count": 48,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# pixelate: first downscale by averaging, then upscale back using the same pattern\n",
+    "averaged = bm.ein_reduce(ims, 'b (h h2) (w w2) c -> b h w c', 'mean', h2=6, w2=8)\n",
+    "bm.ein_repeat(averaged, 'b h w c -> (h h2) (b w w2) c', h2=6, w2=8).shape"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 49,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2024-01-09T03:16:53.783169200Z",
+     "start_time": "2024-01-09T03:16:53.742684900Z"
+    }
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": "(96, 576, 3)"
+     },
+     "execution_count": 49,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "bm.ein_rearrange(ims, 'b h w c -> w (b h) c').shape"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 50,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2024-01-09T03:16:53.827528Z",
+     "start_time": "2024-01-09T03:16:53.765960100Z"
+    }
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": "(96, 576)"
+     },
+     "execution_count": 50,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# let's bring color dimension as part of horizontal axis\n",
+    "# at the same time horizontal axis is downsampled by 2x\n",
+    "bm.ein_reduce(ims, 'b (h h2) (w w2) c -> (h w2) (b w c)', 'mean', h2=3, w2=3).shape"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Ok, numpy is fun, but how do I use einops with some other framework?\n",
+    "\n",
+    "If that's what you've done with `ims` being numpy array:\n",
+    "```python\n",
+    "bm.ein_rearrange(ims, 'b h w c -> w (b h) c')\n",
+    "```\n",
+    "That's how you adapt the code for other frameworks:\n",
+    "\n",
+    "```python\n",
+    "# pytorch:\n",
+    "bm.ein_rearrange(ims, 'b h w c -> w (b h) c')\n",
+    "# tensorflow:\n",
+    "bm.ein_rearrange(ims, 'b h w c -> w (b h) c')\n",
+    "# chainer:\n",
+    "bm.ein_rearrange(ims, 'b h w c -> w (b h) c')\n",
+    "# gluon:\n",
+    "bm.ein_rearrange(ims, 'b h w c -> w (b h) c')\n",
+    "# cupy:\n",
+    "bm.ein_rearrange(ims, 'b h w c -> w (b h) c')\n",
+    "# jax:\n",
+    "bm.ein_rearrange(ims, 'b h w c -> w (b h) c')\n",
+    "\n",
+    "...well, you got the idea.\n",
+    "```\n",
+    "\n",
+    "Einops allows backpropagation as if all operations were native to framework.\n",
+    "Operations do not change when moving to another framework - einops notation is universal"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "pycharm": {
+     "name": "#%% md\n"
+    }
+   },
+   "source": [
+    "# Summary\n",
+    "\n",
+    "- `rearrange` doesn't change number of elements and covers different numpy functions (like `transpose`, `reshape`, `stack`, `concatenate`,  `squeeze` and `expand_dims`)\n",
+    "- `reduce` combines same reordering syntax with reductions (`mean`, `min`, `max`, `sum`, `prod`, and any others)\n",
+    "- `repeat` additionally covers repeating and tiling\n",
+    "- composition and decomposition of axes are a corner stone, they can and should be used together\n"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.9.0"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
diff --git a/docs/tutorial_math/index.rst b/docs/tutorial_math/index.rst
index 6ad69939d..d5b764761 100644
--- a/docs/tutorial_math/index.rst
+++ b/docs/tutorial_math/index.rst
@@ -8,3 +8,4 @@ Math Foundation
    control_flows
    Numpy_like_Operations.ipynb
    Dedicated_Operators.ipynb
+   einops_in_brainpy.ipynb
diff --git a/docs/tutorial_math/test_images.npy b/docs/tutorial_math/test_images.npy
new file mode 100644
index 0000000000000000000000000000000000000000..bbff7bd9b6195476d50b664bdbdda7f6a798774d
GIT binary patch
literal 1327232
zcmeIbN01%YwXVsY_DB>ZQKHpM-4S=dC-F*@o^;o&XYU0?QKCq>xY6&v=s`z6UU|w7
zIRrqMB4LA}U@#aA1>s0|L+^<K$^r$Tyg@liq7}_;?7jBCaZ$xYxvZ)@C%-=exgt*;
zWSz5e?Qj2;xpVIa{`+_T+wFJ%Z&&^Ds^5O;fBf6e{_mfD=^J1C(l`I-ZD0DQFMjEN
z{;!|@{HOo=|NZ2D{q%qQTX_Bd`q%&WXa9Eb^`HGe|9bDgUEKdazWPsJ{CfKLmH+df
zzWDnuE?EEjRZ(x+*`A<@9C_XIK_18h?E*f)0XP5$-~b$e18@KizyUY_2jBo4fCF#<
z4!{9800-az9DoCG01m(bH~<IW033h=Z~zX#0XP5$-~b$e18@KizyUY_2jBo4fCF#<
z4!{9800-az9DoCG01m(bH~<IW033h=Z~zX#0XP5$-~b$e18@KizyUY_2jBo4fCF#<
z4!{9800-az9DoCG01m(bH~<IW033h=Z~zX#0XP5$-~b$e18@KizyUY_2jBo4fCF#<
z4!{9800-az9DoCG01m(bH~<IW033h=Z~zX#0XP5$-~b$e18@KizyUY_2jBo4fCF#<
z4!{9800-az9DoCG01m(bH~<IW033h=Z~zX#0XP5$-~b$e18@KizyUY_2jBo4fCF#<
z4!{9800-az9DoCG01m(bH~<IW033h=Z~zX#0XP5$-~b$e18@KizyUY_2jBo4fCF#<
z4!{9800-az9B>ZY95rn>L`~#)Thyp^v<D64cyCm%ZQB}jHplNqwcB^1PUZOds26SH
zhM?zi{9M$N_L-<XGvc4^(@_uS=*PL=J~=*We!NF!#M#IHY}9jhP1G|vo;n$2S4Xwx
z_=Bi3_Fzj;SB`f^4cZ-1{W)&g95im-T|uKcQg8G-^gHxBZ~zX#fu(bx_k7T_9X=E^
zk>k7X1P$BXJwbyxc6LPBhV?;hIl7;i?H|0IvMnnE53CD%+SaZKI+o-4b5YiFCMrkw
zUE{WYU(jfd@DqN*&!zK|c0@a(9aX0t9XlE{ZCkelP2}j;M+Vb%kp3Lq$2@E5R~5_Z
zuQ|@@X5dp%f3lge=sMfxhk`EHd+!E~<Ot{CJe;o%=V_0$$2zyiFGP*mrj0?vIX?7Y
z&`Wl0)C)PTZZFJlZ<GwUPhOwyBVLL(MV+(jqt4{`A0G!z*@qtljpqoT;d7n&O#M^;
z)c^edUMRY5Q#4NIlgF>zZ;jH=@;RF$1J^~pY%QCDUdZuxQ4@Axf6!Qt^b7P0^b7Ov
z7icfEmu6}&M~?(e*`b3$<2m|$r<2xtA}Yt`=q33Ub<2SJ^L=R?`I9;N^`mk7b<}8%
za30P#6X&UC>bcJKeC}+}q<!#@pbI(rxSUr__o)=?)<5!n>zM)fKL^wEcl&bmaogjj
z`?`wsBlIJ6?nkJ9>c7hMUo`%~$4k7P#!0-AV?BF4zBL)Rat3@H+w<vp!e?{z=h~*y
z^K9ce`q8+LKpyBHtK2_QpUqvLJ`YTLn*Zfgj#s8VKP3ZXU_Ke}>!aOi9Diqy{(td5
z<e|Cy7tSx|w+jB&$LI9h&L@I;bDR&uc@G&N1LbDG=aC$>uSN~$=)VPDL>}ljtI%&!
zhmBQ-8n5&Hs8LhYJi<PY%8oRj$_}2lDz{tXXk=hs8F2seVj5Tde2zu)dHO-|MdYEe
z`XkO4=W8W>?7Fo<t+r!(&>K0<3(CBS43L3U&A>0ChV9m<!5sZV;eW^j{pm{fr_@_D
z>+Q;%TvdeE;+Ok6VUL^7AM;p_e!Z<LtrOOsqt@Zry(_5S6pOCc{U~bEeiL;e$Dc$^
z+c%>oa{PVNq-p+R#iHZ?B5EkDPd1Qa(dYN;hEJq*(H_sS7AJcJzof52_eRR^t2z2d
z!2ghkYW7E*7tYH{`cgf2&OFKGbWTO{F&51mTr@vGvvQ#5`<@NWzNaVoztcI+T(3g=
z<(#eh9p$%CtS8&7?ThafFfx(GW8m|t2d{^HPX&#buXi+*WBq!7a=x*y*J!;{MfZ^-
z*4`F0m?L$#I4*Ux2MyVt-9ZC6YW$gBKiiV7Z){;4+8pIwG><y<xw!fyZ{%$yuUAlS
zi>o(Z&!{b}_jxkMIdEC~AA3K0ZyM))De8F}t*5F;eJ-v(7n(Qup_8paeWqCZf}JA)
zKA+-CX`Q_1b6jZtdH>4PIe8|}_^bI}$MgFvQ)!;T@f`JcvR|iYN%PIN@b~f@{RQ)`
zG!}mtw}0h=*8Rhov4|V=zl~M@oG;E-o&A!oL-=#*{pmVSKR)Nm-NIEx{?NyzOxQo(
z4;sr6&Qy~#zls{Me~uc;agNhp`yczf)5CUW)L@R(c{S^t^TK(lv(G82f9F_7;#XUa
zwFmqh32^SW*&R_Qb1cd;uUC^le~g;8A4W~&sPS2!i+Uu@lk*6F7o>J^OH{Aj7<D#B
z>b#nD&UxXy;D2V{7xU}d#}dB|=a?`@D9C{R?sfmupZt%mo6XKM{An!yT={wMf}*I=
z%RbJ$-PS)4bSg*cys_$>^Tqk9gD>*G(<bd3Q5SNo5wc7Ds%K6IU9j7tMstKG&A}6Y
z?y)`1Pjo8BQj4g852pDqHmnc2wC*bP-yHSNdE`7clP~h?`J?va-k_lz8yLZSbTU8&
zRw)D5MD4b#qju)#ABXFo<N^QCO#Y$K>zQ%RtJG2OYso+pW}s-EhoX7HJq7q8@_;YG
z7p*)8oL9Z3d2<z;5Y2ogGC&4aDFcTO1@Zi&r+?)+&+o-w;jd=5yPh*q=JT0#=eSCk
z<kyk`GSGw>_(9aT?cWzPnxlW!*<9jv`fK`Y+Qo8n;O?ky>+Fciu?Z2(S0V#sV3jg(
zPt*m|`Z<dJao{|8!0+MrR-Oa?yhdkQU#cy~Rl+2{mJE=ACd|NnQDgS*J3+%a`d3|f
z?(uu^8~6>{9sZ{YI~%?d8K_tW@IUws{04r*u~eNrTzOrpib0f*Lk7q|LuNqpi|M|X
zqJNp{@Q}}gPs69xr#&+Mxgom~J~SC<s0@7YkDv?o{ix9#{Y%t`w|rmx75=I^?fu@U
zUfZ@c=xmM+1y?>S86X3d&VbfY^XtuHrsyBKIy~g_;7{-;^=Xe+UazXsNae$k0W#1a
z8Q8upXv(&14jRwVKYo39%lE}!;jgOG-gQ0towtLAOtC?*TDgyT<@MC6bKc=F9EQWR
zM>qfn>db*7hl3_f>vk$u&R8_a2YUa#pi%o})Nqb)2F}!(Gt@uzPyOS!8-)YDj^Q(D
z{leopHi!*bwU2q_^+c%uMyY>3C;c`3HN1iY_2huB3*TiO?NK?dnz3k*U;MeKC(YND
z-IL>wqoz#j04UPF>e;?H@0@r175)kiRFeZw?hP8U^$!FM<k%oKq;el~)5f3<^YP|w
zIrg3pnl`NusYtu1X1n0La9;2&_!c<O+#L9ikAtS{!w-VSbG$C<G4ty~kLFl8!_gog
z=)sntF7xX~9Xb9aYTCXTHIXCjuDRPC=a=(~zrtU^fogJ~yDQ4R7uAttgV>FR`Is%6
zf?lxCMm?9~ZBe7?e&kS&w3lkOm$J^wj_pCy_PeNw9KRMdVsF0{G?=5;Pp5r1sy#=>
zr!{w-lvQsWhdRNB&-USZK1}_vqI30Ln;E?wv0d&A+!Xb?T_5#Ijy`^C$oA|G8p!eJ
zk)SC%bTDW<N4Qrm_x!v}rSmbK<Hij^L-w<%fgHO!qs)Kao})jHusc16urtSn)|KDS
zK4t6H2DRqMc`x_8lTY%Af0*qbnzdbQSRd4853CD1nd3tF%==d+1MX9I*bPxza@6^2
zJrQM!Mb{1d`ofTXC2AnYqV=A(ZV8&O)7?R1ITrbZqU#=ih?=tBMvdpF>mvvD2aTEH
zZ=x=k`-_nr{XXufX}p`F=b_(rJ`vQLW6^T~e*g0MG~WH$9R2<Gr~Aa0;>w)m*}uj8
z&ojyQ9M93;Ssz~Ueeoyw6WV)O9Po9-#_bDHV>y00>fv-<Xh)8-$;>NcfD9~m2HeN4
zN$v1b{A1L#X`Kp1|2%M#Jm8n`OK@PJ99X+1sMpS&jmmMkInVLQ02!En2L5Z*xcPYe
z(H#8<aF0CTfABwWV4)l+x_{xm=!9w9isJl<%{$2e87M0Q{yf95-5NERqaPpckq7(_
z{s#_}#er`{P1-$C7jkqTbR>-vJ(y!zoaZ=XfDFvez`sZJna{g$K1Y8i+#?V8AN&s-
zSZ)qzeQ1Bqe$*8Gxu%0@9O=FsXOo;)$N(8wXa;=TUbktUPDMWs+#?V8AN&s-SS=1{
zy@>t$f<{fz$5%g{#&aIaaiRR@{bZod8Q8up=ymhIKVHeva{~9s1O5m90|%<c0UuXA
zX&?L}=t7PsTceU6I+0_Y`#8Qk87Rs?(Ykg1ez->-@IUw;IM65@C|Xx*_pYG+G_Or>
zjy@00<EHgJ6pQ)`UMB-}$iT-R1x@C88~r}?xcPU7hvWfYgfD^vjn09hdHL>+>bA~~
zs2qKL_w{L>yGyYSeInnnJ{iz_UOpaZIE@b)#{YzFzh{CxkcU_YXUPM3xT-!r(61Lw
znERu#9NkCtq<M|IbG-8Xh5GcTeAAW5fUg_3*R*biV$nF~YokW2y)9@kN9IKg9UDG}
z{@v?<JdlU0YUFdgw_3WdE9ysn88v&}6<;5s+wO|$%+cpfI+o^5I+Ekc^s4-pvNE7`
zM|^#Zmh?PhOCG0oW7KP=cuQ1o+CQ7)ucJoI{q=B;Me}ZzRnHs;U&MW%rS}W8Bia$|
zXzA@p>lHloV9=N;7R_sSZ&a^s+ZuE>NBJ$ke$bxQoo~0I`6`RnlPbF2rEyi>??0OQ
z{YO*3e`!4it<zAnepb<VEkFJN^XnN;<>=Qxo=NM59?#M9^iq1?3pw5y)n<1@oy_sa
zQRmE`M?90G<_-LERKF>HJ!;IJdNOD@$N%^^Xv#kPAZR>CcnA;S;nI0XJE9%Yj%Y`4
z01m(bH~<IW033h=Z~zX#0XP5$-~b$e18@KizyUY_2jBo4fCF#<4!{9800-az9DoCG
z01m(bH~<IW033h=Z~zX#0XP5$-~b$e18@KizyUY_2jBo4fCF#<4!{9800-az9DoCG
z01m(bH~<IW033h=Z~zX#0XP5$-~b$e18@KizyUY_2jBo4fCF#<4!{9800-az9DoCG
z01m(bH~<IW033h=Z~zX#0XP5$-~b$e18@KizyUY_2jBo4fCF#<4!{9800-az9DoCG
z01m(bH~<IW033h=Z~zX#0XP5$-~b$e18@KizyUY_2jBo4fCF#<4!{9800-az9DoCG
z01m(bH~<H#$$_)sWA#iJ=Z8GNX*dlB-~b$e18@KizyUY_2jBo4fCF#<4!{9800&lw
z1Ghzu+L_Zq7fk!7yQ6IDmY^Qf{%50}vzARkFPQdkh<e;U74?K^|MO8TX@84p{|iz3
z)9Z&!`|H*Qwc7dzf=-$CPqqg2*`D1&1E&3NMopM&+O*I0XRjynKpw8TivETEWp(-&
z_yd1dl0WV{ChX>@>EtgatgAE1?uj~`_ua4jDfyO5(Roq*?e$Ry?WU+Bru|<;4W;in
zV%n!Y(jHf`JyLJf8}-IGC^%3q2YSu~jobNiK~tvv>!Mz^8>3z`?N{H)l=}(&I@A5>
zdef!&>8O7DOw^!hpZYAfKFKHfq~E9ChXZh6IXG}<RGV$u7}R0fFV{sh-_hLntWW;u
zQv7PvuszrkG-}$X{nCDyqy17>)D?AwKZFBU=0Imh(6HSdHD=mxK4-GJ-^i~Qy=a>s
z3VJzRH+s>i&nwp_`5`~_<MiWj01hk%2i)&;B>&T${7*-kS7LQNLmhtmT~Q-x{P&n?
zpY~3BUyk-pT~Sxm6@C&9I0t-u^0D+B)bTV<`569Zh0K?aBR_B7j_OO}%FlCMD3s^L
z`Q2!z<bnR0{u&Ozf#u*p(Rh^=>i%p0Yy7@gSGq5DI^7rRqTdhA9siT-0RB@C%fWx@
zin^*pU2WMMG;Z6s1x=auYi~eSIRWlJ+w6N$9ckQjTkiM#<1xPqd2mjX2lBw*1@I6K
zEQ|w(4+TwF@A;r<)4t#TT9fX3U5cJxw#%IXpXaeN%_DF+&HLDi{|V*I^`Be^{GCoc
z)Zy=R>YO@X$vXFWd=BSvKiXgJexKue2K;*5Is5CILA~kv+_~KE`MdaghdekJ$pd-d
z?*e!T2g>1q?z?n%MdkZ2o?o`hoq>1X2^zNhqQ*@7{2j#KK`Z%p5cNjAQE!YhuNwz`
z88vM0zZW!W+F$N|zi!9(c_g1t>yW*a=9PRt--n!eU$Sm>$oHlHrT?Y>g#&P)3LMb-
zpdIZ|dA-oO`Id!#-=gtUp2wlQ{;T$B@3i+Sw0G)|I;0NqpK!oAaOhyrxE(zbG-cZN
zb<|%<&yT&5uV*aOCp7x|i`I43y74pDk*BWcU+7=xU*G^7fCG)o0j+~Nb6wP;d8CVc
zLZkQ}_j~)(I`@Z?-`j8h95rOWiW<rL@D0Afx5nig?U(jT`=$Ms#Q~q+Z`2MR2)dBw
z_Z!XYZQKymllM2Q4|>+Fje5bf@9%ruZi@O#@*|Jubw>U7w9MGwXV*s^Ouu(unrFGy
zu8BI8_t&ir>a(9k4W#+7`_k{3Nb@pJ=l!zs!*Sp(yoI-L01m(bH~<IW033h=Z~zX#
z0XP5$-~b$e18@KizyUY_2jBo4fCF#<4!{9800-az9DoCG01m(bH~<IW033h=Z~zX#
z0XP5$-~b$e18@KizyUY_2jBo4fCC5i2aVb9q9#oHjWRC{pR>6$pmj)Qu1iAyMgK+r
zMgIi{;J{KiAisnEX>K3UsNW9%1K*ayH~L5Fhx)-M7jXdp)2Pm%9-kZkL;pqpMgLXQ
zfAKneSPCETKlShtjruL|Kk#iSe4~G)eyAV(a}fvdKaJ`P>hZbpKlESpU-Vx^{THvp
zho$fV|5Fbi(Wu`N{{!Ea!Z-Ry>WBKlKNoQT|I?_>pdOza|3m*p|3&{*)PM0hd{_z}
z@IUqN5smsS@jviwDSV@Uq<*L${Bscp@IQ^}4C?W@@jvum^k4K}Mg14A!-u8t0sm7E
zAJM4a68{6=mclpsN9u?A!9N#q0RPjd&Y&Kj8~;Q9MgK+rRn&j+I(&c+)#1bDhk`EH
z@1rJ7`}J@u<$p{34}61f)!`d{gZ4;!q&*gKKz;}RQ;rwH|G-;#3vY{f%j@s~K2(Pf
z_@8ookB)->fp74wI()-#&>m@zw8tV2;D6*?%8B?NcnfdgZ4qyI9X`N^>hJ;oQ;zS^
zQSd+T4Zc-}Z}<(`Bkhs)Si}MRkDN<65&r{k;Vrx^;w`Vk2l!AOKHz`K@jW^U{s+Fn
zx9ac>zd?JXJ<=YFIDr3=b15g{f8Z^=g||h#<#qS~AF9I#{7*T)M@PZ`z&H3-9lqf=
zXpgi<+G7z1@IP`c<wX1syoI;$wura94j<q{b@+h)DaZHdDEJ@v2H&c~H~a?ek@iS?
zEaCwEN6w|3i2s4N@D|<{@s`)&1AM3sAMiir_#Pbv{{!FPTXpz`-=ICx9%+w79Kip`
zxs(&}Kkydb!rLO=@;ZEg57prV{-+$@qod$|;2V6a4&U$_v`5+_?Xidh_#Zi!aw7f*
z-ojgWTf|#lhY#?fI()$Yl;eAJ6#Nf-gKyR08-9cKNPDC`7I6UoBj-|1#Q(rscnfcf
zc+2bX0X|fR4+jqfU9fFigC<S;<+y?R9R>db-{4zy_=ex0J<=X&k3}4~K5E=P{wQeD
zv|m5xQqDKV|G-;#3vY{f%j@s~K2(PfTQ&!c+xBfiQ>OiL+(7+~^0}yC`*qZ)X&>&v
zJ-AmL?$I7;kF-bH<L{y-Y{&MXY14lFoJ%?1c<-K|L3{U|pkdQKyoI;$7T&@EH~<H#
z&w<ZJJ!;oQJ!aZ3$0aoHQNA42Z?{Gbn)cxy+=F}7=N|2pc1k;?o!$_2$Zn20YT9pJ
z=Tgo$^XKVK+mE8on)cx;T!pJ}6%N1wI8c2KJoI4DOSW}O(5t5Xa$G|59>uSdy_kH@
z%jvq=iw5`L9^9)w_h_fIQ`#x*wCz+>8o%3N+HYRxvb^8y`%$Cz!9Rj7nD*f<yoI;$
z77oAxI8Yr9{4lD=eiC)ww7<MAp?QzqdLqgmSQm8Cv=8^-9^9)A_h^r_N7^Ip@hedS
z_T8u<(|+?hl-2oWU7b;OPt<AC{x_p0lK-4G?bnUZ?z<=Kk5SXfhfmnsZv_q7*P=%9
zerHFN-5u3!+W&S`t6dXy%Cuj%eDZzi-|64!-{HWOIpF=rSQ_^{k@}M{TOF5Bzu)?)
zsAINnZP0Pke((7x+qE;O&$QpNDQF};uXn+;f97=11v_#$XwtN=c@AdId+^(+ar;Bm
zlxhETchFcG4>n=iFPiW1?|%#G&%Q+aKZrVGKaM(Q+V}a#U$fhzPNezBUrYC2AIZ-j
z`g}@{EVvH!9HWope>^|e2nVjr0rEqB@Jsk5IN%)6^I%s;9na5))z6{ujhi6@@^L=S
z|B|nxz0>d0@6+$Y0XXoVIiULr{&#dwzR%DMeIehWei`^|)R5g3HDcO_dvNbRbB}tU
z9;gTWu=)x8-K_ok^_+a;X2^iwZ|k*BNA;)sZ@tV{5W2gb3GzT5E;jQwa~<FVX&2PP
zY#->?gI-Pd!`?{ahF)FpIE`lLQ~3^!lL7bdooPH%ce)PLnR)Nu1pbo;@At?9dEoB?
zcsP>-8#e?E*?aE>jhOcReQe1<H8bGr{l1dc|9w5J|N9F5hwB~}3l`3Mt0fQQ0pGhA
z-|OQk+VkJl+OKAB$fv4X27Da<&h%XTZqq*3`JH#<fjr;?XW#FjUC=J((=M)$8n>IG
zrcC=j-r{NdOw_Y!{KeDUFRPmYsp$9pP1J?7j@zVZAMU}u`EZXqp-!k1d}9#@-hVG>
z)D9gC8aM6xJS_+8#;C(-o|XgnqKf*2X8l-C?hP7B^A3%g_KWJ7*Xi%*@96K~0328f
z2lQOo%=uWZjd~%k!>ave?SJ`RE0Y1g{<uE*oR)Mwa=qOdHE5eQ1`X$Z_y*sW!Z+%N
z`k{XC$qVIx)_?Hz$S$o1;n#`I+AUGN>AKO`bU*Fk)Svx1@2^a+%5Pb%3=~~IbU)o=
zKHv5EG>>&po-g~O4}&K1{zB`M_tUS@uhFl;0XR?v4)}Un7wlV6lWBde3+efuGwFGs
zb9vv#B^^xTl#ZlvN(Z^`wn9cm_dB)DhxWf8)s@!mIBnY3x^(-V3L44l)akn6v7<rL
zdA$y}1edD7CF+nmqz>`PRpY?nLqQX1J*{ce{=4r44cmQDW2SvwZ~EI`gZlFIrw!|a
z+VZ-a+JET5pqEmA{A#|B_l2naX}+XGrv0Mzs@xx}N$Y`Kin@<8GwOLJU-$FU^G!Y;
z{fMny6Lc*5KL7nM<>$b)ujlqOUyk<AosCMz={4>DENURngP?uS!&qL&RQpw{2R;w|
z6#W$a6dZs9Z~zX#0XP5$-~b$e18@KizyUY_2jBo4fCF#<4!{9800-az9DoCG01m(b
zH~<IW033h=Z~zX#0XP5$-~b$e18@KizyUY_2jBo4fCF#<4!{9800-az9DoCG01m(b
zH~<IW033h=Z~zX#0XP5$-~b$e18@KizyUY_2jBo4fCF#<4!{9800-az9DoCG01m(b
zH~<IW033h=Z~zX#0XP5$-~b$e18@KizyUY_2jBo4fCF#<4!{9800-az9DoCG01m(b
zH~<IW033h=Z~zX#0XP5$-~b$e18@KizyUY_2jBo4fCF#<4!{9800-az9DoCG01m(b
zH~<IW033h=Z~zX#0XP5$;D9)wcvaM!rhP?TcV3VO^3c;Wk=qCLKs`_o^cQdd4!{98
z00-az9DoCG01m(bH~<IW033h=Z~zX#0XP5$-~b$e18@KizyUY_2jBo4fCF#<4!{98
z00-az9DoCG01m(bH~<IW033h=Z~zX#0XP5$-~b$e18@KizyUY_2jBo4fCF#<4!{98
z00-az9DoCG01m(bH~<IW033h=Z~zX#0XP5$-~b$e18@KizyUY_2jBo4fCF#<4!{98
z00-az9DoCG01m(bH~<IW033h=Z~zX#0XP5$-~b$e18@KizyUY_2jBo4fCF#<4!{98
z00-az9DoCG01m(bH~<IW033h=Z~zX#0XP5$-~b$e18@KizyUY_2jBo4fCF#<4!{98
z00-az9DoCG01m(bH~<IW033h=Z~zX#0XP5$-~b$e18@KizyUY_2jBo4fCF#<4!{98
z00-az9DoCG01m(bH~<IW033h=Z~zXh1_%D}e$bdH9zGN_VPB0J&hf#Ppi$esE2uxm
zy?cTNZTq&Mb2<L?&7fXu*%Z{7<8@J|ZS9(%6FJ@ybuzv0g&h6;FWWUy&*b>osORja
zs3ST0agUk*d~c4|Mm=SIyvK69A?k7aOw{^W@l#Q2Y-TJvPSNk(us-NnQ{1#M=z05g
zRBMj@Jx^J8S5%H$HV2)tpGNiP`17d2^uEy?{r%&1HZ-fAi5%e!oS`1@8*l&)tOf_R
zZV8$&#bZZ<rmd?pXe7tqL|w4I{WYjB$F@^J1NQx>t{ktAI%D4Ny^^E*lGjcByrTE#
zPnfp@w#}M>Pe(m%ith8CwHu?(=6Fj~uN~MQG?wG<qAvL-`6>7bKUagFsy{{QkUFFe
zm*3C3zZf^gk3R~Uv>!*E%dz!DP_OOS9`r_z>(&Of+UKJl%~978nzjGsdo@-DW?l#L
zdTq(q?{q!S>-MzW6LlfSZ$(Y!>rM1i-anBC@^IBv%im9_&J^*7_`@>%;k)ky4V&V9
zQDa;`X{_#@Pg`;Zu8!JiihloYxBdNZLH#+t{Z`N*z9{COe$dO1Jdg+cNg4h`epj*V
zx{~iv#Q%7P*^+^UWWevg_N4o<J$7@{d3;e!B7VSYi9C=8d{LErk$jIL{%0Y*9`7Xs
zSI)o>qRyD2?knr>8AU$|{Udpxf2>mfsP?HyeKuEp%J(SZf3Do0@l!HTMh3nd)o+To
zMsXj_bBhlo56#uj%0DS`eyeeQ<$DzIKV|f89D@v8IRkn=KvB;x@|;@88~(?uhdhvn
zYWM{CHpON+FY-N#_#Z#BY{|fCWMKURL8nZyqdkgoMDb+c2fQZ919`v~m2sa`zDE)N
zvl_h~zmp7j27Laid(*s>m*T#sf<{d78&MbVMKNvo0k0nNKpyZ#3-LwrJ&O1rZ$WIy
zK$SAk*%4)mcSm*Of8t4^-}hQ2599&=ga4^g&&6jU153|9(fp!^4hD^z;?W~PQ~08o
za{Pc-9(f=S_@YnpMY`^z$n~G4^;Ne=^_rru&(LpLPfSty*8L1c-Ph2%Wr`Yirg^#*
zHDC8Pqb5wzqNeS`4}!*X{Ev@=rcCEqQUA}!n+%!aZBe7P`JtfGIsPc>Y?|+TUylB{
z4l;k{(#mMw|6A*%DL(aN&@lcdW{`g03yD0C2mBBIXWre)qTlGx^IS;lkY2ElJ`8H!
zb*BDn)VL}BOVm_apS3;5yP`VrMT@SiPyUZZ>&eP5EBgPYA0!X-gP+t7s+}lOC-tZk
zU58NA^?*;NEw9d%fx7u0@j_AD&^ow^S`XLfMSdgAkNienAAYX9^7Bv5fM0KWC9U80
zirp0TI@f<<y6}5mHRORj)WgrHA5xs_yvz3};(t72i{0XXLgtEj`~G`Dqo$~FD@D&e
zE%rRr<^Or!o=DG&Jz;*obr=387C><wijRt3FT9^T;1lX}eN4Vb5&yHao)iBQ^0)AQ
zMAt18TTeuBzig4^$oln*t3&xAMNj|2>skL_k>lVi>f|d__li}mfB7Co{Ewfh#cuIG
zA#)4wFZ?=2Z@P}rYbRTS`tU!CEW{`G1`U~_KR-6gb)lHgh1aJzuIQ(Q<B$h@T_gKC
z`5s05&mudV^xyD5ArTAj7k(KvY>K)+wzRWR-~TQDsJLf$&;b4?mc+vAQ@%&hD~sch
z2mBBIr@sBn%x{kW2^m><ztD3gXxtRfp9`Ae`p?V)o9F&5QN5=4x4#DU;eTTO7G9t7
zJ&JxBI1YKh|KNY-3FR8R4gV8Pa(Vs2)loZ5@tUaJ_@X)5(zdNZXHD_ms9yX}%v*W&
zDL<s>r-kE_2YgW@_#)l+Qsn;E939IXKgR!r1eM4C4eNuRHN|VAUcmp%QHL|<tCSB?
zKdtD8DervgxQZN~ey|bxLDj1wbz6_RmG4o+|IE>`%<*IVPe@RC{i*lcuchl4ucdjl
zPvC#%s6!v`+?D3#?#lCai|>klnDWk-j;qM=>(P(sd@FL^Iq!|&d*yo+@jr8PEOY!A
z{}U2aUVrNAh@8m2$LF^>iT|0S4t@QS_Oy;kd-g%<pA`Kt<()4bSCQk>4>m$SsCreT
zZtGFE@;!?9pE){~Iev`)2?;8%Kb7B6{7h6E{%4Ll+_)jA-4wezqZt1a^HyGcs=rk9
z)57t|Lp}V1>Rz$Y&bxe%BL2tE;9|G<pOCro_`m*vpi`!(>l%waKXv(k@;!?9ANnJ%
z!<5G*_4kU-6^`F%d64f>#Q!X<hs6Jc{FTQ6`5s05&;06BzDE)Nga5(*l*c9c9z}1J
z9KX@>Am5{i|CwL^QuHqTPsmJp9FXr(#Qzjk-sP{$_bB3j@IUyU^0*}5qv)-V<2PC!
z<a-qHKbQ9-3*E#2gp8HP0r?(9{LlO<QoctK|AYU*|CGli`5r}Yl^nm(@*v-%i2s>i
z|5EfW{7=YCc^r`MQN;fgRo>;V%l9bafABx}pYpgQ-=pZQkK;F59^`uz@jsXMBMaTb
z|AdT{#{u~sMf}hFDpI~j5&wh#!T*%UCHWpjZ<QRs(efbQqlo{RU;k3{F8oi(OnDrT
z?@`466jk2kugmu+;(zcz_@DB)B;TXxt&ih3S{~$k6!AZo_ah74!~cYgmB#`39!31m
z{3=quM-l&n|H1#1$0hk5MQ@cHztQp_-=m2CnP2}>^e+5Q$V_=0knd5%{}ff;<*&>4
zDB^$cKlq>WxFp}B=&g_AH(DO#dld0Mm-izJ-NXNcjFra$`5s05&-^M<zDE)Nga5(*
zl*c9c9z}1J9KX@>Am5{i|CwL^QuHqTPsmJp9FXr(#Qzjk-sP{$_bB3j@IUyU^0*}5
zqv)-V<2PC!<a-qHKbQ9-3*E#2gp8HP0r?(9{LlO<QoctK|AYU*|CGli`5r}Yl^nm(
z@*v-%i2s>i|5EfW{7=YCc^r`MQN;fgRo>;V%l9bafABx}pYpgQ-=pZQkK;F59^`uz
z@jsXMBMaTb|AdT{#{u~sMf}hFDpI~j5&wh#!T*%UCHWpjZ<QRs(efbQqlo{RU;k3{
zF8oi(OnDrT?@`466jk2kugmu+;(zcz_@DB)B;TXxt&ih3S{~$k6!AZo_ah74!~cYg
zmB#`39!31m{3=quM-l&n|H1#1$0hk5MQ@cHztQp_-=m2CnP2}>^e+5Q$V_=0knd5%
z{}ff;<*&>4DB^$cKlq>WxFp}B=&g_AH(DO#dld0Mm-izJ-NXNcjFra$`5s05&-^M<
zzDE)Nga5(*l*c9c9z}1J9KX@>Am5{i|CwL^QuHqTPsmJp9FXr(#Qzjk-sP{$_bB3j
z@IUyU^0*}5qv)-V<2PC!<a-qHKbQ9-3*E#2gp8HP0r?(9{LlO<QoctK|AYU*|CGli
z`5r}Yl^nm(@*v-%i2s>i|5EfW{7=YCc^r`MQN;fgRo>;V%l9bafABx}pYpgQ-=pZQ
zkK;F59^`uz@jsXMBMaTb|AdT{#{u~sMf}hFDpI~j5&wh#!T*%UCHWpjZ<QRs(efbQ
zqlo{RU;k3{F8oi(OnDrT?@`466jk2kugmu+;(zcz_@DB)B;TXxt&ih3S{~$k6!AZo
z_ah74!~cYgmB#`39!31m{3=quM-l&n|H1#1$0hk5MQ@cHztQp_-=m2CnP2}>^e+5Q
z$V_=0knd5%{}ff;<*&>4DB^$cKlq>WxFp}B=&g_AH(DO#dld0Mm-izJ-NXNcjFra$
z`5s05&-^M<zDE)Nga5(*l*c9c9z}1J9KX@>Am5{i|CwL^QuHqTPsmJp9FXr(#Qzjk
z-sP{$_bB3j@IUyU^0*}5qv)-V<C6zC183j>9DoCG01m(bH~<IW033h=Z~zX#0XP5$
z-~b$e18@KizyUY_2jBo4fCF#<4!{9800-az9DoCG01m(bH~<IW033h=Z~zX#0XP5$
z-~b$e18@KizyUY_2jBo4fCF#<4!{9800-az9DoCG01m(bH~<IW033h=Z~zX#0XP5$
z-~b$e18@KizyUY_2jBo4fCF#<4!{9800-az9DoCG01m(bH~<IW033h=Z~zX#0XP5$
z-~b$e18@KizyUY_2jBo4fCF#<4!{9800-az9DoCG01m(bH~<IW033h=Z~zX#0XP5$
z-~b$e18@KizyUY_2jBo4fCF#<4!{9800-az9DoCG01m(bH~<IW033h=Z~zW83kMGD
z4;nLlA4OgdnRwI419=E_K)didAP?k${(^o14!{9800-az9DoCG01m(bH~<IW033h=
zZ~zX#0XP5$-~b$e18@KizyUY_2jBo4fCF#<4!{9800-az9DoCG01m(bH~<IW033h=
zZ~zX#0XP5$-~b$e18@KizyUY_2jBo4fCF#<4!{9800-az9DoCG01m(bH~<IW033h=
zZ~zX#0XP5$-~b$e18@KizyUY_2jBo4fCF#<4!{9800-az9DoCG01m(bH~<IW033h=
zZ~zX#0XP5$-~b$e18@KizyUY_2jBo4fCF#<4!{9800-az9DoCG01m(bH~<IW033h=
zZ~zX#0XP5$-~b$e18@KizyUY_2jBo4fCF#<4!{9800-az9DoCG01m(bH~<IW033h=
zZ~zX#0XP5$-~b$e18@KizyUY_2jBo4fCF#<4!{9800-az9DoCG01m(bH~<IW033h=
zZ~zX}n*-mA8n&J@LE}00o{zG&Q$YhcemAP!Ha`?}I>#-WgI=+3N44g7ebhm_I_h|i
zH$**dpN`s<<L9GV%zwTm$4^DAvCl-U&+)pb$IO5JM2^=+J!LmW9nP_3Q_u^xZEMhL
zIesVVl=<)P%5m$KpdLHf8q}BL-aSEsc3^+dSdRBajoGoILDM<veQIxtyzcFfJdlU0
zu6omG7hVVCfjr>Pi~ad8qlWFn4}!*XyeI0k-5%AO<IPb=?V6}(a`bk_wyX@?6!n*s
z-`zQWHmcopeiUzw8niD&jpc|RaQ=}8@_;WI%f8S#RopMoUT80k*j|1SHDm`51YO9{
zuaBIxO&fzca`fvRe@fRW?$6QN2ir!^fM17t(cDkHp5v*LLH%}D)JTqS5AHQ0_o!Ry
zmbz^|fBRrd(5P+S7BrQk-<Rn$-B(e>&osJYS>4aSV|&mWw*G;jQ#p2a1P$8{qb72s
z-=W{3-)X*nhjvZ7UdeV{G!8`LUe>M&I%bNi+e`A>>y`l@FZE~BI4Q;d64htA?x+Yq
z;pa;7lX|1xsJH6-@SRTt^_u*SqQ>j`^_v}O{Oyi(zqW3@E#J4%GvL=PyVCW_b2-va
z(NED&RllF2ozhO1(oTO8bw1rc7|2ogr)I`R@9+5hWMDoS=xC2Jzpi;EM|c~u$9Pn)
zA@V>TygiVIYUiQ#L{P82{Z`Omj_z;nP1kWQML(Zx$p9Jn&l&Lf1Ul_UQD<|+7lkb0
z2fS9u19|ZFKps{$4?m9@G|kJRsQG35`KpK0{4)=8-^vS=Eg4u+25yP!H9g<2h%X9h
z#SeIKkq7eN?SVWjEf2qr8nr(}P37qGs_Zw-f1~L6W=jUhz;a}uXk47%M;l7_&&G4a
z7lpHcAMjcv59Gny19_NV9@^WY>~wcjj-EfZWPl7*Ap@F!+vnweEzQgQTE0HZyxBgV
zxSuQLLnaT@1AiA>;=qRWL4)?tgF$0CGH#~|y&Io_3@k1K-hO-St-l8i<j8fQkWa3M
zcvX=H^5E^Eu6fWpQ@);5Yg$jLm2o>>5X-(T8aJi$q4mKOeLb)~+qfZUC`bSKP+C9i
zLXLjlWYj+XC}=Xr{riGO?dXx9sT@BSHEh3)8qe`7Q3KZ188nik&j;9%<^k-<(buQg
zo94+mkmIr^lwS}t;PVn4NY5P{%TeosSk$x!o6hU4dHOE>eS}{M8Nny0E~p3UA+(1^
zs|WoZ*xeO0YKs0n*e-Vld|v8T(|Vt+Ir{w7BX;;u&_s@n&Y8QT#_W$#(>Z$o+nvUb
zozGG8Ah=&+{gLIWN1xaFyzSW?G>~JX*BhTdR17}N+YNaj58fV>hoW__{O4@T%7CvU
z)tSa&oXrtFhYS~U>5iy=)A|yM_@c5Z(~slV<(^B|<6h2@`n+;|k{|NZh}S*pR=2u8
zv2AP6q$&Em{p-`Zo-H|+OV|1w&F4EhVHzi{sBu)i=YytAv2Jx(*7x;wwt8&c+MvE1
z>r-hf^G&^e2h(#7V>y;}9yv}pkH!9<{E#2~)zaq!FS@?HG8ApT-{RwQ9#8Wb?aJ}C
zsL`~J_*9N?Y-Kp+e*RhW`)M!b*nAbZGT*4^z9jX$GWATJ$us_E>Hf#(<9Wlr6Ll)b
zm7!?!{T5w+*tIjL&lKU+>henGLF2U*o38>_<{MoT^^ECxdPVAdb?cn`lYe~CY+s~#
zmuJrJq4lYK9O<L!KG&l>&oxIkQj|B16F+}0Xv!4fQZ>2cb-yoNk2{j%9H*Zj7i3`H
zQ$Zv4`>4qrssC!$Kj($>QU|~E&rw5mPt=7R3oypZb7Vl{!@n2RVT$miv3TOwfjaC5
zQD<_T<Mi|6f(-mLs^8vyCulfF+C^iv3(gnktC@Te<9rH$y~vk~=GB24jm3@KyMp@d
zYf&RPF0w3Y|KI%YtJ8LGRBw*7i^ggfoG;E-Gx;CC&fk@;@1M)D_P}0b0yb?78a5yA
zJdq>ZXe@60Dr&^O6*ZaTBFm!o|INn*ol5gvcIHUCXsmX@`Qm&vlm98YUsHQv&yfI)
zi@GW5FUdFlh52t|-kSIMT%Ski3G;DYdvctkTx$PgpGW#-yFKbej`XLYG@7@)aDF+z
zRrX5{wgkOs+qVV1o@4ESJx2n3-H7MX`VcSYNPBPY_U_mBo=?~HUd?fia;g1~{qKY$
z_Q1NJr*oveG<SR9{BnM)?0?oj5cG`62PxJb*mER6>z(;{pH@?(y*GDz_xqZ!ruh?F
zbDX1GYX9S+b=+t#&D~x&zntGH`yZd@=cvgCDb^m?b0nZ>-e%f+bGP@R^FBwp)c(hQ
z|Ld`|e(DoB(q5Xoy>Naxzg6}>*GC;p>$M-xvG%~8BLRLL|D^f4qn$a@-kZC<7oGPx
z%BA){_Hk1?(s-#q=SX{L?)Jj@<@{FJ|6Cn)+_avdV(o!FM*@nj1JK@^yS;lqej@eX
zZ8^?SF17!$o)hrrz;>G=?WMWf3+I>fTV?;VVSUiE=Id>~lw<9IJx2nzZV7tTPPRtn
zNPBPY_P+U{pqFj!nxGRo&QUJ4|FK^ueJYI?Ka?ZwrMcS+=a=(aW&Z=u7uPlU{J_tf
z#z8C6-kZC<vtH`r+RM7NLB~zcO)An}n!CMlemTEY_CK1J%-1RGOzReQvd-5Wm@it-
zo%Y@w?cJ}3j#=;dpy?cazO0AM{qG|=&QUJ4|M87cubJi>Q>48#M|<HsavrPfbBgNo
z3sK$X&w-uKvG%}TWCApAu=9S>6pQlD>(%6j)`9<R)VL`wvMg%<-~2v%Te`pAog?*M
z&HCrOa9$eG&n&cVcRmr+YkT(u4dz&TU@tNOcSQADXGhR*jtf06ynl81L%UdH11R`!
zK3=XTjg#xkkvd=9Iw$|+zYad9sLngugNDqXgCEba0Asv7M+Uwd)t=^qKAmGx{&~Hc
z{LwlMcSg0D;vA=+9~Wd`&+ecB`*qZ4j?{lO>!0(&d8vahDjMfdw7yjV-gr4>K+n1O
zeXqmm{@3Ajy?^$25;(KEoKb)L<ESZ9^!jf}^Yrb=kwU+620U-~r+lr+(dU&PP3!AS
z<w*UnZvB&g@?Qs^v(P&Cd4XEfyg=<aUYUOUv?v4q{7awddYj@x&kOG_msj^holfg<
z_U2f09(n!B8PGf$nkPw-x-GYE$tU^5w~VFdwBF3ItaforRIl9{HJIa->BmotGN5%u
zw7!~RS?7!6T$zXao(dYV-$za6sOPx!JfC9GdF1sgXW)lXJ+^CSP+yMJ=auV|{E#2~
z%<}n}g|@Hv-VGYDb7zAlbL76jmHEHWpZ1)IO6#Zg<+#xE#rvI~_eG7_v7<rLIr{v=
zuiL4UQ8_O3eDQwIfM3`7Q@X!(DQ?;rG@Rx)oXC;-^!g+Z<N=>k9iOw%_Ne~juBc8^
zq+eKE2Goye+=L>}kzJXSMdQ}|dQOM^AnHtxi#t2~uc8cW+ZyzmX+1Q>h4O^=U%Af7
z5Bb3_RmU&+b1;MU^{BBNi#W;ai_1WfFDja^>F%g8`(xB}j*ZSw%_H!;s0mZ_=M`Tu
zt;43cxNP&k7Mg)CNA;Vo_b4`cz47@&CA?|mfjr>5=I^_T{vKQ?S9$;9GT`gh?zQWp
z4&><mp~EyTL-8w71J>0UG?L?2qlVME9^*N7cLj~w`Ex;2Io=jEYJPudD#u?$4cWm1
zK^JoTc2r+l2Vp2jjknkQ0gBvjU0gks^}n<pve(tisjgn8e~bynCwQHa2l5ayQ2sdN
zM)Be8Z9#+fkN1Pda@2K<nbGToEg2vKRnCCl*XlQ)*I^{bM&S~lGgKG;#M>BoAP?Ri
zRw@tQi|R0q^HucU&z203fhuG`*Wdhl&!zb;@DJha&~JK;k_YnO?SVYZFAqiYp80de
zhfM22D|%kpk^wTX92wB_{l5O_uCyNME@R!gkZt;RuQ~ES9=tt}hvmt`+iwL8rtt(9
za%6pNuQ;}3U?~~U{boH+sJMGqP``aGY9vRl^Mq{E?|aRW2lC+Ufjq2i9{xS5&z{^H
zG?b%wJ~R6HVM_+cz<<tw#<?>O;r@%7|8pGpz)(f_A8$+KfjoG7SnWKlT@%!68`cL6
z<~Y+2Y2QCDTQbnx893b?mF6Mr&5?dQoOk*!uYU4C9=tt}hvv#d(eu>YNApT&+gKTp
zf1eq*Yz{hO|0Sv~M|_U+j69Hs;CH>X<M+sew}-ms!JlIpvTsCP$WiOEyFYq5jsJW)
zwHwbCTQacR8PGUnf9`6hDKf8M$UOa**B^Nx58fWg19>RpKK}lFL8Er`NYGS{KEFw)
zZQT;olcW9)*Z5IIuSd3I;K~{B`?e?2{o1x1echN5`{=`<i5$6(5%NjD?^Q(}$b+{B
z@<1NS@IMRfU#>ilaci2txIIVrOZTR67?+~gFWW}XfcxJU(>&L&=U6nK_Co88_dDmv
z19{+n7UTBeKx1*h*WVnr-$sq+=<^_T+4gNg=W_J^bYJR6kL1|sKA~!#zvzCG&s*MR
z8b_>nYt*3mdTV1j(jFVDJ#xM{U-&3|6dY(Q4(R#+5B?E!!4&16u8A5jMSm`?&3s+-
z?i{sl)9q0wOi}BE`Sq$}>3Y>MEAk1SiCS;&*IH)8`_unllcT@y(d6g;lB4^>C(YL{
zJ(#24e|#a`cYHNR|9+>^`edCs>iKS;H?uFzquFOaj~cZ5qQ+*$9ovJZO_6rdSnY!I
z#reX|;Ah|f9DoCG01m(bH~<IW033h=Z~zX#0XP5$-~b$e18@KizyUY_2jBo4fCF#<
z4!{9800-az9DoCG01m(bH~<IW033h=Z~zX#0XP5$-~b$e18@KizyUY_2jBo4fCF#<
z4!{9800-az9DoCG01m(bH~<IW033h=Z~zX#0XP5$-~b$e18@KizyUY_2jBo4fCF#<
z4!{9800-az9DoCG01m(bH~<IW033h=Z~zX#0XP5$-~b$e18@KizyUY_2jBo4fCF#<
z4!{9800-az9DoCG01m(bH~<IW033h=Z~zX#0XP5$-~b$e18@KizyUY_2jBo4fCF#<
z4!{9800-az9DoCG01m(bIMBEp2-O?TAMKbt(2g6o9n*ekzqDW4FC2gaZ~zX#0XP5$
z-~b$e18@KizyUY_2jBo4fCF#<4!{9800-az9DoCG01m(bH~<IW033h=Z~zX#0XP5$
z-~b$e18@KizyUY_2jBo4fCF#<4!{9800-az9DoCG01m(bH~<IW033h=Z~zX#0XP5$
z-~b$e18@KizyUY_2jBo4fCF#<4!{9800-az9DoCG01m(bH~<IW033h=Z~zX#0XP5$
z-~b$e18@KizyUY_2jBo4fCF#<4!{9800-az9DoCG01m(bH~<IW033h=Z~zX#0XP5$
z-~b$e18@KizyUY_2jBo4fCF#<4!{9800-az9DoCG01m(bH~<IW033h=Z~zX#0XP5$
z-~b$e18@KizyUY_2jBo4fCF#<4!{9800-az9DoCG01m(bH~<IW033h=Z~zX#0XP5$
z-~b$e18@KizyUY_2jBo4Xig5CKNmD*-CaSWrv2-pF4)^|1r3_^zZx}cCtHL1O#8cc
z1@+raQRhtiKaM(Q-;e4t?YC?S>a=U3PMP-Ct_eC}n;!~#*|fiDW6<;V>8Ph|=KhUQ
zhx5;`k2+|di+a+ue|6MOyCLdv)BdNTcG%~mT1@*M|CIJGML+H%>GwUHe$ONJ*{Hqw
zxY}R;K+rRGbJS7O{<^h6t?Bc;YTCa$s@v|2YBTLCuQT&}x;tph+S`H#P5Xb0noj?J
z%(Ne>IdGkJO&;jq>EGc199T^bxG$KpUq_85pD<-#iyFy(Li@Kz_1ZmAr%n6nFMa<+
z?mxYsKR&|;9G}&%tG^G`fEqso?vwtKeBGnTpZ+EL&6)nwedVxyIjTST%waorG-%rH
ziyF)Oo*{e`c?cPJ)5rsP*uNh?Sr<R~@kc?Esec(Y?cW=9$)9N7uRFY&uU}~2uRGkI
z{fPFRCv3^U>SjRK-DX~AbN@78+qVUsOV{%T>?cvv_VuW-yzk!vzeFBF2JlO&1I4=Z
zn|$ApfyS?Y{XP8isKNaAv7cYIWT2iI*tj9+xpY78xM}~ks8Rb_)PQN<vxF}q4<Q4M
z-!E_;IgfR?4u=2nmd=(8%qIhWfA!vUo%2%E^`i${g1YkcB<{O~tl|f}w#Wl{*nd&n
zU#rS>Li~@nNVa5PIWq9Tx}cM`ZEMhLrhSdmpE+*dKLgj1$V12geo1wpSe1U0&l57x
z`1OhZ@k(Y(239fyKAykTwrmc1B|UG~nx9|97lpHd5A>QO59Gny19>Pn5BMLiD7Ivv
z3K`h2KBz6vo1lIFZ1jWVA!MNO|L*2IavrO29}xfJEs-r5sB#AW`esnC{UGX$Y2QCH
zzKA@84B(4Y2a4P;^SU4ppPUE$kGB-IWT46!(0me_-%9(MCvxVzk^WuqMdTr50AHj!
zP{a>-UDQ1f_#bZ(Y{@_)X29n)e>u(D{al*o{AJs*J?PSWkJ|U|h%X`!Ap`g#)q!H&
z{XpIFfdBC`&Xx=`QwBa0HE6d)^_uqmJK}%HL&yOBM|Gfx|M9x0OCIn)UioEh|Mu6Q
zK0DhJW!f)Vm!T;ChYkjf+xMb|P5a-9nzUa<jhOb29toPV{riGOP5W9e%h%t!w0_o+
z!$Fhw-n&5~ru{EOjoGbHgQop&NA+3HnJCl#wNWpm=enOZ?U$8Lj<dK7_`KGa)@#$g
zt^=v=wC}$KUql{42Jl6y1I4=dfx3NPagh0+{md_Q`_YF%6LweBh-trWb>;IV4`rXE
zeXZx}>rw2=>$@&B@BIJyWuT)yXeiBV-DldbTYd6<Lk93UsslxQj@Lz9@__&G$}ejx
zzf*VrQ^b=WMom~}N6@fozi54~vd#;~@eF9aSf8)@^|apEDc1dr>A??pm5>MWP>*@4
z@jqVaWo_|4A(IQ|r|$2~yx&)}{(D*Ht8T~9b7M1~Bl}U*q_u1c8ZqsAez-429zq5d
z-mb|Td8@)aKSgzj|MAR}wZ;FG+upwuHDEVKoj2{5b-p-`XTaxK98B{r_L}znr}!W8
z5Hf)OQ5`7O<vvlfeqa2LpV_ju_@9u;@;F&Eep};t-4`8A_i>&{_i@TPoAo)4U-#+E
z_rtXB8NmOLhme8t+BNwkpH;Xng#YoDSk@N*voh^k&!_46Zta(K9_x0Tb!&r;+YM2N
z()j4(j0=i|K)>(pfjp3h<{B4-|0%bp!~cYQE<Y##5H)4Ljv6)X`?!`J=HpyiO#9`Y
z<>ot@=jq;*uS@aAQJ3<pea`^bg~&t5!1A|u@=TtW<9ZSP$6HoeTl~-Jw$nSK+U)kI
z6Q=#L&Trk0qjl)t`g_oTY2WjJ{~-?{1Na}+fg<CAye_Jq2mFs$Y*}0UkJr8DbLIQT
zJEHp2^Gjz<`(>Tqx*g~3w}J+(^+ZswX@9lOFTXcr0KcI+P^@}C%jXOks9xRTf4nlw
z+Twp!x6U=6wVo5xep%<YZpZO)L49eSrk*@+qUT}d{Vl&YWT1Ne2<L_KvK;rPmb<R-
zKc3;Tw)mgb<$%`lknhocS?9NI$GI`;Y?{x!)3on-!2ghkkOBOd>Ohg}QC=5S&jbF)
zE4HjH{>SUy^I4t#(bxZ6pT_sDN$Y{G=Xv=#3jdy{(`kMAcGG^<&M%)QWB?zcI#8^7
z|I6nL8Cbo##s7E(mbJzItZqH~d7SOLwBDFM-~G60zpV4#tjF=^=(^Iplc&;ibzM9c
zA5Zw|al5F_c|Bx+^P)OX#0Pp^%*+G+$19|)E&j*r%=0<hhpg7?Yu5xFv+JV{rg@H!
zas6kG!r!zps3ZBslcxREI*<I`kOBI4)q&#7evQ{d2B;I&f#S?M!T)$>%G%<8LME%%
z?tFfmV|HECfpi}X|1($JbN#29_3Y=ZI)7jAd8iY{KYJ^v^A`V8O%C|@pFPR<Jdwu#
z?BV*)9F6sKca*iaMVa>fztR7ahibNW&I{*dK7NerKi*Qx+Twqz+3tM2&%<ebhb`%T
z&cpbhISPO4mY^QHHmcLK@Ba<|LmsNx-Z?Lv7yOU6f|cCje^$5MX?+H*i=q9Mtbu0w
zT~|l-*$<-5nD+fU;(y3P$Uyb(cXM7iFZdsC1uMD5|EzAi<9cSf)!ue0Xux{UN168h
zx8Q%s!|Jwo@=yNf<8zAYg#N-aQ`Q#$v%2kZ@1CGR+xbLLuW7%m^INy$?0YI`#IA`N
zFztID@IT}sWPs~GssqKD&l&N0$Uyb#1pnieS=JW+v$}Pz@wuNoUbn2XT%Y6kc%KP-
z_nn~OH2!CT`F~;wRFCIgnfOxj5Hdi0sty!q`cnLlXQr$z{wHK|_1c}U&u}=6*L^Cj
z*Kn9|&T|yF&kOKadM@tIX`X<`@HsJWt5?r{+VC^vA!L9$R2?YJ^fL>69slDQD{G7Y
z37K5}cDjFG&}bfyqW!YYW8IG9^FyAn%?|~=oaTo-f&Ynxu>AGx^@h(O4<Q58o9aMu
zsXnLZIQSpWOj%p}Psrr*x8DcW1)a2QTZ3LR?U!{P>vo)*qt4qGqPk7{o(KA6@(?nx
z{Oy`Nljo(3XI^Oi;D0=0Wo_|4<+k5nMUB|0lR^Ea{j$y%$MFnkTr=aGy?&SL?@iA(
zev&+Rosoyt%>(|&E3m9B{%3LZ{NV>d<Mzv_Vbi|%$Gh#OsK2D=kajbkXORV~Jo>s8
zkLJAoBx>5e9yMm#_g_PQO&%84uE`sDTgmG=3#~W&k7ul`?RTTv(>!slrv2YUT}b;w
zru~JUuSS2r&m%mN?pO7vampiUUGaU{-+2GL51&(3VW#7Bc0}17Q729Njb2xL{*VFg
zPpJ+RxsK;`K^{Ii5BMK%DXX>pY}9kfFC91S|0wFLwQLINH0{6rR?wieo(SqS?fbaZ
zOa4mxzmJ;C^R{STf8Wpid;jB)f+mw688z)6JP>rjei1cf+V^#}2JM!p-t_#?VD^7o
zHV3_u_kA7Gd(%2xm*Q&Wziz*m`&{9S;}g8@$pd-Ve^F!kIS1vVE`ACB<7b>L8EB>q
z_;~QA)BUX@>3ZAKX&&ckd+5QSvApl!o%>wmA!Gnwq&iTniyx@l_r?GC8D~odnkfT&
zb_WgEj_pC`P5b^G@j2unWB{L|I#9&tcwN*b5BMLie70nu5i{WH+CQD=dDQ+Nqo(cd
zs4>&Ne@FZec?cQ6|ELZW>(cM@eM1J`H2jZO0$VcBh#4pvPww9ZpF<u(2Jktm14Vp}
z*G1j)fdBCp!IlhEIRke^b*6ccPMY@pGvjl}L&yL=M|GfB_x_zcgbWn73;d5)2wO5x
zg$!)k81%g9K9=^kZws2T2U~(hP5b_t@jv7tWT2Sy<cIv!`8ltm`p5rxX4sN}mCS(8
zLwF?5d!&8M`#5tRNdIj39P$t{fX`7KC~}?0>w-KiED!h}Z$WIyz;a|@-P)kzc0<%5
z(>^{YWSf56YmPjShkDP$qw~Mg{tf@*nP*D|mXd+*M4d|OhQ49i*Zh$_pY^5rB>i`B
zJ%>Dm4B&rM2Z}4*5Au6M2I^Ob_#dzQvbNujYBkM!tbL#Ncul(Baetcscn#}2l$H3B
zkK^m0?@#kt?@8;T@3)^v4cfg?y?Nh%J${BfgbdWLAK|=n-Wzdz5&p+pZdu!+^-m8U
z3Yy64U1|UJs9w7#>a=O!pKE?0&)cj0vOZ_sj-z!&u8z92&dBavLH%i-<a4HdjW1@r
zai|o0owpbAKpwn3kcZ~T1OCS=xvXuG|0%1z{4i=Fty?*4+SmNPGw1s~b2`d4ZU|~O
z?R)?JdU~Gd#q>PU>*@Z@(R{t*=BT52o$#6KhWoy*XZ<N(-}+qCBk4J&hfVuh*Hr7e
zYyY~am-G5(+Sha5r%p!Y=f3^<sQxtX-`Vtg`tACtaqDOg8p``+<)7n(3~+rzb)d-o
z8m|lTKpq;25BMJ^SXtZU_CNkP`8$+6gbd(AR0oRmzg`#Qfjq#6xz}~Y|2R3z+Twql
zJDyK`8hM}|8gX0%bxYmWp>FX%-qOq3;(y>q9k}7`1%FE(Xz#T5#^nJ1#|c)}7XJee
z;bG(Q&^f^8Cl7F}4jjY(cuOy9i~oTeb>N1#7kmPFpuN-H8<zw4A17E@Tl^0^golmG
zL+1dWpFF^^I&cjC<1M|cE&c~?)PWn`UhoOzf%Z;&Z(I)Gf1F@tZSg<w5FR!z51j*i
ze)0gv>cBDlkGJ%)w)h{oQ3q~#d%-7=2iiOBy>U5!|8atqwZ;FyLwMM@Jai85`N;zu
zs{_aIKi<;I+TwrUMjg1}?FFAe9%%2h_r~P_{>KSc))xN*58+|s^3XZJ=O+(vtPUK*
z|9DF;Ym5JZ8+G7@w-<Z@d7!=1-W!(#_#Y=&SzG)MJcNgh%R}b?pPxLyu{v-J|Klya
ztS$ZrZq$Jr-d^wt<bn21dv9D0;D4N8Wo_|4@DLt0E)Sgpe17r($Lhc_{ExTvvbOji
zxKRggczeMokO$g3?Hvxl0XP5$-~b$e18@KizyUY_2jBo4fCF%#m;<yA+6V2UxP2Vh
zA2eoaCyKo8?S?#%hpVo7Gq(@wfqI}G=r7;^9DoCG01m(bH~<IW033h=Z~zX#0XP5$
z-~b$e18@KizyUY_2jBo4fCF#<4!{9800-az9DoCG01m(bH~<IW033h=Z~zX#0XP5$
z-~b$e18@KizyUY_2jBo4fCF#<4!{9800-az9DoCG01m(bH~<IW033h=Z~zX#0XP5$
z-~b$e18@KizyUY_2jBo4fCF#<4!{9800-az9DoCG01m(bH~<IW033h=Z~zX#0XP5$
z-~b$e18@KizyUY_2jBo4fCF#<4!{9800-az9DoCG01m(bH~<IW033h=Z~zX#0XP5$
z-~b$e18@KizyUY_2jBo4fCF#<4!{9800-az9DoCG01m(bH~<IW033h=Z~zX#0XP5$
z-~b$e18@KizyUY_2jBo4fCF#<4!{9800-az9DoCG01m(bH~<IW033h=Z~zX#0XWdS
z9QbY2xG8=vYS?}rHJIbQQN8x#sB<}fKdQ^_jB3mA_NWtfL)5D|ZrT|1yscdmbS%f~
zqYm2DQ9E<=|NppsI_lvZ{q-&Osi;5Aiu$?Y=c8KE|KF9Pj(=m+VN<*|>M8q7)U!GI
z?>lDKMZKKkXQQ68_O_@TTTevUwyi;DbL{Gjvadu9<aqvEl(lRM8p-kRe+%kQpKmh9
zUqy{rpcL)2dD|)HcmG9ke(__@Ir31CJoKIqnl{CU9t;|@_udT}$?-c;1F1hgoul{5
z?a6-}%2EBhbCqqQXF$G8(f!$h<exfoye8_D?cEbJnB&eTf_m-T*`UcBRcD@`H?z-z
zKOrM;8hIcO{M~^sg#+{Bz)zy4P4VkdW9I$RK#tqD1)a0?4+K4v<IPb=Q$G}bCn5vo
zWkA=-6#e?w`gHxPEl2lRC(M2HNRAJz3+m5)eWq{sPl}Hs4_wc|Kf{4$<$&%>D4seQ
z)SrEfd<}l4yuO<wRVM>&r-BAd(S6YQeEmr6p1+$yCb*77J$N0E2l7zG`@F}F22Gpd
zeNkiSKEc@>{qN?kG#=&NS<&m8Eg5L;4E#E3)D-vZ4jRD!gv%a`W8uDGZ2tIOryzN#
zN*?|gHEoJ_M~zuWdsL3=)&{k5AF8>#alYLmGw|_8L6fHFe;176fABvSD~kTkn-zH=
z4<Q59xexeYOHh|7`t|!$_#7`(wq)S44EXbNMdRwL)2{tP^LfYv*Wd7;a6r#XDf&1Z
zjZe4?@$8TRGEk5KpO@fdnor<ldM^JpuJdr6=VHZheZ!j#c_0rV11s;lz7;iTiW-mX
z&j&x2o)dp8J=gMQuJd@IvLyra%7D+8bT;{(vvyZhC;kWjbFqHt_q`F32l5ayu=4(A
z>z1G%Q}pL`#b-Y&Y{>u_SX>7D`S`)~_rsvw64i_U!T($=DEfVGRpfy@gbXa7|MB}*
z)28_TdqJal9%0QZtoeTZe6S@0WMC;7(0rVV8`cL6;(zcz7t4x%-`f^>AP*q}%jbVO
zJA#Hy@s6l|?q_-7u_Xgb%7D*zID4HDU$0}z6uY~Ec+QJ?jY7H6?|Y*o59Gny!*b=J
zXniIBeQe9hfRA^by?)l_hk{O<;&-Fk?d`XM26OcFuzF3`3lw$Tz{iu!Uf=2uQB$U<
zd1ZEN51KYb&C8_uwG}nbwy$3{mDVqt%Jb74J`^-zikjD0$5Yh0Xz#ooG-Qgqb_VsC
zuY1;;<2_NQ&A(q;j=uiciFDs+Pmb<0o~WE}_w)XA8t43U8t-$M>pWcNxri3}cjqg4
zAP*q}%Xghe`Bl{4S$@6eq$#=|(sL;OZnpE!fYvot^!ch!W&bhLhv<AL>bz)Qk=H{e
zE^$)pLnz7*&W!Sh?&r>%;uoU2P3y`k%CC6dm%jCJKC{=q@FTB$d*}BC<Hfp5@FcW^
zOFW@n)U#dq-|H`=b%<UtUx(=F%KaVf^;w;*e2t=CM|p$4>zm~7`YoG-#!d07QNyO|
zRf<21>Ph40-pJAY)k|p}^WGd!cSlhtO;RWHoAvBBIq%JJ-hDpi<EE&2*{b6-s`*KL
zefS|${Og-RJio)=m7!kwyV6@Nc_0to9>@cEsN!{3@4uf*{rGd~xv%3~|M3P}xb<-;
zy{4$~J8I{O_<@iK`a$Z!>wr9vhbs1is!K(<?RDDhdGK+wujTn?TQ&u~fdBD^P`K5+
z0E(I)pyy1`xG8=wYMA>S-0!$pOZY%<P~?F;gbXnLO?9By?0!SvQ*pJwuh#2QWZf=5
zGmG8&c&ew;c&ew;b3eQ6Yf&RP?%oxoc?RHs+JoY1w+np^MLti+1ipoO@H!w5<e?d^
z<9#J+z!ZHxl5<@D@dm!^Tlq~ze~$Zv`S{tf9B++cems1LbC*1j2d<y+93LE58V5f5
zFlfRQeI25a^t{UKdE1uV5bz7Gm;ryj`=vAw&`W7voM-Lte+%l*(dX+qXFi|7WR6{(
zK_lFc<9^)5I^%kdH#YJ>9zq7V4zD^;tgFBBaVs6BsP(EnFKinn13v%E`ZQ1A`s8cZ
z<aNh&UNsI(*8vrEJ+NqgMYR*fy0#a6e?{^TGJ&t49=s0719@nq`(+<~5HxOzzl<7A
z{pS-oYCU+bZ??^vfzL<XZ;D#i;oDKIrs(tKykc+tJ!l|De}Aw2OVm`3K2Id;KH!(U
ze;^O!f$K2b*Q+WAG!LBO$=0C0<cIE^6`R!$@Vyq70rx>WQl0EbeyYXV+oE#x^*&D6
zlY4`Pa{TwGzPw)V%yoPHQ*j-eJT%ArY@A2VBiB)@&vRdg;JhjNIF)Y3Q7_IW@W1BB
zfX`dCH_cnMH_cnMGp%QKKF7PGy3_m)v)4zXf1!V&f2n@|LOZ3MHex&dX;i-{YQ0t6
zM^V)MlsURReoO`yl>wihX^Sa(-L=@6(@{BY+z`}m=gtOA=J?*bL5v%(>UC)^3i@&K
zfKO;dpHR2DefOQ9VN;Yp^5?$yr{})+=k>i8)#>q{$iN~qP;`CC=P4L4MP2XoUqk;(
z9_WARf8hWeSR4oZx$P-a{B_i*ZQT;olcUcw*qNVe(sNFJ9cc~s-4@yJ^WPSi0iQRb
z(-cpgj561Bj+)mAJ}iz8<c+-HukcrJz&Y^SsBu&LT-0#BpP~C3y078)JswTtZXZql
z>CgD3#dZ7quQ@WXV|&mWX*|*!dETavKMI;OMXyVqgX8{MC@AI)@Mc9G$b+{B@=#75
z)ZZwo-}!P>zbXDG>TF&YP4nYueKfz$#Pe=mmuwfB0nNLlsQU8zY18RBg6ZU!cpesC
z<oqKK<RSQD=MDalJa~Ip={)%TxQTQ>ZX%Ck^8Tzp%`e!W{L|Sq|MF`&;-mbGvAsM4
zzMgJhdY-M%eh_sg^KPZ<T)$6|{x@WTeusMSIv@|^p^@fwRJ|&$HqQ?q3YsuQAJ;Tu
z53CF7&(WVpI-B~%)*Ln8pvFHbG7iemaHDN?-W2`%<o$L})P)@XIcf-B6mEm>=la-c
zb4)x{<nx3~;FGBbuLJV1ym`1OYRVLEiyBSy2@dAC^NFBd`)*Wwj(**tHC?xQI7dA%
z?&pu~YGuHmM{7^>9kplw>z{|~;p71yh~Iz%Z~zW0g#&u-PEpU_{VZz06n_$R-gF;N
z@tUY-7*AfZQPOoiMO`0ao+8@&QrbKH1@%My;2Yt<;y6$=UU%=Fpg~jAbtmt)SwF5U
zd*s)JhSGJRA^Yo_LA~4;3+2pxG2=c|XldMsQvXIhP!FL!P!C=QE1id;aV$P=ocnHm
zKIXsm>vZSRd|Btr*OTqa^}W*PSHDjY-iA!j|56WL2jrnzc{p}7XxbF-i()(>-1a(M
zsXX{Nqn@-rMUQRT7}UXiGzVtk*5|wGFhyUFxE24y^`DD1gYWW&LmtRO$iVFD9uKwz
zjhdq7_4mKWpHAaopH9zlAGX`0dUN!@=TE18Ya&PQzlPy39FBE2myhysErV%3%0as&
zO80TR^|GCv0e?>LWO`2UWS(Ea-#NRz^E&s%@I|x>IN%)k``?24O;Oj=XOs5IulTr*
z+4BMV>nBXn*M)y6|Gt^O*XMy_o$k8T=lOF{rr6yT#r2;nH@mCwKlmB^41T6={Y6>d
zck@F*r%h4MMbClVW&5%Fq(3FUbSd5(^@b_>I!Z6)`I=|W+pP646tyl!S$X3)i|cnf
zJA#Hy@s6nebpPg&95o(j*$S;m4y^lxik`QVTz3uS$@pi+KgV{?{V(-bb6tnm`xKk~
z{XSpJ>uKJY*VDRKFE;Jp-R_@`CLeV)uQRXyLs9+6x1;(@vFA)w8lN(pqx{!5qAr-C
zp3nJx)TAlydn#xoUGJF6(XV@q=JgNse3f4}8%yg*P3NfbTL1BJ(3B~@|6b5&x*sx`
zV^?R;i21yU<2imcYB=8qoq0c0&v7a0d9H0+gEYVOS@Yv}<~s8E`(Cv54+QaCymO_z
zt=IGURL|#AJ$K`O@IM!ejDEi^{kXoTVzc`eeNV;3eqX=d!t=<UN48|3u`(dvs_5gN
z`|v;bpNoZz|MAvL9_pTlA4W}>qV7L57VY`8WT2iI@at$T>AG4=8YjL!&-XoZ-fzzm
z<08ld_W|*NaNwyYgN98}<1XsioA9m4Kx1W~y)7!OW71}=C!)B{!*!mEh0OI0Z_VU^
zJcJC)_Ip=eSK712mJE=AYG&Z>s4-K#Ginh3ga5f$u=M-hYRLn62pO2|e>ARFQR90(
zUu?+$8K^=AG=H<Aj~hRQ&%x(hEKmA(Z=>XaJcJC)_BqOnqQ*tL?|D8wuk?JL*U;a`
zmJE=AmCQg<d;CMxlqqWcS^sSFZ|}UHaWD8K+DqLy;OlZvo8s4_#`3sHjib~$(ksb)
zeis?2TLygogXhxvq|ce(hdyo}eHb*6<6TiBT<76B&&4Lf^$q6)c_0rV1M~N>M-B%~
zn&Nx!292co4PVZ2<A$K;5^p`rY{>u_SX>7D?}@hbTvl72@BCk)rcBZAYmMS_@HrO?
zjsD$R8+jlPAp`UGIm(x!?sI8hajw^W9_6t#Z`zpM8a0^Lk)AmY!RJwCo?g!{+h)yx
zKTm!nJzst#T|eEIp1&H*QS(5my(!Ms?({xI-XAhSKS(`z9gqj|P{sM%eSNSgQ`Gnb
zAE!Bx=J_2+_un)x?d)|ceBQ;`>xFsUvn2!NW<c{$D{4OKuS5-)qSha?sA&&2o&P@Y
zPqp%Yr|^3j&x1dt-NAw8;J`PdCQMQD60KVs)Mtu5Ut>?cpD=SBRgGg(^q;TczFxWZ
zVcC!7zU2ONo#p<#{+rgzRot*X=-IUX$N3!nx$tiL$NNEJIr_XgJSWcerBK3LU-E`b
z9>{~Y2lCL^dC+wzMO}}&Eo#&he-<@RxL&IFuid8T{cM~0Jdr1J^!2V@vumSX$kEr~
zemdVL@&0_DT@$t2Zisq3N8LyC*SDnCx8!l^Mc2pV?-b=5-2eSKjX(Kw@((-GbAG#W
z)aTGRT*d374w}bTa@07n??jz4MPGO7jkJE)nH>Flb)<at=lG+jv!?NAihjHS+rBMm
zDo3q9ruL~wKkn_7JdlUg7%!oAs7M`Bht%O1<2f16c@YHsT}i*-{Q-F(4^<q$p}JJW
zzj>XK2l9ZA#mB+{H~<IW033h=Z~zX#0XP5$-~b$e18@KizyUY_2jBo4fCF#<4!{98
z00-az9DoCG01m(bH~<IW033h=Z~zX#0XP5$-~b$e18@KizyUY_2jBo4fCF#<4!{98
z00-az9DoCG01m(bH~<IW033h=Z~zX#0XP5$-~b$e18@KizyUY_2jBo4fCF#<4!{98
z00-az9DoCG01m(bH~<IW033h=Z~zX#0XP5$-~b$e18@KizyUY_2jBo4fCF#<4!{98
z00-az9DoCG01m(bH~<IW033h=Z~zX#0XP5$-~b$e18@KizyUY_2jBo4fCF#<4!{98
K00)XW@c#q2Wx>Dz

literal 0
HcmV?d00001


From ba36bfecaa1d695f20dd279d21860752e9d62696 Mon Sep 17 00:00:00 2001
From: chaoming <adaduo@outlook.com>
Date: Tue, 9 Jan 2024 11:21:07 +0800
Subject: [PATCH 2/4] updater version

---
 brainpy/__init__.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/brainpy/__init__.py b/brainpy/__init__.py
index c8f834c6d..98b69c2ef 100644
--- a/brainpy/__init__.py
+++ b/brainpy/__init__.py
@@ -1,6 +1,6 @@
 # -*- coding: utf-8 -*-
 
-__version__ = "2.4.6.post5"
+__version__ = "2.5.0.post1"
 
 # fundamental supporting modules
 from brainpy import errors, check, tools

From 15ffba6751093f923da0ffb19b864e5927e3b55e Mon Sep 17 00:00:00 2001
From: chaoming <adaduo@outlook.com>
Date: Tue, 9 Jan 2024 11:25:09 +0800
Subject: [PATCH 3/4] update version

---
 brainpy/__init__.py | 3 ++-
 setup.py            | 4 +++-
 2 files changed, 5 insertions(+), 2 deletions(-)

diff --git a/brainpy/__init__.py b/brainpy/__init__.py
index 98b69c2ef..a3a1de694 100644
--- a/brainpy/__init__.py
+++ b/brainpy/__init__.py
@@ -1,6 +1,7 @@
 # -*- coding: utf-8 -*-
 
-__version__ = "2.5.0.post1"
+
+__version__ = "2.5.0"
 
 # fundamental supporting modules
 from brainpy import errors, check, tools
diff --git a/setup.py b/setup.py
index b9f51dd6b..d03fd91fd 100644
--- a/setup.py
+++ b/setup.py
@@ -3,6 +3,7 @@
 import io
 import os
 import re
+import time
 
 from setuptools import find_packages
 from setuptools import setup
@@ -26,6 +27,7 @@
 except ModuleNotFoundError:
   pass
 
+
 # version
 here = os.path.abspath(os.path.dirname(__file__))
 with open(os.path.join(here, 'brainpy', '__init__.py'), 'r') as f:
@@ -42,7 +44,7 @@
 # setup
 setup(
   name='brainpy',
-  version=version,
+  version=version + '.post{}'.format(time.strftime("%Y%m%d", time.localtime())),
   description='BrainPy: Brain Dynamics Programming in Python',
   long_description=README,
   long_description_content_type="text/markdown",

From 07559d1795d815150c5cb42d765a94e37def78c2 Mon Sep 17 00:00:00 2001
From: chaoming <adaduo@outlook.com>
Date: Tue, 9 Jan 2024 12:07:20 +0800
Subject: [PATCH 4/4] fix bug

---
 brainpy/_src/math/tests/test_einops.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/brainpy/_src/math/tests/test_einops.py b/brainpy/_src/math/tests/test_einops.py
index e6738009e..2f018d973 100644
--- a/brainpy/_src/math/tests/test_einops.py
+++ b/brainpy/_src/math/tests/test_einops.py
@@ -77,7 +77,6 @@ def test_rearrange_consistency_numpy():
   ]:
     result = ein_rearrange(x, pattern)
     assert len(numpy.setdiff1d(x, result)) == 0
-    assert result.dtype == x.dtype
 
   result = ein_rearrange(x, "a b c d e f -> a (b) (c d e) f")
   assert numpy.array_equal(x.flatten(), result.flatten())