Merge pull request #38 from Nonannet/dev

Dev

.github/workflows/ci.yml

@@ -299,7 +299,7 @@ jobs:
 
   build-windows:
     needs: [build_stencils]
-    runs-on: windows-latest
+    runs-on: windows-2022
 
     strategy:
       matrix:

README.md

@@ -39,7 +39,7 @@ Despite missing SIMD-optimization, benchmark performance shows promising numbers
 
 For the benchmark (`tests/benchmark.py`) the timing of 30000 iterations for calculating the therm `sum((v1 + i) @ v2 for i in range(10))` where measured on an Ryzen 5 3400G. Where the vectors `v1` and `v2` both have a lengths of `v_size` which was varied according to the chart from 10 to 500. For the NumPy case the "i in range(10)" loop was vectorized like this: `np.sum((v1 + i) @ v2)` with i being here a `NDArray` with a dimension of `[10, 1]`. The number of calculated scalar operations is the same for both contenders. Obviously Copapy profits from less overheat by calling a single function from python per iteration, where the NumPy variant requires 3. Interestingly there is no indication visible in the chart that for increasing `v_size` the calling overhead for NumPy will be compensated by using faster SIMD instructions. It is to note that in this benchmark the Copapy case does not move any data between python and the compiled code.
 
-Furthermore for many applications copypy will benefit by reducing the actual number of operations significantly compared to a NumPy implementation, by precompute constant values know at compile time and benefiting from sparcity. Multiplying by zero (e.g. in a diagonal matrix) eliminate a hole branch in the computation graph. Operations without effect, like multiplications by 1 oder additions with zero gets eliminated at compile time.
+Furthermore for many applications Copapy performance will benefit by reducing the actual number of operations significantly compared to a NumPy implementation, by precompute constant values know at compile time and benefiting from sparcity. Multiplying by zero (e.g. in a diagonal matrix) eliminate a hole branch in the computation graph. Operations without effect, like multiplications by 1 oder additions with zero gets eliminated at compile time.
 
 For Testing and using Copapy to speed up computations in conventional Python programs there is also the `@cp.jit` decorator available, to compile functions on first use and cache the compiled version for later calls:
 
@@ -188,7 +188,7 @@ For more complex operations - where inlining is less useful - stencils call a no
 
 Unlike stencils, non-stencil functions like `sinf` are not stripped and do not need to be tail-call-optimizable. These functions can be provided as C code and compiled together with the stencils or can be object files like in the case of `sinf` compiled from C and assembly code and merged into the stencil object files. Math functions like `sinf` are currently provided by the MUSL C library, with architecture-specific optimizations. 
 
-Non-stencil functions and constants are stored together with the stencils in an ELF object file for each supported CPU architecture. The required non-stencil functions and constants are bundled during compilation. The compiler includes only the data and code required for a specific Copapy program.
+Non-stencil functions and constants used by them are stored together with the stencils in ELF object files for each supported CPU architecture. The required non-stencil functions and constants are bundled during compilation. The compiler includes only the data and code required for a specific Copapy program.
 
 The Copapy compilation process is independent of the actual instruction set. It relies purely on relocation entries and symbol metadata from the ELF file generated by the C compiler.
 

docs/source/generate_class_list.py

@@ -99,6 +99,8 @@ if __name__ == "__main__":
 
         write_functions(f, ['*'], 'copapy', title='Tensor/Matrix functions', path_patterns=['*_tensors*'], api_dir=api_dir)
 
+        write_functions(f, ['*'], 'copapy', title='ANN functions', path_patterns=['*_nn*'], api_dir=api_dir)
+
         #write_manual(f, ['NumLike'], title='Types')
 
     with open(f'{api_dir}/backend.md', 'w') as f:

pyproject.toml

@@ -30,7 +30,7 @@ where = ["src"]
 copapy = ["obj/*.o", "py.typed"]
 
 [tool.setuptools_scm]
-version_scheme = "post-release"
+version_scheme = "no-guess-dev"
 local_scheme = "no-local-version"
 tag_regex = "^v(?P<version>\\d+\\.\\d+\\.\\d+(?:-beta)?)$"
 fallback_version = "0.0.0"

src/copapy/__init__.py

@@ -37,8 +37,9 @@ from ._target import Target, jit
 from ._basic_types import NumLike, value, generic_sdb, iif
 from ._vectors import vector, distance, scalar_projection, angle_between, rotate_vector, vector_projection
 from ._quaternion import quaternion
-from ._tensors import tensor, zeros, ones, arange, eye, identity, diagonal
+from ._tensors import tensor, zeros, ones, arange, eye, identity, diagonal, concat
-from ._math import sqrt, abs, sign, sin, cos, tan, asin, acos, atan, atan2, log, exp, pow, get_42, clamp, min, max, relu
+from ._math import sqrt, abs, sign, sin, cos, tan, asin, acos, atan, atan2, log, exp, pow, get_42, clamp, minimum, maximum
+from ._nn import relu, sigmoid
 from ._autograd import grad
 from ._tensors import tensor as matrix
 from ._version import __version__  # Run "pip install -e ." to generate _version.py
@@ -74,16 +75,18 @@ __all__ = [
     "pow",
     "get_42",
     "clamp",
-    "min",
+    "minimum",
-    "max",
+    "maximum",
     "relu",
     "distance",
     "scalar_projection",
     "angle_between",
     "rotate_vector",
-"vector_projection",
+    "vector_projection",
     "quaternion",
     "grad",
     "eye",
+    "concat",
+    "sigmoid",
     "jit"
 ]

src/copapy/_math.py

@@ -393,12 +393,20 @@ def clamp(x: U | value[U] | vector[U], min_value: U | value[U], max_value:  U |
 
 
 @overload
-def min(x: value[U], y: U | value[U]) -> value[U]: ...
+def minimum(x: U, y: U) -> U: ...
 @overload
-def min(x: U | value[U], y: value[U]) -> value[U]: ...
+def minimum(x: value[U], y: U | value[U]) -> value[U]: ...
 @overload
-def min(x: U, y: U) -> U: ...
+def minimum(x: U | value[U], y: value[U]) -> value[U]: ...
-def min(x: U | value[U], y: U | value[U]) -> Any:
+@overload
+def minimum(x: vector[U], y: U | value[U] | vector[U]) -> vector[U]: ...
+@overload
+def minimum(x: U | value[U] | vector[U], y: vector[U]) -> vector[U]: ...
+@overload
+def minimum(x: tensor[U], y: U | value[U] | tensor[U]) -> tensor[U]: ...
+@overload
+def minimum(x: U | value[U] | tensor[U], y: tensor[U]) -> tensor[U]: ...
+def minimum(x: U | value[U] | vector[U] | tensor[U], y: U | value[U] | vector[U] | tensor[U]) -> Any:
     """Minimum function to get the smaller of two values.
 
     Arguments:
@@ -408,22 +416,30 @@ def min(x: U | value[U], y: U | value[U]) -> Any:
     Returns:
         Minimum of x and y
     """
-    if isinstance(x, value):
+    if isinstance(x, tensor) or isinstance(y, tensor):
+        return _map2_tensor(x, y, minimum)
+    if isinstance(x, vector) or isinstance(y, vector):
+        return _map2_vector(x, y, minimum)
+    if isinstance(x, value) or isinstance(y, value):
         return add_op('min', [x, y])
-    if isinstance(x, tensor):
-        return _map2_tensor(x, y, min)
-    if isinstance(x, vector):
-        return _map2_vector(x, y, min)
     return x if x < y else y
 
 
 @overload
-def max(x: value[U], y: U | value[U]) -> value[U]: ...
+def maximum(x: U, y: U) -> U: ...
 @overload
-def max(x: U | value[U], y: value[U]) -> value[U]: ...
+def maximum(x: value[U], y: U | value[U]) -> value[U]: ...
 @overload
-def max(x: U, y: U) -> U: ...
+def maximum(x: U | value[U], y: value[U]) -> value[U]: ...
-def max(x: U | value[U], y: U | value[U]) -> Any:
+@overload
+def maximum(x: vector[U], y: U | value[U] | vector[U]) -> vector[U]: ...
+@overload
+def maximum(x: U | value[U] | vector[U], y: vector[U]) -> vector[U]: ...
+@overload
+def maximum(x: tensor[U], y: U | value[U] | tensor[U]) -> tensor[U]: ...
+@overload
+def maximum(x: U | value[U] | tensor[U], y: tensor[U]) -> tensor[U]: ...
+def maximum(x: U | value[U] | vector[U] | tensor[U], y: U | value[U] | vector[U] | tensor[U]) -> Any:
     """Maximum function to get the larger of two values.
 
     Arguments:
@@ -433,12 +449,12 @@ def max(x: U | value[U], y: U | value[U]) -> Any:
     Returns:
         Maximum of x and y
     """
-    if isinstance(x, value):
+    if isinstance(x, tensor) or isinstance(y, tensor):
+        return _map2_tensor(x, y, maximum)
+    if isinstance(x, vector) or isinstance(y, vector):
+        return _map2_vector(x, y, maximum)
+    if isinstance(x, value) or isinstance(y, value):
         return add_op('max', [x, y])
-    if isinstance(x, tensor):
-        return _map2_tensor(x, y, max)
-    if isinstance(x, vector):
-        return _map2_vector(x, y, max)
     return x if x > y else y
 
 
@@ -470,20 +486,6 @@ def lerp(v1: U | value[U] | vector[U], v2: U | value[U] | vector[U], t:  unifloa
     return v1 * (1 - t) + v2 * t
 
 
-@overload
-def relu(x: U) -> U: ...
-@overload
-def relu(x: value[U]) -> value[U]: ...
-@overload
-def relu(x: vector[U]) -> vector[U]: ...
-@overload
-def relu(x: tensor[U]) -> tensor[U]: ...
-def relu(x: U | value[U] | vector[U] | tensor[U]) -> Any:
-    """Returns x for x > 0 and otherwise 0."""
-    ret = x * (x > 0)
-    return ret
-
-
 def _map2_vector(self: VecNumLike, other: VecNumLike, func: Callable[[Any, Any], value[U] | U]) -> vector[U]:
     """Applies a function to each element of the vector and a second vector or scalar."""
     if isinstance(self, vector) and isinstance(other, vector):
@@ -499,9 +501,9 @@ def _map2_vector(self: VecNumLike, other: VecNumLike, func: Callable[[Any, Any],
 def _map2_tensor(self: TensorNumLike, other: TensorNumLike, func: Callable[[Any, Any], value[U] | U]) -> tensor[U]:
     """Applies a function to each element of the vector and a second vector or scalar."""
     if isinstance(self, vector):
-        self = tensor(self.values, (len(self.values),))
+        self = tensor(self)
     if isinstance(other, vector):
-        other = tensor(other.values, (len(other.values),))
+        other = tensor(other)
     if isinstance(self, tensor) and isinstance(other, tensor):
         assert self.shape == other.shape, "Tensors must have the same shape"
         return tensor([func(x, y) for x, y in zip(self.values, other.values)], self.shape)

src/copapy/_nn.py

@@ -0,0 +1,34 @@
+from . import vector
+from . import tensor 
+from . import value
+from typing import TypeVar, Any, overload
+import copapy as cp
+
+U = TypeVar("U", int, float)
+
+
+@overload
+def relu(x: U) -> U: ...
+@overload
+def relu(x: value[U]) -> value[U]: ...
+@overload
+def relu(x: vector[U]) -> vector[U]: ...
+@overload
+def relu(x: tensor[U]) -> tensor[U]: ...
+def relu(x: U | value[U] | vector[U] | tensor[U]) -> Any:
+    """Returns x for x > 0 and otherwise 0."""
+    ret = x * (x > 0)
+    return ret
+
+
+@overload
+def sigmoid(x: U) -> float: ...
+@overload
+def sigmoid(x: value[U]) -> value[float]: ...
+@overload
+def sigmoid(x: vector[U]) -> vector[float]: ...
+@overload
+def sigmoid(x: tensor[U]) -> tensor[float]: ...
+def sigmoid(x: U | value[U] | vector[U] | tensor[U]) -> Any:
+    """Sigmoid function to map any value to the range (0, 1)."""
+    return 1 / (1 + cp.exp(-x))

src/copapy/_tensors.py

@@ -32,7 +32,7 @@ class tensor(ArrayType[TNum]):
             self.shape: tuple[int, ...] = tuple(shape)
             assert (isinstance(values, Sequence) and
                     any(isinstance(v, (value, int, float)) for v in values)), \
-                    "Values must be a sequence of values if shape is provided"
+                    "Values must be a sequence of scalars if shape is provided"
             self.values: tuple[TNum | value[TNum], ...] = tuple(v for v in values if not isinstance(v, Sequence))
             self.ndim: int = len(shape)
         elif isinstance(values, (int, float)):
@@ -128,7 +128,14 @@ class tensor(ArrayType[TNum]):
         return self.shape[0]
 
     def get_scalar(self: 'tensor[TNum]', *key: int) -> TNum | value[TNum]:
-        """Get a single scalar value from the tensor given multi-dimensional indices."""
+        """Get a single scalar value from the tensor given multi-dimensional indices.
+
+        Arguments:
+            *key: Variable number of indices for each dimension.
+
+        Returns:
+            The scalar value at the specified indices.
+        """
         assert len(key) == self.ndim, f"Expected {self.ndim} indices, got {len(key)}"
         flat_idx = self._get_flat_index(key)
         return self.values[flat_idx]
@@ -567,7 +574,11 @@ class tensor(ArrayType[TNum]):
     @overload
     def trace(self: 'tensor[float]') -> float | value[float]: ...
     def trace(self) -> Any:
-        """Calculate the trace (sum of diagonal elements)."""
+        """Calculate the trace (sum of diagonal elements).
+
+        Returns:
+            The sum of diagonal elements.
+        """
         assert self.ndim == 2 and self.shape[0] == self.shape[1], "Trace is only defined for square matrices"
         return mixed_sum(self.get_scalar(i, i) for i in range(self.shape[0]))
 
@@ -615,7 +626,11 @@ class tensor(ArrayType[TNum]):
         return self.reshape(-1)
 
     def size(self) -> int:
-        """Return total number of elements."""
+        """Count number of elements over all dimensions.
+
+        Returns:
+            Total number of elements.
+        """
         size = 1
         for dim in self.shape:
             size *= dim
@@ -820,7 +835,11 @@ class tensor(ArrayType[TNum]):
         return tensor(result_vals, self.shape)
 
     def homogenize(self) -> 'tensor[TNum]':
-        """Convert all elements to copapy values if any element is a copapy value."""
+        """Convert all elements to Copapy values if any element is a Copapy value.
+
+        Returns:
+            Tensor with all elements converted to Copapy values, or the input tensor if no value is a Copapy value.
+        """
         if any(isinstance(val, value) for val in self.values):
             homogenized: tuple[value[Any], ...] = tuple(value(val) if not isinstance(val, value) else val for val in self.values)
             return tensor(homogenized, self.shape)
@@ -833,7 +852,13 @@ class tensor(ArrayType[TNum]):
 
 
 def zeros(shape: Sequence[int] | int) -> tensor[int]:
-    """Create a zero tensor of given shape."""
+    """Create a zero tensor of given shape.
+    
+    Arguments:
+        shape: shape of the tensor to create.
+
+    Returns:
+        New tensor of given shape, initialized with zeros."""
     if isinstance(shape, int):
         shape = (shape,)
 
@@ -845,7 +870,13 @@ def zeros(shape: Sequence[int] | int) -> tensor[int]:
 
 
 def ones(shape: Sequence[int] | int) -> tensor[int]:
-    """Create a tensor of ones with given shape."""
+    """Create a tensor of ones with given shape.
+    
+    Arguments:
+        shape: shape of the tensor to create.
+
+    Returns:
+        New tensor of given shape, initialized with ones."""
     if isinstance(shape, int):
         shape = (shape,)
 
@@ -856,6 +887,91 @@ def ones(shape: Sequence[int] | int) -> tensor[int]:
     return tensor([1] * size, tuple(shape))
 
 
+@overload
+def concat(tensors: Sequence[vector[U]]) -> vector[U]: ...
+@overload
+def concat(tensors: Sequence[tensor[U]], axis: int = 0) -> tensor[U]: ...
+def concat(tensors: Sequence[tensor[U] | vector[U]], axis: int = 0) -> tensor[U] | vector[U]:
+    """Concatenate tensors or vectors along a specified axis.
+
+    Arguments:
+        tensors: Tensors or vectors to concatenate. Must have the same shape except for the specified axis.
+        axis: Axis along which to concatenate (default 0).
+
+    Returns:
+        A new tensor or vector resulting from concatenation.
+    """
+    assert tensors, "At least one tensor must be provided"
+
+    # Check if all inputs are vectors
+    all_vectors = all(isinstance(item, vector) for item in tensors)
+
+    # Convert vectors to tensors for uniform processing
+    tensor_list: list[tensor[U]] = []
+    for item in tensors:
+        if isinstance(item, vector):
+            tensor_list.append(tensor(item.values, item.shape))
+        else:
+            tensor_list.append(item)
+
+    first_shape = tensor_list[0].shape
+    ndim = len(first_shape)
+
+    if axis < 0:
+        axis += ndim
+
+    assert 0 <= axis < ndim
+
+    # Shape checks
+    for t in tensor_list:
+        assert len(t.shape) == ndim
+        for i in range(ndim):
+            if i != axis:
+                assert t.shape[i] == first_shape[i]
+
+    # Output shape
+    new_shape = list(first_shape)
+    new_shape[axis] = sum(t.shape[axis] for t in tensor_list)
+
+    # Compute block sizes
+    inner_block: int = 1
+    for s in first_shape[axis + 1:]:
+        inner_block *= s
+
+    outer_block: int = 1
+    for s in first_shape[:axis]:
+        outer_block *= s
+
+    new_values: list[value[U] | U] = []
+
+    for outer in range(outer_block):
+        for t in tensor_list:
+            axis_size = t.shape[axis]
+            start = outer * axis_size * inner_block
+            end = start + axis_size * inner_block
+            new_values.extend(t.values[start:end])
+
+    result_tensor = tensor(new_values, tuple(new_shape))
+
+    # If all inputs were vectors and result is 1D, return as vector
+    if all_vectors and result_tensor.ndim == 1:
+        return vector(result_tensor.values)
+
+    return result_tensor
+
+
+def flatten(t: tensor[U]) -> tensor[U]:
+    """Flatten a tensor to a 1D tensor.
+
+    Arguments:
+        t: n-dimensional tensor.
+
+    Returns:
+        A 1D tensor containing all elements from the input tensor.
+    """
+    return t.flatten()
+
+
 def arange(start: int | float, stop: int | float | None = None,
            step: int | float = 1) -> tensor[int] | tensor[float]:
     """Create a tensor with evenly spaced values.
@@ -927,10 +1043,10 @@ def identity(size: int) -> tensor[int]:
 
 
 @overload
-def diagonal(vec: 'tensor[int] | vector[int]') -> tensor[int]: ...
+def diagonal(vec: tensor[int] | vector[int]) -> tensor[int]: ...
 @overload
-def diagonal(vec: 'tensor[float] | vector[float]') -> tensor[float]: ...
+def diagonal(vec: tensor[float] | vector[float]) -> tensor[float]: ...
-def diagonal(vec: 'tensor[Any] | vector[Any]') -> 'tensor[Any]':
+def diagonal(vec: tensor[Any] | vector[Any]) -> tensor[Any]:
     """Create a diagonal tensor from a 1D tensor.
 
     Arguments:

src/copapy/_vectors.py

@@ -54,7 +54,14 @@ class vector(ArrayType[TNum]):
         return iter(self.values)
 
     def get_scalar(self, index: int) -> TNum | value[TNum]:
-        """Get a single scalar value from the vector."""
+        """Get a single scalar value from the vector.
+
+        Arguments:
+            index: The index of the element to retrieve.
+
+        Returns:
+            The scalar value at the specified index.
+        """
         return self.values[index]
 
     @overload
@@ -200,6 +207,14 @@ class vector(ArrayType[TNum]):
     @overload
     def dot(self, other: 'vector[int] | vector[float]') -> float | int | value[float] | value[int]: ...
     def dot(self, other: 'vector[int] | vector[float]') -> Any:
+        """Calculate the dot product of this vector with another.
+
+        Arguments:
+            other: Another vector of the same length.
+
+        Returns:
+            The dot product as a scalar value.
+        """
         assert len(self.values) == len(other.values), "Vectors must be of same length."
         return mixed_sum(a * b for a, b in zip(self.values, other.values))
 
@@ -216,7 +231,14 @@ class vector(ArrayType[TNum]):
         return self.dot(other)
 
     def cross(self: 'vector[float]', other: 'vector[float]') -> 'vector[float]':
-        """3D cross product"""
+        """Calculate the cross product of this vector with another 3D vector.
+
+        Arguments:
+            other: Another 3D vector.
+
+        Returns:
+            A new vector perpendicular to both vectors.
+        """
         assert len(self.values) == 3 and len(other.values) == 3, "Both vectors must be 3-dimensional."
         a1, a2, a3 = self.values
         b1, b2, b3 = other.values
@@ -286,20 +308,37 @@ class vector(ArrayType[TNum]):
     @overload
     def sum(self: 'vector[float]') -> float | value[float]: ...
     def sum(self) -> Any:
-        """Sum of all vector elements."""
+        """Sum of all vector elements.
+
+        Returns:
+            The sum of all vector elements as a scalar value.
+        """
         return mixed_sum(self.values)
 
     def magnitude(self) -> 'float | value[float]':
-        """Magnitude (length) of the vector."""
+        """Magnitude (length) of the vector.
+
+        Returns:
+            The magnitude of the vector as a scalar value.
+        """
         s = mixed_sum(a * a for a in self.values)
         return cp.sqrt(s)
 
     def normalize(self) -> 'vector[float]':
-        """Returns a normalized (unit length) version of the vector."""
+        """Calculate a normalized (unit length) version of the vector.
+
+        Returns:
+            A normalized (unit length) vector.
+        """
         mag = self.magnitude() + epsilon
         return self / mag
 
     def homogenize(self) -> 'vector[TNum]':
+        """Convert all elements to Copapy values if any element is a Copapy value.
+
+        Returns:
+            Vector with homogeneous element type.
+        """
         if any(isinstance(val, value) for val in self.values):
             return vector(mixed_homogenize(self))
         else:

tests/test_math.py

@@ -23,8 +23,8 @@ def test_fine():
                 cp.abs(-c_f),
                 cp.sign(c_i),
                 cp.sign(-c_f),
-                cp.min(c_i, 5),
+                cp.minimum(c_i, 5),
-                cp.max(c_f, 5))
+                cp.maximum(c_f, 5))
 
     re2_test = (a_f ** 2,
                 a_i ** -1,
@@ -39,8 +39,8 @@ def test_fine():
                 cp.abs(-a_f),
                 cp.sign(a_i),
                 cp.sign(-a_f),
-                cp.min(a_i, 5),
+                cp.minimum(a_i, 5),
-                cp.max(a_f, 5))
+                cp.maximum(a_f, 5))
 
     ret_refe = (a_f ** 2,
                 a_i ** -1,

tests/test_tensor_basic.py

@@ -3,6 +3,7 @@
 
 import copapy as cp
 
+
 def test_tensor_basic():
     # Test 1: Create a scalar tensor
     print("Test 1: Scalar tensor")
@@ -225,6 +226,7 @@ def test_tensor_basic():
     assert t.ndim == 3
     print()
 
+
 def test_tensor_slicing():
     print("Test Numpy-style slicing")
     t = cp.tensor([[10, 20, 30], [40, 50, 60], [70, 80, 90]])
@@ -251,6 +253,25 @@ def test_tensor_slicing():
     assert slice4[2] == 90
     print()
 
+
+def test_tensor_concat():
+    print("Test tensor concatenation")
+    t1 = cp.tensor([[1, 2], [3, 4]])
+    t2 = cp.tensor([[5, 6], [7, 8]])
+    t3 = cp.tensor([[9, 10], [11, 12]])
+
+    concat_0 = cp.concat([t1, t2, t3], axis=0)
+    print(f"Concatenate along axis 0: shape={concat_0.shape}")
+    assert concat_0.shape == (6, 2)
+    assert concat_0[4, 1] == 10
+
+    concat_1 = cp.concat([t1, t2, t3], axis=1)
+    print(f"Concatenate along axis 1: shape={concat_1.shape}")
+    assert concat_1.shape == (2, 6)
+    assert concat_1[1, 4] == 11
+    print()
+
+
 if __name__ == "__main__":
     test_tensor_basic()
     print("All tests completed!")

tests/test_vector.py

@@ -114,6 +114,19 @@ def test_sort_vector():
     assert ref == result
 
 
+def test_vector_concat():
+    print("Test vector concatenation")
+    v1 = cp.vector([1, 2, 3])
+    v2 = cp.vector([4, 5])
+    v3 = cp.vector([6])
+
+    concat_vec = cp.concat([v1, v2, v3])
+    print(f"Concatenate vectors: shape={concat_vec.shape}")
+    assert concat_vec.shape == (6,)
+    assert concat_vec[4] == 5
+    print()
+
+
 if __name__ == "__main__":
     #test_vectors_init()
     #test_compiled_vectors()