Merge pull request #13 from Nonannet/dev

Dev

.github/workflows/build_wheels.yml

@@ -65,21 +65,21 @@ jobs:
           name: wheels-${{ matrix.os }}
           path: wheelhouse/*.whl
 
-#  publish:
+  publish:
-#    if: contains(github.ref, '-beta') == false
+    if: contains(github.ref, '-beta') == false
-#    needs: [build_wheels]
+    needs: [build_wheels]
-#    runs-on: ubuntu-latest
+    runs-on: ubuntu-latest
-#    steps:
+    steps:
-#      - name: Install Twine
+      - name: Install Twine
-#        run: pip install twine
+        run: pip install twine
 
-#      - uses: actions/download-artifact@v4
+      - uses: actions/download-artifact@v4
-#        with:
+        with:
-#          path: wheelhouse
+          path: wheelhouse
 
-#      - name: Publish to PyPI
+      - name: Publish to PyPI
-#        env:
+        env:
-#          TWINE_USERNAME: __token__
+          TWINE_USERNAME: __token__
-#          TWINE_PASSWORD: ${{ secrets.PYPI_API_TOKEN }}
+          TWINE_PASSWORD: ${{ secrets.PYPI_API_TOKEN }}
-#        
+        
-#        run: python -m twine upload wheelhouse/*
+        run: python -m twine upload wheelhouse/*

.github/workflows/ci.yml

@@ -30,11 +30,6 @@ jobs:
           name: musl-object-files
           path: /object_files/musl_objects_*.*o
 
-      - uses: actions/upload-artifact@v4
-        with:
-          name: cross-runner
-          path: build/runner/coparun-*
-
   build-ubuntu:
     needs: [build_stencils]
     runs-on: ubuntu-latest
@@ -52,11 +47,6 @@ jobs:
         name: stencil-object-files
         path: src/copapy/obj
 
-    #- uses: actions/download-artifact@v4
-    #  with:
-    #    name: cross-runner
-    #    path: build/runner
-
     - name: Set up Python ${{ matrix.python-version }}
       uses: actions/setup-python@v5
       with:
@@ -65,46 +55,23 @@ jobs:
     - name: Install Python dependencies
       run: python -m pip install -e .[dev]
 
-    #- name: Install ARM binutils and qemu
-    #  if: strategy.job-index == 0
-    #  run: |
-    #    echo "set man-db/auto-update false" | sudo debconf-communicate
-    #    sudo dpkg-reconfigure man-db
-    #    sudo apt-get update
-    #    sudo apt-get install --no-install-recommends --no-install-suggests binutils-aarch64-linux-gnu qemu-user gcc-aarch64-linux-gnu libc6-dev-arm64-cross
-
     - name: Compile coparun
       run: |
         mkdir -p build/runner
         gcc -O3 -DENABLE_BASIC_LOGGING -o build/runner/coparun src/coparun/runmem.c src/coparun/coparun.c src/coparun/mem_man.c
-     #   aarch64-linux-gnu-gcc -O3 -static -DENABLE_BASIC_LOGGING src/coparun/runmem.c src/coparun/coparun.c src/coparun/mem_man.c -o build/runner/coparun-aarch64
 
     - name: Generate debug asm files
       if: strategy.job-index == 0
       run: |
         set -e
         set -v
-        python tools/make_example.py
+        bash tools/create_asm.sh
 
-        echo "- Patch code..."
+        echo '<p>example</p>' >> $GITHUB_STEP_SUMMARY
-        build/runner/coparun build/runner/test.copapy build/runner/test.copapy.bin
+        python tools/clean_asm.py build/runner/example.asm >> $GITHUB_STEP_SUMMARY
-        #qemu-aarch64 build/runner/coparun-aarch64 build/runner/test-arm64.copapy build/runner/test-arm64.copapy.bin
 
-        objdump -D -b binary -m i386:x86-64 --adjust-vma=0x1000 build/runner/test.copapy.bin > build/runner/test.copapy.asm
+        echo '<p>stencils_x86_64_O3.o</p>' >> $GITHUB_STEP_SUMMARY
-        echo '<p>test.copapy.asm</p>' >> $GITHUB_STEP_SUMMARY
+        python tools/clean_asm.py build/runner/stencils.asm >> $GITHUB_STEP_SUMMARY
-        python tools/clean_asm.py build/runner/test.copapy.asm >> $GITHUB_STEP_SUMMARY
-
-        #aarch64-linux-gnu-objdump -D -b binary -m aarch64 --adjust-vma=0x1000 build/runner/test-arm64.copapy.bin > build/runner/test-arm64.copapy.asm
-        #echo '<p>test-arm64.copapy.asm</p>' >> $GITHUB_STEP_SUMMARY
-        #python tools/clean_asm.py build/runner/test-arm64.copapy.asm >> $GITHUB_STEP_SUMMARY
-
-        objdump -d -x src/copapy/obj/stencils_x86_64_O3.o > build/runner/stencils_x86_64_O3.asm
-        echo '<p>stencils_x86_64_O3.asm</p>' >> $GITHUB_STEP_SUMMARY
-        python tools/clean_asm.py build/runner/stencils_x86_64_O3.asm >> $GITHUB_STEP_SUMMARY
-
-        #aarch64-linux-gnu-objdump -d -x src/copapy/obj/stencils_arm64_O3.o > build/runner/stencils_arm64_O3.asm
-        #echo '<p>stencils_arm64_O3.asm</p>' >> $GITHUB_STEP_SUMMARY
-        #python tools/clean_asm.py build/runner/stencils_arm64_O3.asm >> $GITHUB_STEP_SUMMARY
 
     - name: Run tests with pytest
       run: pytest
@@ -142,7 +109,8 @@ jobs:
                   mkdir -p build/runner && \
                   gcc -O3 -DENABLE_LOGGING -o build/runner/coparun src/coparun/runmem.c \
                   src/coparun/coparun.c src/coparun/mem_man.c && \
-                  pytest"
+                  pytest && \
+                  bash tools/create_asm.sh"
 
     - uses: actions/upload-artifact@v4
       with:
@@ -170,7 +138,8 @@ jobs:
                   mkdir -p build/runner && \
                   gcc -O3 -DENABLE_LOGGING -o build/runner/coparun src/coparun/runmem.c \
                   src/coparun/coparun.c src/coparun/mem_man.c && \
-                  pytest"
+                  pytest && \
+                  bash tools/create_asm.sh"
 
     - uses: actions/upload-artifact@v4
       with:
@@ -219,12 +188,6 @@ jobs:
     - name: Run tests with pytest
       run: pytest
 
-    - name: Type checking with mypy
-      run: mypy
-
-    #- name: Lint code with flake8
-    #  run: flake8
-
     - uses: actions/upload-artifact@v4
       if: strategy.job-index == 0
       with:
@@ -234,7 +197,7 @@ jobs:
   release-stencils:
     needs: [build_stencils, build-ubuntu, build-windows, build-arm64, build-armv7]
     runs-on: ubuntu-latest
-    if: github.event_name == 'push'
+    if: github.ref == 'refs/heads/main' && github.event_name == 'push'
     permissions:
       contents: write
     steps:
@@ -262,7 +225,6 @@ jobs:
           mkdir -p release
           cp tmp/stencil-object-files/* release/
           cp tmp/musl-object-files/* release/
-          cp tmp/cross-runner/coparun-* release/
           cp tmp/runner-linux/coparun release/
           cp tmp/runner-linux-arm64/coparun release/coparun-aarch64
           cp tmp/runner-linux-armv7/coparun release/coparun-armv7

.gitignore

@@ -20,3 +20,6 @@ build/*
 /*.obj
 /src/*.pyd
 vc140.pdb
+benchmark_results*
+docs/build
+docs/source/api

README.md

@@ -1,35 +1,46 @@
 # Copapy
+Copapy is a python framework for deterministic low latency realtime computations, targeting hardware applications - for example in the field of robotics, aerospace, embedded systems and control systems in general.
 
-Copapy is a python based embedded domain specific language (eDSL) with copy & patch compiler. It uses the python interpreter for compilation. It generates a directed graph of variables and operations. The compiler generates machine code by composing pre-compiled stencils derived from compiled C code.
+GPU frameworks like PyTorch, JAX and TensorFlow jump started the development in the field of AI. With the right balance of flexibility and performance they allow for fast iterations of new ideas while being performant enough to test them or even use them in production.
 
-The Project targets applications that profit from fast implementation (e.g. prototyping) and require low latency realtime execution as well as minimizing risk of implementation errors not caught during compile time. This applies primarily to applications interfacing hardware, where runtime errors might lead to physical damage - for example in the field of robotics, aerospace, embedded systems and control systems in general.
+This is exactly what Copapy is aiming for - but in the field of embedded realtime computation. While making use of the ergonomics of Python, the tooling and the general Python ecosystem, Copapy runs seamlessly optimized machine code. Despite being highly portable, the **copy and patch** compiler allows for effortless and fast deployment, without any dependencies beyond Python. It's designed to feel like writing python scripts, with a flat learning curve. But under the hood it produces high performance static typed and memory save code with a minimized set of possible runtime errors[^1]. To maximize productivity the framework provides detailed type hints to catch most errors even before compilation.
 
-The language aims to be:
+Embedded systems comes with a variety of CPU architectures. The **copy and patch** compiler already supports the most common ones [^3] and porting it to new architectures is effortless if a C compiler for the target architecture is available [^2]. The generated code depends only on the CPU architecture. The actual generated code does neither do system calls nor calling external libraries like libc. This allows Copapy for one to be highly deterministic and for the other it makes targeting different realtime operating systems or bare metal straight forward. 
+
+The summarized main features are:
 - Fast to write & easy to read
-- Type safe
+- Memory and type safety, minimal set of runtime errors
-- Having a predictable runtime
+- deterministic execution
-- No runtime errors
+- Auto grad for efficient realtime optimizations
+- Optimized machine code for the target architectures x68_64, Aarch64 and ARMv7
+- Very portable to new architectures
+- Small python package, minimal dependencies, no cross compile toolchain required
 
-Because the language is an embedded language, it can relay heavily on **python tooling**. While copapy is static typed, it uses Python to derive types during compile time wherever possible. It can get full benefit from python type checkers, to catch type errors even before compilation to improve ergonomics.
+## Current state
+While obviously hardware IO is a core aspect, this is not yet available. Therefore this package is at the moment a proof of concept with limited direct use. However the computation part is fully working and available for testing and playing with it by simply installing the package. At this point the project is quite close to being ready for integration into the first demonstration hardware platform.
 
-## How it works
+Currently worked on:
-The **Compilation** step starts with tracing the python code to generate an acyclic directed graph (DAG) of variables and operations. The DAG can be optimized and gets than linearized to a sequence of operations. Each operation gets mapped to a pre-compiled stencil, which is a piece of machine code with placeholders for memory addresses. The compiler generates patch instructions to fill the placeholders with the correct memory addresses. The binary code build from the stencils, data for constants and the patch instructions are than passed to the Runner for execution. The runner allocates memory for the code and data, applies the patch instructions to correct memory addresses and finally executes the code.
+- Array stencils for handling very large arrays and generate SIMD optimized code - e.g. for machine vision and neural network applications.
+- For targeting Crossover‑MCUs, support for Thumb instructions required by ARM*-M is on the way.
+- Constant-regrouping for symbolic optimization of the computation graph.
 
-## Getting started
+## Install
-To install copapy, you can use pip:
+To install copapy, you can use pip. Precompiled wheels are available for Linux (x86_64, Aarch64 and ARMv7), Windows (x86_64) and Mac OS (x86_64, Aarch64):
 
 ```bash
 pip install copapy
 ```
 
+## Examples
+### Basic example
 A very simple example program using copapy can look like this:
 
 ```python
 import copapy as cp
 
 # Define variables
-a = cp.variable(0.25)
+a = cp.value(0.25)
-b = cp.variable(0.87)
+b = cp.value(0.87)
 
 # Define computations
 c = a + b * 2.0
@@ -47,10 +58,63 @@ print("Result d:", tg.read_value(d))
 print("Result e:", tg.read_value(e))
 ```
 
+### Inverse kinematics
+An other example using autograd in copapy. Here for for implementing
+gradient descent to solve a reverse kinematic problem for
+a two joint 2D arm:
+
+```python
+import copapy as cp
+
+# Arm lengths
+l1, l2 = 1.8, 2.0
+
+# Target position
+target = cp.vector([0.7, 0.7])
+
+# Learning rate for iterative adjustment
+alpha = 0.1
+
+def forward_kinematics(theta1, theta2):
+    """Return positions of joint and end-effector."""
+    joint = cp.vector([l1 * cp.cos(theta1), l1 * cp.sin(theta1)])
+    end_effector = joint + cp.vector([l2 * cp.cos(theta1 + theta2),
+                                    l2 * cp.sin(theta1 + theta2)])
+    return joint, end_effector
+
+# Start values
+theta = cp.vector([cp.value(0.0), cp.value(0.0)])
+
+# Iterative inverse kinematic
+for _ in range(48):
+    joint, effector = forward_kinematics(theta[0], theta[1])
+    error = ((target - effector) ** 2).sum()
+
+    theta -= alpha * cp.grad(error, theta)
+
+tg = cp.Target()
+tg.compile(error, theta, joint)
+tg.run()
+
+print(f"Joint angles: {tg.read_value(theta)}")
+print(f"Joint position: {tg.read_value(joint)}")
+print(f"End-effector position: {tg.read_value(effector)}")
+print(f"quadratic error = {tg.read_value(error)}")
+```
+```
+Joint angles: [-0.7221821546554565, 2.6245293617248535]
+Joint position: [1.3509329557418823, -1.189529299736023]       
+End-effector position: [0.6995794177055359, 0.7014330625534058]
+quadratic error = 2.2305819129542215e-06
+```
+
+## How it works
+The **Compilation** step starts with tracing the python code to generate an acyclic directed graph (DAG) of variables and operations. The DAG can be optimized and gets than linearized to a sequence of operations. Each operation gets mapped to a pre-compiled stencil, which is a piece of machine code with placeholders for memory addresses. The compiler generates patch instructions to fill the placeholders with the correct memory addresses. The binary code build from the stencils, data for constants and the patch instructions are than passed to the runner for execution. The runner allocates memory for the code and data, applies the patch instructions to correct memory addresses and finally executes the code.
+
 ## Developer Guide
 Contributions are welcome, please open an issue or submit a pull request on GitHub.
 
-To get started with developing the package, first clone the repository to your local machine using Git:
+To get started with developing the package, first clone the repository using Git:
 
 ```bash
 git clone https://github.com/Nonannet/copapy.git
@@ -70,13 +134,13 @@ Build and install the package and dev dependencies:
 pip install -e .[dev]
 ```
 
-If the build fails because you have no suitable c compiler installed, you can either install a compiler or use the binary from pypi:
+If the build fails because you have no suitable c compiler installed, you can either install a compiler (obviously) or use the binary from pypi:
 
 ```bash
 pip install copapy[dev]
 ```
 
-When running pytest it will use the binary from pypi but the local python code gets executed from the local repo.
+When running pytest it will use the binary part from pypi but all the python code gets executed from the local repo.
 
 For running all tests you need the stencil object files and the compiled runner. You can download the stencils and binary runner from GitHub or build them with gcc yourself.
 
@@ -92,12 +156,6 @@ To build the binaries from source on Linux run:
 bash tools/build.sh
 ```
 
-The runner (without the stencils) can be build on windows with:
-
-```
-tools\build
-```
-
 Ensure that everything is set up correctly by running the tests:
 
 ```bash
@@ -106,3 +164,7 @@ pytest
 
 ## License
 This project is licensed under GPL - see the [LICENSE](LICENSE) file for details.
+
+[^1]: Currently errors like divide by zero are possible. The feasibility of tacking value ranges in the type system is under investigation to be able to do checks at compile time.
+[^2]: The compiler must support TCO (tail call optimization). Currently gcc as C compiler is supported. Porting to a new architecture requires to implement a subset of relocation types used by the architecture.
+[^3]: Supported are x68_64, Aarch64, ARMv7 (non-Thumb); ARMv6/7-M (Thumb) is under development; code for x68 32 Bit is present but has open issues (low priority).

docs/Makefile

@@ -0,0 +1,20 @@
+# Minimal makefile for Sphinx documentation
+#
+
+# You can set these variables from the command line, and also
+# from the environment for the first two.
+SPHINXOPTS    ?=
+SPHINXBUILD   ?= sphinx-build
+SOURCEDIR     = source
+BUILDDIR      = build
+
+# Put it first so that "make" without argument is like "make help".
+help:
+	@$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)
+
+.PHONY: help Makefile
+
+# Catch-all target: route all unknown targets to Sphinx using the new
+# "make mode" option.  $(O) is meant as a shortcut for $(SPHINXOPTS).
+%: Makefile
+	@$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)

docs/make.bat

@@ -0,0 +1,35 @@
+@ECHO OFF
+
+pushd %~dp0
+
+REM Command file for Sphinx documentation
+
+if "%SPHINXBUILD%" == "" (
+	set SPHINXBUILD=sphinx-build
+)
+set SOURCEDIR=source
+set BUILDDIR=build
+
+%SPHINXBUILD% >NUL 2>NUL
+if errorlevel 9009 (
+	echo.
+	echo.The 'sphinx-build' command was not found. Make sure you have Sphinx
+	echo.installed, then set the SPHINXBUILD environment variable to point
+	echo.to the full path of the 'sphinx-build' executable. Alternatively you
+	echo.may add the Sphinx directory to PATH.
+	echo.
+	echo.If you don't have Sphinx installed, grab it from
+	echo.https://www.sphinx-doc.org/
+	exit /b 1
+)
+
+if "%1" == "" goto help
+
+%SPHINXBUILD% -M %1 %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O%
+goto end
+
+:help
+%SPHINXBUILD% -M help %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O%
+
+:end
+popd

docs/source/conf.py

@@ -0,0 +1,37 @@
+# Configuration file for the Sphinx documentation builder.
+#
+# For the full list of built-in configuration values, see the documentation:
+# https://www.sphinx-doc.org/en/master/usage/configuration.html
+
+# -- Project information -----------------------------------------------------
+# https://www.sphinx-doc.org/en/master/usage/configuration.html#project-information
+
+import os
+import sys
+sys.path.insert(0, os.path.abspath("../src/"))
+
+project = 'copapy'
+copyright = '2025, Nicolas Kruse'
+author = 'Nicolas Kruse'
+
+# -- General configuration ---------------------------------------------------
+# https://www.sphinx-doc.org/en/master/usage/configuration.html#general-configuration
+
+extensions = ["sphinx.ext.autodoc", "sphinx.ext.napoleon", "sphinx_autodoc_typehints", "myst_parser"]
+
+templates_path = ['_templates']
+exclude_patterns = []
+
+# -- Options for HTML output -------------------------------------------------
+# https://www.sphinx-doc.org/en/master/usage/configuration.html#options-for-html-output
+
+# html_theme = 'alabaster'
+html_theme = 'pydata_sphinx_theme'
+html_static_path = ['_static']
+html_theme_options = {
+  "secondary_sidebar_items": ["page-toc"],
+  "footer_start": ["copyright"]
+}
+html_theme_options["footer_end"] = []
+
+autodoc_inherit_docstrings = True

docs/source/extract_section.py

@@ -0,0 +1,29 @@
+import re
+
+def extract_sections(md_text: str) -> dict[str, str]:
+    """
+    Extracts sections based on headings (#...).
+    Returns {heading_text: section_content}
+    Works for simple Markdown, not fully strict.
+    """
+
+    # regex captures: heading marks (###...), heading text, and the following content
+    pattern = re.compile(
+        r'^(#{1,6})\s+(.*?)\s*$'          # heading level + heading text
+        r'(.*?)'                          # section content (lazy)
+        r'(?=^#{1,6}\s+|\Z)',             # stop at next heading or end of file
+        re.MULTILINE | re.DOTALL
+    )
+
+    sections: dict[str, str] = {}
+    for _, title, content in pattern.findall(md_text):
+        sections[title] = content.strip()
+
+    return sections
+
+if __name__ == '__main__':
+    with open('README.md', 'rt') as f:
+        readme = extract_sections(f.read())
+
+    with open('docs/source/start.md', 'wt') as f:
+        f.write('\n'.join(readme[s] for s in ['Copapy', 'Current state']))

docs/source/generate_class_list.py

@@ -0,0 +1,86 @@
+# This script generates the source md-files for all classes and functions for the docs
+
+import importlib
+import inspect
+import fnmatch
+from io import TextIOWrapper
+import os
+
+
+def write_manual(f: TextIOWrapper, doc_files: list[str], title: str) -> None:
+    write_dochtree(f, title, doc_files)
+
+
+def write_classes(f: TextIOWrapper, patterns: list[str], module_name: str, title: str, description: str = '', exclude: list[str] = []) -> None:
+    """Write the classes to the file."""
+    module = importlib.import_module(module_name)
+
+    classes = [
+        name for name, obj in inspect.getmembers(module, inspect.isclass)
+        if (any(fnmatch.fnmatch(name, pat) for pat in patterns if pat not in exclude) and
+            obj.__doc__ and '(Automatic generated stub)' not in obj.__doc__)
+    ]
+
+    if description:
+        f.write(f'{description}\n\n')
+
+    write_dochtree(f, title, classes)
+
+    for cls in classes:
+        with open(f'docs/source/api/{cls}.md', 'w') as f2:
+            f2.write(f'# {module_name}.{cls}\n')
+            f2.write('```{eval-rst}\n')
+            f2.write(f'.. autoclass:: {module_name}.{cls}\n')
+            f2.write('   :members:\n')
+            f2.write('   :undoc-members:\n')
+            f2.write('   :show-inheritance:\n')
+            f2.write('   :inherited-members:\n')
+            f2.write('```\n\n')
+
+
+def write_functions(f: TextIOWrapper, patterns: list[str], module_name: str, title: str, description: str = '', exclude: list[str] = []) -> None:
+    """Write the classes to the file."""
+    module = importlib.import_module(module_name)
+
+    functions = [
+        name for name, _ in inspect.getmembers(module, inspect.isfunction)
+        if (any(fnmatch.fnmatch(name, pat) for pat in patterns if pat not in exclude))
+    ]
+
+    if description:
+        f.write(f'{description}\n\n')
+
+    write_dochtree(f, title, functions)
+
+    for func in functions:
+        if not func.startswith('_'):
+            with open(f'docs/source/api/{func}.md', 'w') as f2:
+                f2.write(f'# {module_name}.{func}\n')
+                f2.write('```{eval-rst}\n')
+                f2.write(f'.. autofunction:: {module_name}.{func}\n')
+                f2.write('```\n\n')
+
+
+def write_dochtree(f: TextIOWrapper, title: str, items: list[str]):
+    f.write('```{toctree}\n')
+    f.write(':maxdepth: 1\n')
+    f.write(f':caption: {title}:\n')
+    #f.write(':hidden:\n')
+    for text in items:
+        if not text.startswith('_'):
+            f.write(f"{text}\n")
+    f.write('```\n\n')
+
+
+if __name__ == "__main__":
+    # Ensure the output directory exists
+    os.makedirs('docs/source/api', exist_ok=True)
+
+    with open('docs/source/api/index.md', 'w') as f:
+        f.write('# Classes and functions\n\n')
+
+        write_classes(f, ['*'], 'copapy', title='Classes')
+
+        write_functions(f, ['*'], 'copapy', title='Functions')
+
+        #write_manual(f, ['../ndfloat', '../floatarray'], title='Types')

docs/source/index.md

@@ -0,0 +1,9 @@
+```{toctree}
+:maxdepth: 1
+:hidden:
+api/index
+repo
+```
+
+```{include} ../../README.md
+```

docs/source/repo.md

@@ -0,0 +1,3 @@
+# Code repository
+
+Code repository is on GitHub: [github.com/Nonannet/copapy](https://github.com/Nonannet/copapy).

pyproject.toml

@@ -36,6 +36,12 @@ dev = [
     "mypy",
     "pytest"
 ]
+doc_build = [
+  "sphinx",
+  "pydata_sphinx_theme",
+  "sphinx-autodoc-typehints",
+  "myst-parser"
+]
 
 [tool.mypy]
 files = ["src", "tools", "stencils"]

src/copapy/__init__.py

@@ -1,18 +1,26 @@
 from ._target import Target
-from ._basic_types import NumLike, variable, generic_sdb, iif
+from ._basic_types import NumLike, value, generic_sdb, iif
 from ._vectors import vector, distance, scalar_projection, angle_between, rotate_vector, vector_projection
-from ._math import sqrt, abs, sin, cos, tan, asin, acos, atan, atan2, log, exp, pow, get_42, clamp, min, max
+from ._matrices import matrix, identity, zeros, ones, diagonal, eye
+from ._math import sqrt, abs, sign, sin, cos, tan, asin, acos, atan, atan2, log, exp, pow, get_42, clamp, min, max, relu
+from ._autograd import grad
 
 __all__ = [
     "Target",
     "NumLike",
-    "variable",
+    "value",
     "generic_sdb",
     "iif",
     "vector",
+    "matrix",
+    "identity",
+    "zeros",
+    "ones",
+    "diagonal",
     "sqrt",
     "abs",
     "sin",
+    "sign",
     "cos",
     "tan",
     "asin",
@@ -26,9 +34,12 @@ __all__ = [
     "clamp",
     "min",
     "max",
+    "relu",
     "distance",
     "scalar_projection",
     "angle_between",
     "rotate_vector",
     "vector_projection",
+    "grad",
+    "eye"
 ]

src/copapy/_autograd.py

@@ -0,0 +1,126 @@
+from . import value, vector, matrix
+import copapy.backend as cpb
+from typing import Any, Sequence, overload
+import copapy as cp
+from ._basic_types import Net, unifloat
+
+
+@overload
+def grad(x: Any, y: value[Any]) -> unifloat: ...
+@overload
+def grad(x: Any, y: vector[Any]) -> vector[float]: ...
+@overload
+def grad(x: Any, y: Sequence[value[Any]]) -> list[unifloat]: ...
+@overload
+def grad(x: Any, y: matrix[Any]) -> matrix[float]: ...
+def grad(x: Any, y: value[Any] | Sequence[value[Any]] | vector[Any] | matrix[Any]) -> Any:
+    """Returns the partial derivative dx/dy where x needs to be a scalar
+    and y might be a scalar, a list of scalars, a vector or matrix.
+
+    Arguments:
+        x: Value to return derivative of
+        y: Value(s) to derive in respect to
+
+    Returns:
+        Derivative of x with the type and dimensions of y
+    """
+    assert isinstance(x, value), f"Argument x for grad function must be a copapy value but is {type(x)}."
+
+    if isinstance(y, value):
+        y_set = {y}
+    if isinstance(y, matrix):
+        y_set = {v for row in y for v in row}
+    else:
+        assert isinstance(y, Sequence) or isinstance(y, vector)
+        y_set = {v for v in y}
+
+    edges = cpb.get_all_dag_edges_between([x.source], (net.source for net in y_set if isinstance(net, Net)))
+    ordered_ops = cpb.stable_toposort(edges)
+
+    net_lookup = {net.source: net for node in ordered_ops for net in node.args}
+    grad_dict: dict[Net, unifloat] = dict()
+
+    def add_grad(val: value[Any], gradient_value: unifloat) -> None:
+        grad_dict[val] = grad_dict.get(val, 0.0) + gradient_value
+
+    for node in reversed(ordered_ops):
+        #print(f"-->   {'x' if node in net_lookup else ' '}", node, f"{net_lookup.get(node)}")
+        if node.args:
+            args: Sequence[Any] = list(node.args)
+            g = 1.0 if node is x.source else grad_dict[net_lookup[node]]
+            opn = node.name.split('_')[0]
+            a: value[Any] = args[0]
+            b: value[Any] = args[1] if len(args) > 1 else a
+
+            if opn in ['ge', 'gt', 'eq', 'ne', 'floordiv', 'bwand', 'bwor', 'bwxor']:
+                pass  # Derivative is 0 for all ops returning integers
+
+            elif opn == 'add':
+                add_grad(a, g)
+                add_grad(b, g)
+
+            elif opn == 'sub':
+                add_grad(a, g)
+                add_grad(b, -g)
+
+            elif opn == 'mul':
+                add_grad(a, b * g)
+                add_grad(b, a * g)
+
+            elif opn == 'div':
+                add_grad(a, g / b)
+                add_grad(b, -a * g / (b**2))
+
+            elif opn == 'mod':
+                add_grad(a, g)
+                add_grad(b, -a * g / b)
+
+            elif opn == 'log':
+                add_grad(a, g / a)
+
+            elif opn == 'exp':
+                add_grad(a, g * cp.exp(a))
+
+            elif opn == 'pow':
+                add_grad(a, (b * (a ** (b - 1))) * g)
+                add_grad(b, (a ** b * cp.log(a)) * g)
+
+            elif opn == 'sqrt':
+                add_grad(a, g * (0.5 / cp.sqrt(a)))
+
+            #elif opn == 'abs':
+            #    add_grad(x, g * cp.sign(x))
+
+            elif opn == 'sin':
+                add_grad(a, g * cp.cos(a))
+
+            elif opn == 'cos':
+                add_grad(a, g * -cp.sin(a))
+
+            elif opn == 'tan':
+                add_grad(a, g * (1 / cp.cos(a) ** 2))
+
+            elif opn == 'asin':
+                add_grad(a, g * (1 / cp.sqrt(1 - a**2)))
+
+            elif opn == 'acos':
+                add_grad(a, g * (-1 / cp.sqrt(1 - a**2)))
+
+            elif opn == 'atan':
+                add_grad(a, g * (1 / (1 + a**2)))
+
+            elif opn == 'atan2':
+                denom = a**2 + b**2
+                add_grad(a, g * (-b / denom))
+                add_grad(b, g * ( a / denom))
+
+            else:
+                raise ValueError(f"Operation {opn} not yet supported for auto diff.")
+
+    if isinstance(y, value):
+        return grad_dict[y]
+    if isinstance(y, vector):
+        return vector(grad_dict[yi] if isinstance(yi, value) else 0.0 for yi in y)
+    if isinstance(y, matrix):
+        return matrix((grad_dict[yi] if isinstance(yi, value) else 0.0 for yi in row) for row in y)
+    return [grad_dict[yi] for yi in y]

src/copapy/_basic_types.py

@@ -1,17 +1,15 @@
 import pkgutil
-from typing import Any, TypeVar, overload, TypeAlias, Generic, cast
+from typing import Any, Sequence, TypeVar, overload, TypeAlias, Generic, cast
 from ._stencils import stencil_database, detect_process_arch
 import copapy as cp
+from ._helper_types import TNum
 
-NumLike: TypeAlias = 'variable[int] | variable[float] | variable[bool] | int | float | bool'
+NumLike: TypeAlias = 'value[int] | value[float] | int | float'
-unifloat: TypeAlias = 'variable[float] | float'
+unifloat: TypeAlias = 'value[float] | float'
-uniint: TypeAlias = 'variable[int] | int'
+uniint: TypeAlias = 'value[int] | int'
-unibool: TypeAlias = 'variable[bool] | bool'
-uniboolint: TypeAlias = 'variable[bool] | bool | variable[int] | int'
 
-TCPNum = TypeVar("TCPNum", bound='variable[Any]')
+TCPNum = TypeVar("TCPNum", bound='value[Any]')
-TNum = TypeVar("TNum", int, float, bool)
+TVarNumb: TypeAlias = 'value[Any] | int | float'
-TVarNumb: TypeAlias = 'variable[Any] | int | float | bool'
 
 stencil_cache: dict[tuple[str, str], stencil_database] = {}
 
@@ -51,15 +49,25 @@ class Node:
         name (str): The name of the operation this Node represents.
     """
     def __init__(self) -> None:
-        self.args: list[Net] = []
+        self.args: tuple[Net, ...] = tuple()
         self.name: str = ''
+        self.node_hash = 0
 
     def __repr__(self) -> str:
         return f"Node:{self.name}({', '.join(str(a) for a in self.args) if self.args else (self.value if isinstance(self, CPConstant) else '')})"
 
+    def get_node_hash(self, commutative: bool = False) -> int:
+        if commutative:
+            return hash(self.name) ^ hash(frozenset(a.source.node_hash for a in self.args))
+        return hash(self.name) ^ hash(tuple(a.source.node_hash for a in self.args))
+
+
+    def __hash__(self) -> int:
+        return self.node_hash
+
 
 class Net:
-    """A Net represents a variable in the computation graph - or more generally it
+    """A Net represents a scalar type in the computation graph - or more generally it
     connects Nodes together.
 
     Attributes:
@@ -72,25 +80,25 @@ class Net:
 
     def __repr__(self) -> str:
         names = get_var_name(self)
-        return f"{'name:' + names[0] if names else 'id:' + str(id(self))[-5:]}"
+        return f"{'name:' + names[0] if names else 'id:' + str(hash(self))[-5:]}"
 
     def __hash__(self) -> int:
-        return id(self)
+        return self.source.node_hash
 
 
-class variable(Generic[TNum], Net):
+class value(Generic[TNum], Net):
-    """A "variable" represents a typed variable. It supports arithmetic and
+    """A "value" represents a typed scalar variable. It supports arithmetic and
     comparison operations.
 
     Attributes:
-        dtype (str): Data type of this variable.
+        dtype (str): Data type of this value.
     """
     def __init__(self, source: TNum | Node, dtype: str | None = None):
-        """Instance a variable.
+        """Instance a value.
 
         Args:
             source: A numeric value or Node object.
-            dtype: Data type of this variable. Required if source is a Node.
+            dtype: Data type of this value. Required if source is a Node.
         """
         if isinstance(source, Node):
             self.source = source
@@ -107,147 +115,180 @@ class variable(Generic[TNum], Net):
             self.dtype = 'int'
 
     @overload
-    def __add__(self, other: TCPNum) -> TCPNum: ...
+    def __add__(self: 'value[TNum]', other: 'value[TNum] | TNum') -> 'value[TNum]': ...
     @overload
-    def __add__(self: TCPNum, other: uniint) -> TCPNum: ...
+    def __add__(self: 'value[int]', other: uniint) -> 'value[int]': ...
     @overload
-    def __add__(self, other: unifloat) -> 'variable[float]': ...
+    def __add__(self, other: unifloat) -> 'value[float]': ...
     @overload
-    def __add__(self, other: NumLike) -> 'variable[float] | variable[int]': ...
+    def __add__(self: 'value[float]', other: NumLike) -> 'value[float]': ...
-    def __add__(self, other: NumLike) -> Any:
+    @overload
-        if isinstance(other, int | float) and other == 0:
+    def __add__(self, other: TVarNumb) -> 'value[float] | value[int]': ...
+    def __add__(self, other: TVarNumb) -> Any:
+        if not isinstance(other, value) and other == 0:
             return self
         return add_op('add', [self, other], True)
 
     @overload
-    def __radd__(self: TCPNum, other: int) -> TCPNum: ...
+    def __radd__(self: 'value[TNum]', other: TNum) -> 'value[TNum]': ...
     @overload
-    def __radd__(self, other: float) -> 'variable[float]': ...
+    def __radd__(self: 'value[int]', other: int) -> 'value[int]': ...
+    @overload
+    def __radd__(self, other: float) -> 'value[float]': ...
     def __radd__(self, other: NumLike) -> Any:
-        if isinstance(other, int | float) and other == 0:
+        return self + other
-            return self
-        return add_op('add', [self, other], True)
 
     @overload
-    def __sub__(self, other: TCPNum) -> TCPNum: ...
+    def __sub__(self: 'value[TNum]', other: 'value[TNum] | TNum') -> 'value[TNum]': ...
     @overload
-    def __sub__(self: TCPNum, other: uniint) -> TCPNum: ...
+    def __sub__(self: 'value[int]', other: uniint) -> 'value[int]': ...
     @overload
-    def __sub__(self, other: unifloat) -> 'variable[float]': ...
+    def __sub__(self, other: unifloat) -> 'value[float]': ...
     @overload
-    def __sub__(self, other: NumLike) -> 'variable[float] | variable[int]': ...
+    def __sub__(self: 'value[float]', other: NumLike) -> 'value[float]': ...
-    def __sub__(self, other: NumLike) -> Any:
+    @overload
+    def __sub__(self, other: TVarNumb) -> 'value[float] | value[int]': ...
+    def __sub__(self, other: TVarNumb) -> Any:
+        if isinstance(other, int | float) and other == 0:
+            return self
         return add_op('sub', [self, other])
 
     @overload
-    def __rsub__(self: TCPNum, other: int) -> TCPNum: ...
+    def __rsub__(self: 'value[TNum]', other: TNum) -> 'value[TNum]': ...
     @overload
-    def __rsub__(self, other: float) -> 'variable[float]': ...
+    def __rsub__(self: 'value[int]', other: int) -> 'value[int]': ...
+    @overload
+    def __rsub__(self, other: float) -> 'value[float]': ...
     def __rsub__(self, other: NumLike) -> Any:
         return add_op('sub', [other, self])
 
     @overload
-    def __mul__(self, other: TCPNum) -> TCPNum: ...
+    def __mul__(self: 'value[TNum]', other: 'value[TNum] | TNum') -> 'value[TNum]': ...
     @overload
-    def __mul__(self: TCPNum, other: uniint) -> TCPNum: ...
+    def __mul__(self: 'value[int]', other: uniint) -> 'value[int]': ...
     @overload
-    def __mul__(self, other: unifloat) -> 'variable[float]': ...
+    def __mul__(self, other: unifloat) -> 'value[float]': ...
     @overload
-    def __mul__(self, other: NumLike) -> 'variable[float] | variable[int]': ...
+    def __mul__(self: 'value[float]', other: NumLike) -> 'value[float]': ...
-    def __mul__(self, other: NumLike) -> Any:
+    @overload
+    def __mul__(self, other: TVarNumb) -> 'value[float] | value[int]': ...
+    def __mul__(self, other: TVarNumb) -> Any:
+        if self.dtype == 'float' and isinstance(other, int):
+            other = float(other)  # Prevent runtime conversion of consts; TODO: add this for other operations
+        if not isinstance(other, value):
+            if other == 1:
+                return self
+            elif other == 0:
+                return 0
         return add_op('mul', [self, other], True)
 
     @overload
-    def __rmul__(self: TCPNum, other: int) -> TCPNum: ...
+    def __rmul__(self: 'value[TNum]', other: TNum) -> 'value[TNum]': ...
     @overload
-    def __rmul__(self, other: float) -> 'variable[float]': ...
+    def __rmul__(self: 'value[int]', other: int) -> 'value[int]': ...
+    @overload
+    def __rmul__(self, other: float) -> 'value[float]': ...
     def __rmul__(self, other: NumLike) -> Any:
-        return add_op('mul', [self, other], True)
+        return self * other
 
-    def __truediv__(self, other: NumLike) -> 'variable[float]':
+    def __truediv__(self, other: NumLike) -> 'value[float]':
         return add_op('div', [self, other])
 
-    def __rtruediv__(self, other: NumLike) -> 'variable[float]':
+    def __rtruediv__(self, other: NumLike) -> 'value[float]':
         return add_op('div', [other, self])
 
     @overload
-    def __floordiv__(self, other: TCPNum) -> TCPNum: ...
+    def __floordiv__(self: 'value[TNum]', other: 'value[TNum] | TNum') -> 'value[TNum]': ...
     @overload
-    def __floordiv__(self: TCPNum, other: uniint) -> TCPNum: ...
+    def __floordiv__(self: 'value[int]', other: uniint) -> 'value[int]': ...
     @overload
-    def __floordiv__(self, other: unifloat) -> 'variable[float]': ...
+    def __floordiv__(self, other: unifloat) -> 'value[float]': ...
     @overload
-    def __floordiv__(self, other: NumLike) -> 'variable[float] | variable[int]': ...
+    def __floordiv__(self: 'value[float]', other: NumLike) -> 'value[float]': ...
-    def __floordiv__(self, other: NumLike) -> Any:
+    @overload
+    def __floordiv__(self, other: TVarNumb) -> 'value[float] | value[int]': ...
+    def __floordiv__(self, other: TVarNumb) -> Any:
         return add_op('floordiv', [self, other])
 
     @overload
-    def __rfloordiv__(self: TCPNum, other: int) -> TCPNum: ...
+    def __rfloordiv__(self: 'value[TNum]', other: TNum) -> 'value[TNum]': ...
     @overload
-    def __rfloordiv__(self, other: float) -> 'variable[float]': ...
+    def __rfloordiv__(self: 'value[int]', other: int) -> 'value[int]': ...
+    @overload
+    def __rfloordiv__(self, other: float) -> 'value[float]': ...
     def __rfloordiv__(self, other: NumLike) -> Any:
         return add_op('floordiv', [other, self])
 
     def __neg__(self: TCPNum) -> TCPNum:
-        return cast(TCPNum, add_op('sub', [variable(0), self]))
+        if self.dtype == 'int':
+            return cast(TCPNum, add_op('sub', [value(0), self]))
+        return cast(TCPNum, add_op('sub', [value(0.0), self]))
 
-    def __gt__(self, other: TVarNumb) -> 'variable[bool]':
+    def __gt__(self, other: TVarNumb) -> 'value[int]':
         ret = add_op('gt', [self, other])
-        return variable(ret.source, dtype='bool')
+        return value(ret.source, dtype='bool')
 
-    def __lt__(self, other: TVarNumb) -> 'variable[bool]':
+    def __lt__(self, other: TVarNumb) -> 'value[int]':
         ret = add_op('gt', [other, self])
-        return variable(ret.source, dtype='bool')
+        return value(ret.source, dtype='bool')
 
-    def __ge__(self, other: TVarNumb) -> 'variable[bool]':
+    def __ge__(self, other: TVarNumb) -> 'value[int]':
         ret = add_op('ge', [self, other])
-        return variable(ret.source, dtype='bool')
+        return value(ret.source, dtype='bool')
 
-    def __le__(self, other: TVarNumb) -> 'variable[bool]':
+    def __le__(self, other: TVarNumb) -> 'value[int]':
         ret = add_op('ge', [other, self])
-        return variable(ret.source, dtype='bool')
+        return value(ret.source, dtype='bool')
 
-    def __eq__(self, other: TVarNumb) -> 'variable[bool]':  # type: ignore
+    def __eq__(self, other: TVarNumb) -> 'value[int]':  # type: ignore
         ret = add_op('eq', [self, other], True)
-        return variable(ret.source, dtype='bool')
+        return value(ret.source, dtype='bool')
 
-    def __ne__(self, other: TVarNumb) -> 'variable[bool]':  # type: ignore
+    def __ne__(self, other: TVarNumb) -> 'value[int]':  # type: ignore
         ret = add_op('ne', [self, other], True)
-        return variable(ret.source, dtype='bool')
+        return value(ret.source, dtype='bool')
 
     @overload
-    def __mod__(self, other: TCPNum) -> TCPNum: ...
+    def __mod__(self: 'value[TNum]', other: 'value[TNum] | TNum') -> 'value[TNum]': ...
     @overload
-    def __mod__(self: TCPNum, other: uniint) -> TCPNum: ...
+    def __mod__(self: 'value[int]', other: uniint) -> 'value[int]': ...
     @overload
-    def __mod__(self, other: unifloat) -> 'variable[float]': ...
+    def __mod__(self, other: unifloat) -> 'value[float]': ...
     @overload
-    def __mod__(self, other: NumLike) -> 'variable[float] | variable[int]': ...
+    def __mod__(self: 'value[float]', other: NumLike) -> 'value[float]': ...
-    def __mod__(self, other: NumLike) -> Any:
+    @overload
+    def __mod__(self, other: TVarNumb) -> 'value[float] | value[int]': ...
+    def __mod__(self, other: TVarNumb) -> Any:
         return add_op('mod', [self, other])
 
     @overload
-    def __rmod__(self: TCPNum, other: int) -> TCPNum: ...
+    def __rmod__(self: 'value[TNum]', other: TNum) -> 'value[TNum]': ...
     @overload
-    def __rmod__(self, other: float) -> 'variable[float]': ...
+    def __rmod__(self: 'value[int]', other: int) -> 'value[int]': ...
+    @overload
+    def __rmod__(self, other: float) -> 'value[float]': ...
     def __rmod__(self, other: NumLike) -> Any:
         return add_op('mod', [other, self])
 
     @overload
-    def __pow__(self, other: TCPNum) -> TCPNum: ...
+    def __pow__(self: 'value[TNum]', other: 'value[TNum] | TNum') -> 'value[TNum]': ...
     @overload
-    def __pow__(self: TCPNum, other: uniint) -> TCPNum: ...
+    def __pow__(self: 'value[int]', other: uniint) -> 'value[int]': ...
     @overload
-    def __pow__(self, other: unifloat) -> 'variable[float]': ...
+    def __pow__(self, other: unifloat) -> 'value[float]': ...
     @overload
-    def __pow__(self, other: NumLike) -> 'variable[float] | variable[int]': ...
+    def __pow__(self: 'value[float]', other: NumLike) -> 'value[float]': ...
-    def __pow__(self, other: NumLike) -> Any:
+    @overload
+    def __pow__(self, other: TVarNumb) -> 'value[float] | value[int]': ...
+    def __pow__(self, other: TVarNumb) -> Any:
         return cp.pow(self, other)
 
     @overload
-    def __rpow__(self: TCPNum, other: int) -> TCPNum: ...
+    def __rpow__(self: 'value[TNum]', other: TNum) -> 'value[TNum]': ...
     @overload
-    def __rpow__(self, other: float) -> 'variable[float]': ...
+    def __rpow__(self: 'value[int]', other: int) -> 'value[int]': ...
+    @overload
+    def __rpow__(self, other: float) -> 'value[float]': ...
     def __rpow__(self, other: NumLike) -> Any:
         return cp.pow(other, self)
 
@@ -255,34 +296,34 @@ class variable(Generic[TNum], Net):
         return super().__hash__()
 
     # Bitwise and shift operations for cp[int]
-    def __lshift__(self, other: uniboolint) -> 'variable[int]':
+    def __lshift__(self, other: uniint) -> 'value[int]':
         return add_op('lshift', [self, other])
 
-    def __rlshift__(self, other: uniboolint) -> 'variable[int]':
+    def __rlshift__(self, other: uniint) -> 'value[int]':
         return add_op('lshift', [other, self])
 
-    def __rshift__(self, other: uniboolint) -> 'variable[int]':
+    def __rshift__(self, other: uniint) -> 'value[int]':
         return add_op('rshift', [self, other])
 
-    def __rrshift__(self, other: uniboolint) -> 'variable[int]':
+    def __rrshift__(self, other: uniint) -> 'value[int]':
         return add_op('rshift', [other, self])
 
-    def __and__(self, other: uniboolint) -> 'variable[int]':
+    def __and__(self, other: uniint) -> 'value[int]':
         return add_op('bwand', [self, other], True)
 
-    def __rand__(self, other: uniboolint) -> 'variable[int]':
+    def __rand__(self, other: uniint) -> 'value[int]':
-        return add_op('rwand', [other, self], True)
+        return add_op('bwand', [other, self], True)
 
-    def __or__(self, other: uniboolint) -> 'variable[int]':
+    def __or__(self, other: uniint) -> 'value[int]':
         return add_op('bwor', [self, other], True)
 
-    def __ror__(self, other: uniboolint) -> 'variable[int]':
+    def __ror__(self, other: uniint) -> 'value[int]':
         return add_op('bwor', [other, self], True)
 
-    def __xor__(self, other: uniboolint) -> 'variable[int]':
+    def __xor__(self, other: uniint) -> 'value[int]':
         return add_op('bwxor', [self, other], True)
 
-    def __rxor__(self, other: uniboolint) -> 'variable[int]':
+    def __rxor__(self, other: uniint) -> 'value[int]':
         return add_op('bwxor', [other, self], True)
 
 
@@ -290,7 +331,8 @@ class CPConstant(Node):
     def __init__(self, value: int | float):
         self.dtype, self.value = _get_data_and_dtype(value)
         self.name = 'const_' + self.dtype
-        self.args = []
+        self.args = tuple()
+        self.node_hash = id(self)
 
 
 class Write(Node):
@@ -302,40 +344,42 @@ class Write(Node):
             net = Net(node.dtype, node)
 
         self.name = 'write_' + transl_type(net.dtype)
-        self.args = [net]
+        self.args = (net,)
+        self.node_hash = hash(self.name) ^ hash(net.source.node_hash)
 
 
 class Op(Node):
-    def __init__(self, typed_op_name: str, args: list[Net]):
+    def __init__(self, typed_op_name: str, args: Sequence[Net], commutative: bool = False):
         assert not args or any(isinstance(t, Net) for t in args), 'args parameter must be of type list[Net]'
         self.name: str = typed_op_name
-        self.args: list[Net] = args
+        self.args: tuple[Net, ...] = tuple(args)
+        self.node_hash = self.get_node_hash(commutative)
 
 
-def net_from_value(value: Any) -> Net:
+def net_from_value(val: Any) -> value[Any]:
-    vi = CPConstant(value)
+    vi = CPConstant(val)
-    return Net(vi.dtype, vi)
+    return value(vi, vi.dtype)
 
 
 @overload
-def iif(expression: variable[Any], true_result: unibool, false_result: unibool) -> variable[bool]:  ... # pyright: ignore[reportOverlappingOverload]
+def iif(expression: value[Any], true_result: uniint, false_result: uniint) -> value[int]: ...  # pyright: ignore[reportOverlappingOverload]
 @overload
-def iif(expression: variable[Any], true_result: uniint, false_result: uniint) -> variable[int]: ...
+def iif(expression: value[Any], true_result: unifloat, false_result: unifloat) -> value[float]: ...
-@overload
-def iif(expression: variable[Any], true_result: unifloat, false_result: unifloat) -> variable[float]: ...
 @overload
 def iif(expression: float | int, true_result: TNum, false_result: TNum) -> TNum: ...
 @overload
-def iif(expression: float | int, true_result: TNum, false_result: variable[TNum]) -> variable[TNum]: ...
+def iif(expression: float | int, true_result: TNum | value[TNum], false_result: value[TNum]) -> value[TNum]: ...
 @overload
-def iif(expression: float | int, true_result: variable[TNum], false_result: TNum | variable[TNum]) -> variable[TNum]: ...
+def iif(expression: float | int, true_result: value[TNum], false_result: TNum | value[TNum]) -> value[TNum]: ...
+@overload
+def iif(expression: float | int | value[Any], true_result: TNum | value[TNum], false_result: TNum | value[TNum]) -> value[TNum] | TNum: ...
 def iif(expression: Any, true_result: Any, false_result: Any) -> Any:
-    allowed_type = (variable, int, float, bool)
+    allowed_type = (value, int, float)
     assert isinstance(true_result, allowed_type) and isinstance(false_result, allowed_type), "Result type not supported"
     return (expression != 0) * true_result + (expression == 0) * false_result
 
 
-def add_op(op: str, args: list[variable[Any] | int | float], commutative: bool = False) -> variable[Any]:
+def add_op(op: str, args: list[value[Any] | int | float], commutative: bool = False) -> value[Any]:
     arg_nets = [a if isinstance(a, Net) else net_from_value(a) for a in args]
 
     if commutative:
@@ -349,15 +393,13 @@ def add_op(op: str, args: list[variable[Any] | int | float], commutative: bool =
     result_type = generic_sdb.stencil_definitions[typed_op].split('_')[0]
 
     if result_type == 'float':
-        return variable[float](Op(typed_op, arg_nets), result_type)
+        return value[float](Op(typed_op, arg_nets, commutative), result_type)
     else:
-        return variable[int](Op(typed_op, arg_nets), result_type)
+        return value[int](Op(typed_op, arg_nets, commutative), result_type)
 
 
 def _get_data_and_dtype(value: Any) -> tuple[str, float | int]:
-    if isinstance(value, bool):
+    if isinstance(value, int):
-        return ('bool', int(value))
-    elif isinstance(value, int):
         return ('int', int(value))
     elif isinstance(value, float):
         return ('float', float(value))

src/copapy/_compiler.py

@@ -55,6 +55,43 @@ def stable_toposort(edges: Iterable[tuple[Node, Node]]) -> list[Node]:
     return result
 
 
+def get_all_dag_edges_between(roots: Iterable[Node], leaves: Iterable[Node]) -> Generator[tuple[Node, Node], None, None]:
+    """Get all edges in the DAG connecting given roots with given leaves
+
+    Arguments:
+        nodes: Iterable of nodes to start the traversal from
+
+    Yields:
+        Tuples of (source_node, target_node) representing edges in the DAG
+    """
+    # Walk the full DAG starting from given roots to final leaves
+    parent_lookup: dict[Node, set[Node]] = dict()
+    node_list: list[Node] = [n for n in roots]
+    while(node_list):
+        node = node_list.pop()
+        for net in node.args:
+            if net.source in parent_lookup:
+                parent_lookup[net.source].add(node)
+            else:
+                parent_lookup[net.source] = {node}
+                node_list.append(net.source)
+
+    # Walk the DAG in reverse direction starting from given leaves to given roots
+    emitted_edges: set[tuple[Node, Node]] = set()
+    node_list = [n for n in leaves]
+    while(node_list):
+        child_node = node_list.pop()
+        if child_node in parent_lookup:
+            for node in parent_lookup[child_node]:
+                edge = (child_node, node)
+                if edge not in emitted_edges:
+                    yield edge
+                    node_list.append(node)
+                    emitted_edges.add(edge)
+
+    assert all(r in {e[0] for e in emitted_edges} for r in leaves)
+
+
 def get_all_dag_edges(nodes: Iterable[Node]) -> Generator[tuple[Node, Node], None, None]:
     """Get all edges in the DAG by traversing from the given nodes
 
@@ -64,9 +101,17 @@ def get_all_dag_edges(nodes: Iterable[Node]) -> Generator[tuple[Node, Node], Non
     Yields:
         Tuples of (source_node, target_node) representing edges in the DAG
     """
-    for node in nodes:
+    emitted_edges: set[tuple[Node, Node]] = set()
-        yield from get_all_dag_edges(net.source for net in node.args)
+    node_list: list[Node] = [n for n in nodes]
-        yield from ((net.source, node) for net in node.args)
+
+    while(node_list):
+        node = node_list.pop()
+        for net in node.args:
+            edge = (net.source, node)
+            if edge not in emitted_edges:
+                yield edge
+                node_list.append(net.source)
+                emitted_edges.add(edge)
 
 
 def get_const_nets(nodes: list[Node]) -> list[Net]:

src/copapy/_helper_types.py

@@ -0,0 +1,4 @@
+from typing import TypeVar
+
+TNum = TypeVar("TNum", int, float)
+U = TypeVar("U", int, float)

src/copapy/_math.py

@@ -1,17 +1,18 @@
 from . import vector
 from ._vectors import VecNumLike
-from . import variable, NumLike
+from . import value, NumLike
 from typing import TypeVar, Any, overload, Callable
-from ._basic_types import add_op
+from ._basic_types import add_op, unifloat
 import math
 
-T = TypeVar("T", int, float, variable[int], variable[float])
+T = TypeVar("T", int, float, value[int], value[float])
 U = TypeVar("U", int, float)
 
+
 @overload
 def exp(x: float | int) -> float: ...
 @overload
-def exp(x: variable[Any]) -> variable[float]: ...
+def exp(x: value[Any]) -> value[float]: ...
 @overload
 def exp(x: vector[Any]) -> vector[float]: ...
 def exp(x: Any) -> Any:
@@ -23,7 +24,7 @@ def exp(x: Any) -> Any:
     Returns:
         result of e**x
     """
-    if isinstance(x, variable):
+    if isinstance(x, value):
         return add_op('exp', [x])
     if isinstance(x, vector):
         return x.map(exp)
@@ -33,7 +34,7 @@ def exp(x: Any) -> Any:
 @overload
 def log(x: float | int) -> float: ...
 @overload
-def log(x: variable[Any]) -> variable[float]: ...
+def log(x: value[Any]) -> value[float]: ...
 @overload
 def log(x: vector[Any]) -> vector[float]: ...
 def log(x: Any) -> Any:
@@ -45,7 +46,7 @@ def log(x: Any) -> Any:
     Returns:
         result of ln(x)
     """
-    if isinstance(x, variable):
+    if isinstance(x, value):
         return add_op('log', [x])
     if isinstance(x, vector):
         return x.map(log)
@@ -55,9 +56,9 @@ def log(x: Any) -> Any:
 @overload
 def pow(x: float | int, y: float | int) -> float: ...
 @overload
-def pow(x: variable[Any], y: NumLike) -> variable[float]: ...
+def pow(x: value[Any], y: NumLike) -> value[float]: ...
 @overload
-def pow(x: NumLike, y: variable[Any]) -> variable[float]: ...
+def pow(x: NumLike, y: value[Any]) -> value[float]: ...
 @overload
 def pow(x: vector[Any], y: Any) -> vector[float]: ...
 def pow(x: VecNumLike, y: VecNumLike) -> Any:
@@ -80,7 +81,7 @@ def pow(x: VecNumLike, y: VecNumLike) -> Any:
         return m
     if y == -1:
         return 1 / x
-    if isinstance(x, variable) or isinstance(y, variable):
+    if isinstance(x, value) or isinstance(y, value):
         return add_op('pow', [x, y])
     else:
         return float(x ** y)
@@ -89,7 +90,7 @@ def pow(x: VecNumLike, y: VecNumLike) -> Any:
 @overload
 def sqrt(x: float | int) -> float: ...
 @overload
-def sqrt(x: variable[Any]) -> variable[float]: ...
+def sqrt(x: value[Any]) -> value[float]: ...
 @overload
 def sqrt(x: vector[Any]) -> vector[float]: ...
 def sqrt(x: Any) -> Any:
@@ -101,7 +102,7 @@ def sqrt(x: Any) -> Any:
     Returns:
         Square root of x
     """
-    if isinstance(x, variable):
+    if isinstance(x, value):
         return add_op('sqrt', [x])
     if isinstance(x, vector):
         return x.map(sqrt)
@@ -111,7 +112,7 @@ def sqrt(x: Any) -> Any:
 @overload
 def sin(x: float | int) -> float: ...
 @overload
-def sin(x: variable[Any]) -> variable[float]: ...
+def sin(x: value[Any]) -> value[float]: ...
 @overload
 def sin(x: vector[Any]) -> vector[float]: ...
 def sin(x: Any) -> Any:
@@ -123,7 +124,7 @@ def sin(x: Any) -> Any:
     Returns:
         Square root of x
     """
-    if isinstance(x, variable):
+    if isinstance(x, value):
         return add_op('sin', [x])
     if isinstance(x, vector):
         return x.map(sin)
@@ -133,7 +134,7 @@ def sin(x: Any) -> Any:
 @overload
 def cos(x: float | int) -> float: ...
 @overload
-def cos(x: variable[Any]) -> variable[float]: ...
+def cos(x: value[Any]) -> value[float]: ...
 @overload
 def cos(x: vector[Any]) -> vector[float]: ...
 def cos(x: Any) -> Any:
@@ -145,7 +146,7 @@ def cos(x: Any) -> Any:
     Returns:
         Cosine of x
     """
-    if isinstance(x, variable):
+    if isinstance(x, value):
         return add_op('cos', [x])
     if isinstance(x, vector):
         return x.map(cos)
@@ -155,7 +156,7 @@ def cos(x: Any) -> Any:
 @overload
 def tan(x: float | int) -> float: ...
 @overload
-def tan(x: variable[Any]) -> variable[float]: ...
+def tan(x: value[Any]) -> value[float]: ...
 @overload
 def tan(x: vector[Any]) -> vector[float]: ...
 def tan(x: Any) -> Any:
@@ -167,7 +168,7 @@ def tan(x: Any) -> Any:
     Returns:
         Tangent of x
     """
-    if isinstance(x, variable):
+    if isinstance(x, value):
         return add_op('tan', [x])
     if isinstance(x, vector):
         #return x.map(tan)
@@ -178,7 +179,7 @@ def tan(x: Any) -> Any:
 @overload
 def atan(x: float | int) -> float: ...
 @overload
-def atan(x: variable[Any]) -> variable[float]: ...
+def atan(x: value[Any]) -> value[float]: ...
 @overload
 def atan(x: vector[Any]) -> vector[float]: ...
 def atan(x: Any) -> Any:
@@ -190,7 +191,7 @@ def atan(x: Any) -> Any:
     Returns:
         Inverse tangent of x
     """
-    if isinstance(x, variable):
+    if isinstance(x, value):
         return add_op('atan', [x])
     if isinstance(x, vector):
         return x.map(atan)
@@ -200,9 +201,9 @@ def atan(x: Any) -> Any:
 @overload
 def atan2(x: float | int, y: float | int) -> float: ...
 @overload
-def atan2(x: variable[Any], y: NumLike) -> variable[float]: ...
+def atan2(x: value[Any], y: NumLike) -> value[float]: ...
 @overload
-def atan2(x: NumLike, y: variable[Any]) -> variable[float]: ...
+def atan2(x: NumLike, y: value[Any]) -> value[float]: ...
 @overload
 def atan2(x: vector[float], y: VecNumLike) -> vector[float]: ...
 @overload
@@ -219,7 +220,7 @@ def atan2(x: VecNumLike, y: VecNumLike) -> Any:
     """
     if isinstance(x, vector) or isinstance(y, vector):
         return _map2(x, y, atan2)
-    if isinstance(x, variable) or isinstance(y, variable):
+    if isinstance(x, value) or isinstance(y, value):
         return add_op('atan2', [x, y])
     return math.atan2(x, y)
 
@@ -227,7 +228,7 @@ def atan2(x: VecNumLike, y: VecNumLike) -> Any:
 @overload
 def asin(x: float | int) -> float: ...
 @overload
-def asin(x: variable[Any]) -> variable[float]: ...
+def asin(x: value[Any]) -> value[float]: ...
 @overload
 def asin(x: vector[Any]) -> vector[float]: ...
 def asin(x: Any) -> Any:
@@ -239,7 +240,7 @@ def asin(x: Any) -> Any:
     Returns:
         Inverse sine of x
     """
-    if isinstance(x, variable):
+    if isinstance(x, value):
         return add_op('asin', [x])
     if isinstance(x, vector):
         return x.map(asin)
@@ -249,7 +250,7 @@ def asin(x: Any) -> Any:
 @overload
 def acos(x: float | int) -> float: ...
 @overload
-def acos(x: variable[Any]) -> variable[float]: ...
+def acos(x: value[Any]) -> value[float]: ...
 @overload
 def acos(x: vector[Any]) -> vector[float]: ...
 def acos(x: Any) -> Any:
@@ -261,7 +262,7 @@ def acos(x: Any) -> Any:
     Returns:
         Inverse cosine of x
     """
-    if isinstance(x, variable):
+    if isinstance(x, value):
         return add_op('acos', [x])
     if isinstance(x, vector):
         return x.map(acos)
@@ -271,18 +272,22 @@ def acos(x: Any) -> Any:
 @overload
 def get_42(x: float | int) -> float: ...
 @overload
-def get_42(x: variable[Any]) -> variable[float]: ...
+def get_42(x: value[Any]) -> value[float]: ...
-def get_42(x: NumLike) -> variable[float] | float:
+def get_42(x: NumLike) -> value[float] | float:
-    """Returns the variable representing the constant 42"""
+    """Returns the value representing the constant 42"""
-    if isinstance(x, variable):
+    if isinstance(x, value):
         return add_op('get_42', [x, x])
     return float((int(x) * 3.0 + 42.0) * 5.0 + 21.0)
 
+
+#TODO: Add vector support
 @overload
 def abs(x: U) -> U: ...
 @overload
-def abs(x: variable[U]) -> variable[U]: ...
+def abs(x: value[U]) -> value[U]: ...
-def abs(x: U | variable[U]) -> Any:
+@overload
+def abs(x: vector[U]) -> vector[U]: ...
+def abs(x: U | value[U] | vector[U]) -> Any:
     """Absolute value function
 
     Arguments:
@@ -291,21 +296,42 @@ def abs(x: U | variable[U]) -> Any:
     Returns:
         Absolute value of x
     """
+    #tt = -x * (x < 0)
     ret = (x < 0) * -x + (x >= 0) * x
-    return ret  # pyright: ignore[reportReturnType]
+    return ret  # REMpyright: ignore[reportReturnType]
 
 
 @overload
-def clamp(x: variable[U], min_value: U | variable[U], max_value: U | variable[U]) -> variable[U]: ...
+def sign(x: U) -> U: ...
 @overload
-def clamp(x: U | variable[U], min_value: variable[U], max_value: U | variable[U]) -> variable[U]: ...
+def sign(x: value[U]) -> value[U]: ...
 @overload
-def clamp(x: U | variable[U], min_value: U | variable[U], max_value: variable[U]) -> variable[U]: ...
+def sign(x: vector[U]) -> vector[U]: ...
+def sign(x: U | value[U] | vector[U]) -> Any:
+    """Return 1 for positive numbers and -1 for negative numbers.
+    For an input of 0 the return value is 0.
+
+    Arguments:
+        x: Input value
+
+    Returns:
+        -1, 0 or 1
+    """
+    ret = (x > 0) - (x < 0)
+    return ret
+
+
+@overload
+def clamp(x: value[U], min_value: U | value[U], max_value: U | value[U]) -> value[U]: ...
+@overload
+def clamp(x: U | value[U], min_value: value[U], max_value: U | value[U]) -> value[U]: ...
+@overload
+def clamp(x: U | value[U], min_value: U | value[U], max_value: value[U]) -> value[U]: ...
 @overload
 def clamp(x: U, min_value: U, max_value: U) -> U: ...
 @overload
-def clamp(x: vector[U], min_value: 'U | variable[U]', max_value: 'U | variable[U]') -> vector[U]: ...
+def clamp(x: vector[U], min_value: 'U | value[U]', max_value: 'U | value[U]') -> vector[U]: ...
-def clamp(x: U | variable[U] | vector[U], min_value: U | variable[U], max_value:  U | variable[U]) -> Any:
+def clamp(x: U | value[U] | vector[U], min_value: U | value[U], max_value:  U | value[U]) -> Any:
     """Clamp function to limit a value between a minimum and maximum.
 
     Arguments:
@@ -325,12 +351,12 @@ def clamp(x: U | variable[U] | vector[U], min_value: U | variable[U], max_value:
 
 
 @overload
-def min(x: variable[U], y: U | variable[U]) -> variable[U]: ...
+def min(x: value[U], y: U | value[U]) -> value[U]: ...
 @overload
-def min(x: U | variable[U], y: variable[U]) -> variable[U]: ...
+def min(x: U | value[U], y: value[U]) -> value[U]: ...
 @overload
 def min(x: U, y: U) -> U: ...
-def min(x: U | variable[U], y: U | variable[U]) -> Any:
+def min(x: U | value[U], y: U | value[U]) -> Any:
     """Minimum function to get the smaller of two values.
 
     Arguments:
@@ -344,12 +370,12 @@ def min(x: U | variable[U], y: U | variable[U]) -> Any:
 
 
 @overload
-def max(x: variable[U], y: U | variable[U]) -> variable[U]: ...
+def max(x: value[U], y: U | value[U]) -> value[U]: ...
 @overload
-def max(x: U | variable[U], y: variable[U]) -> variable[U]: ...
+def max(x: U | value[U], y: value[U]) -> value[U]: ...
 @overload
 def max(x: U, y: U) -> U: ...
-def max(x: U | variable[U], y: U | variable[U]) -> Any:
+def max(x: U | value[U], y: U | value[U]) -> Any:
     """Maximum function to get the larger of two values.
 
     Arguments:
@@ -363,16 +389,16 @@ def max(x: U | variable[U], y: U | variable[U]) -> Any:
 
 
 @overload
-def lerp(v1: variable[U], v2: U | variable[U], t: U | variable[U]) -> variable[U]: ...
+def lerp(v1: value[U], v2: U | value[U], t: unifloat) -> value[U]: ...
 @overload
-def lerp(v1: U | variable[U], v2: variable[U], t: U | variable[U]) -> variable[U]: ...
+def lerp(v1: U | value[U], v2: value[U], t: unifloat) -> value[U]: ...
 @overload
-def lerp(v1: U | variable[U], v2: U | variable[U], t: variable[U]) -> variable[U]: ...
+def lerp(v1: U | value[U], v2: U | value[U], t: value[float]) -> value[U]: ...
 @overload
-def lerp(v1: U, v2: U, t: U) -> U: ...
+def lerp(v1: U, v2: U, t: float) -> U: ...
 @overload
-def lerp(v1: vector[U], v2: vector[U], t: 'U | variable[U]') -> vector[U]: ...
+def lerp(v1: vector[U], v2: vector[U], t: unifloat) -> vector[U]: ...
-def lerp(v1: U | variable[U] | vector[U], v2: U | variable[U] | vector[U], t:  U | variable[U]) -> Any:
+def lerp(v1: U | value[U] | vector[U], v2: U | value[U] | vector[U], t:  unifloat) -> Any:
     """Linearly interpolate between two values or vectors v1 and v2 by a factor t."""
     if isinstance(v1, vector) or isinstance(v2, vector):
         assert isinstance(v1, vector) and isinstance(v2, vector), "None or both v1 and v2 must be vectors."
@@ -381,7 +407,19 @@ def lerp(v1: U | variable[U] | vector[U], v2: U | variable[U] | vector[U], t:  U
     return v1 * (1 - t) + v2 * t
 
 
-def _map2(self: VecNumLike, other: VecNumLike, func: Callable[[Any, Any], variable[U] | U]) -> vector[U]:
+@overload
+def relu(x: U) -> U: ...
+@overload
+def relu(x: value[U]) -> value[U]: ...
+@overload
+def relu(x: vector[U]) -> vector[U]: ...
+def relu(x: U | value[U] | vector[U]) -> Any:
+    """Returns x for x > 0 and otherwise 0."""
+    ret = (x > 0) * x
+    return ret
+
+
+def _map2(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):
         return vector(func(x, y) for x, y in zip(self.values, other.values))

src/copapy/_matrices.py

@@ -0,0 +1,323 @@
+from . import value
+from ._vectors import vector
+from ._mixed import mixed_sum
+from typing import TypeVar, Iterable, Any, overload, TypeAlias, Callable, Iterator, Generic
+from ._helper_types import TNum
+
+MatNumLike: TypeAlias = 'matrix[int] | matrix[float] | value[int] | value[float] | int | float'
+MatIntLike: TypeAlias = 'matrix[int] | value[int] | int'
+MatFloatLike: TypeAlias = 'matrix[float] | value[float] | float'
+U = TypeVar("U", int, float)
+
+
+class matrix(Generic[TNum]):
+    """Mathematical matrix class supporting basic operations and interactions with values.
+    """
+    def __init__(self, values: Iterable[Iterable[TNum | value[TNum]]] | vector[TNum]):
+        """Create a matrix with given values.
+
+        Args:
+            values: iterable of iterable of constant values
+        """
+        if isinstance(values, vector):
+            rows = [values.values]
+        else:
+            rows = [tuple(row) for row in values]
+
+        if rows:
+            row_len = len(rows[0])
+            assert all(len(row) == row_len for row in rows), "All rows must have the same length"
+        self.values: tuple[tuple[value[TNum] | TNum, ...], ...] = tuple(rows)
+        self.rows = len(self.values)
+        self.cols = len(self.values[0]) if self.values else 0
+
+    def __repr__(self) -> str:
+        return f"matrix({self.values})"
+
+    def __len__(self) -> int:
+        """Return the number of rows in the matrix."""
+        return self.rows
+
+    @overload
+    def __getitem__(self, key: int) -> vector[TNum]: ...
+    @overload
+    def __getitem__(self, key: tuple[int, int]) -> value[TNum] | TNum: ...
+    def __getitem__(self, key: int | tuple[int, int]) -> Any:
+        """Get a row as a vector or a specific element.
+            Args:
+                key: row index or (row, col) tuple
+            
+            Returns:
+                vector if row index is given, else the element at (row, col)
+        """
+        if isinstance(key, tuple):
+            assert len(key) == 2
+            return self.values[key[0]][key[1]]
+        else:
+            return vector(self.values[key])
+
+    def __iter__(self) -> Iterator[tuple[value[TNum] | TNum, ...]]:
+        return iter(self.values)
+
+    def __neg__(self) -> 'matrix[TNum]':
+        return matrix((-a for a in row) for row in self.values)
+
+    @overload
+    def __add__(self: 'matrix[int]', other: MatFloatLike) -> 'matrix[float]': ...
+    @overload
+    def __add__(self: 'matrix[int]', other: MatIntLike) -> 'matrix[int]': ...
+    @overload
+    def __add__(self: 'matrix[float]', other: MatNumLike) -> 'matrix[float]': ...
+    @overload
+    def __add__(self, other: MatNumLike) -> 'matrix[int] | matrix[float]': ...
+    def __add__(self, other: MatNumLike) -> Any:
+        if isinstance(other, matrix):
+            assert self.rows == other.rows and self.cols == other.cols, \
+                "Matrices must have the same dimensions"
+            return matrix(
+                tuple(a + b for a, b in zip(row1, row2))
+                for row1, row2 in zip(self.values, other.values)
+            )
+        return matrix(
+            tuple(a + other for a in row)
+            for row in self.values
+        )
+
+    @overload
+    def __radd__(self: 'matrix[float]', other: MatNumLike) -> 'matrix[float]': ...
+    @overload
+    def __radd__(self: 'matrix[int]', other: value[int] | int) -> 'matrix[int]': ...
+    def __radd__(self, other: Any) -> Any:
+        return self + other
+
+    @overload
+    def __sub__(self: 'matrix[int]', other: MatFloatLike) -> 'matrix[float]': ...
+    @overload
+    def __sub__(self: 'matrix[int]', other: MatIntLike) -> 'matrix[int]': ...
+    @overload
+    def __sub__(self: 'matrix[float]', other: MatNumLike) -> 'matrix[float]': ...
+    @overload
+    def __sub__(self, other: MatNumLike) -> 'matrix[int] | matrix[float]': ...
+    def __sub__(self, other: MatNumLike) -> Any:
+        if isinstance(other, matrix):
+            assert self.rows == other.rows and self.cols == other.cols, \
+                "Matrices must have the same dimensions"
+            return matrix(
+                tuple(a - b for a, b in zip(row1, row2))
+                for row1, row2 in zip(self.values, other.values)
+            )
+        return matrix(
+            tuple(a - other for a in row)
+            for row in self.values
+        )
+
+    @overload
+    def __rsub__(self: 'matrix[float]', other: MatNumLike) -> 'matrix[float]': ...
+    @overload
+    def __rsub__(self: 'matrix[int]', other: value[int] | int) -> 'matrix[int]': ...
+    def __rsub__(self, other: MatNumLike) -> Any:
+        if isinstance(other, matrix):
+            assert self.rows == other.rows and self.cols == other.cols, \
+                "Matrices must have the same dimensions"
+            return matrix(
+                tuple(b - a for a, b in zip(row1, row2))
+                for row1, row2 in zip(self.values, other.values)
+            )
+        return matrix(
+            tuple(other - a for a in row)
+            for row in self.values
+        )
+
+    @overload
+    def __mul__(self: 'matrix[int]', other: MatFloatLike) -> 'matrix[float]': ...
+    @overload
+    def __mul__(self: 'matrix[int]', other: MatIntLike) -> 'matrix[int]': ...
+    @overload
+    def __mul__(self: 'matrix[float]', other: MatNumLike) -> 'matrix[float]': ...
+    @overload
+    def __mul__(self, other: MatNumLike) -> 'matrix[int] | matrix[float]': ...
+    def __mul__(self, other: MatNumLike) -> Any:
+        """Element-wise multiplication"""
+        if isinstance(other, matrix):
+            assert self.rows == other.rows and self.cols == other.cols, \
+                "Matrices must have the same dimensions"
+            return matrix(
+                tuple(a * b for a, b in zip(row1, row2))
+                for row1, row2 in zip(self.values, other.values)
+            )
+        return matrix(
+            tuple(a * other for a in row)
+            for row in self.values
+        )
+
+    @overload
+    def __rmul__(self: 'matrix[float]', other: MatNumLike) -> 'matrix[float]': ...
+    @overload
+    def __rmul__(self: 'matrix[int]', other: value[int] | int) -> 'matrix[int]': ...
+    def __rmul__(self, other: MatNumLike) -> Any:
+        return self * other
+
+    def __truediv__(self, other: MatNumLike) -> 'matrix[float]':
+        """Element-wise division"""
+        if isinstance(other, matrix):
+            assert self.rows == other.rows and self.cols == other.cols, \
+                "Matrices must have the same dimensions"
+            return matrix(
+                tuple(a / b for a, b in zip(row1, row2))
+                for row1, row2 in zip(self.values, other.values)
+            )
+        return matrix(
+            tuple(a / other for a in row)
+            for row in self.values
+        )
+
+    def __rtruediv__(self, other: MatNumLike) -> 'matrix[float]':
+        if isinstance(other, matrix):
+            assert self.rows == other.rows and self.cols == other.cols, \
+                "Matrices must have the same dimensions"
+            return matrix(
+                tuple(b / a for a, b in zip(row1, row2))
+                for row1, row2 in zip(self.values, other.values)
+            )
+        return matrix(
+            tuple(other / a for a in row)
+            for row in self.values
+        )
+
+    @overload
+    def __matmul__(self: 'matrix[TNum]', other: 'vector[TNum]') -> 'vector[TNum]': ...
+    @overload
+    def __matmul__(self: 'matrix[TNum]', other: 'matrix[TNum]') -> 'matrix[TNum]': ...
+    def __matmul__(self: 'matrix[TNum]', other: 'matrix[TNum] | vector[TNum]') -> 'matrix[TNum] | vector[TNum]':
+        """Matrix multiplication using @ operator"""
+        if isinstance(other, vector):
+            assert self.cols == len(other.values), \
+                f"Matrix columns ({self.cols}) must match vector length ({len(other.values)})"
+            vec_result = (mixed_sum(a * b for a, b in zip(row, other.values)) for row in self.values)
+            return vector(vec_result)
+        else:
+            assert isinstance(other, matrix), "Cannot multiply matrix with {type(other)}"
+            assert self.cols == other.rows, \
+                f"Matrix columns ({self.cols}) must match other matrix rows ({other.rows})"
+            result: list[list[TNum | value[TNum]]] = []
+            for row in self.values:
+                new_row: list[TNum | value[TNum]] = []
+                for col_idx in range(other.cols):
+                    col = tuple(other.values[i][col_idx] for i in range(other.rows))
+                    element = sum(a * b for a, b in zip(row, col))
+                    new_row.append(element)
+                result.append(new_row)
+            return matrix(result)
+
+    def transpose(self) -> 'matrix[TNum]':
+        """Return the transpose of the matrix."""
+        if not self.values:
+            return matrix([])
+        return matrix(
+            tuple(self.values[i][j] for i in range(self.rows))
+            for j in range(self.cols)
+        )
+    
+    @property
+    def shape(self) -> tuple[int, int]:
+        """Return the shape of the matrix as (rows, cols)."""
+        return (self.rows, self.cols)
+
+    @property
+    def T(self) ->  'matrix[TNum]':
+        return self.transpose()
+
+    def row(self, index: int) -> vector[TNum]:
+        """Get a row as a vector."""
+        assert 0 <= index < self.rows, f"Row index {index} out of bounds"
+        return vector(self.values[index])
+
+    def col(self, index: int) -> vector[TNum]:
+        """Get a column as a vector."""
+        assert 0 <= index < self.cols, f"Column index {index} out of bounds"
+        return vector(self.values[i][index] for i in range(self.rows))
+
+    @overload
+    def trace(self: 'matrix[TNum]') -> TNum | value[TNum]: ...
+    @overload
+    def trace(self: 'matrix[int]') -> int | value[int]: ...
+    @overload
+    def trace(self: 'matrix[float]') -> float | value[float]: ...
+    def trace(self) -> Any:
+        """Calculate the trace (sum of diagonal elements)."""
+        assert self.rows == self.cols, "Trace is only defined for square matrices"
+        return mixed_sum(self.values[i][i] for i in range(self.rows))
+
+    @overload
+    def sum(self: 'matrix[TNum]') -> TNum | value[TNum]: ...
+    @overload
+    def sum(self: 'matrix[int]') -> int | value[int]: ...
+    @overload
+    def sum(self: 'matrix[float]') -> float | value[float]: ...
+    def sum(self) -> Any:
+        """Calculate the sum of all elements."""
+        return mixed_sum(a for row in self.values for a in row)
+
+    def map(self, func: Callable[[Any], value[U] | U]) -> 'matrix[U]':
+        """Applies a function to each element of the matrix and returns a new matrix."""
+        return matrix(
+            tuple(func(a) for a in row)
+            for row in self.values
+        )
+
+    def homogenize(self) -> 'matrix[TNum]':
+        """Convert all elements to copapy values if any element is a copapy value."""
+        if any(isinstance(val, value) for row in self.values for val in row):
+            return matrix(
+                tuple(value(val) if not isinstance(val, value) else val for val in row)
+                for row in self.values
+            )
+        else:
+            return self
+
+
+def identity(size: int) -> matrix[int]:
+    """Create an identity matrix of given size."""
+    return matrix(
+        tuple(1 if i == j else 0 for j in range(size))
+        for i in range(size)
+    )
+
+
+def zeros(rows: int, cols: int) -> matrix[int]:
+    """Create a zero matrix of given dimensions."""
+    return matrix(
+        tuple(0 for _ in range(cols))
+        for _ in range(rows)
+    )
+
+
+def ones(rows: int, cols: int) -> matrix[int]:
+    """Create a matrix of ones with given dimensions."""
+    return matrix(
+        tuple(1 for _ in range(cols))
+        for _ in range(rows)
+    )
+
+
+def eye(rows: int, cols: int | None = None) -> matrix[int]:
+    """Create a matrix with ones on the diagonal and zeros elsewhere."""
+    cols = cols if cols else rows
+    return matrix(
+        tuple(1 if i == j else 0 for j in range(cols))
+        for i in range(rows)
+    )
+
+
+@overload
+def diagonal(vec: 'vector[int]') -> matrix[int]: ...
+@overload
+def diagonal(vec: 'vector[float]') -> matrix[float]: ...
+def diagonal(vec: vector[Any]) -> matrix[Any]:
+    """Create a diagonal matrix from a vector."""
+    size = len(vec)
+
+    return matrix(
+        tuple(vec[i] if i == j else 0 for j in range(size))
+        for i in range(size)
+    )

src/copapy/_mixed.py

@@ -0,0 +1,24 @@
+
+from . import value
+from typing import TypeVar, Iterable, Any, overload
+
+T = TypeVar("T", int, float)
+
+
+@overload
+def mixed_sum(scalars: Iterable[float | value[float]]) -> float | value[float]: ...
+@overload
+def mixed_sum(scalars: Iterable[int | value[int]]) -> int | value[int]: ...
+@overload
+def mixed_sum(scalars: Iterable[T | value[T]]) -> T | value[T]: ...
+def mixed_sum(scalars: Iterable[int | float | value[Any]]) -> Any:
+    sl = list(scalars)
+    return sum(a for a in sl if not isinstance(a, value)) +\
+           sum(a for a in sl if isinstance(a, value))
+
+
+def mixed_homogenize(scalars: Iterable[T | value[T]]) -> Iterable[T] | Iterable[value[T]]:
+    if any(isinstance(val, value) for val in scalars):
+        return (value(val) if not isinstance(val, value) else val for val in scalars)
+    else:
+        return (val for val in scalars if not isinstance(val, value))

src/copapy/_target.py

@@ -1,11 +1,13 @@
-from typing import Iterable, overload
+from typing import Iterable, overload, TypeVar, Any
 from . import _binwrite as binw
 from coparun_module import coparun, read_data_mem
 import struct
 from ._basic_types import stencil_db_from_package
-from ._basic_types import variable, Net, Node, Write, NumLike
+from ._basic_types import value, Net, Node, Write, NumLike
 from ._compiler import compile_to_dag
 
+T = TypeVar("T", int, float)
+
 
 def add_read_command(dw: binw.data_writer, variables: dict[Net, tuple[int, int, str]], net: Net) -> None:
     assert net in variables, f"Variable {net} not found in data writer variables"
@@ -26,24 +28,25 @@ class Target():
             optimization: Optimization level
         """
         self.sdb = stencil_db_from_package(arch, optimization)
-        self._variables: dict[Net, tuple[int, int, str]] = {}
+        self._values: dict[Net, tuple[int, int, str]] = {}
 
-    def compile(self, *variables: int | float | variable[int] | variable[float] | variable[bool] | Iterable[int | float | variable[int] | variable[float] | variable[bool]]) -> None:
+    def compile(self, *values: int | float | value[int] | value[float] | Iterable[int | float | value[int] | value[float]]) -> None:
-        """Compiles the code to compute the given variables.
+        """Compiles the code to compute the given values.
 
         Arguments:
-            variables: Variables to compute
+            values: Values to compute
         """
         nodes: list[Node] = []
-        for s in variables:
+        for s in values:
             if isinstance(s, Iterable):
                 for net in s:
-                    assert isinstance(net, Net), f"The folowing element is not a Net: {net}"
+                    if isinstance(net, Net):
                         nodes.append(Write(net))
             else:
+                if isinstance(s, Net):
                     nodes.append(Write(s))
 
-        dw, self._variables = compile_to_dag(nodes, self.sdb)
+        dw, self._values = compile_to_dag(nodes, self.sdb)
         dw.write_com(binw.Command.END_COM)
         assert coparun(dw.get_data()) > 0
 
@@ -56,59 +59,56 @@ class Target():
         assert coparun(dw.get_data()) > 0
 
     @overload
-    def read_value(self, net: variable[bool]) -> bool:
+    def read_value(self, net: value[T]) -> T: ...
-        ...
-
     @overload
-    def read_value(self, net: variable[float]) -> float:
+    def read_value(self, net: NumLike) -> float | int | bool: ...
-        ...
-
     @overload
-    def read_value(self, net: variable[int]) -> int:
+    def read_value(self, net: Iterable[T | value[T]]) -> list[T]: ...
-        ...
+    def read_value(self, net: NumLike | value[T] | Iterable[T | value[T]]) -> Any:
-
+        """Reads the numeric value of a copapy type.
-    @overload
-    def read_value(self, net: NumLike) -> float | int | bool:
-        ...
-
-    def read_value(self, net: NumLike) -> float | int | bool:
-        """Reads the value of a variable.
 
         Arguments:
-            net: Variable to read
+            net: Values to read
 
         Returns:
-            Value of the variable
+            Numeric value
         """
-        assert isinstance(net, Net), "Variable must be a copapy variable object"
+        if isinstance(net, Iterable):
-        assert net in self._variables, f"Variable {net} not found. It might not have been compiled for the target."
+            return [self.read_value(ni) if isinstance(ni, value) else ni for ni in net]
-        addr, lengths, var_type = self._variables[net]
+
+        if isinstance(net, float | int):
+            print("Warning: value is not a copypy value")
+            return net
+
+        assert isinstance(net, Net), "Argument must be a copapy value"
+        assert net in self._values, f"Value {net} not found. It might not have been compiled for the target."
+        addr, lengths, var_type = self._values[net]
         assert lengths > 0
         data = read_data_mem(addr, lengths)
-        assert data is not None and len(data) == lengths, f"Failed to read variable {net}"
+        assert data is not None and len(data) == lengths, f"Failed to read value {net}"
         en = {'little': '<', 'big': '>'}[self.sdb.byteorder]
         if var_type == 'float':
             if lengths == 4:
-                value = struct.unpack(en + 'f', data)[0]
+                val = struct.unpack(en + 'f', data)[0]
             elif lengths == 8:
-                value = struct.unpack(en + 'd', data)[0]
+                val = struct.unpack(en + 'd', data)[0]
             else:
                 raise ValueError(f"Unsupported float length: {lengths} bytes")
-            assert isinstance(value, float)
+            assert isinstance(val, float)
-            return value
+            return val
         elif var_type == 'int':
             assert lengths in (1, 2, 4, 8), f"Unsupported int length: {lengths} bytes"
-            value = int.from_bytes(data, byteorder=self.sdb.byteorder, signed=True)
+            val = int.from_bytes(data, byteorder=self.sdb.byteorder, signed=True)
-            return value
+            return val
         elif var_type == 'bool':
             assert lengths in (1, 2, 4, 8), f"Unsupported int length: {lengths} bytes"
-            value = bool.from_bytes(data, byteorder=self.sdb.byteorder, signed=True)
+            val = bool.from_bytes(data, byteorder=self.sdb.byteorder, signed=True)
-            return value
+            return val
         else:
-            raise ValueError(f"Unsupported variable type: {var_type}")
+            raise ValueError(f"Unsupported value type: {var_type}")
 
     def read_value_remote(self, net: Net) -> None:
-        """Reads the raw data of a variable by the runner."""
+        """Reads the raw data of a value by the runner."""
         dw = binw.data_writer(self.sdb.byteorder)
-        add_read_command(dw, self._variables, net)
+        add_read_command(dw, self._values, net)
         assert coparun(dw.get_data()) > 0

src/copapy/_vectors.py

@@ -1,26 +1,28 @@
-from . import variable
+from . import value
-from typing import Generic, TypeVar, Iterable, Any, overload, TypeAlias, Callable
+from ._mixed import mixed_sum, mixed_homogenize
+from typing import TypeVar, Iterable, Any, overload, TypeAlias, Callable, Iterator, Generic
 import copapy as cp
+from ._helper_types import TNum
 
-VecNumLike: TypeAlias = 'vector[int] | vector[float] | variable[int] | variable[float] | variable[bool] | int | float | bool'
+#VecNumLike: TypeAlias = 'vector[int] | vector[float] | value[int] | value[float] | int | float | bool'
-VecIntLike: TypeAlias = 'vector[int] | variable[int] | int'
+VecNumLike: TypeAlias = 'vector[Any] | value[Any] | int | float | bool'
-VecFloatLike: TypeAlias = 'vector[float] | variable[float] | float'
+VecIntLike: TypeAlias = 'vector[int] | value[int] | int'
-T = TypeVar("T", int, float)
+VecFloatLike: TypeAlias = 'vector[float] | value[float] | float'
 U = TypeVar("U", int, float)
 
 epsilon = 1e-20
 
 
-class vector(Generic[T]):
+class vector(Generic[TNum]):
-    """Mathematical vector class supporting basic operations and interactions with variables.
+    """Mathematical vector class supporting basic operations and interactions with values.
     """
-    def __init__(self, values: Iterable[T | variable[T]]):
+    def __init__(self, values: Iterable[TNum | value[TNum]]):
-        """Create a vector with given values and variables.
+        """Create a vector with given values.
 
         Args:
-            values: iterable of constant values and variables
+            values: iterable of constant values
         """
-        self.values: tuple[variable[T] | T, ...] = tuple(values)
+        self.values: tuple[value[TNum] | TNum, ...] = tuple(values)
 
     def __repr__(self) -> str:
         return f"vector({self.values})"
@@ -28,9 +30,21 @@ class vector(Generic[T]):
     def __len__(self) -> int:
         return len(self.values)
 
-    def __getitem__(self, index: int) -> variable[T] | T:
+    @overload
+    def __getitem__(self, index: int) -> value[TNum] | TNum: ...
+    @overload
+    def __getitem__(self, index: slice) -> 'vector[TNum]': ...
+    def __getitem__(self, index: int | slice) -> 'vector[TNum] | value[TNum] | TNum':
+        if isinstance(index, slice):
+            return vector(self.values[index])
         return self.values[index]
 
+    def __neg__(self) -> 'vector[TNum]':
+        return vector(-a for a in self.values)
+
+    def __iter__(self) -> Iterator[value[TNum] | TNum]:
+        return iter(self.values)
+
     @overload
     def __add__(self: 'vector[int]', other: VecFloatLike) -> 'vector[float]': ...
     @overload
@@ -48,7 +62,9 @@ class vector(Generic[T]):
     @overload
     def __radd__(self: 'vector[float]', other: VecNumLike) -> 'vector[float]': ...
     @overload
-    def __radd__(self: 'vector[int]', other: variable[int] | int) -> 'vector[int]': ...
+    def __radd__(self: 'vector[int]', other: value[int] | int) -> 'vector[int]': ...
+    @overload
+    def __radd__(self, other: VecNumLike) -> 'vector[Any]': ...
     def __radd__(self, other: Any) -> Any:
         return self + other
 
@@ -69,7 +85,9 @@ class vector(Generic[T]):
     @overload
     def __rsub__(self: 'vector[float]', other: VecNumLike) -> 'vector[float]': ...
     @overload
-    def __rsub__(self: 'vector[int]', other: variable[int] | int) -> 'vector[int]': ...
+    def __rsub__(self: 'vector[int]', other: value[int] | int) -> 'vector[int]': ...
+    @overload
+    def __rsub__(self, other: VecNumLike) -> 'vector[Any]': ...
     def __rsub__(self, other: VecNumLike) -> Any:
         if isinstance(other, vector):
             assert len(self.values) == len(other.values)
@@ -93,10 +111,35 @@ class vector(Generic[T]):
     @overload
     def __rmul__(self: 'vector[float]', other: VecNumLike) -> 'vector[float]': ...
     @overload
-    def __rmul__(self: 'vector[int]', other: variable[int] | int) -> 'vector[int]': ...
+    def __rmul__(self: 'vector[int]', other: value[int] | int) -> 'vector[int]': ...
+    @overload
+    def __rmul__(self, other: VecNumLike) -> 'vector[Any]': ...
     def __rmul__(self, other: VecNumLike) -> Any:
         return self * other
 
+    @overload
+    def __pow__(self: 'vector[int]', other: VecFloatLike) -> 'vector[float]': ...
+    @overload
+    def __pow__(self: 'vector[int]', other: VecIntLike) -> 'vector[int]': ...
+    @overload
+    def __pow__(self: 'vector[float]', other: VecNumLike) -> 'vector[float]': ...
+    @overload
+    def __pow__(self, other: VecNumLike) -> 'vector[int] | vector[float]': ...
+    def __pow__(self, other: VecNumLike) -> Any:
+        if isinstance(other, vector):
+            assert len(self.values) == len(other.values)
+            return vector(a ** b for a, b in zip(self.values, other.values))
+        return vector(a ** other for a in self.values)
+
+    @overload
+    def __rpow__(self: 'vector[float]', other: VecNumLike) -> 'vector[float]': ...
+    @overload
+    def __rpow__(self: 'vector[int]', other: value[int] | int) -> 'vector[int]': ...
+    @overload
+    def __rpow__(self, other: VecNumLike) -> 'vector[Any]': ...
+    def __rpow__(self, other: VecNumLike) -> Any:
+        return self ** other
+
     def __truediv__(self, other: VecNumLike) -> 'vector[float]':
         if isinstance(other, vector):
             assert len(self.values) == len(other.values)
@@ -110,26 +153,26 @@ class vector(Generic[T]):
         return vector(other / a for a in self.values)
 
     @overload
-    def dot(self: 'vector[int]', other: 'vector[int]') -> int | variable[int]: ...
+    def dot(self: 'vector[int]', other: 'vector[int]') -> int | value[int]: ...
     @overload
-    def dot(self, other: 'vector[float]') -> float | variable[float]: ...
+    def dot(self, other: 'vector[float]') -> float | value[float]: ...
     @overload
-    def dot(self: 'vector[float]', other: 'vector[int] | vector[float]') -> float | variable[float]: ...
+    def dot(self: 'vector[float]', other: 'vector[int] | vector[float]') -> float | value[float]: ...
     @overload
-    def dot(self, other: 'vector[int] | vector[float]') -> float | int | variable[float] | variable[int]: ...
+    def dot(self, other: 'vector[int] | vector[float]') -> float | int | value[float] | value[int]: ...
     def dot(self, other: 'vector[int] | vector[float]') -> Any:
         assert len(self.values) == len(other.values), "Vectors must be of same length."
-        return sum(a * b for a, b in zip(self.values, other.values))
+        return mixed_sum(a * b for a, b in zip(self.values, other.values))
 
     # @ operator
     @overload
-    def __matmul__(self: 'vector[int]', other: 'vector[int]') -> int | variable[int]: ...
+    def __matmul__(self: 'vector[int]', other: 'vector[int]') -> int | value[int]: ...
     @overload
-    def __matmul__(self, other: 'vector[float]') -> float | variable[float]: ...
+    def __matmul__(self, other: 'vector[float]') -> float | value[float]: ...
     @overload
-    def __matmul__(self: 'vector[float]', other: 'vector[int] | vector[float]') -> float | variable[float]: ...
+    def __matmul__(self: 'vector[float]', other: 'vector[int] | vector[float]') -> float | value[float]: ...
     @overload
-    def __matmul__(self, other: 'vector[int] | vector[float]') -> float | int | variable[float] | variable[int]: ...
+    def __matmul__(self, other: 'vector[int] | vector[float]') -> float | int | value[float] | value[int]: ...
     def __matmul__(self, other: 'vector[int] | vector[float]') -> Any:
         return self.dot(other)
 
@@ -144,55 +187,93 @@ class vector(Generic[T]):
             a1 * b2 - a2 * b1
         ])
     
+    def __gt__(self, other: VecNumLike) -> 'vector[int]':
+        if isinstance(other, vector):
+            assert len(self.values) == len(other.values)
+            return vector(a > b for a, b in zip(self.values, other.values))
+        return vector(a > other for a in self.values)
+
+    def __lt__(self, other: VecNumLike) -> 'vector[int]':
+        if isinstance(other, vector):
+            assert len(self.values) == len(other.values)
+            return vector(a < b for a, b in zip(self.values, other.values))
+        return vector(a < other for a in self.values)
+
+    def __ge__(self, other: VecNumLike) -> 'vector[int]':
+        if isinstance(other, vector):
+            assert len(self.values) == len(other.values)
+            return vector(a >= b for a, b in zip(self.values, other.values))
+        return vector(a >= other for a in self.values)
+
+    def __le__(self, other: VecNumLike) -> 'vector[int]':
+        if isinstance(other, vector):
+            assert len(self.values) == len(other.values)
+            return vector(a <= b for a, b in zip(self.values, other.values))
+        return vector(a <= other for a in self.values)
+
+    def __eq__(self, other: VecNumLike) -> 'vector[int]':  # type: ignore
+        if isinstance(other, vector):
+            assert len(self.values) == len(other.values)
+            return vector(a == b for a, b in zip(self.values, other.values))
+        return vector(a == other for a in self.values)
+
+    def __ne__(self, other: VecNumLike) -> 'vector[int]':  # type: ignore
+        if isinstance(other, vector):
+            assert len(self.values) == len(other.values)
+            return vector(a != b for a, b in zip(self.values, other.values))
+        return vector(a != other for a in self.values)
+    
+    @property
+    def shape(self) -> tuple[int]:
+        """Return the shape of the vector as (length,)."""
+        return (len(self.values),)
+
     @overload
-    def sum(self: 'vector[int]') -> int | variable[int]: ...
+    def sum(self: 'vector[int]') -> int | value[int]: ...
     @overload
-    def sum(self: 'vector[float]') -> float | variable[float]: ...
+    def sum(self: 'vector[float]') -> float | value[float]: ...
     def sum(self) -> Any:
         """Sum of all vector elements."""
-        return sum(a for a in self.values if isinstance(a, variable)) +\
+        return mixed_sum(self.values)
-               sum(a for a in self.values if not isinstance(a, variable))
 
-    def magnitude(self) -> 'float | variable[float]':
+    def magnitude(self) -> 'float | value[float]':
         """Magnitude (length) of the vector."""
-        s = sum(a * a for a in self.values)
+        s = mixed_sum(a * a for a in self.values)
-        return cp.sqrt(s) if isinstance(s, variable) else cp.sqrt(s)
+        return cp.sqrt(s)
 
     def normalize(self) -> 'vector[float]':
         """Returns a normalized (unit length) version of the vector."""
         mag = self.magnitude() + epsilon
         return self / mag
 
-    def __neg__(self) -> 'vector[float] | vector[int]':
+    def homogenize(self) -> 'vector[TNum]':
-        return vector(-a for a in self.values)
+        if any(isinstance(val, value) for val in self.values):
+            return vector(mixed_homogenize(self))
+        else:
+            return self
 
-    def __iter__(self) -> Iterable[variable[T] | T]:
+    def map(self, func: Callable[[Any], value[U] | U]) -> 'vector[U]':
-        return iter(self.values)
-
-    def map(self, func: Callable[[Any], variable[U] | U]) -> 'vector[U]':
         """Applies a function to each element of the vector and returns a new vector."""
         return vector(func(x) for x in self.values)
 
 
-# Utility functions for 3D vectors with two arguments
-
 def cross_product(v1: vector[float], v2: vector[float]) -> vector[float]:
     """Calculate the cross product of two 3D vectors."""
     return v1.cross(v2)
 
 
-def dot_product(v1: vector[float], v2: vector[float]) -> 'float | variable[float]':
+def dot_product(v1: vector[float], v2: vector[float]) -> 'float | value[float]':
     """Calculate the dot product of two vectors."""
     return v1.dot(v2)
 
 
-def distance(v1: vector[float], v2: vector[float]) -> 'float | variable[float]':
+def distance(v1: vector[float], v2: vector[float]) -> 'float | value[float]':
     """Calculate the Euclidean distance between two vectors."""
     diff = v1 - v2
     return diff.magnitude()
 
 
-def scalar_projection(v1: vector[float], v2: vector[float]) -> 'float | variable[float]':
+def scalar_projection(v1: vector[float], v2: vector[float]) -> 'float | value[float]':
     """Calculate the scalar projection of v1 onto v2."""
     dot_prod = v1.dot(v2)
     mag_v2 = v2.magnitude() + epsilon
@@ -207,7 +288,7 @@ def vector_projection(v1: vector[float], v2: vector[float]) -> vector[float]:
     return v2 * scalar_proj
 
 
-def angle_between(v1: vector[float], v2: vector[float]) -> 'float | variable[float]':
+def angle_between(v1: vector[float], v2: vector[float]) -> 'float | value[float]':
     """Calculate the angle in radians between two vectors."""
     dot_prod = v1.dot(v2)
     mag_v1 = v1.magnitude()
@@ -216,7 +297,7 @@ def angle_between(v1: vector[float], v2: vector[float]) -> 'float | variable[flo
     return cp.acos(cos_angle)
 
 
-def rotate_vector(v: vector[float], axis: vector[float], angle: 'float | variable[float]') -> vector[float]:
+def rotate_vector(v: vector[float], axis: vector[float], angle: 'float | value[float]') -> vector[float]:
     """Rotate vector v around a given axis by a specified angle using Rodrigues' rotation formula."""
     k = axis.normalize()
     cos_angle = cp.cos(angle)

src/copapy/backend.py

@@ -1,7 +1,7 @@
 from ._target import add_read_command
 from ._basic_types import Net, Op, Node, CPConstant, Write, stencil_db_from_package
 from ._compiler import compile_to_dag, \
-    stable_toposort, get_const_nets, get_all_dag_edges, add_read_ops, \
+    stable_toposort, get_const_nets, get_all_dag_edges, add_read_ops, get_all_dag_edges_between, \
     add_write_ops
 
 __all__ = [
@@ -15,6 +15,7 @@ __all__ = [
     "stable_toposort",
     "get_const_nets",
     "get_all_dag_edges",
+    "get_all_dag_edges_between",
     "add_read_ops",
     "add_write_ops",
     "stencil_db_from_package"

src/copapy/filters.py

@@ -0,0 +1,65 @@
+from . import value, vector
+from ._basic_types import iif, unifloat
+from._helper_types import TNum
+from typing import Any, Iterable
+
+
+def homogenize_vector(input_values: Iterable[TNum | value[TNum]]) -> Iterable[TNum] | Iterable[value[TNum]]:
+    input_list = list(input_values)
+    if any(isinstance(val, value) for val in input_list):
+        return (v if isinstance(v, value) else value(v) for v in input_list)
+    else:
+        return (v for v in input_list if not isinstance(v, value))
+
+
+def _inv_argsort(input_vector: vector[TNum]) -> vector[int]:
+    positions = (sum((v1 > v2) for v2 in input_vector) for v1 in input_vector)
+    return vector(positions)
+
+
+def argsort(input_vector: vector[TNum]) -> vector[int]:
+    """
+    Perform an indirect sort. It returns an array of indices that index data
+    in sorted order.
+
+    Args:
+        input_vector: The input vector containing numerical values.
+
+    Returns:
+        Index array.
+    """
+    return _inv_argsort(_inv_argsort(input_vector))
+
+
+def median(input_vector: vector[TNum]) -> TNum | value[TNum]:
+    """
+    Applies a median filter to the input vector and returns the median as a unifloat.
+
+    Args:
+        input_vector: The input vector containing numerical values.
+
+    Returns:
+        The median value of the input vector.
+    """
+    vec = input_vector
+    ret = vec[0]
+    for v1 in vec:
+        n2 = len(vec) // 2 + 1
+        lt = sum(v1 < v2 for v2 in vec)
+        gt = sum(v1 > v2 for v2 in vec)
+        ret = iif((lt < n2) & (gt < n2), v1, ret)
+
+    return ret
+
+
+def mean(input_vector: vector[Any]) -> unifloat:
+    """
+    Applies a mean filter to the input vector and returns the mean as a unifloat.
+
+    Args:
+        input_vector (vector): The input vector containing numerical values.
+
+    Returns:
+        unifloat: The mean value of the input vector.
+    """
+    return input_vector.sum() / len(input_vector)

tests/benchmark.py

@@ -0,0 +1,224 @@
+import copapy as cp
+import time
+import json
+import os
+import subprocess
+import sys
+import numpy as np
+from numpy.core._multiarray_umath import __cpu_features__
+
+from copapy._matrices import diagonal
+
+CPU_SIMD_FEATURES = "SSE SSE2 SSE3 SSSE3 SSE41 SSE42 AVX AVX2 AVX512F FMA3"
+
+
+def cp_vs_python(path: str):
+    os.environ.get("NPY_DISABLE_CPU_FEATURES")
+    cpu_f = CPU_SIMD_FEATURES.split(' ')
+    print('\n'.join(f"> {k}: {v}" for k, v in __cpu_features__.items() if k in cpu_f))
+
+
+    results: list[dict[str, str | float | int]] = []
+
+    for _ in range(15):
+        for v_size in [10, 30, 60] + list(range(100, 600, 100)):
+
+            sum_size = 10
+            #v_size = 400
+            iter_size = 30000
+
+            v1 = cp.vector(cp.value(float(v)) for v in range(v_size))
+            v2 = cp.vector(cp.value(float(v)) for v in [5]*v_size)
+
+            v3 = sum((v1 + i) @ v2 for i in range(sum_size))
+
+            tg = cp.Target()
+            tg.compile(v3)
+
+            time.sleep(0.1)
+            t0 = time.perf_counter()
+            for _ in range(iter_size):
+                tg.run()
+            elapsed_cp = time.perf_counter() - t0
+
+            #print(f"Copapy: {elapsed_cp:.4f} s")
+            results.append({'benchmark': 'Copapy', 'iter_size': iter_size, 'elapsed_time': elapsed_cp, 'sum_size': sum_size, 'v_size': v_size})
+
+
+
+            v1 = cp.vector(float(v) for v in range(v_size))
+            v2 = cp.vector(float(v) for v in [5]*v_size)
+
+            time.sleep(0.1)
+            t0 = time.perf_counter()
+            for _ in range(iter_size//100):
+                v3 = sum((v1 + i) @ v2 for i in range(sum_size))
+
+            elapsed_python = time.perf_counter() - t0
+
+            #print(f"Python: {elapsed_python:.4f} s")
+            results.append({'benchmark': 'Python','iter_size': iter_size//10, 'elapsed_time': elapsed_python, 'sum_size': sum_size, 'v_size': v_size})
+
+            v1 = np.array(list(range(v_size)), dtype=np.float32)
+            v2 = np.array([5]*v_size, dtype=np.float32)
+            i = np.array(list(range(sum_size)), dtype=np.int32).reshape([sum_size, 1])
+
+            time.sleep(0.1)
+            t0 = time.perf_counter()
+            for _ in range(iter_size):
+                v3 = np.sum((v1 + i) @ v2)
+
+            elapsed_np = time.perf_counter() - t0
+
+            #print(f"Numpy 2: {elapsed_np2:.4f} s")
+            results.append({'benchmark': 'NumPy', 'iter_size': iter_size, 'elapsed_time': elapsed_np, 'sum_size': sum_size, 'v_size': v_size})
+
+
+            print(f"{v_size} {elapsed_cp}, {elapsed_python}, {elapsed_np}")
+
+    with open(path, 'w') as f:
+        json.dump(results, f)
+
+
+def cp_vs_python_sparse(path: str = 'benchmark_results_001_sparse.json'):
+    results: list[dict[str, str | float | int]] = []
+
+    for _ in range(7):
+        for v_size in [8, 8, 16, 20, 24, 32]:
+
+            n_ones = int((v_size ** 2) * 0.5)
+            n_zeros = (v_size ** 2) - n_ones
+            mask = np.array([1] * n_ones + [0] * n_zeros).reshape((v_size, v_size))
+            np.random.shuffle(mask)
+
+            sum_size = 10
+            #v_size = 400
+            iter_size = 3000
+
+            v1 = cp.vector(cp.value(float(v)) for v in range(v_size))
+            v2 = cp.vector(cp.value(float(v)) for v in [5]*v_size)
+
+            test = cp.vector(np.linspace(0, 1, v_size))
+
+            assert False, test * v2
+
+            v3 = sum(((cp.diagonal(v1) + i) * cp.matrix(mask)) @ v2 for i in range(sum_size))
+
+            tg = cp.Target()
+            tg.compile(v3)
+
+            time.sleep(0.1)
+            t0 = time.perf_counter()
+            for _ in range(iter_size):
+                tg.run()
+            elapsed_cp = time.perf_counter() - t0
+
+            #print(f"Copapy: {elapsed_cp:.4f} s")
+            results.append({'benchmark': 'Copapy', 'iter_size': iter_size, 'elapsed_time': elapsed_cp, 'sum_size': sum_size, 'v_size': v_size})
+
+
+
+            v1 = cp.vector(float(v) for v in range(v_size))
+            v2 = cp.vector(float(v) for v in [5]*v_size)
+
+            time.sleep(0.1)
+            t0 = time.perf_counter()
+            for _ in range(iter_size//1000):
+                v3 = sum(((cp.diagonal(v1) + i) * cp.matrix(mask)) @ v2 for i in range(sum_size))
+
+            elapsed_python = time.perf_counter() - t0
+
+            #print(f"Python: {elapsed_python:.4f} s")
+            results.append({'benchmark': 'Python','iter_size': iter_size//10, 'elapsed_time': elapsed_python, 'sum_size': sum_size, 'v_size': v_size})
+
+            v1 = np.array(list(range(v_size)), dtype=np.float32)
+            v2 = np.array([5]*v_size, dtype=np.float32)
+            i_arr = np.array(list(range(sum_size)), dtype=np.int32).reshape([sum_size, 1, 1])
+            tmp1 = v1 * np.eye(v_size) + i_arr
+
+            time.sleep(0.1)
+            t0 = time.perf_counter()
+            for _ in range(iter_size):
+                v3 = np.sum(((tmp1) * mask) @ v2)
+
+            elapsed_np = time.perf_counter() - t0
+
+            #print(f"Numpy 2: {elapsed_np2:.4f} s")
+            results.append({'benchmark': 'NumPy', 'iter_size': iter_size, 'elapsed_time': elapsed_np, 'sum_size': sum_size, 'v_size': v_size})
+
+
+            print(f"{v_size} {elapsed_cp}, {elapsed_python}, {elapsed_np}")
+
+    with open(path, 'w') as f:
+        json.dump(results, f)
+
+
+def plot_results(path: str):
+    import json
+    import matplotlib.pyplot as plt
+    import numpy as np
+    from collections import defaultdict
+
+    # Load the benchmark results
+    with open(path, 'r') as f:
+        results = json.load(f)
+
+    # Group data by benchmark and v_size, then calculate medians
+    data_by_benchmark = defaultdict(lambda: defaultdict(list))
+
+    for entry in results:
+        benchmark = entry['benchmark']
+        v_size = entry['v_size']
+        elapsed_time = entry['elapsed_time']
+        data_by_benchmark[benchmark][v_size].append(elapsed_time)
+
+    # Calculate medians
+    medians_by_benchmark = {}
+    for benchmark, v_sizes in data_by_benchmark.items():
+        medians_by_benchmark[benchmark] = {
+            v_size: np.median(times)
+            for v_size, times in v_sizes.items()
+        }
+
+    # Sort by v_size for plotting
+    benchmarks = sorted(medians_by_benchmark.keys())
+    v_sizes_set = sorted(set(v for benchmark_data in medians_by_benchmark.values() for v in benchmark_data.keys()))
+
+    # Create the plot
+    plt.figure(figsize=(10, 6))
+
+    for benchmark in benchmarks:
+        if benchmark != 'Python':
+            v_sizes = sorted(medians_by_benchmark[benchmark].keys())
+            elapsed_times = [medians_by_benchmark[benchmark][v] for v in v_sizes]
+            plt.plot(v_sizes, elapsed_times, '.', label=benchmark)
+
+    plt.xlabel('Vector Size (v_size)')
+    plt.ylabel('Elapsed Time (seconds)')
+    #plt.title('Benchmark Results: Elapsed Time vs Vector Size')
+    plt.legend()
+    #plt.grid(True, alpha=0.3)
+    plt.ylim(bottom=0)
+    plt.tight_layout()
+
+    # Save to PNG
+    plt.savefig(path.replace('.json', '') + '.png', dpi=300)
+    print("Plot saved")
+
+
+if __name__ == "__main__":
+    path1 = 'benchmark_results_001.json'
+    path2 = 'benchmark_results_001_sparse.json'
+
+    if 'no_simd' in sys.argv[1:]:
+        os.environ["NPY_DISABLE_CPU_FEATURES"] = CPU_SIMD_FEATURES
+        subprocess.run([sys.executable, "tests/benchmark.py"])
+    elif 'plot' in sys.argv[1:]:
+        plot_results(path1)
+        #plot_results(path2)
+    else:
+        cp_vs_python(path1)
+        plot_results(path1)
+
+        #cp_vs_python_sparse(path2)
+        #plot_results(path2)

tests/test_ast_gen.py

@@ -1,6 +1,6 @@
-from copapy import variable
+from copapy import value
 from copapy.backend import Write
-import copapy.backend as cpbe
+import copapy.backend as cpb
 
 
 def test_ast_generation():
@@ -21,8 +21,8 @@ def test_ast_generation():
     #r2 = i1 + 9
     #out = [Write(r1), Write(r2)]
 
-    c1 = variable(4)
+    c1 = value(4)
-    c2 = variable(2)
+    c2 = value(2)
     #i1 = c1 * 2
     #r1 = i1 + 7 + (c2 + 7 * 9)
     #r2 = i1 + 9
@@ -33,27 +33,27 @@ def test_ast_generation():
     print(out)
     print('-- get_edges:')
 
-    edges = list(cpbe.get_all_dag_edges(out))
+    edges = list(cpb.get_all_dag_edges(out))
     for p in edges:
         print('#', p)
 
     print('-- get_ordered_ops:')
-    ordered_ops = list(cpbe.stable_toposort(edges))
+    ordered_ops = cpb.stable_toposort(edges)
     for p in ordered_ops:
         print('#', p)
 
     print('-- get_consts:')
-    const_list = cpbe.get_const_nets(ordered_ops)
+    const_list = cpb.get_const_nets(ordered_ops)
     for p in const_list:
         print('#', p)
 
     print('-- add_read_ops:')
-    output_ops = list(cpbe.add_read_ops(ordered_ops))
+    output_ops = list(cpb.add_read_ops(ordered_ops))
     for p in output_ops:
         print('#', p)
 
     print('-- add_write_ops:')
-    extended_output_ops = list(cpbe.add_write_ops(output_ops, const_list))
+    extended_output_ops = list(cpb.add_write_ops(output_ops, const_list))
     for p in extended_output_ops:
         print('#', p)
     print('--')

tests/test_autograd.py

@@ -0,0 +1,38 @@
+from copapy import value, grad
+import copapy as cp
+import pytest
+
+
+def test_autograd():
+    # Validate against micrograd results from Andrej Karpathy
+    # https://github.com/karpathy/micrograd/blob/master/test/test_engine.py
+    a = value(-4.0)
+    b = value(2.0)
+    c = a + b
+    d = a * b + b**3
+    c += c + 1
+    c += 1 + c + (-a)
+    d += d * 2 + cp.relu(b + a)
+    d += 3 * d + cp.relu(b - a)
+    e = c - d
+    f = e**2
+    g = f / 2.0
+    g += 10.0 / f
+
+    dg = grad(g, (a, b))
+
+    tg = cp.Target()
+    tg.compile(g, dg)
+    tg.run()
+
+
+    print(f"g = {tg.read_value(g)}")
+    print(f"dg/da = {tg.read_value(dg[0])}   grad:{dg[0]}    val:{a} = {tg.read_value(a)}")
+    print(f"dg/db = {tg.read_value(dg[1])}   grad:{dg[1]}    val:{b} = {tg.read_value(b)}")
+
+    assert pytest.approx(dg[0], abs=1e-4) == 138.83381  # pyright: ignore[reportUnknownMemberType]
+    assert pytest.approx(dg[1], abs=1e-4) == 645.57725  # pyright: ignore[reportUnknownMemberType]
+
+
+if __name__ == "__main__":
+    test_autograd()

tests/test_branching_stencils.py

@@ -1,4 +1,4 @@
-from copapy import variable
+from copapy import value
 from copapy.backend import Write, compile_to_dag, add_read_command
 import copapy as cp
 import subprocess
@@ -20,7 +20,7 @@ def test_compile():
     test_vals = [0.0, -1.5, -2.0, -2.5, -3.0]
 
     # Function with no passing-on-jump as last instruction:
-    ret_test = [r for v in test_vals for r in (cp.tan(variable(v)),)]
+    ret_test = [r for v in test_vals for r in (cp.tan(value(v)),)]
 
     out = [Write(r) for r in ret_test]
 

tests/test_comp_timing.py

@@ -1,17 +1,17 @@
 import time
-from copapy import variable
 from copapy import backend
 from copapy.backend import Write, stencil_db_from_package
-import copapy.backend as cpbe
+import copapy.backend as cpb
 import copapy as cp
 import copapy._binwrite as binw
 from copapy._compiler import get_nets, get_section_layout, get_data_layout
 from copapy._compiler import patch_entry, CPConstant, get_aux_func_layout
 
+
 def test_timing_compiler():
-    t1 = cp.vector([10, 11]*128) + cp.vector(cp.variable(v) for v in range(256))
+    t1 = cp.vector([10, 11]*128) + cp.vector(cp.value(v) for v in range(256))
-    t2 = t1.sum()
+    #t2 = t1.sum()
-    t3 = cp.vector(cp.variable(1 / (v + 1)) for v in range(256))
+    t3 = cp.vector(cp.value(1 / (v + 1)) for v in range(256))
     t5 = ((t3 * t1) * 2).magnitude()
     out = [Write(t5)]
 
@@ -19,7 +19,7 @@ def test_timing_compiler():
 
     print('-- get_edges:')
     t0 = time.time()
-    edges = list(cpbe.get_all_dag_edges(out))
+    edges = list(cpb.get_all_dag_edges(out))
     t1 = time.time()
     print(f'   found {len(edges)} edges')
     #for p in edges:
@@ -28,7 +28,7 @@ def test_timing_compiler():
 
     print('-- get_ordered_ops:')
     t0 = time.time()
-    ordered_ops = list(cpbe.stable_toposort(edges))
+    ordered_ops = cpb.stable_toposort(edges)
     t1 = time.time()
     print(f'   found {len(ordered_ops)} ops')
     #for p in ordered_ops:
@@ -37,7 +37,7 @@ def test_timing_compiler():
 
     print('-- get_consts:')
     t0 = time.time()
-    const_net_list = cpbe.get_const_nets(ordered_ops)
+    const_net_list = cpb.get_const_nets(ordered_ops)
     t1 = time.time()
     #for p in const_list:
     #    print('#', p)
@@ -45,7 +45,7 @@ def test_timing_compiler():
 
     print('-- add_read_ops:')
     t0 = time.time()
-    output_ops = list(cpbe.add_read_ops(ordered_ops))
+    output_ops = list(cpb.add_read_ops(ordered_ops))
     t1 = time.time()
     #for p in output_ops:
     #    print('#', p)
@@ -53,7 +53,7 @@ def test_timing_compiler():
 
     print('-- add_write_ops:')
     t0 = time.time()
-    extended_output_ops = list(cpbe.add_write_ops(output_ops, const_net_list))
+    extended_output_ops = list(cpb.add_write_ops(output_ops, const_net_list))
     t1 = time.time()
     #for p in extended_output_ops:
     #    print('#', p)

tests/test_compile.py

@@ -49,10 +49,10 @@ def test_compile():
     #ret = function(c1, c2)
     #ret = [c1 // 3.3 + 5]
 
-    t1 = cp.vector([10, 11, 12]) + cp.vector(cp.variable(v) for v in range(3))
+    t1 = cp.vector([10, 11, 12]) + cp.vector(cp.value(v) for v in range(3))
     t2 = t1.sum()
 
-    t3 = cp.vector(cp.variable(1 / (v + 1)) for v in range(3))
+    t3 = cp.vector(cp.value(1 / (v + 1)) for v in range(3))
     t4 = ((t3 * t1) * 2).sum()
     t5 = ((t3 * t1) * 2).magnitude()
 

tests/test_compile_aarch64.py

@@ -26,7 +26,7 @@ def run_command(command: list[str]) -> str:
 def check_for_qemu() -> bool:
     command = qemu_command + ['--version']
     try:
-        result = subprocess.run(command, stdout=subprocess.PIPE, stderr=subprocess.PIPE, encoding='utf8', check=False)
+        result = subprocess.run(command, stdout=subprocess.PIPE, stderr=subprocess.PIPE, check=False)
     except Exception:
         return False
     return result.returncode == 0
@@ -43,10 +43,10 @@ def function(c1: NumLike, c2: NumLike) -> tuple[NumLike, ...]:
 
 @pytest.mark.runner
 def test_compile():
-    t1 = cp.vector([10, 11, 12]) + cp.vector(cp.variable(v) for v in range(3))
+    t1 = cp.vector([10, 11, 12]) + cp.vector(cp.value(v) for v in range(3))
     t2 = t1.sum()
 
-    t3 = cp.vector(cp.variable(1 / (v + 1)) for v in range(3))
+    t3 = cp.vector(cp.value(1 / (v + 1)) for v in range(3))
     t4 = ((t3 * t1) * 2).sum()
     t5 = ((t3 * t1) * 2).magnitude()
 

tests/test_compile_armv7.py

@@ -26,7 +26,7 @@ def run_command(command: list[str]) -> str:
 def check_for_qemu() -> bool:
     command = qemu_command + ['--version']
     try:
-        result = subprocess.run(command, stdout=subprocess.PIPE, stderr=subprocess.PIPE, encoding='utf8', check=False)
+        result = subprocess.run(command, stdout=subprocess.PIPE, stderr=subprocess.PIPE, check=False)
     except Exception:
         return False
     return result.returncode == 0
@@ -43,10 +43,10 @@ def function(c1: NumLike, c2: NumLike) -> tuple[NumLike, ...]:
 
 @pytest.mark.runner
 def test_compile():
-    t1 = cp.vector([10, 11, 12]) + cp.vector(cp.variable(v) for v in range(3))
+    t1 = cp.vector([10, 11, 12]) + cp.vector(cp.value(v) for v in range(3))
     t2 = t1.sum()
 
-    t3 = cp.vector(cp.variable(1 / (v + 1)) for v in range(3))
+    t3 = cp.vector(cp.value(1 / (v + 1)) for v in range(3))
     t4 = ((t3 * t1) * 2).sum()
     t5 = ((t3 * t1) * 2).magnitude()
 

tests/test_compile_div.py

@@ -1,4 +1,4 @@
-from copapy import variable, NumLike
+from copapy import value, NumLike
 from copapy.backend import Write, compile_to_dag
 import copapy
 import subprocess
@@ -22,7 +22,7 @@ def function(c1: NumLike) -> list[NumLike]:
 @pytest.mark.runner
 def test_compile():
 
-    c1 = variable(16)
+    c1 = value(16)
 
     ret = function(c1)
 

tests/test_compile_math.py

@@ -1,4 +1,4 @@
-from copapy import variable
+from copapy import value
 from copapy.backend import Write, compile_to_dag, add_read_command
 import copapy as cp
 import subprocess
@@ -18,7 +18,7 @@ def run_command(command: list[str]) -> str:
 def test_compile_sqrt():
     test_vals = [0.0, 0.0001, 0.1, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0, 6.28318530718, 100.0, 1000.0, 100000.0]
 
-    ret = [r for v in test_vals for r in (cp.sqrt(variable(v)),)]
+    ret = [r for v in test_vals for r in (cp.sqrt(value(v)),)]
 
 
     out = [Write(r) for r in ret]
@@ -52,7 +52,7 @@ def test_compile_sqrt():
 def test_compile_log():
     test_vals = [0.0, 0.0001, 0.1, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0, 6.28318530718, 100.0, 1000.0, 100000.0]
 
-    ret = [r for v in test_vals for r in (cp.log(variable(v)),)]
+    ret = [r for v in test_vals for r in (cp.log(value(v)),)]
 
 
     out = [Write(r) for r in ret]
@@ -86,7 +86,7 @@ def test_compile_log():
 def test_compile_sin():
     test_vals = [0.0, 0.0001, 0.1, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0, 6.28318530718, 100.0, 1000.0, 100000.0]
 
-    ret = [r for v in test_vals for r in (cp.sin(variable(v)),)]
+    ret = [r for v in test_vals for r in (cp.sin(value(v)),)]
 
 
     out = [Write(r) for r in ret]

tests/test_coparun_module.py

@@ -1,4 +1,4 @@
-from copapy import variable, Target, NumLike
+from copapy import value, Target, NumLike
 import pytest
 
 
@@ -14,7 +14,7 @@ def function(c1: NumLike) -> list[NumLike]:
 
 def test_compile():
 
-    c1 = variable(16)
+    c1 = value(16)
 
     ret = function(c1)
 

tests/test_math.py

@@ -1,4 +1,4 @@
-from copapy import variable, Target
+from copapy import value, Target
 import pytest
 import copapy as cp
 import math as ma
@@ -7,8 +7,8 @@ import warnings
 def test_fine():
     a_i = 9
     a_f = 2.5
-    c_i = variable(a_i)
+    c_i = value(a_i)
-    c_f = variable(a_f)
+    c_f = value(a_f)
     # c_b = variable(True)
 
     ret_test = (c_f ** 2,
@@ -49,7 +49,7 @@ def test_fine():
     print('* finished')
 
     for test, val2, ref, name in zip(ret_test, re2_test, ret_refe, ('^2', '**-1', 'sqrt_int', 'sqrt_float', 'sin', 'cos', 'tan')):
-        assert isinstance(test, cp.variable)
+        assert isinstance(test, cp.value)
         val = tg.read_value(test)
         print('+', val, ref, type(val), test.dtype)
         #for t in (int, float, bool):
@@ -63,7 +63,7 @@ def test_trig_precision():
     test_vals = [0.0, 0.0001, 0.1, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0, 6.28318530718, 100.0, 1000.0, 100000.0,
                  -0.0001, -0.1, -0.5, -1.0, -1.5, -2.0, -2.5, -3.0, -3.5, -4.0, -4.5, -5.0, -5.5, -6.0, -6.28318530718, -100.0, -1000.0, -100000.0]
 
-    ret_test = [r for v in test_vals for r in (cp.sin(variable(v)), cp.cos(variable(v)), cp.tan(variable(v)))]
+    ret_test = [r for v in test_vals for r in (cp.sin(value(v)), cp.cos(value(v)), cp.tan(value(v)))]
     ret_refe = [r for v in test_vals for r in (ma.sin(v), ma.cos(v), ma.tan(v))]
 
     tg = Target()
@@ -72,7 +72,7 @@ def test_trig_precision():
 
     for i, (v, test, ref) in enumerate(zip(test_vals, ret_test, ret_refe)):
         func_name = ['sin', 'cos', 'tan'][i % 3]
-        assert isinstance(test, cp.variable)
+        assert isinstance(test, cp.value)
         val = tg.read_value(test)
         print(f"+ Result of {func_name}: {val}; reference: {ref}")
         assert val == pytest.approx(ref, abs=1e-3), f"Result of {func_name} for input {test_vals[i // 3]} does not match: {val} and reference: {ref} (value: {v})"  # pyright: ignore[reportUnknownMemberType]
@@ -85,11 +85,11 @@ def test_arcus_trig_precision():
     test_vals = [0.0, 0.01, 0.1, 0.5, 0.7, 0.9, 0.95,
                  -0.01, -0.1, -0.5, -0.7, -0.9, 0.95]
 
-    ret_test = [r for v in test_vals for r in (cp.asin(variable(v)),
+    ret_test = [r for v in test_vals for r in (cp.asin(value(v)),
-                                               cp.acos(variable(v)),
+                                               cp.acos(value(v)),
-                                               cp.atan(variable(v)),
+                                               cp.atan(value(v)),
-                                               cp.atan2(variable(v), variable(3)),
+                                               cp.atan2(value(v), value(3)),
-                                               cp.atan2(variable(v), variable(-3)),)]
+                                               cp.atan2(value(v), value(-3)),)]
     ret_refe = [r for v in test_vals for r in (ma.asin(v),
                                                ma.acos(v),
                                                ma.atan(v),
@@ -102,7 +102,7 @@ def test_arcus_trig_precision():
 
     for i, (test, ref) in enumerate(zip(ret_test, ret_refe)):
         func_name = ['asin', 'acos', 'atan', 'atan2[1]', 'atan2[2]'][i % 5]
-        assert isinstance(test, cp.variable)
+        assert isinstance(test, cp.value)
         val = tg.read_value(test)
         print(f"+ Result of {func_name}: {val}; reference: {ref}")
         #assert val == pytest.approx(ref, abs=1e-5), f"Result of {func_name} for input {test_vals[i // 5]} does not match: {val} and reference: {ref}"  # pyright: ignore[reportUnknownMemberType]
@@ -113,7 +113,7 @@ def test_arcus_trig_precision():
 def test_sqrt_precision():
     test_vals = [0.0, 0.0001, 0.1, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0, 6.28318530718, 100.0, 1000.0, 100000.0]
 
-    ret_test = [r for v in test_vals for r in (cp.sqrt(variable(v)),)]
+    ret_test = [r for v in test_vals for r in (cp.sqrt(value(v)),)]
     ret_refe = [r for v in test_vals for r in (cp.sqrt(v),)]
 
     tg = Target()
@@ -122,7 +122,7 @@ def test_sqrt_precision():
 
     for i, (test, ref) in enumerate(zip(ret_test, ret_refe)):
         func_name = 'sqrt'
-        assert isinstance(test, cp.variable)
+        assert isinstance(test, cp.value)
         val = tg.read_value(test)
         print(f"+ Result of {func_name}: {val}; reference: {ref}")
         assert val == pytest.approx(ref, rel=1e-5), f"Result of {func_name} for input {test_vals[i]} does not match: {val} and reference: {ref}"  # pyright: ignore[reportUnknownMemberType]
@@ -135,8 +135,8 @@ def test_log_exp_precision():
     test_vals = [0.1, 0.5, 0.9, 0.999, 1.0, 2.5,
                  -0.1, -0.5, -0.9, -0.999, -1.0, 2.5]
 
-    ret_test = [r for v in test_vals for r in (cp.log(variable(abs(v))),
+    ret_test = [r for v in test_vals for r in (cp.log(value(abs(v))),
-                                               cp.exp(variable(v)))]
+                                               cp.exp(value(v)))]
     ret_refe = [r for v in test_vals for r in (ma.log(abs(v)),
                                                ma.exp(v))]
 
@@ -146,7 +146,7 @@ def test_log_exp_precision():
 
     for i, (test, ref) in enumerate(zip(ret_test, ret_refe)):
         func_name = ['log', 'exp'][i % 2]
-        assert isinstance(test, cp.variable)
+        assert isinstance(test, cp.value)
         val = tg.read_value(test)
         print(f"+ Result of {func_name}: {val}; reference: {ref}")
         assert val == pytest.approx(ref, abs=1e-5), f"Result of {func_name} for input {test_vals[i // 2]} does not match: {val} and reference: {ref}"  # pyright: ignore[reportUnknownMemberType]

tests/test_matrix.py

@@ -0,0 +1,273 @@
+import copapy as cp
+import pytest
+
+
+def test_matrix_init():
+    """Test basic matrix initialization"""
+    m1 = cp.matrix([[1, 2, 3], [4, 5, 6]])
+    assert m1.rows == 2
+    assert m1.cols == 3
+    assert m1[0] == (1, 2, 3)
+    assert m1[1] == (4, 5, 6)
+
+
+def test_matrix_with_variables():
+    """Test matrix initialization with variables"""
+    m1 = cp.matrix([[cp.value(1), 2], [3, cp.value(4)]])
+    assert m1.rows == 2
+    assert m1.cols == 2
+    assert isinstance(m1[0][0], cp.value)
+    assert isinstance(m1[1][1], cp.value)
+
+
+def test_matrix_addition():
+    """Test matrix addition"""
+    m1 = cp.matrix([[1, 2], [3, 4]])
+    m2 = cp.matrix([[5, 6], [7, 8]])
+    m3 = m1 + m2
+
+    assert m3[0] == (6, 8)
+    assert m3[1] == (10, 12)
+
+
+def test_matrix_scalar_addition():
+    """Test matrix addition with scalar"""
+    m1 = cp.matrix([[1, 2], [3, 4]])
+    m2 = m1 + 5
+
+    assert m2[0] == (6, 7)
+    assert m2[1] == (8, 9)
+
+
+def test_matrix_subtraction():
+    """Test matrix subtraction"""
+    m1 = cp.matrix([[5, 6], [7, 8]])
+    m2 = cp.matrix([[1, 2], [3, 4]])
+    m3 = m1 - m2
+
+    assert m3[0] == (4, 4)
+    assert m3[1] == (4, 4)
+
+
+def test_matrix_scalar_subtraction():
+    """Test matrix subtraction with scalar"""
+    m1 = cp.matrix([[5, 6], [7, 8]])
+    m2 = m1 - 2
+
+    assert m2[0] == (3, 4)
+    assert m2[1] == (5, 6)
+
+
+def test_matrix_negation():
+    """Test matrix negation"""
+    m1 = cp.matrix([[1, 2], [3, 4]])
+    m2 = -m1
+
+    assert m2[0] == (-1, -2)
+    assert m2[1] == (-3, -4)
+
+
+def test_matrix_element_wise_multiplication():
+    """Test element-wise matrix multiplication"""
+    m1 = cp.matrix([[1, 2], [3, 4]])
+    m2 = cp.matrix([[5, 6], [7, 8]])
+    m3 = m1 * m2
+
+    assert m3[0] == (5, 12)
+    assert m3[1] == (21, 32)
+
+
+def test_matrix_scalar_multiplication():
+    """Test matrix multiplication with scalar"""
+    m1 = cp.matrix([[1, 2], [3, 4]])
+    m2 = m1 * 3
+
+    assert m2[0] == (3, 6)
+    assert m2[1] == (9, 12)
+
+
+def test_matrix_element_wise_division():
+    """Test element-wise matrix division"""
+    m1 = cp.matrix([[6.0, 8.0], [12.0, 16.0]])
+    m2 = cp.matrix([[2.0, 2.0], [3.0, 4.0]])
+    m3 = m1 / m2
+
+    assert m3[0][0] == pytest.approx(3.0)  # pyright: ignore[reportUnknownMemberType]
+    assert m3[0][1] == pytest.approx(4.0)  # pyright: ignore[reportUnknownMemberType]
+    assert m3[1][0] == pytest.approx(4.0)  # pyright: ignore[reportUnknownMemberType]
+    assert m3[1][1] == pytest.approx(4.0)  # pyright: ignore[reportUnknownMemberType]
+
+
+def test_matrix_scalar_division():
+    """Test matrix division by scalar"""
+    m1 = cp.matrix([[6.0, 8.0], [12.0, 16.0]])
+    m2 = m1 / 2.0
+
+    assert m2[0] == pytest.approx((3.0, 4.0))  # pyright: ignore[reportUnknownMemberType]
+    assert m2[1] == pytest.approx((6.0, 8.0))  # pyright: ignore[reportUnknownMemberType]
+
+
+def test_matrix_vector_multiplication():
+    """Test matrix-vector multiplication using @ operator"""
+    m = cp.matrix([[1, 2, 3], [4, 5, 6]])
+    v = cp.vector([7, 8, 9])
+    result = m @ v
+
+    assert isinstance(result, cp.vector)
+    assert len(result.values) == 2
+    assert result.values[0] == 1*7 + 2*8 + 3*9
+    assert result.values[1] == 4*7 + 5*8 + 6*9
+
+
+def test_matrix_matrix_multiplication():
+    """Test matrix-matrix multiplication using @ operator"""
+    m1 = cp.matrix([[1, 2], [3, 4]])
+    m2 = cp.matrix([[5, 6], [7, 8]])
+    result = m1 @ m2
+
+    assert isinstance(result, cp.matrix)
+    assert result.rows == 2
+    assert result.cols == 2
+    assert result[0][0] == 1*5 + 2*7
+    assert result[0][1] == 1*6 + 2*8
+    assert result[1][0] == 3*5 + 4*7
+    assert result[1][1] == 3*6 + 4*8
+
+
+def test_matrix_transpose():
+    """Test matrix transpose"""
+    m = cp.matrix([[1, 2, 3], [4, 5, 6]])
+    mt = m.transpose()
+
+    assert mt.rows == 3
+    assert mt.cols == 2
+    assert mt[0] == (1, 4)
+    assert mt[1] == (2, 5)
+    assert mt[2] == (3, 6)
+
+
+def test_matrix_transpose_property():
+    """Test matrix transpose using .T property"""
+    m = cp.matrix([[1, 2, 3], [4, 5, 6]])
+    mt = m.T
+
+    assert mt.rows == 3
+    assert mt.cols == 2
+    assert mt[0] == (1, 4)
+
+
+def test_matrix_row_access():
+    """Test getting a row as a vector"""
+    m = cp.matrix([[1, 2, 3], [4, 5, 6]])
+    row0 = m.row(0)
+
+    assert isinstance(row0, cp.vector)
+    assert row0.values == (1, 2, 3)
+
+
+def test_matrix_col_access():
+    """Test getting a column as a vector"""
+    m = cp.matrix([[1, 2, 3], [4, 5, 6]])
+    col1 = m.col(1)
+
+    assert isinstance(col1, cp.vector)
+    assert col1.values == (2, 5)
+
+
+def test_matrix_trace():
+    """Test matrix trace (sum of diagonal elements)"""
+    m = cp.matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
+    trace = m.trace()
+
+    assert trace == 1 + 5 + 9
+
+
+def test_matrix_sum():
+    """Test sum of all matrix elements"""
+    m = cp.matrix([[1, 2, 3], [4, 5, 6]])
+    total = m.sum()
+
+    assert total == 1 + 2 + 3 + 4 + 5 + 6
+
+
+def test_matrix_map():
+    """Test mapping a function over matrix elements"""
+    m = cp.matrix([[1, 2], [3, 4]])
+    m_doubled = m.map(lambda x: x * 2)
+
+    assert m_doubled[0] == (2, 4)
+    assert m_doubled[1] == (6, 8)
+
+
+def test_matrix_homogenize():
+    """Test homogenizing matrix (converting to all variables)"""
+    m = cp.matrix([[1, cp.value(2)], [3, 4]])
+    m_homo = m.homogenize()
+
+    for row in m_homo:
+        for elem in row:
+            assert isinstance(elem, cp.value)
+
+
+def test_identity_matrix():
+    """Test identity matrix creation"""
+    m = cp.identity(3)
+
+    assert m.rows == 3
+    assert m.cols == 3
+    assert m[0] == (1, 0, 0)
+    assert m[1] == (0, 1, 0)
+    assert m[2] == (0, 0, 1)
+
+
+def test_zeros_matrix():
+    """Test zeros matrix creation"""
+    m = cp.zeros(2, 3)
+
+    assert m.rows == 2
+    assert m.cols == 3
+    assert m[0] == (0, 0, 0)
+    assert m[1] == (0, 0, 0)
+
+
+def test_ones_matrix():
+    """Test ones matrix creation"""
+    m = cp.ones(2, 3)
+
+    assert m.rows == 2
+    assert m.cols == 3
+    assert m[0] == (1, 1, 1)
+    assert m[1] == (1, 1, 1)
+
+
+def test_diagonal_matrix():
+    """Test diagonal matrix creation from vector"""
+    v = cp.vector([1, 2, 3])
+    m = cp.diagonal(v)
+
+    assert m.rows == 3
+    assert m.cols == 3
+    assert m[0] == (1, 0, 0)
+    assert m[1] == (0, 2, 0)
+    assert m[2] == (0, 0, 3)
+
+
+def test_matrix_with_variables_compiled():
+    """Test matrix operations with variables in compilation"""
+    m = cp.matrix([[cp.value(1), 2], [3, cp.value(4)]])
+    v = cp.vector([cp.value(5), 6])
+    result = m @ v
+
+    # result[0] = 1*5 + 2*6 = 17
+    # result[1] = 3*5 + 4*6 = 39
+
+    tg = cp.Target()
+    tg.compile(result)
+    tg.run()
+
+    assert tg.read_value(result.values[0]) == pytest.approx(17)  # pyright: ignore[reportUnknownMemberType]
+    assert tg.read_value(result.values[1]) == pytest.approx(39)  # pyright: ignore[reportUnknownMemberType]
+
+
+if __name__ == "__main__":
+    pytest.main([__file__, "-v"])

tests/test_ops.py

@@ -1,4 +1,4 @@
-from copapy import variable, Target, NumLike, iif
+from copapy import value, Target, NumLike, iif
 import pytest
 import copapy
 
@@ -42,11 +42,11 @@ def iiftests(c1: NumLike) -> list[NumLike]:
 
 
 def test_compile():
-    c_i = variable(9)
+    c_i = value(9)
-    c_f = variable(1.111)
+    c_f = value(1.111)
-    c_b = variable(True)
+    c_b = value(True)
 
-    ret_test = function1(c_i) + function1(c_f) + function2(c_i) + function2(c_f) + function3(c_i) + function4(c_i) + function5(c_b) + [variable(9) % 2] + iiftests(c_i) + iiftests(c_f)
+    ret_test = function1(c_i) + function1(c_f) + function2(c_i) + function2(c_f) + function3(c_i) + function4(c_i) + function5(c_b) + [value(9) % 2] + iiftests(c_i) + iiftests(c_f)
     ret_ref = function1(9) + function1(1.111) + function2(9) + function2(1.111) + function3(9) + function4(9) + function5(True) + [9 % 2] + iiftests(9) + iiftests(1.111)
 
     tg = Target()
@@ -57,7 +57,7 @@ def test_compile():
     print('* finished')
 
     for test, ref in zip(ret_test, ret_ref):
-        assert isinstance(test, copapy.variable)
+        assert isinstance(test, copapy.value)
         val = tg.read_value(test)
         print('+', val, ref, test.dtype)
         for t in (int, float, bool):

tests/test_ops_aarch64.py

@@ -1,4 +1,4 @@
-from copapy import NumLike, iif, variable
+from copapy import NumLike, iif, value
 from copapy.backend import Write, compile_to_dag, add_read_command
 import subprocess
 from copapy import _binwrite
@@ -41,8 +41,8 @@ def run_command(command: list[str]) -> str:
 def check_for_qemu() -> bool:
     command = qemu_command + ['--version']
     try:
-        result = subprocess.run(command, stdout=subprocess.PIPE, stderr=subprocess.PIPE, encoding='utf8', check=False)
+        result = subprocess.run(command, stdout=subprocess.PIPE, stderr=subprocess.PIPE, check=False)
-    except:
+    except Exception:
         return False
     return result.returncode == 0
 
@@ -84,11 +84,11 @@ def iiftests(c1: NumLike) -> list[NumLike]:
 
 @pytest.mark.runner
 def test_compile():
-    c_i = variable(9)
+    c_i = value(9)
-    c_f = variable(1.111)
+    c_f = value(1.111)
-    c_b = variable(True)
+    c_b = value(True)
 
-    ret_test = function1(c_i) + function1(c_f) + function2(c_i) + function2(c_f) + function3(c_i) + function4(c_i) + function5(c_b) + [variable(9) % 2] + iiftests(c_i) + iiftests(c_f) + [cp.asin(c_i/10)]
+    ret_test = function1(c_i) + function1(c_f) + function2(c_i) + function2(c_f) + function3(c_i) + function4(c_i) + function5(c_b) + [value(9) % 2] + iiftests(c_i) + iiftests(c_f) + [cp.asin(c_i/10)]
     ret_ref = function1(9) + function1(1.111) + function2(9) + function2(1.111) + function3(9) + function4(9) + function5(True) + [9 % 2] + iiftests(9) + iiftests(1.111) + [cp.asin(9/10)]
 
     out = [Write(r) for r in ret_test]
@@ -145,7 +145,7 @@ def test_compile():
     result_data = parse_results(result)
 
     for test, ref in zip(ret_test, ret_ref):
-        assert isinstance(test, variable)
+        assert isinstance(test, value)
         address = variables[test][0]
         data = result_data[address]
         if test.dtype == 'int':

tests/test_ops_armv7.py

@@ -1,4 +1,4 @@
-from copapy import NumLike, iif, variable
+from copapy import NumLike, iif, value
 from copapy.backend import Write, compile_to_dag, add_read_command
 import subprocess
 from copapy import _binwrite
@@ -17,6 +17,7 @@ if os.name == 'nt':
 else:
     qemu_command = ['qemu-arm']
 
+
 def parse_results(log_text: str) -> dict[int, bytes]:
     regex = r"^READ_DATA offs=(\d*) size=(\d*) data=(.*)$"
     matches = re.finditer(regex, log_text, re.MULTILINE)
@@ -31,6 +32,7 @@ def parse_results(log_text: str) -> dict[int, bytes]:
 
     return var_dict
 
+
 def run_command(command: list[str]) -> str:
     result = subprocess.run(command, stdout=subprocess.PIPE, stderr=subprocess.PIPE, encoding='utf8', check=False)
     assert result.returncode != 11, f"SIGSEGV (segmentation fault)\n -Error occurred: {result.stderr}\n -Output: {result.stdout}"
@@ -41,8 +43,8 @@ def run_command(command: list[str]) -> str:
 def check_for_qemu() -> bool:
     command = qemu_command + ['--version']
     try:
-        result = subprocess.run(command, stdout=subprocess.PIPE, stderr=subprocess.PIPE, encoding='utf8', check=False)
+        result = subprocess.run(command, stdout=subprocess.PIPE, stderr=subprocess.PIPE, check=False)
-    except:
+    except Exception:
         return False
     return result.returncode == 0
 
@@ -84,11 +86,11 @@ def iiftests(c1: NumLike) -> list[NumLike]:
 
 @pytest.mark.runner
 def test_compile():
-    c_i = variable(9)
+    c_i = value(9)
-    c_f = variable(1.111)
+    c_f = value(1.111)
-    c_b = variable(True)
+    c_b = value(True)
 
-    ret_test = function1(c_i) + function1(c_f) + function2(c_i) + function2(c_f) + function3(c_i) + function4(c_i) + function5(c_b) + [variable(9) % 2] + iiftests(c_i) + iiftests(c_f) + [cp.asin(c_i/10)]
+    ret_test = function1(c_i) + function1(c_f) + function2(c_i) + function2(c_f) + function3(c_i) + function4(c_i) + function5(c_b) + [value(9) % 2] + iiftests(c_i) + iiftests(c_f) + [cp.asin(c_i/10)]
     ret_ref = function1(9) + function1(1.111) + function2(9) + function2(1.111) + function3(9) + function4(9) + function5(True) + [9 % 2] + iiftests(9) + iiftests(1.111) + [cp.asin(9/10)]
 
     #ret_test = (c_i * 100 // 5, c_f * 10 // 5)
@@ -145,7 +147,7 @@ def test_compile():
     result_data = parse_results(result)
 
     for test, ref in zip(ret_test, ret_ref):
-        assert isinstance(test, variable)
+        assert isinstance(test, value)
         address = variables[test][0]
         data = result_data[address]
         if test.dtype == 'int':

tests/test_ops_x86.py

@@ -1,4 +1,4 @@
-from copapy import NumLike, iif, variable, sin
+from copapy import NumLike, iif, value
 from copapy.backend import Write, compile_to_dag, add_read_command
 import subprocess
 from copapy import _binwrite
@@ -70,16 +70,16 @@ def iiftests(c1: NumLike) -> list[NumLike]:
 
 @pytest.mark.runner
 def test_compile():
-    t1 = cp.vector([10, 11, 12]) + cp.vector(cp.variable(v) for v in range(3))
+    #t1 = cp.vector([10, 11, 12]) + cp.vector(cp.variable(v) for v in range(3))
-    t2 = t1.sum()
+    #t2 = t1.sum()
 
-    t3 = cp.vector(cp.variable(1 / (v + 1)) for v in range(3))
+    #t3 = cp.vector(cp.variable(1 / (v + 1)) for v in range(3))
-    t4 = ((t3 * t1) * 2).sum()
+    #t4 = ((t3 * t1) * 2).sum()
-    t5 = ((t3 * t1) * 2).magnitude()
+    #t5 = ((t3 * t1) * 2).magnitude()
 
-    c_i = variable(9)
+    c_i = value(9)
-    c_f = variable(1.111)
+    #c_f = variable(1.111)
-    c_b = variable(True)
+    #c_b = variable(True)
 
     #ret_test = function1(c_i) + function1(c_f) + function2(c_i) + function2(c_f) + function3(c_i) + function4(c_i) + function5(c_b) + [c_i % 2, sin(c_f)] + iiftests(c_i) + iiftests(c_f)
     #ret_ref = function1(9) + function1(1.111) + function2(9) + function2(1.111) + function3(9) + function4(9) + function5(True) + [9 % 2, sin(1.111)] + iiftests(9) + iiftests(1.111)
@@ -87,7 +87,7 @@ def test_compile():
     #ret_test = [cp.sin(c_i), cp.asin(variable(0.0))]
     #ret_ref = [cp.sin(9), cp.asin(0.0)]
 
-    ret_test: list[variable[float]] = []
+    ret_test: list[value[float]] = []
     ret_ref: list[float] = []
     #sval = variable(8.0)
     #tval = 8.0
@@ -143,7 +143,7 @@ def test_compile():
         try:
             result = run_command(command)
         except FileNotFoundError:
-            warnings.warn(f"Test skipped, executable not found.", UserWarning)
+            warnings.warn("Test skipped, executable not found.", UserWarning)
             return
 
         print('* Output from runner:\n--')
@@ -155,7 +155,7 @@ def test_compile():
         result_data = parse_results(result)
 
         for test, ref in zip(ret_test, ret_ref):
-            assert isinstance(test, variable)
+            assert isinstance(test, value)
             address = variables[test][0]
             data = result_data[address]
             if test.dtype == 'int':

tests/test_prog_flow.py

@@ -1,11 +1,11 @@
-from copapy import variable, Target, iif
+from copapy import value, Target, iif
 import pytest
 import copapy
 
 
 def test_compile():
-    c_i = variable(9)
+    c_i = value(9)
-    c_f = variable(2.5)
+    c_f = value(2.5)
     # c_b = variable(True)
 
     ret_test = (iif(c_f > 5, c_f, -1), iif(c_i > 5, c_f, 8.8), iif(c_i > 2, c_i, 1))
@@ -19,7 +19,7 @@ def test_compile():
     print('* finished')
 
     for test, ref in zip(ret_test, ret_ref):
-        assert isinstance(test, copapy.variable)
+        assert isinstance(test, copapy.value)
         val = tg.read_value(test)
         print('+', val, ref, type(val), test.dtype)
         #for t in (int, float, bool):

tests/test_readme_example.py

@@ -3,8 +3,8 @@ import pytest
 
 def test_readme_example():
     # Define variables
-    a = cp.variable(0.25)
+    a = cp.value(0.25)
-    b = cp.variable(0.87)
+    b = cp.value(0.87)
 
     # Define computations
     c = a + b * 2.0

tests/test_rev_kinematics.py

@@ -0,0 +1,48 @@
+import copapy as cp
+
+# Arm lengths
+l1, l2 = 1.8, 2.0
+
+# Target position
+target = cp.vector([0.7, 0.7])
+
+# Learning rate for iterative adjustment
+alpha = 0.1
+
+
+def forward_kinematics(theta1: cp.value[float] | float, theta2: cp.value[float] | float) -> tuple[cp.vector[float], cp.vector[float]]:
+    """Return positions of joint and end-effector."""
+    joint = cp.vector([l1 * cp.cos(theta1), l1 * cp.sin(theta1)])
+    end_effector = joint + cp.vector([l2 * cp.cos(theta1 + theta2),
+                                    l2 * cp.sin(theta1 + theta2)])
+    return joint, end_effector
+
+
+def test_two_arms():    
+    target_vec = cp.vector(target)
+    theta = cp.vector([cp.value(0.0), cp.value(0.0)])
+
+    joint = cp.vector([0.0, 0.0])
+    effector = cp.vector([0.0, 0.0])
+    error = 0.0
+
+    # Iterative IK
+    for _ in range(48):
+        joint, effector = forward_kinematics(theta[0], theta[1])
+        error = ((target_vec - effector) ** 2).sum()
+
+        grad_vec = cp.grad(error, theta)
+        theta -= alpha * grad_vec
+
+    tg = cp.Target()
+    tg.compile(error, theta, joint)
+    tg.run()
+
+    print(f"Joint angles: {tg.read_value(theta)}")
+    print(f"Joint position: {tg.read_value(joint)}")
+    print(f"End-effector position: {tg.read_value(effector)}")
+    print(f"quadratic error = {tg.read_value(error)}")
+
+
+if __name__ == '__main__':
+    test_two_arms()

tests/test_vector.py

@@ -1,42 +1,42 @@
 import math
 import copapy as cp
 import pytest
-
+from copapy import filters
 
 def test_vectors_init():
     tt1 = cp.vector(range(3)) + cp.vector([1.1, 2.2, 3.3])
-    tt2 = cp.vector([1.1, 2, cp.variable(5)]) + cp.vector(range(3))
+    tt2 = cp.vector([1.1, 2, cp.value(5)]) + cp.vector(range(3))
     tt3 = (cp.vector(range(3)) + 5.6)
-    tt4 = cp.vector([1.1, 2, 3]) + cp.vector(cp.variable(v) for v in range(3))
+    tt4 = cp.vector([1.1, 2, 3]) + cp.vector(cp.value(v) for v in range(3))
     tt5 = cp.vector([1, 2, 3]).dot(tt4)
 
     print(tt1, tt2, tt3, tt4, tt5)
 
 
 def test_compiled_vectors():
-    t1 = cp.vector([10, 11, 12]) + cp.vector(cp.variable(v) for v in range(3))
+    t1 = cp.vector([10, 11, 12]) + cp.vector(cp.value(v) for v in range(3))
     t2 = t1.sum()
 
-    t3 = cp.vector(cp.variable(1 / (v + 1)) for v in range(3))
+    t3 = cp.vector(cp.value(1 / (v + 1)) for v in range(3))
     t4 = ((t3 * t1) * 2).sum()
     t5 = ((t3 * t1) * 2).magnitude()
 
-    t6 = cp.angle_between(cp.vector([cp.variable(5.0), 0.0, 0.0]), cp.vector([5.0, 5.0, 0.0]))
+    t6 = cp.angle_between(cp.vector([cp.value(5.0), 0.0, 0.0]), cp.vector([5.0, 5.0, 0.0]))
 
     tg = cp.Target()
     tg.compile(t2, t4, t5, t6)
     tg.run()
 
-    assert isinstance(t2, cp.variable)
+    assert isinstance(t2, cp.value)
     assert tg.read_value(t2) == 10 + 11 + 12 + 0 + 1 + 2
 
-    assert isinstance(t4, cp.variable)
+    assert isinstance(t4, cp.value)
     assert tg.read_value(t4) == pytest.approx(((10/1*2) + (12/2*2) + (14/3*2)), 0.001)  # pyright: ignore[reportUnknownMemberType]
 
-    assert isinstance(t5, cp.variable)
+    assert isinstance(t5, cp.value)
     assert tg.read_value(t5) == pytest.approx(((10/1*2)**2 + (12/2*2)**2 + (14/3*2)**2) ** 0.5, 0.001)  # pyright: ignore[reportUnknownMemberType]
 
-    assert isinstance(t6, cp.variable)
+    assert isinstance(t6, cp.value)
     assert tg.read_value(t6) == pytest.approx(math.pi / 4, 0.001), tg.read_value(t6)  # pyright: ignore[reportUnknownMemberType]
 
 
@@ -95,12 +95,27 @@ def test_non_compiled_vector_operations():
     assert rotated.values[2] == pytest.approx(3.0, abs=1e-6)  # pyright: ignore[reportUnknownMemberType]
 
 
-if __name__ == "__main__":
+def test_sort_vector():
-    test_vectors_init()
+    vlist = [50, 21, 20, 10, 22, 1, 80, 70, 90]
-    test_compiled_vectors()
+    t1 = cp.vector(cp.value(v) for v in vlist)
-    test_vector_operations()
+    #t1 = cp.vector(v for v in vlist)
-    print('Finished!')
+
+    t2 = filters.median(t1)
+
+    tg = cp.Target()
+    tg.compile(t2)
+    tg.run()
+
+    result = tg.read_value(t2)
+
+    ref = sorted(vlist)[len(vlist) // 2]
+    print(sorted(vlist))
+
+    assert ref == result
+
 
 if __name__ == "__main__":
-    test_compiled_vectors()
+    #test_vectors_init()
+    #test_compiled_vectors()
+    test_sort_vector()
     print('Finished!')

tools/build.bat

@@ -46,8 +46,8 @@ wsl aarch64-linux-gnu-gcc-12 -static -Wall -Wextra -Wconversion -Wsign-conversio
 
 echo --------------arm-v6  32 bit----------------
 REM  sh ../copapy/tools/cross_compiler_unix/packobjs.sh arm-none-eabi-gcc arm-none-eabi-ld ../copapy/build/musl/musl_objects_armv6.o "-march=armv6 -mfpu=vfp -marm"
-wsl arm-none-eabi-gcc -fno-pic -ffunction-sections -march=armv6 -mfpu=vfp -marm -c build/stencils/stencils.c -O3 -o build/stencils/stencils.o
+wsl arm-none-eabi-gcc -fno-pic -ffunction-sections -march=armv6 -mfpu=vfp -mfloat-abi=hard -marm -c build/stencils/stencils.c -O3 -o build/stencils/stencils.o
-wsl arm-none-eabi-ld -r build/stencils/stencils.o build/musl/musl_objects_armv6.o -o src/copapy/obj/stencils_armv6_O3.o
+wsl arm-none-eabi-ld -r build/stencils/stencils.o build/musl/musl_objects_armv6.o $(arm-none-eabi-gcc -print-libgcc-file-name) -o src/copapy/obj/stencils_armv6_O3.o
 wsl arm-none-eabi-objdump -d -x src/copapy/obj/stencils_armv6_O3.o > build/stencils/stencils_armv6_O3.asm
 echo ------------------------------
 REM echo - Build runner
@@ -57,9 +57,11 @@ REM wsl arm-linux-gnueabihf-gcc -march=armv6 -mfpu=vfp -marm -static -Wall -Wext
 
 echo --------------arm-v7  32 bit----------------
 REM  sh ../copapy/tools/cross_compiler_unix/packobjs.sh arm-none-eabi-gcc arm-none-eabi-ld ../copapy/build/musl/musl_objects_armv7.o "-march=armv7-a -mfpu=neon-vfpv3 -marm"
-wsl arm-none-eabi-gcc -fno-pic -ffunction-sections -march=armv7-a -mfpu=neon-vfpv3 -marm -c build/stencils/stencils.c -O3 -o build/stencils/stencils.o
+wsl arm-none-eabi-gcc -fno-pic -ffunction-sections -march=armv7-a -mfpu=neon-vfpv3 -mfloat-abi=hard -marm -c build/stencils/stencils.c -O3 -o build/stencils/stencils.o
-wsl arm-none-eabi-ld -r build/stencils/stencils.o build/musl/musl_objects_armv7.o -o src/copapy/obj/stencils_armv7_O3.o
+wsl arm-none-eabi-ld -r build/stencils/stencils.o build/musl/musl_objects_armv7.o $(arm-none-eabi-gcc -print-libgcc-file-name) -o src/copapy/obj/stencils_armv7_O3.o
 wsl arm-none-eabi-objdump -d -x src/copapy/obj/stencils_armv7_O3.o > build/stencils/stencils_armv7_O3.asm
+
+
 echo ------------------------------
 echo - Build runner
 wsl arm-linux-gnueabihf-gcc -march=armv7-a -mfpu=neon-vfpv3 -marm -static -Wall -Wextra -Wconversion -Wsign-conversion -Wshadow -Wstrict-overflow -O3 -DENABLE_LOGGING src/coparun/runmem.c src/coparun/coparun.c src/coparun/mem_man.c -o build/runner/coparun-armv7

tools/clean_asm.py

@@ -28,7 +28,7 @@ def main() -> None:
         if outp_flag:
             print(line + '<br>')
 
-        if "Disassembly of section .text:" in line:
+        if "Disassembly of section .text" in line:
             outp_flag = True
 
     print('</code>')

tools/create_asm.sh

@@ -0,0 +1,33 @@
+#!/bin/bash
+
+set -e
+set -v
+
+mkdir -p build/runner
+
+cparch=$(python3 -c "import copapy; print(copapy._stencils.detect_process_arch())")
+
+# Disassemble stencil object file
+objdump -d -x src/copapy/obj/stencils_${cparch}_O3.o > build/runner/stencils.asm
+
+# Create example code disassembly
+python3 tools/make_example.py
+build/runner/coparun build/runner/test.copapy build/runner/test.copapy.bin
+
+if [ "$cparch" = 'x86_64' ]; then
+	cparch="i386:x86-64"
+elif [ "$cparch" = 'x86' ]; then
+	cparch="i386"
+elif [ "$cparch" = 'arm64' ]; then
+	cparch="aarch64"
+elif [ "$cparch" = 'armv6' ]; then
+	cparch="arm"
+elif [ "$cparch" = 'armv7' ]; then
+	cparch="arm"
+fi
+
+echo "Archtitecture: '$cparch'"
+
+objdump -D -b binary -m $cparch --adjust-vma=0x10000 build/runner/test.copapy.bin > build/runner/example.asm
+
+rm build/runner/test.copapy.bin

tools/crosscompile.sh

@@ -58,10 +58,3 @@ arm-none-eabi-ld -r $STMP /object_files/musl_objects_armv7.o $LIBGCC -o $DEST/st
 
 # RISCV 64 Bit
 #riscv64-linux-gnu-gcc-13 $FLAGS -$OPT -c $SRC -o $DEST/stencils_riscv64_$OPT.o
-
-
-# -------------- Cross compile runner --------------
-mkdir -p build/runner
-
-# Aarch64
-aarch64-linux-gnu-gcc-13 -static -O3 -DENABLE_LOGGING -o build/runner/coparun-aarch64 src/coparun/runmem.c src/coparun/coparun.c src/coparun/mem_man.c

tools/make_example.py

@@ -1,12 +1,12 @@
-from copapy import variable
+from copapy import value
 from copapy.backend import Write, compile_to_dag, stencil_db_from_package
 from copapy._binwrite import Command
 import copapy as cp
 
 
-def compile_to_x86_64() -> None:
+def compile_example(arch: str = 'native') -> None:
     """Test compilation of a simple program for x86_64."""
-    c1 = variable(9.0)
+    c1 = value(9.0)
 
     #ret = [c1 / 4, c1 / -4, c1 // 4, c1 // -4, (c1 * -1) // 4]
     ret = [c1 // 3.3 + 5]
@@ -16,65 +16,16 @@ def compile_to_x86_64() -> None:
 
     out = [Write(r) for r in ret]
 
-    sdb = stencil_db_from_package('x86_64')
+    sdb = stencil_db_from_package(arch)
     dw, _ = compile_to_dag(out, sdb)
 
     dw.write_com(Command.DUMP_CODE)
 
-    print('* Data to runner:')
+    #print('* Data to runner:')
-    dw.print()
+    #dw.print()
 
     dw.to_file('build/runner/test.copapy')
 
 
-def compile_to_x86() -> None:
-    """Test compilation of a simple program for x86 32 bit."""
-    c1 = variable(9.0)
-
-    #ret = [c1 / 4, c1 / -4, c1 // 4, c1 // -4, (c1 * -1) // 4]
-    ret = [c1 // 3.3 + 5]
-    #ret = [cp.sqrt(c1)]
-    #c2 = cp._math.get_42()
-    #ret = [c2]
-    ret = [cp.sin(variable(2.5))]
-
-    out = [Write(r) for r in ret]
-
-    sdb = stencil_db_from_package('x86')
-    dw, _ = compile_to_dag(out, sdb)
-
-    dw.write_com(Command.DUMP_CODE)
-
-    print('* Data to runner:')
-    dw.print()
-
-    dw.to_file('build/runner/test-x86.copapy')
-
-
-def compile_to_aarch64() -> None:
-    """Test compilation of a simple program for arm64."""
-    c1 = variable(9.0)
-
-    #ret = [c1 / 4, c1 / -4, c1 // 4, c1 // -4, (c1 * -1) // 4]
-    #ret = [cp.sin(c1), cp.sqrt(c1) + 5]
-    ret = [c1 // 3.3 + 5]
-    #c2 = cp._math.get_42()
-    #ret = [c2]
-
-    out = [Write(r) for r in ret]
-
-    sdb = stencil_db_from_package('arm64')
-    dw, _ = compile_to_dag(out, sdb)
-
-    dw.write_com(Command.DUMP_CODE)
-
-    print('* Data to runner:')
-    dw.print()
-
-    dw.to_file('build/runner/test-arm64.copapy')
-
-
 if __name__ == "__main__":
-    compile_to_x86_64()
+    compile_example()
-    compile_to_x86()
-    compile_to_aarch64()