diff --git a/.github/workflows/release.yml b/.github/workflows/release.yml index 7a78f71..b8e2c22 100644 --- a/.github/workflows/release.yml +++ b/.github/workflows/release.yml @@ -133,7 +133,7 @@ jobs: name: python-package-distributions path: dist/ - name: Sign the dists with Sigstore - uses: sigstore/gh-action-sigstore-python@v2.1.1 + uses: sigstore/gh-action-sigstore-python@v3.0.0 with: inputs: | ./dist/*.tar.gz diff --git a/CITATION.bib b/CITATION.bib deleted file mode 100644 index e6d3ef5..0000000 --- a/CITATION.bib +++ /dev/null @@ -1,6 +0,0 @@ -@article{georgescu2024qlbm, - title={qlbm -- A Quantum Lattice Boltzmann Software Framework}, - author={Georgescu, C{\u{a}}lin A. and Schalkers, Merel A. and M{\"o}ller, Matthias}, - journal={arXiv preprint arXiv:2411.19439}, - year={2024} -} diff --git a/CITATION.cff b/CITATION.cff new file mode 100644 index 0000000..17b415c --- /dev/null +++ b/CITATION.cff @@ -0,0 +1,35 @@ +cff-version: 1.2.0 +message: "If you use this software, please cite it as below." +authors: +- family-names: "Georgescu" + given-names: "Calin A." + orcid: "https://orcid.org/0000-0002-8102-6389" +- family-names: "Schalkers" + given-names: "Merel A." + orcid: "https://orcid.org/0000-0001-7751-9060" +- family-names: "Möller" + given-names: "Matthias" + orcid: "https://orcid.org/0000-0003-0802-945X" +title: "qlbm – A Qantum Lattice Boltzmann Software Framework" +version: 0.0.5 +doi: "10.1016/j.cpc.2025.109699" +date-released: 2025-03-24 +url: "https://doi.org/10.1016/j.cpc.2025.109699" +preferred-citation: + type: article + authors: + - family-names: "Georgescu" + given-names: "Calin A." + orcid: "https://orcid.org/0000-0002-8102-6389" + - family-names: "Schalkers" + given-names: "Merel A." + orcid: "https://orcid.org/0000-0001-7751-9060" + - family-names: "Möller" + given-names: "Matthias" + orcid: "https://orcid.org/0000-0003-0802-945X" + doi: "10.1016/j.cpc.2025.109699" + journal: "Computer Physics Communications" + title: "qlbm – A Qantum Lattice Boltzmann Software Framework" + volume: 315 + year: 2025 + issue: 109699 \ No newline at end of file diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md new file mode 100644 index 0000000..979ea1d --- /dev/null +++ b/CONTRIBUTING.md @@ -0,0 +1,31 @@ +# Contributing + +Pull requests are welcome. To make a PR, first fork the repo, make your proposed changes on the main branch, and open a PR from your fork. If it passes tests and is accepted after review, it will be merged. For discussions about feature requests and contributions, please contact `c.a.georgescu@tudelft.nl`. + +## Code style + +### Formatting and Linting + +All code should be formatted using `ruff` with default options. You should verify that your changes adhere to the `ruff` standards before submitting a PR. The exact version used in the CI pipeline can be found in the `pyproject.toml` file. Both the source code `ruff check qlbm` and the demos `ruff check demos` should adhere to these standards. + +### Type annotation + +`mypy` is used as a static type checker and all submissions must pass its checks. YYou should verify that your changes adhere to the `mypy` standards before submitting a PR. The exact version used in the CI pipeline can be found in the `pyproject.toml` file. There are some custom rules used for type checking: you can check whether your submission adheres to them by running `mypy qlbm test --config-file pyproject.toml`. +Linting + + +### Tests + +We encourage each contribution to include unit tests for the added functionality. Tests reside in the `test` directory and can be executed using `pytest test/unit`. +Please sure all tests are passing before submitting a PR, and verify that simple flow simulations work as intended. +When fixing a bug, please add a test that demonstrates the fix. + +### Documentation + +We encourage all contributions to contain an adequate and thorough documentation of the changes. All contributed classes, methods, and attributes should be documented using the established style. Where appropriate, please include code-block examples and references to the scientific literature. + +We use `sphinx` to automatically generate our documentation website. To make sure your contribution fits the standards of the code base, you can execute `make docs` in the root directory. To check the documentation website locally, you can first `cd docs` and then `make clean html` to build the website locally. You can view the updated web pages by opening `docs/build/html/index.html`. + +### Easy verification + +To verify that your changes comply with our checks, you can simply run `make check-ci` to run all checks. \ No newline at end of file diff --git a/Makefile b/Makefile index 07eece2..19b4430 100644 --- a/Makefile +++ b/Makefile @@ -13,8 +13,8 @@ check-python-version: @PYTHON_VERSION=$$($(PYTHON) --version 2>&1 | awk '{print $$2}'); \ MAJOR_VERSION=$$(echo $$PYTHON_VERSION | cut -d. -f1); \ MINOR_VERSION=$$(echo $$PYTHON_VERSION | cut -d. -f2); \ - if [ "$$MAJOR_VERSION" -ne 3 ] || [ "$$MINOR_VERSION" -lt 8 ] || [ "$$MINOR_VERSION" -gt 13 ]; then \ - echo "Python version must be between 3.8 and 3.13"; \ + if [ "$$MAJOR_VERSION" -ne 3 ] || [ "$$MINOR_VERSION" -lt 12 ] || [ "$$MINOR_VERSION" -gt 13 ]; then \ + echo "Python version must be between 3.12 and 3.13"; \ exit 1; \ fi diff --git a/README.md b/README.md index 622e710..f6b98fd 100644 --- a/README.md +++ b/README.md @@ -62,8 +62,9 @@ source qlbm-cpu-venv/bin/activate Currently, `qlbm` supports two algorithms: - - The Collisionless QLBM described in [Efficient and fail-safe quantum algorithm for the transport equation](https://doi.org/10.1016/j.jcp.2024.112816) ([arXiv:2211.14269](https://arxiv.org/abs/2211.14269)) by M.A. Schalkers and M. Möller. - - The Space-Time QLBM described in [On the importance of data encoding in quantum Boltzmann methods](https://link.springer.com/article/10.1007/s11128-023-04216-6) by M.A. Schalkers and M. Möller. + - The Quantum Transport Method (Collisionless QLBM) described in [Efficient and fail-safe quantum algorithm for the transport equation](https://doi.org/10.1016/j.jcp.2024.112816) ([arXiv:2211.14269](https://arxiv.org/abs/2211.14269)) by M.A. Schalkers and M. Möller. + - The Space-Time QLBM/QLGA described in [On the importance of data encoding in quantum Boltzmann methods](https://link.springer.com/article/10.1007/s11128-023-04216-6) by M.A. Schalkers and M. Möller and expanded in [Fully Quantum Lattice Gas Automata Building Blocks for Computational Basis State Encodings](https://arxiv.org/abs/2506.12662). + - The Linear-encoding Quantum Lattice Gas Automata (LQLGA) described in [On quantum extensions of hydrodynamic lattice gas automata](https://www.mdpi.com/2410-3896/4/2/48) by P. Love and [Fully Quantum Lattice Gas Automata Building Blocks for Computational Basis State Encodings](https://arxiv.org/abs/2506.12662). The `demos` directory contains several use cases for simulating and analyzing these algorithms. Each demo requires minimal setup once the virtual environment has been configured. Consult the `README.md` file in the `demos` directory for further details. @@ -110,4 +111,18 @@ The `demos` directory contains several use cases for simulating and analyzing th ## Citation -A preprint describing `qlbm` in detail is currently available on [arXiv](https://arxiv.org/abs/2411.19439). If you use `qlbm`, you can cite it as per the [CITATION.bib file](CITATION.bib). \ No newline at end of file +An open access peer-reviewed article describing `qlbm` is available [here](https://doi.org/10.1016/j.cpc.2025.109699). If you use `qlbm`, you can cite it as: + +``` +@article{georgescu2025qlbm, +title = {qlbm – A quantum lattice Boltzmann software framework}, +journal = {Computer Physics Communications}, +volume = {315}, +pages = {109699}, +year = {2025}, +issn = {0010-4655}, +doi = {https://doi.org/10.1016/j.cpc.2025.109699}, +url = {https://www.sciencedirect.com/science/article/pii/S0010465525002012}, +author = {C\u{{a}}lin A. Georgescu and Merel A. Schalkers and Matthias M\"{o}ller}, +} +``` \ No newline at end of file diff --git a/demos/simulation/lqlga_simulation.ipynb b/demos/simulation/lqlga_simulation.ipynb new file mode 100644 index 0000000..d92178b --- /dev/null +++ b/demos/simulation/lqlga_simulation.ipynb @@ -0,0 +1,268 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "70305ef3", + "metadata": {}, + "outputs": [], + "source": [ + "from qiskit_aer import AerSimulator\n", + "\n", + "from qlbm.components import EmptyPrimitive\n", + "from qlbm.components.lqlga.initial import LQGLAInitialConditions\n", + "from qlbm.components.lqlga.lqlga import LQLGA\n", + "from qlbm.components.lqlga.measurement import LQLGAGridVelocityMeasurement\n", + "from qlbm.infra import QiskitRunner, SimulationConfig\n", + "from qlbm.lattice.lattices.lqlga_lattice import LQLGALattice\n", + "from qlbm.tools.utils import create_directory_and_parents\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ec2f4062", + "metadata": {}, + "outputs": [], + "source": [ + "from qlbm.components.lqlga import LQGLAInitialConditions\n", + "from qlbm.lattice import LQLGALattice\n", + "\n", + "lattice = LQLGALattice(\n", + " {\n", + " \"lattice\": {\n", + " \"dim\": {\"x\": 7},\n", + " \"velocities\": \"D1Q3\",\n", + " },\n", + " \"geometry\": [{\"shape\": \"cuboid\", \"x\": [3, 5], \"boundary\": \"bounceback\"}],\n", + " },\n", + ")\n", + "initial_conditions = LQGLAInitialConditions(lattice, [(tuple([2]), (True, True, True))])\n", + "initial_conditions.draw(\"mpl\")\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "05e422b4", + "metadata": {}, + "outputs": [], + "source": [ + "from qlbm.components.lqlga import LQLGAStreamingOperator\n", + "from qlbm.lattice import LQLGALattice\n", + "\n", + "lattice = LQLGALattice(\n", + " {\n", + " \"lattice\": {\n", + " \"dim\": {\"x\": 4},\n", + " \"velocities\": \"D1Q3\",\n", + " },\n", + " \"geometry\": [],\n", + " # \"geometry\": [{\"shape\": \"cuboid\", \"x\": [3, 5], \"boundary\": \"bounceback\"}],\n", + " },\n", + ")\n", + "streaming_operator = LQLGAStreamingOperator(lattice)\n", + "streaming_operator.draw(\"mpl\")\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "020bb111", + "metadata": {}, + "outputs": [], + "source": [ + "from qlbm.components.lqlga import LQLGAReflectionOperator\n", + "from qlbm.lattice import LQLGALattice\n", + "\n", + "lattice = LQLGALattice(\n", + " {\n", + " \"lattice\": {\n", + " \"dim\": {\"x\": 7},\n", + " \"velocities\": \"D1Q3\",\n", + " },\n", + " \"geometry\": [{\"shape\": \"cuboid\", \"x\": [3, 5], \"boundary\": \"bounceback\"}],\n", + " },\n", + ")\n", + "reflection_operator = LQLGAReflectionOperator(\n", + " lattice, shapes=lattice.shapes[\"bounceback\"]\n", + ")\n", + "reflection_operator.draw(\"mpl\")\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "db7e5c99", + "metadata": {}, + "outputs": [], + "source": [ + "from qlbm.components.lqlga import LQLGA\n", + "from qlbm.lattice import LQLGALattice\n", + "\n", + "lattice = LQLGALattice(\n", + " {\n", + " \"lattice\": {\n", + " \"dim\": {\"x\": 7},\n", + " \"velocities\": \"D1Q3\",\n", + " },\n", + " \"geometry\": [{\"shape\": \"cuboid\", \"x\": [3, 5], \"boundary\": \"bounceback\"}],\n", + " },\n", + ")\n", + "\n", + "LQLGA(lattice=lattice).draw(\"mpl\")\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bb1f01f6", + "metadata": {}, + "outputs": [], + "source": [ + "from qlbm.components.lqlga import LQLGAGridVelocityMeasurement\n", + "from qlbm.lattice import LQLGALattice\n", + "\n", + "lattice = LQLGALattice(\n", + " {\n", + " \"lattice\": {\n", + " \"dim\": {\"x\": 5},\n", + " \"velocities\": \"D1Q3\",\n", + " },\n", + " \"geometry\": [],\n", + " },\n", + ")\n", + "\n", + "LQLGAGridVelocityMeasurement(lattice=lattice).draw(\"mpl\")\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "34715b68", + "metadata": {}, + "outputs": [], + "source": [ + "lattice = LQLGALattice(\n", + " {\n", + " \"lattice\": {\n", + " \"dim\": {\"x\": 4},\n", + " \"velocities\": \"D1Q3\",\n", + " },\n", + " \"geometry\": [],\n", + " # \"geometry\": [{\"shape\": \"cuboid\", \"x\": [3, 5], \"boundary\": \"bounceback\"}],\n", + " },\n", + ")\n", + "\n", + "output_dir = f\"qlbm-output/lqlga-{lattice.logger_name()}-qiskit\"\n", + "create_directory_and_parents(output_dir)\n", + "\n", + "lattice.circuit.draw(\"mpl\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b084dda3", + "metadata": {}, + "outputs": [], + "source": [ + "LQLGA(lattice).circuit.draw(\"mpl\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "86ec87a5", + "metadata": {}, + "outputs": [], + "source": [ + "cfg = SimulationConfig(\n", + " initial_conditions=LQGLAInitialConditions(\n", + " lattice, [(tuple([2]), (True, True, True))]\n", + " ),\n", + " algorithm=LQLGA(lattice),\n", + " postprocessing=EmptyPrimitive(lattice),\n", + " measurement=LQLGAGridVelocityMeasurement(lattice),\n", + " target_platform=\"QISKIT\",\n", + " compiler_platform=\"QISKIT\",\n", + " optimization_level=0,\n", + " statevector_sampling=True,\n", + " execution_backend=AerSimulator(method=\"statevector\"),\n", + " sampling_backend=AerSimulator(method=\"statevector\"),\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4d54ef63", + "metadata": {}, + "outputs": [], + "source": [ + "cfg.prepare_for_simulation()\n", + "\n", + "# Number of shots to simulate for each timestep when running the circuit\n", + "# NUM_SHOTS = 2**10\n", + "\n", + "NUM_SHOTS = 2**10\n", + "\n", + "# Number of timesteps to simulate\n", + "NUM_STEPS = 100\n", + "\n", + "# Create a runner object to simulate the circuit\n", + "runner = QiskitRunner(\n", + " cfg,\n", + " lattice,\n", + ")\n", + "\n", + "# Simulate the circuits using both snapshots\n", + "runner.run(\n", + " NUM_STEPS, # Number of time steps\n", + " NUM_SHOTS, # Number of shots per time step\n", + " output_dir,\n", + " statevector_snapshots=True,\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "76799b9e", + "metadata": {}, + "outputs": [], + "source": [ + "cfg.initial_conditions.draw(\"mpl\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "244ed487", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.10" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/demos/simulation/spacetime_simulation.ipynb b/demos/simulation/spacetime_simulation.ipynb index 3359f2d..fdf50f7 100644 --- a/demos/simulation/spacetime_simulation.ipynb +++ b/demos/simulation/spacetime_simulation.ipynb @@ -38,7 +38,7 @@ "lattice = SpaceTimeLattice(\n", " num_timesteps=1,\n", " lattice_data={\n", - " \"lattice\": {\"dim\": {\"x\": 4, \"y\": 8}, \"velocities\": {\"x\": 2, \"y\": 2}},\n", + " \"lattice\": {\"dim\": {\"x\": 4, \"y\": 8}, \"velocities\": \"D2Q4\"},\n", " \"geometry\": [],\n", " },\n", ")\n", @@ -116,6 +116,7 @@ " NUM_SHOTS, # Number of shots per time step\n", " output_dir,\n", " statevector_snapshots=True,\n", + " recompile_each_step=True, # This is necessary due to Qiskit transpiler settings. WIP.\n", ")" ] }, @@ -143,7 +144,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.12" + "version": "3.12.10" } }, "nbformat": 4, diff --git a/demos/visualization/lqlga_components.ipynb b/demos/visualization/lqlga_components.ipynb new file mode 100644 index 0000000..4786403 --- /dev/null +++ b/demos/visualization/lqlga_components.ipynb @@ -0,0 +1,105 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "ddea9d2e", + "metadata": {}, + "outputs": [], + "source": [ + "from qlbm.components.lqlga.collision import GenericLQLGACollisionOperator\n", + "from qlbm.components.lqlga.initial import LQGLAInitialConditions\n", + "from qlbm.components.lqlga.lqlga import LQLGA\n", + "from qlbm.components.lqlga.streaming import LQLGAStreamingOperator\n", + "from qlbm.lattice.lattices.lqlga_lattice import LQLGALattice" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "311e7286", + "metadata": {}, + "outputs": [], + "source": [ + "lattice = LQLGALattice(\n", + " {\n", + " \"lattice\": {\n", + " \"dim\": {\"x\": 8},\n", + " \"velocities\": \"d1q3\",\n", + " },\n", + " \"geometry\": [{\"shape\": \"cuboid\", \"x\": [3, 5], \"boundary\": \"bounceback\"}],\n", + " },\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0f787330", + "metadata": {}, + "outputs": [], + "source": [ + "LQGLAInitialConditions(lattice, [(tuple([2]), (False, False))]).draw(\"mpl\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a8967c83", + "metadata": {}, + "outputs": [], + "source": [ + "LQLGAStreamingOperator(lattice).draw(\"mpl\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c32cccb3", + "metadata": {}, + "outputs": [], + "source": [ + "GenericLQLGACollisionOperator(lattice).draw(\"mpl\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fcc13530", + "metadata": {}, + "outputs": [], + "source": [ + "LQLGA(lattice).circuit.draw(\"mpl\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d31eb2aa", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.10" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/demos/visualization/spacetime_components.ipynb b/demos/visualization/spacetime_components.ipynb index 68f64ec..79fa0a4 100644 --- a/demos/visualization/spacetime_components.ipynb +++ b/demos/visualization/spacetime_components.ipynb @@ -17,7 +17,7 @@ "source": [ "from qlbm.components.spacetime import (\n", " PointWiseSpaceTimeInitialConditions,\n", - " SpaceTimeCollisionOperator,\n", + " SpaceTimeD2Q4CollisionOperator,\n", " SpaceTimeQLBM,\n", " SpaceTimeStreamingOperator,\n", ")\n", @@ -35,7 +35,7 @@ "example_lattice = SpaceTimeLattice(\n", " 1,\n", " {\n", - " \"lattice\": {\"dim\": {\"x\": 16, \"y\": 16}, \"velocities\": {\"x\": 2, \"y\": 2}},\n", + " \"lattice\": {\"dim\": {\"x\": 16, \"y\": 16}, \"velocities\": \"d2q4\"},\n", " \"geometry\": [],\n", " },\n", ")" @@ -68,7 +68,7 @@ "metadata": {}, "outputs": [], "source": [ - "SpaceTimeCollisionOperator(example_lattice, 1).draw(\"mpl\")" + "SpaceTimeD2Q4CollisionOperator(example_lattice, 1).draw(\"mpl\")" ] }, { @@ -92,7 +92,7 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 3 (ipykernel)", + "display_name": "qiskit2-venv", "language": "python", "name": "python3" }, @@ -106,7 +106,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.12.5" + "version": "3.13.0" } }, "nbformat": 4, diff --git a/docs/source/code/base_comps.rst b/docs/source/code/comps_base.rst similarity index 94% rename from docs/source/code/base_comps.rst rename to docs/source/code/comps_base.rst index 196441e..370c4f5 100644 --- a/docs/source/code/base_comps.rst +++ b/docs/source/code/comps_base.rst @@ -20,6 +20,9 @@ Lattice Base .. autoclass:: qlbm.lattice.lattices.base.Lattice :members: +.. autoclass:: qlbm.lattice.geometry.shapes.base.Shape + :members: + .. _components_bases: Components Base diff --git a/docs/source/code/comps_collision.rst b/docs/source/code/comps_collision.rst new file mode 100644 index 0000000..1eb7a13 --- /dev/null +++ b/docs/source/code/comps_collision.rst @@ -0,0 +1,65 @@ +.. _comps_collision: + +==================================== +Collision Logic Classes +==================================== + +.. testcode:: + :hide: + + from qlbm.components import ( + EQCPermutation, + EQCRedistribution, + ) + from qlbm.lattice.spacetime.properties_base import ( + LatticeDiscretization, + LatticeDiscretizationProperties, + ) + + from qlbm.lattice.eqc import ( + EquivalenceClass, + EquivalenceClassGenerator, + ) + print("ok") + +.. testoutput:: + :hide: + + ok + + +This page contain collision-related components shared +between the :ref:`stqlbm_components` and :ref:`lqlga_components` algorithms. +Collision in these algorithms is based on the concept of equivalence classes +described in Section 4 of :cite:p:`spacetime2`, and follows a PRP (permute-redistribute-unpermute) approach. +All components of this module may be used for different variations of the Computational Basis State Encoding (CBSE) +of the velocity register. +The components of this module consist of: + +#. :ref:`collision_logic` classes that provide information about the lattice discretization and equivalence classes. +#. :ref:`collision_components` that implement the small parts of the collision operator. + +.. _collision_logic: + +Collision Logic Classes +----------------------------------- + +.. autoclass:: qlbm.lattice.spacetime.properties_base.LatticeDiscretization + +.. autoclass:: qlbm.lattice.spacetime.properties_base.LatticeDiscretizationProperties + +.. autoclass:: qlbm.lattice.eqc.EquivalenceClass + +.. autoclass:: qlbm.lattice.eqc.EquivalenceClassGenerator + +.. _collision_components: + +Collision Components +----------------------------------- + +.. autoclass:: qlbm.components.common.EQCPermutation + +.. autoclass:: qlbm.components.common.EQCRedistribution + +.. autoclass:: qlbm.components.common.EQCCollisionOperator + diff --git a/docs/source/code/cqlbm_comps.rst b/docs/source/code/comps_cqlbm.rst similarity index 100% rename from docs/source/code/cqlbm_comps.rst rename to docs/source/code/comps_cqlbm.rst diff --git a/docs/source/code/comps_lqlga.rst b/docs/source/code/comps_lqlga.rst new file mode 100644 index 0000000..da2f7b4 --- /dev/null +++ b/docs/source/code/comps_lqlga.rst @@ -0,0 +1,84 @@ +.. _lqlga_components: + +==================================== +LQLGA Circuits +==================================== + +.. testcode:: + :hide: + + from qlbm.components import ( + LQLGA, + LQGLAInitialConditions, + LQLGAGridVelocityMeasurement, + LQLGAReflectionOperator, + LQLGAStreamingOperator, + ) + from qlbm.lattice import LQLGALattice + print("ok") + +.. testoutput:: + :hide: + + ok + + +This page contains documentation about the quantum circuits that make up the +**L**\ inear **Q**\ uantum **L**\ attice **G**\ as **A**\ utomata (LQLGA). +For a more in-depth depth description of the LQLGA algorithm, +we suggest :cite:p:`lqlga1`, :cite:p:`lqlga2`, and :cite:p:`spacetime2`. +The LQGLA encodes a lattice of :math:`N_g` gridpoints with :math:`q` discrete velocities +each into :math:`N_g \cdot q` qubits. +The time-evolution of the system consists of the following steps: + +#. :ref:`lqlga_initial` prepare the starting state of the flow field. +#. :ref:`lqlga_streaming` move particles across gridpoints according to the velocity discretization. +#. :ref:`lqlga_reflection` circuits apply boundary conditions that affect particles that come in contact with solid obstacles. Reflection places those particles back in the appropriate position of the fluid domain. +#. :ref:`lqlga_collision` operators create superposed local configurations of velocity profiles. +#. :ref:`lqlga_measurement` operations extract information out of the quantum state, which can later be post-processed classically. + +This page documents the individual components that make up the CQLBM algorithm. +Subsections follow a top-down approach, where end-to-end operators are introduced first, +before being broken down into their constituent parts. + +.. _lqlga_e2e: + +End-to-end algorithms +---------------------------------- + +.. autoclass:: qlbm.components.LQLGA + +.. _lqlga_initial: + +Initial Conditions +---------------------------------- + +.. autoclass:: qlbm.components.LQGLAInitialConditions + +.. _lqlga_streaming: + +Streaming +---------------------------------- + +.. autoclass:: qlbm.components.LQLGAStreamingOperator + +.. _lqlga_reflection: + +Reflection +---------------------------------- + +.. autoclass:: qlbm.components.LQLGAReflectionOperator + +.. _lqlga_collision: + +Collision +----------------------------------- + +.. autoclass:: qlbm.components.GenericLQLGACollisionOperator + +.. _lqlga_measurement: + +Measurement +---------------------------------- + +.. autoclass:: qlbm.components.LQLGAGridVelocityMeasurement \ No newline at end of file diff --git a/docs/source/code/stqlbm_comps.rst b/docs/source/code/comps_stqbm.rst similarity index 85% rename from docs/source/code/stqlbm_comps.rst rename to docs/source/code/comps_stqbm.rst index 9851a60..a37faa3 100644 --- a/docs/source/code/stqlbm_comps.rst +++ b/docs/source/code/comps_stqbm.rst @@ -65,8 +65,13 @@ Streaming Collision ---------------------------------- -.. autoclass:: qlbm.components.spacetime.collision.SpaceTimeCollisionOperator +The collision module contains collision operators and adjacent logic classes. +The former implements the circuits that perform collision in computational basis state encodings, +while the latter contains useful abstractions that circuits build on top of. +.. autoclass:: qlbm.components.spacetime.collision.GenericSpaceTimeCollisionOperator + +.. autoclass:: qlbm.components.spacetime.collision.SpaceTimeD2Q4CollisionOperator .. _stqlbm_others: diff --git a/docs/source/code/index.rst b/docs/source/code/index.rst index ca25609..f869ba8 100644 --- a/docs/source/code/index.rst +++ b/docs/source/code/index.rst @@ -14,9 +14,11 @@ The :ref:`tools` module contains miscellaneous utilities. .. toctree:: - base_comps - cqlbm_comps - stqlbm_comps + comps_base + comps_cqlbm + comps_collision + comps_stqbm + comps_lqlga lattice infra tools diff --git a/docs/source/code/infra.rst b/docs/source/code/infra.rst index 1e4e119..a25f96f 100644 --- a/docs/source/code/infra.rst +++ b/docs/source/code/infra.rst @@ -40,7 +40,7 @@ Performance .. autoclass:: qlbm.infra.reinitialize.base.Reinitializer :members: -.. autoclass:: qlbm.infra.reinitialize.collisionless_reinitializer.CollisionlessReinitializer +.. autoclass:: qlbm.infra.reinitialize.identity_reinitializer.IdentityReinitializer :members: .. autoclass:: qlbm.infra.reinitialize.spacetime_reinitializer.SpaceTimeReinitializer diff --git a/docs/source/index.rst b/docs/source/index.rst index b80fd1d..5d0dfa1 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -6,22 +6,24 @@ quantum computers running on a heterogeneous quantum-classical high-performance Currently, the primary aim of ``qlbm`` is to accelerate and improve the research surrounding **Q**\ uantum **L**\ attice **B**\ oltzmann **M**\ ethods (QLBMs). On this website, you can find the :ref:`internal_docs` of the source code components that make up ``qlbm``. -A preprint describing `qlbm` in detail is available on `arXiv `_ :cite:p:`qlbm`. +A paper describing `qlbm` in detail is available on `here `_ :cite:p:`qlbm`. ``qlbm`` is made up of 4 main modules. -Together, the :ref:`base_components`, :ref:`cqlbm_components`, and :ref:`stqlbm_components` -module handle the parameterized creation of quantum circuits that compose QBMs. +Together, the :ref:`base_components`, :ref:`cqlbm_components`, :ref:`stqlbm_components`, and :ref:`lqlga_components` +modules handle the parameterized creation of quantum circuits that compose QBMs. The :ref:`lattice` module parses external information into quantum registers and provides uniform interfaces for underlying algorithms. The :ref:`infra` module integrates the quantum components with Tket, Qiskit, and Qulacs transpilers and runners. The :ref:`tools` module contains miscellaneous utilities. -``qlbm`` currently supports two algorithms: +``qlbm`` currently supports three algorithms: #. The **C**\ ollisionless **QLBM** (CQLBM) first described in :cite:p:`collisionless` and later expanded in :cite:p:`qmem`. -#. **S**\ pace-\ **T**\ ime **QLBM** (STQLBM) described in :cite:p:`spacetime`. +#. **S**\ pace-\ **T**\ ime **QLBM** (STQLBM) described in :cite:p:`spacetime` and :cite:p:`spacetime2`. + +#. **L**\ inear \ **Q**\ uantum **L**\ attice **G**\ as **A**\ utomata (LQLGA) described in :cite:p:`spacetime2`, :cite:p:`lqlga1`, and :cite:p:`lqlga2`. .. card:: Internal Documentation :link: internal_docs diff --git a/docs/source/refs.bib b/docs/source/refs.bib index 61a4313..ab56d73 100644 --- a/docs/source/refs.bib +++ b/docs/source/refs.bib @@ -30,12 +30,25 @@ @article{qmem } @article{qlbm, - title={qlbm -- A Quantum Lattice Boltzmann Software Framework}, - author={Georgescu, C{\u{a}}lin A. and Schalkers, Merel A. and M{\"o}ller, Matthias}, - journal={arXiv preprint arXiv:2411.19439}, - year={2024} +title = {qlbm – A quantum lattice Boltzmann software framework}, +journal = {Computer Physics Communications}, +volume = {315}, +pages = {109699}, +year = {2025}, +issn = {0010-4655}, +doi = {https://doi.org/10.1016/j.cpc.2025.109699}, +url = {https://www.sciencedirect.com/science/article/pii/S0010465525002012}, +author = {C{\u{a}}lin A. Georgescu and Merel A. Schalkers and Matthias M{\"o}ller}, } +@article{spacetime2, + title={Fully Quantum Lattice Gas Automata Building Blocks for Computational Basis State Encodings}, + author={Georgescu, C{\u{a}}lin A and Schalkers, Merel A and M{\"o}ller, Matthias}, + journal={arXiv preprint arXiv:2506.12662}, + year={2025} +} + + @article{mcswap, title={Representation of binary classification trees with binary features by quantum circuits}, author={Heese, Raoul and Bickert, Patricia and Niederle, Astrid Elisa}, @@ -45,3 +58,24 @@ @article{mcswap year={2022}, publisher={Verein zur F{\"o}rderung des Open Access Publizierens in den Quantenwissenschaften} } + +@article{lqlga2, + title={Fully quantum algorithm for mesoscale fluid simulations with application to partial differential equations}, + author={Kocherla, Sriharsha and Song, Zhixin and Chrit, Fatima Ezahra and Gard, Bryan and Dumitrescu, Eugene F and Alexeev, Alexander and Bryngelson, Spencer H}, + journal={AVS Quantum Science}, + volume={6}, + number={3}, + year={2024}, + publisher={AIP Publishing} +} + +@article{lqlga1, + title={On quantum extensions of hydrodynamic lattice gas automata}, + author={Love, Peter}, + journal={Condensed Matter}, + volume={4}, + number={2}, + pages={48}, + year={2019}, + publisher={MDPI} +} diff --git a/examples/index.rst b/examples/index.rst new file mode 100644 index 0000000..4a34a27 --- /dev/null +++ b/examples/index.rst @@ -0,0 +1,20 @@ +.. _tutorials: + +Tutorials +================================ + +We are working on integrating more current tutorials into the web documentation. +More Jupyter notebooks that demonstrate the features of ``qlbm`` are available +`in the demos directory of the GitHub repository `_. +Check the ``README`` of the directory for an explanation of the different available demos. + +Currently, the following visualization tutorials are available online: + +.. toctree:: + :caption: Tutorials + :maxdepth: 2 + + notebooks/collisionless_vis + notebooks/spacetime_vis + notebooks/geometry_vis + notebooks/flowfield_vis diff --git a/examples/notebooks/collisionless_vis.nblink b/examples/notebooks/collisionless_vis.nblink new file mode 100644 index 0000000..35cee09 --- /dev/null +++ b/examples/notebooks/collisionless_vis.nblink @@ -0,0 +1,3 @@ +{ + "path": "../../../../demos/visualization/collisionless_components.ipynb" +} \ No newline at end of file diff --git a/examples/notebooks/flowfield_vis.nblink b/examples/notebooks/flowfield_vis.nblink new file mode 100644 index 0000000..6e10e1e --- /dev/null +++ b/examples/notebooks/flowfield_vis.nblink @@ -0,0 +1,3 @@ +{ + "path": "../../../../demos/visualization/flowfields.ipynb" +} \ No newline at end of file diff --git a/examples/notebooks/geometry_vis.nblink b/examples/notebooks/geometry_vis.nblink new file mode 100644 index 0000000..0da5790 --- /dev/null +++ b/examples/notebooks/geometry_vis.nblink @@ -0,0 +1,3 @@ +{ + "path": "../../../../demos/visualization/geometry.ipynb" +} \ No newline at end of file diff --git a/examples/notebooks/lqlga_vis.nblink b/examples/notebooks/lqlga_vis.nblink new file mode 100644 index 0000000..a8ca849 --- /dev/null +++ b/examples/notebooks/lqlga_vis.nblink @@ -0,0 +1,3 @@ +{ + "path": "../../../../demos/visualization/lqlga_components.ipynb" +} \ No newline at end of file diff --git a/examples/notebooks/spacetime_vis.nblink b/examples/notebooks/spacetime_vis.nblink new file mode 100644 index 0000000..68905b1 --- /dev/null +++ b/examples/notebooks/spacetime_vis.nblink @@ -0,0 +1,3 @@ +{ + "path": "../../../../demos/visualization/spacetime_components.ipynb" +} \ No newline at end of file diff --git a/pyproject.toml b/pyproject.toml index 4aba192..9f8de33 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -29,7 +29,7 @@ dependencies = [ "pytket>=1.29.2", "pytket-qiskit", "pytket-qulacs>=0.33", - "qiskit>=1.2, <1.3", + "qiskit>=2.0", "qiskit_qasm3_import>=0.4.2", "qiskit-qulacs>=0.1.0", "tqdm>=4.66", @@ -45,7 +45,7 @@ readme = "README.md" [project.optional-dependencies] dev = [ - "mypy", + "mypy==1.16.1", "pytest", "pytest-cov", "coverage", @@ -104,6 +104,7 @@ module = [ "qiskit.qasm3", "qiskit.quantum_info", "qiskit.result", + "qiskit.synthesis", "qiskit.transpiler.preset_passmanagers", "qiskit_aer", "qiskit_aer.backends.aerbackend", diff --git a/qlbm/__init__.py b/qlbm/__init__.py index 25fa576..f2d7a69 100644 --- a/qlbm/__init__.py +++ b/qlbm/__init__.py @@ -9,7 +9,7 @@ CollisionlessStreamingOperator, SpecularReflectionOperator, ) -from .infra import CircuitCompiler, CollisionlessResult, QiskitRunner, QulacsRunner +from .infra import CircuitCompiler, CollisionlessResult, QiskitRunner from .lattice import CollisionlessLattice, Lattice from .lattice.lattices.spacetime_lattice import SpaceTimeLattice diff --git a/qlbm/components/__init__.py b/qlbm/components/__init__.py index 56a2846..3c79cd7 100644 --- a/qlbm/components/__init__.py +++ b/qlbm/components/__init__.py @@ -25,11 +25,23 @@ ) from .common import ( EmptyPrimitive, + EQCCollisionOperator, + EQCPermutation, + EQCRedistribution, +) +from .lqlga import ( + LQLGA, + GenericLQLGACollisionOperator, + LQGLAInitialConditions, + LQLGAGridVelocityMeasurement, + LQLGAReflectionOperator, + LQLGAStreamingOperator, ) __all__ = [ "QuantumComponent", "LBMPrimitive", + "GenericLQLGACollisionOperator", "LBMOperator", "CQLBMOperator", "SpaceTimeOperator", @@ -48,4 +60,13 @@ "SpecularReflectionOperator", "BounceBackReflectionOperator", "CQLBM", + "GenericLQLGACollisionOperator", + "LQGLAInitialConditions", + "LQLGA", + "LQLGAGridVelocityMeasurement", + "LQLGAReflectionOperator", + "LQLGAStreamingOperator", + "EQCCollisionOperator", + "EQCPermutation", + "EQCRedistribution", ] diff --git a/qlbm/components/base.py b/qlbm/components/base.py index 4d70670..fc3cf73 100644 --- a/qlbm/components/base.py +++ b/qlbm/components/base.py @@ -10,6 +10,7 @@ from typing_extensions import override from qlbm.lattice import CollisionlessLattice, Lattice +from qlbm.lattice.lattices.lqlga_lattice import LQLGALattice from qlbm.lattice.lattices.spacetime_lattice import SpaceTimeLattice @@ -237,6 +238,33 @@ def __init__( self.lattice = lattice +class LQLGAOperator(LBMOperator): + """ + Specialization of the :class:`.LBMOperator` operator class for the LQLGA algorithm. + + Specializaitons of this class infer their properties + based on a :class:`.LQLGALattice`. + + ========================= ====================================================================== + Attribute Summary + ========================= ====================================================================== + :attr:`circuit` The :class:`.qiskit.QuantumCircuit` of the operator. + :attr:`lattice` The :class:`.LQLGALattice` based on which the properties of the operator are inferred. + :attr:`logger` The performance logger, by default ``getLogger("qlbm")`` + ========================= ====================================================================== + """ + + lattice: LQLGALattice + + def __init__( + self, + lattice: LQLGALattice, + logger: Logger = getLogger("qlbm"), + ) -> None: + super().__init__(lattice, logger) + self.lattice = lattice + + class LBMAlgorithm(QuantumComponent): """ Base class for all end-to-end Quantum Boltzmann Methods. diff --git a/qlbm/components/collisionless/bounceback_reflection.py b/qlbm/components/collisionless/bounceback_reflection.py index ffe9254..7ca3571 100644 --- a/qlbm/components/collisionless/bounceback_reflection.py +++ b/qlbm/components/collisionless/bounceback_reflection.py @@ -5,7 +5,7 @@ from typing import List from qiskit import QuantumCircuit -from qiskit.circuit.library import MCMT, XGate +from qiskit.circuit.library import MCMTGate, XGate from typing_extensions import override from qlbm.components.base import CQLBMOperator, LBMPrimitive @@ -61,7 +61,7 @@ class BounceBackWallComparator(LBMPrimitive): # Comparing on the indices of the inside x-wall on the lower-bound of the obstacle BounceBackWallComparator( - lattice=lattice, wall=lattice.block_list[0].walls_inside[0][0] + lattice=lattice, wall=lattice.shape_list[0].walls_inside[0][0] ).draw("mpl") """ @@ -171,7 +171,7 @@ class BounceBackReflectionOperator(CQLBMOperator): } ) - BounceBackReflectionOperator(lattice=lattice, blocks=lattice.block_list) + BounceBackReflectionOperator(lattice=lattice, blocks=lattice.shape_list) """ def __init__( @@ -321,7 +321,7 @@ def reflect_wall( target_qubits = self.lattice.ancillae_obstacle_index(0) circuit.compose( - MCMT( + MCMTGate( XGate(), len(control_qubits), len(target_qubits), @@ -397,7 +397,7 @@ def reset_edge_state( circuit.x(self.lattice.velocity_dir_index(dim)) circuit.compose( - MCMT( + MCMTGate( XGate(), len(control_qubits), len(target_qubits), @@ -467,7 +467,7 @@ def reset_point_state( target_qubits = self.lattice.ancillae_obstacle_index(0) circuit.compose( - MCMT( + MCMTGate( XGate(), len(control_qubits), len(target_qubits), @@ -509,7 +509,7 @@ def flip_and_stream( target_qubits = self.lattice.velocity_dir_index() circuit.compose( - MCMT( + MCMTGate( XGate(), len(control_qubits), len(target_qubits), @@ -529,4 +529,4 @@ def flip_and_stream( @override def __str__(self) -> str: - return f"[Operator BounceBackReflectionOperator against block {self.lattice.blocks}]" + return f"[Operator BounceBackReflectionOperator against shapes {self.lattice.shapes}]" diff --git a/qlbm/components/collisionless/cqlbm.py b/qlbm/components/collisionless/cqlbm.py index 579fad8..b9b54b5 100644 --- a/qlbm/components/collisionless/cqlbm.py +++ b/qlbm/components/collisionless/cqlbm.py @@ -8,6 +8,8 @@ from qlbm.components.base import LBMAlgorithm from qlbm.lattice import CollisionlessLattice +from qlbm.lattice.geometry.shapes.block import Block +from qlbm.tools.exceptions import LatticeException from qlbm.tools.utils import get_time_series from .bounceback_reflection import BounceBackReflectionOperator @@ -68,21 +70,36 @@ def create_circuit(self): ).circuit, inplace=True, ) - if self.lattice.blocks["specular"]: + if self.lattice.shapes["specular"]: + if not all( + isinstance(shape, Block) + for shape in self.lattice.shapes["specular"] + ): + raise LatticeException( + "All shapes with the 'specular' boundary condition must be of type Block for the CQLBM algorithm. " + ) circuit.compose( SpecularReflectionOperator( self.lattice, - self.lattice.blocks["specular"], + self.lattice.shapes["specular"], # type: ignore logger=self.logger, ).circuit, inplace=True, ) - if self.lattice.blocks["bounceback"]: + if self.lattice.shapes["bounceback"]: + if self.lattice.shapes["specular"]: + if not all( + isinstance(shape, Block) + for shape in self.lattice.shapes["specular"] + ): + raise LatticeException( + "All shapes with the 'bounceback' boundary condition must be of type Block for the CQLBM algorithm. " + ) circuit.compose( BounceBackReflectionOperator( self.lattice, - self.lattice.blocks["bounceback"], + self.lattice.shapes["bounceback"], # type: ignore logger=self.logger, ).circuit, inplace=True, diff --git a/qlbm/components/collisionless/primitives.py b/qlbm/components/collisionless/primitives.py index 81f1b61..d9147ee 100644 --- a/qlbm/components/collisionless/primitives.py +++ b/qlbm/components/collisionless/primitives.py @@ -6,7 +6,7 @@ from typing import List from qiskit import ClassicalRegister, QuantumCircuit -from qiskit.circuit.library import QFT +from qiskit.synthesis import synth_qft_full as QFT from typing_extensions import override from qlbm.components.base import LBMPrimitive @@ -49,7 +49,7 @@ class GridMeasurement(LBMPrimitive): } }, "geometry": [ - { + { "shape": "cuboid", "x": [5, 6], "y": [1, 2], @@ -488,7 +488,7 @@ class EdgeComparator(LBMPrimitive): ) # Draw the edge comparator circuit for one specific corner edge - EdgeComparator(lattice, lattice.block_list[0].corner_edges_3d[0]).draw("mpl") + EdgeComparator(lattice, lattice.shape_list[0].corner_edges_3d[0]).draw("mpl") """ def __init__( diff --git a/qlbm/components/collisionless/specular_reflection.py b/qlbm/components/collisionless/specular_reflection.py index d652baa..b0bfb07 100644 --- a/qlbm/components/collisionless/specular_reflection.py +++ b/qlbm/components/collisionless/specular_reflection.py @@ -5,7 +5,7 @@ from typing import List from qiskit import QuantumCircuit -from qiskit.circuit.library import MCMT, XGate +from qiskit.circuit.library import MCMTGate, XGate from typing_extensions import override from qlbm.components.base import CQLBMOperator, LBMPrimitive @@ -62,7 +62,7 @@ class SpecularWallComparator(LBMPrimitive): # Comparing on the indices of the inside x-wall on the lower-bound of the obstacle SpecularWallComparator( - lattice=lattice, wall=lattice.block_list[0].walls_inside[0][0] + lattice=lattice, wall=lattice.shape_list[0].walls_inside[0][0] ).draw("mpl") """ @@ -166,7 +166,7 @@ class SpecularReflectionOperator(CQLBMOperator): } ) - SpecularReflectionOperator(lattice=lattice, blocks=lattice.block_list) + SpecularReflectionOperator(lattice=lattice, blocks=lattice.shape_list) """ def __init__( @@ -307,7 +307,7 @@ def reflect_wall( target_qubits = self.lattice.ancillae_obstacle_index(wall.dim) circuit.compose( - MCMT( + MCMTGate( XGate(), len(control_qubits), len(target_qubits), @@ -388,7 +388,7 @@ def reset_edge_state( circuit.x(self.lattice.velocity_dir_index(dim)) circuit.compose( - MCMT( + MCMTGate( XGate(), len(control_qubits), len(target_qubits), @@ -463,7 +463,7 @@ def reset_point_state( ) circuit.compose( - MCMT( + MCMTGate( XGate(), len(control_qubits), len(target_qubits), @@ -519,5 +519,5 @@ def flip_and_stream( @override def __str__(self) -> str: return ( - f"[Operator SpecularReflectionOperator against block {self.lattice.blocks}]" + f"[Operator SpecularReflectionOperator against shapes {self.lattice.shapes}]" ) diff --git a/qlbm/components/collisionless/streaming.py b/qlbm/components/collisionless/streaming.py index 8406d4e..8e58210 100644 --- a/qlbm/components/collisionless/streaming.py +++ b/qlbm/components/collisionless/streaming.py @@ -7,7 +7,8 @@ import numpy as np from qiskit import QuantumCircuit -from qiskit.circuit.library import MCMT, QFT, XGate +from qiskit.circuit.library import MCXGate +from qiskit.synthesis import synth_qft_full as QFT from typing_extensions import override from qlbm.components.base import CQLBMOperator, LBMPrimitive @@ -93,7 +94,7 @@ def create_circuit(self) -> QuantumCircuit: circuit.x(velocity_qubit_indices_to_invert) circuit.compose( - MCMT(XGate(), num_velocity_qubits, 1), + MCXGate(num_velocity_qubits), qubits=self.lattice.velocity_index(self.dim) + self.lattice.ancillae_velocity_index(self.dim), inplace=True, diff --git a/qlbm/components/common/__init__.py b/qlbm/components/common/__init__.py index 5ca19c9..75ca482 100644 --- a/qlbm/components/common/__init__.py +++ b/qlbm/components/common/__init__.py @@ -1,9 +1,13 @@ """Common primitives used for multiple encodings.""" +from .cbse_collision import EQCCollisionOperator, EQCPermutation, EQCRedistribution from .primitives import ( EmptyPrimitive, ) __all__ = [ "EmptyPrimitive", + "EQCCollisionOperator", + "EQCPermutation", + "EQCRedistribution", ] diff --git a/qlbm/components/common/cbse_collision/__init__.py b/qlbm/components/common/cbse_collision/__init__.py new file mode 100644 index 0000000..d835c7c --- /dev/null +++ b/qlbm/components/common/cbse_collision/__init__.py @@ -0,0 +1,11 @@ +"""Quantum circuits used for computational basis state encodings.""" + +from .cbse_collision import EQCCollisionOperator +from .cbse_permutation import EQCPermutation +from .cbse_redistribution import EQCRedistribution + +__all__ = [ + "EQCCollisionOperator", + "EQCPermutation", + "EQCRedistribution", +] diff --git a/qlbm/components/common/cbse_collision/cbse_collision.py b/qlbm/components/common/cbse_collision/cbse_collision.py new file mode 100644 index 0000000..0235660 --- /dev/null +++ b/qlbm/components/common/cbse_collision/cbse_collision.py @@ -0,0 +1,147 @@ +"""Collision operators for the :class:`.SpaceTimeQLBM` algorithm :cite:`spacetime`.""" + +from time import perf_counter_ns + +from qiskit import QuantumCircuit +from typing_extensions import override + +from qlbm.components.base import LBMPrimitive +from qlbm.components.common.cbse_collision.cbse_permutation import ( + EQCPermutation, +) +from qlbm.components.common.cbse_collision.cbse_redistribution import ( + EQCRedistribution, +) +from qlbm.lattice.eqc.eqc import EquivalenceClass +from qlbm.lattice.eqc.eqc_generator import ( + EquivalenceClassGenerator, +) +from qlbm.lattice.spacetime.properties_base import ( + LatticeDiscretization, + LatticeDiscretizationProperties, +) + + +class EQCCollisionOperator(LBMPrimitive): + """ + Collision operator based on the equivalence class abstraction described in section 5 of :cite:`spacetime2`. + + Consists of a permutation, redistribution, and inverse permutation of the velocity qubits. + This operator is designed to be applied to a single velocity register, which can be repeated depending on the encoding. + Used in the :class:`.GenericSpaceTimeCollisionOperator` and :class:`.GenericLQLGACollisionOperator`. + + ========================= ====================================================================== + Attribute Summary + ========================= ====================================================================== + :attr:`discretization` The discretization for which this collision operator is defined. + :attr:`num_velocities` The number of velocities in the discretization. + ========================= ====================================================================== + + Simple D2Q4 example usage: + + .. plot:: + :include-source: + + from qlbm.components.common import EQCCollisionOperator + from qlbm.lattice import LatticeDiscretization + + # Select a discretization and draw its circuit + EQCCollisionOperator( + LatticeDiscretization.D2Q4 + ).draw("mpl") + + More complex D3Q6 example usage: + + .. plot:: + :include-source: + + from qlbm.components.common import EQCCollisionOperator + from qlbm.lattice import LatticeDiscretization + + # Select a discretization and draw its circuit + EQCCollisionOperator( + LatticeDiscretization.D3Q6 + ).draw("mpl") + """ + + discretization: LatticeDiscretization + """ + The discretization of the lattice for which this collision operator is defined. + """ + + num_velocities: int + """ + The number of velocities in the discretization.""" + + def __init__( + self, + discretization: LatticeDiscretization, + ) -> None: + super().__init__() + self.discretization = discretization + self.num_velocities = LatticeDiscretizationProperties.get_num_velocities( + self.discretization + ) + + self.logger.info(f"Creating circuit {str(self)}...") + circuit_creation_start_time = perf_counter_ns() + self.circuit = self.create_circuit() + self.logger.info( + f"Creating circuit {str(self)} took {perf_counter_ns() - circuit_creation_start_time} (ns)" + ) + + @override + def create_circuit(self) -> QuantumCircuit: + """ + Applies the collision operator of all equivalence classes onto one velocity register. + + Returns + ------- + QuantumCircuit + The circuit performing the complete collision operator. + """ + circuit = QuantumCircuit(self.num_velocities) + + for eqc in EquivalenceClassGenerator( + self.discretization + ).generate_equivalence_classes(): + circuit.compose(self.create_circuit_one_eqc(eqc), inplace=True) + + return circuit + + def create_circuit_one_eqc( + self, equivalence_class: EquivalenceClass + ) -> QuantumCircuit: + """ + Creates the PRP-based collision operator for one equivalence class. + + Parameters + ---------- + equivalence_class : EquivalenceClass + The equivalence class to collide. + + Returns + ------- + QuantumCircuit + The circuit performing the collision. + """ + circuit = QuantumCircuit(self.num_velocities) + circuit.compose( + EQCPermutation(equivalence_class, logger=self.logger).circuit, + inplace=True, + ) + + circuit.compose( + EQCRedistribution(equivalence_class, logger=self.logger).circuit, + inplace=True, + ) + circuit.compose( + EQCPermutation(equivalence_class, inverse=True, logger=self.logger).circuit, + inplace=True, + ) + + return circuit + + @override + def __str__(self) -> str: + return f"[EQCCollisionOperator for equi {self.discretization}]" diff --git a/qlbm/components/common/cbse_collision/cbse_permutation.py b/qlbm/components/common/cbse_collision/cbse_permutation.py new file mode 100644 index 0000000..c50ee53 --- /dev/null +++ b/qlbm/components/common/cbse_collision/cbse_permutation.py @@ -0,0 +1,176 @@ +"""Permutations of states belonging to equivalence classes, based on the computational basis state encoding.""" + +from logging import Logger, getLogger +from time import perf_counter_ns +from typing import override + +from qiskit import QuantumCircuit + +from qlbm.components.base import LBMPrimitive +from qlbm.lattice.eqc.eqc import EquivalenceClass +from qlbm.lattice.spacetime.properties_base import LatticeDiscretization +from qlbm.tools.exceptions import CircuitException + + +class EQCPermutation(LBMPrimitive): + """Applies a permutation to the velocity qubits of an equivalence Sclass in the CBSE encoding. + + This is used as part of the PRP collision operator described in section 5 of :cite:`spacetime2`. + Utilized in the :class:`.EQCCollisionOperator`. + + ============================== ================================================================= + Attribute Summary + ============================== ================================================================= + :attr:`equivalence_class` The equivalence class of the operator. + :attr:`inverse` Whether to apply the inverse permutation. + ============================== ================================================================= + + Example usage: + + .. plot:: + :include-source: + + from qlbm.components.common import EQCPermutation + from qlbm.lattice import LatticeDiscretization + from qlbm.lattice.eqc import EquivalenceClassGenerator + + # Generate some equivalence classes + eqcs = EquivalenceClassGenerator( + LatticeDiscretization.D3Q6 + ).generate_equivalence_classes() + + # Select one at random and draw its circuit + EQCPermutation(eqcs.pop(), inverse=False).circuit.draw("mpl") + + """ + + equivalence_class: EquivalenceClass + """ + The equivalence class for which the permutation is defined. + """ + + inverse: bool + """ + Whether to apply the inverse permutation. + """ + + def __init__( + self, + equivalence_class: EquivalenceClass, + inverse: bool = False, + logger: Logger = getLogger("qlbm"), + ): + super().__init__(logger) + self.equivalence_class = equivalence_class + self.inverse = inverse + + self.logger.info(f"Creating circuit {str(self)}...") + circuit_creation_start_time = perf_counter_ns() + self.circuit = self.create_circuit() + self.logger.info( + f"Creating circuit {str(self)} took {perf_counter_ns() - circuit_creation_start_time} (ns)" + ) + + @override + def create_circuit(self): + if self.equivalence_class.discretization == LatticeDiscretization.D1Q3: + return self.__create_circuit_d1q3() + elif self.equivalence_class.discretization == LatticeDiscretization.D2Q4: + return self.__create_circuit_d2q4() + elif self.equivalence_class.discretization == LatticeDiscretization.D3Q6: + return self.__create_circuit_d3q6() + else: + raise CircuitException( + f"Collision not yet supported for discretization {self.equivalence_class.discretization}." + ) + + def __create_circuit_d1q3(self): + circuit = QuantumCircuit(3) + + if not self.inverse: + circuit.cx(0, 1) + circuit.cx(0, 2) + else: + circuit.cx(0, 2) + circuit.cx(0, 1) + + return circuit + + def __create_circuit_d2q4(self): + circuit = QuantumCircuit(4) + + if not self.inverse: + circuit.cx(1, 2) + circuit.cx(0, 1) + circuit.cx(0, 3) + else: + circuit.cx(0, 3) + circuit.cx(0, 1) + circuit.cx(1, 2) + + return circuit + + def __create_circuit_d3q6(self): + circuit = QuantumCircuit(6) + + match self.equivalence_class.id(): + case (2, [0, 0, 0]): + circuit.cx(2, 3) + circuit.cx(5, 4) + + circuit.cx(1, 2) + circuit.cx(1, 3) + circuit.cx(1, 5) + + circuit.cx(0, 2) + circuit.cx(0, 4) + circuit.cx(0, 5) + case (4, [0, 0, 0]): + circuit.ccx(4, 5, 3) + circuit.ccx(0, 2, 4) + circuit.ccx(0, 1, 2) + circuit.ccx(0, 1, 5) + circuit.x(1) + circuit.cx(1, 0) + case (3, [1, 0, 0]): + circuit.cx(0, 3) + circuit.cx(1, 2) + circuit.ccx(0, 5, 4) + circuit.cx(1, 5) + circuit.swap(0, 1) + case (3, [-1, 0, 0]): + circuit.cx(1, 2) + circuit.cx(5, 4) + circuit.cx(1, 5) + circuit.x(0) + circuit.swap(0, 1) + case (3, [0, 1, 0]): + circuit.cx(0, 2) + circuit.cx(5, 3) + circuit.cx(0, 5) + circuit.x(4) + case (3, [0, -1, 0]): + circuit.cx(5, 3) + circuit.cx(0, 2) + circuit.cx(0, 5) + circuit.x(1) + case (3, [0, 0, 1]): + circuit.cx(1, 3) + circuit.cx(0, 1) + circuit.cx(0, 4) + circuit.x(5) + case (3, [0, 0, -1]): + circuit.cx(0, 1) + circuit.cx(4, 3) + circuit.cx(0, 4) + circuit.x(2) + case _: + raise CircuitException( + f"Collision not yet supported for discretization {self.equivalence_class.discretization} and equivalence class {self.equivalence_class.id()}." + ) + + return circuit if not self.inverse else circuit.inverse() + + @override + def __str__(self) -> str: + return f"[EQCPermutation for eqc {self.equivalence_class}]" diff --git a/qlbm/components/common/cbse_collision/cbse_redistribution.py b/qlbm/components/common/cbse_collision/cbse_redistribution.py new file mode 100644 index 0000000..25d66fa --- /dev/null +++ b/qlbm/components/common/cbse_collision/cbse_redistribution.py @@ -0,0 +1,141 @@ +"""Permutations of states belonging to equivalence classes, based on the computational basis state encoding.""" + +from logging import Logger, getLogger +from time import perf_counter_ns +from typing import override + +import numpy as np +from qiskit import QuantumCircuit +from qiskit.quantum_info import Operator + +from qlbm.components.base import LBMPrimitive +from qlbm.lattice.eqc.eqc import EquivalenceClass +from qlbm.lattice.spacetime.properties_base import LatticeDiscretizationProperties +from qlbm.tools.utils import is_two_pow + + +class EQCRedistribution(LBMPrimitive): + """ + Redistribution operator for equivalence classes in the CBSE encoding. + + The operator is mathematically described in section 4 of :cite:`spacetime2`. + Redistribution is applied before and after permutations, and consists of a controlled + unitary operator composed of a discrete Fourier transform (DFT)-block matrix. + + + ========================= ====================================================================== + Attribute Summary + ========================= ====================================================================== + :attr:`equivalence_class` The equivalence class of the operator. + :attr:`decompose_block` Whether to decompose the DFT block into a circuit. + ========================= ====================================================================== + + Example usage: + + .. plot:: + :include-source: + + from qlbm.components.common import EQCRedistribution + from qlbm.lattice import LatticeDiscretization + from qlbm.lattice.eqc import EquivalenceClassGenerator + + # Generate some equivalence classes + eqcs = EquivalenceClassGenerator( + LatticeDiscretization.D3Q6 + ).generate_equivalence_classes() + + # Select one at random and draw its circuit in the schematic form + EQCRedistribution(eqcs.pop(), decompose_block=False).circuit.draw("mpl") + + The `decompose_block` parameter can be set to ``True`` to decompose the DFT block into a circuit: + + .. plot:: + :include-source: + + from qlbm.components.common import EQCRedistribution + from qlbm.lattice import LatticeDiscretization + from qlbm.lattice.eqc import EquivalenceClassGenerator + + # Generate some equivalence classes + eqcs = EquivalenceClassGenerator( + LatticeDiscretization.D3Q6 + ).generate_equivalence_classes() + + # Select one at random and draw its decomposed circuit + EQCRedistribution(eqcs.pop(), decompose_block=True).circuit.draw("mpl") + + """ + + equivalence_class: EquivalenceClass + """ + The equivalence class for which the redistribution is defined. + """ + + decompose_block: bool + """ + Whether to decompose the DFT block into a circuit. + If set to ``False``, the block is returned as a matrix. Otherwise, it is decomposed into a circuit. + Defaults to ``True``. + """ + + def __init__( + self, + equivalence_class: EquivalenceClass, + decompose_block: bool = True, + logger: Logger = getLogger("qlbm"), + ): + super().__init__(logger) + self.equivalence_class = equivalence_class + self.decompose_block = decompose_block + + self.logger.info(f"Creating circuit {str(self)}...") + circuit_creation_start_time = perf_counter_ns() + self.circuit = self.create_circuit() + self.logger.info( + f"Creating circuit {str(self)} took {perf_counter_ns() - circuit_creation_start_time} (ns)" + ) + + @override + def create_circuit(self): + nv = LatticeDiscretizationProperties.get_num_velocities( + self.equivalence_class.discretization + ) + circuit = QuantumCircuit(nv) + n = self.equivalence_class.size() + nq = np.ceil(np.log2(n)).astype(int) + + redistribution_circuit = QuantumCircuit(nq) + if is_two_pow(n): + redistribution_circuit.ry(np.pi / 2, list(range(nq)), label="RY(π/2)") + else: + QFT = np.array( + [ + [np.exp(2j * np.pi * i * j / n) / np.sqrt(n) for j in range(n)] + for i in range(n) + ] + ) + + U = np.eye(2**nq, dtype=complex) + U[:n, :n] = QFT + op = Operator(U) + assert op.is_unitary() + + redistribution_circuit.append(op, list(range(nq))) + + circuit.compose( + redistribution_circuit.control( + nv - int(nq), + label=rf"MCRY(π/2, {nq})" if is_two_pow(n) else rf"Coll({n}, {nq})", + ), + qubits=list(range(nv - 1, -1, -1)), + inplace=True, + ) + + if not is_two_pow(n): + return circuit.decompose() if self.decompose_block else circuit + else: + return circuit + + @override + def __str__(self): + return f"[EQCRedistribution({self.equivalence_class})]" diff --git a/qlbm/components/common/primitives.py b/qlbm/components/common/primitives.py index cd37d92..ef629fc 100644 --- a/qlbm/components/common/primitives.py +++ b/qlbm/components/common/primitives.py @@ -5,7 +5,7 @@ from typing import List, Tuple from qiskit import QuantumCircuit -from qiskit.circuit.library import MCMT, XGate +from qiskit.circuit.library import MCMTGate, XGate from typing_extensions import override from qlbm.components.base import LBMPrimitive @@ -44,7 +44,7 @@ def __init__( @override def create_circuit(self) -> QuantumCircuit: - return QuantumCircuit(*self.lattice.registers) + return self.lattice.circuit.copy() @override def __str__(self) -> str: @@ -93,7 +93,9 @@ def create_circuit(self) -> QuantumCircuit: circuit.cx(self.target_qubits[1], self.target_qubits[0]) circuit.compose( - MCMT(XGate(), len(self.control_qubits) + 1, len(self.target_qubits) - 1), + MCMTGate( + XGate(), len(self.control_qubits) + 1, len(self.target_qubits) - 1 + ), qubits=self.control_qubits + list(self.target_qubits), inplace=True, ) diff --git a/qlbm/components/lqlga/__init__.py b/qlbm/components/lqlga/__init__.py new file mode 100644 index 0000000..11a630a --- /dev/null +++ b/qlbm/components/lqlga/__init__.py @@ -0,0 +1,17 @@ +"""Implementation of the Linear Quantum Lattice Gas Algorithm (LQLGA) algorithm.""" + +from .collision import GenericLQLGACollisionOperator +from .initial import LQGLAInitialConditions +from .lqlga import LQLGA +from .measurement import LQLGAGridVelocityMeasurement +from .reflection import LQLGAReflectionOperator +from .streaming import LQLGAStreamingOperator + +__all__ = [ + "GenericLQLGACollisionOperator", + "LQGLAInitialConditions", + "LQLGA", + "LQLGAGridVelocityMeasurement", + "LQLGAReflectionOperator", + "LQLGAStreamingOperator", +] diff --git a/qlbm/components/lqlga/collision.py b/qlbm/components/lqlga/collision.py new file mode 100644 index 0000000..3efbcf6 --- /dev/null +++ b/qlbm/components/lqlga/collision.py @@ -0,0 +1,61 @@ +"""Collision operators for the :class:`.LQLGA` algorithm.""" + +from logging import Logger, getLogger +from time import perf_counter_ns + +from qiskit import QuantumCircuit +from typing_extensions import override + +from qlbm.components.base import LQLGAOperator +from qlbm.components.common.cbse_collision.cbse_collision import EQCCollisionOperator +from qlbm.lattice.lattices.lqlga_lattice import LQLGALattice + + +class GenericLQLGACollisionOperator(LQLGAOperator): + """ + Equivalence class-based LGA collision operator for the :class:`.LQLGA` algorithm. + + This operator applies the :class:`.EQCCollisionOperator` operator to all velocity qubits at each grid point. + """ + + def __init__( + self, + lattice: LQLGALattice, + logger: Logger = getLogger("qlbm"), + ) -> None: + super().__init__(lattice, logger) + self.lattice = lattice + + self.logger.info(f"Creating circuit {str(self)}...") + circuit_creation_start_time = perf_counter_ns() + self.circuit = self.create_circuit() + self.logger.info( + f"Creating circuit {str(self)} took {perf_counter_ns() - circuit_creation_start_time} (ns)" + ) + + @override + def create_circuit(self) -> QuantumCircuit: + local_collision_circuit = EQCCollisionOperator( + self.lattice.discretization + ).circuit + circuit = self.lattice.circuit.copy() + + for velocity_qubit_indices in range( + 0, + self.lattice.num_base_qubits, + self.lattice.num_velocities_per_point, + ): + circuit.compose( + local_collision_circuit, + inplace=True, + qubits=range( + velocity_qubit_indices, + velocity_qubit_indices + self.lattice.num_velocities_per_point, + ), + ) + + return circuit + + @override + def __str__(self) -> str: + return f"[GenericSpaceTimeCollisionOperator for discretization {self.lattice.discretization}]" diff --git a/qlbm/components/lqlga/initial.py b/qlbm/components/lqlga/initial.py new file mode 100644 index 0000000..e20c1da --- /dev/null +++ b/qlbm/components/lqlga/initial.py @@ -0,0 +1,89 @@ +"""Initial conditions for the :class:`.LQLGA` algorithm.""" + +from logging import Logger, getLogger +from time import perf_counter_ns +from typing import List, Tuple + +from typing_extensions import override + +from qlbm.components.base import LBMPrimitive +from qlbm.lattice.lattices.lqlga_lattice import LQLGALattice + + +class LQGLAInitialConditions(LBMPrimitive): + """ + Primitive for setting initial conditions in the :class:`.LQLGA` algorithm. + + This operator allows the construction of arbitrary deterministic initial conditions for the LQLGA algorithm. + The number of gates required by this operator is equal to the number of enabled velocity qubits across all grid points. + The depth of the circuit is 1, as all gates are applied in parallel at each grid point. + + Example usage: + + .. plot:: + :include-source: + + from qlbm.lattice import LQLGALattice + from qlbm.components.lqlga import LQGLAInitialConditions + + lattice = LQLGALattice( + { + "lattice": { + "dim": {"x": 4}, + "velocities": "D1Q3", + }, + "geometry": [], + }, + ) + initial_conditions = LQGLAInitialConditions(lattice, [(tuple([2]), (True, True, True))]) + initial_conditions.draw("mpl") + """ + + grid_data: List[Tuple[Tuple[int, ...], Tuple[bool, ...]]] + """ + Grid data for the initial conditions, where each tuple contains: + #. A tuple of grid point indices (e.g., `(x, y, z)`). + #. A tuple of booleans indicating which velocity qubits are enabled at that grid point. + """ + + def __init__( + self, + lattice: LQLGALattice, + grid_data: List[Tuple[Tuple[int, ...], Tuple[bool, ...]]], + logger: Logger = getLogger("qlbm"), + ): + super().__init__(logger) + + self.lattice = lattice + self.grid_data = grid_data + + self.logger.info(f"Creating circuit {str(self)}...") + circuit_creation_start_time = perf_counter_ns() + self.circuit = self.create_circuit() + self.logger.info( + f"Creating circuit {str(self)} took {perf_counter_ns() - circuit_creation_start_time} (ns)" + ) + + @override + def create_circuit(self): + circuit = self.lattice.circuit.copy() + + for tup in self.grid_data: + gp, vel_profile = tup[0], tup[1] + + if not any(vel_profile): + continue + + circuit.x( + self.lattice.gridpoint_index_tuple(gp) + * self.lattice.num_velocities_per_point + + v + for v, is_enabled in enumerate(vel_profile) + if is_enabled + ) + + return circuit + + @override + def __str__(self): + return f"[Primitive LQGLAInitialConditions on lattice={self.lattice}, grid_data={self.grid_data})]" diff --git a/qlbm/components/lqlga/lqlga.py b/qlbm/components/lqlga/lqlga.py new file mode 100644 index 0000000..d4e6dc5 --- /dev/null +++ b/qlbm/components/lqlga/lqlga.py @@ -0,0 +1,86 @@ +"""The end-to-end algorithm of the Space-Time Quantum Lattice Boltzmann Algorithm described in :cite:`spacetime`.""" + +from logging import Logger, getLogger + +from qiskit import QuantumCircuit +from typing_extensions import override + +from qlbm.components.base import LBMAlgorithm +from qlbm.components.lqlga.collision import GenericLQLGACollisionOperator +from qlbm.components.lqlga.reflection import LQLGAReflectionOperator +from qlbm.components.lqlga.streaming import LQLGAStreamingOperator +from qlbm.lattice.lattices.lqlga_lattice import LQLGALattice + + +class LQLGA(LBMAlgorithm): + r""" + Implementation of the Linear Quantum Lattice Gas Algorithm (LQLGA). + + For a lattice with :math:`N_g` gridpoints and :math:`q` discrete velocities, + LQLGA requires exactly :math:`N_g \cdot q` qubits. + + That is exactly equal to the number of classical bits required for one + deterministic run of the classical LGA algorithm. + + More information about this algorithm can be found in :cite:t:`lqlga1`, :cite:t:`lqlga2`, and :cite:t:`spacetime2`. + + Example usage: + + .. plot:: + :include-source: + + from qlbm.components.lqlga import LQLGA + from qlbm.lattice import LQLGALattice + + lattice = LQLGALattice( + { + "lattice": { + "dim": {"x": 7}, + "velocities": "D1Q3", + }, + "geometry": [{"shape": "cuboid", "x": [3, 5], "boundary": "bounceback"}], + }, + ) + + LQLGA(lattice=lattice).draw("mpl") + + """ + + lattice: LQLGALattice + + def __init__( + self, + lattice: LQLGALattice, + logger: Logger = getLogger("qlbm"), + ): + super().__init__(lattice, logger) + + self.lattice = lattice + + self.circuit = self.create_circuit() + + @override + def create_circuit(self) -> QuantumCircuit: + circuit = self.lattice.circuit.copy() + + circuit.compose( + GenericLQLGACollisionOperator(self.lattice, self.logger).circuit, + inplace=True, + ) + + circuit.compose( + LQLGAStreamingOperator(self.lattice, self.logger).circuit, inplace=True + ) + + circuit.compose( + LQLGAReflectionOperator( + self.lattice, self.lattice.shapes["bounceback"], self.logger + ).circuit, + inplace=True, + ) + + return circuit + + @override + def __str__(self) -> str: + return f"[LQLGA on lattice={self.lattice}]" diff --git a/qlbm/components/lqlga/measurement.py b/qlbm/components/lqlga/measurement.py new file mode 100644 index 0000000..f9f1568 --- /dev/null +++ b/qlbm/components/lqlga/measurement.py @@ -0,0 +1,66 @@ +"""Measurement operator for the :class:`.SpaceTimeQLBM` algorithm :cite:`spacetime`.""" + +from logging import Logger, getLogger + +from qiskit import ClassicalRegister +from typing_extensions import override + +from qlbm.components.base import LQLGAOperator +from qlbm.lattice.lattices.lqlga_lattice import LQLGALattice + + +class LQLGAGridVelocityMeasurement(LQLGAOperator): + # TODO: Improve documentation + """ + Measurement operator for the :class:`.LQLGA` algorithm. + + This operator measures the velocity qubits at each grid point in the LQLGA lattice. + + Example usage: + + .. plot:: + :include-source: + + from qlbm.components.lqlga import LQLGAGridVelocityMeasurement + from qlbm.lattice import LQLGALattice + + lattice = LQLGALattice( + { + "lattice": { + "dim": {"x": 5}, + "velocities": "D1Q3", + }, + "geometry": [], + }, + ) + + LQLGAGridVelocityMeasurement(lattice=lattice).draw("mpl") + """ + + def __init__( + self, + lattice: LQLGALattice, + logger: Logger = getLogger("qlbm"), + ) -> None: + super().__init__(lattice, logger) + self.lattice = lattice + + self.circuit = self.create_circuit() + + @override + def create_circuit(self): + circuit = self.lattice.circuit.copy() + + qubits_to_measure = list(range(self.lattice.num_base_qubits)) + circuit.add_register(ClassicalRegister(self.lattice.num_base_qubits)) + + circuit.measure( + qubits_to_measure, + list(range(len(qubits_to_measure))), + ) + + return circuit + + @override + def __str__(self) -> str: + return f"[LQLGAGridVelocityMeasurement for lattice {self.lattice}]" diff --git a/qlbm/components/lqlga/reflection.py b/qlbm/components/lqlga/reflection.py new file mode 100644 index 0000000..d990d36 --- /dev/null +++ b/qlbm/components/lqlga/reflection.py @@ -0,0 +1,125 @@ +"""Reflection operator for the :class:`.LQLGA` algorithm :cite:`spacetime` that swaps particles one gridpoint at a time.""" + +from logging import Logger, getLogger +from time import perf_counter_ns +from typing import List, cast + +from qiskit import QuantumCircuit +from typing_extensions import override + +from qlbm.components.base import LQLGAOperator +from qlbm.lattice.geometry.shapes.base import LQLGAShape, Shape +from qlbm.lattice.lattices.lqlga_lattice import LQLGALattice +from qlbm.lattice.spacetime.properties_base import LatticeDiscretization +from qlbm.tools.exceptions import CircuitException + + +class LQLGAReflectionOperator(LQLGAOperator): + """ + Operator implementing reflection in the :class:`.LQLGA` algorithm. + + Reflections in this algorithm can be entirely implemented by swap gates. + The number of gates scales with the number of grridpoints of the solid geometry. + The depth of the operator is 1. + + ============================ ====================================================================== + Attribute Summary + ============================ ====================================================================== + :attr:`lattice` The lattice the operator acts on. + :attr:`shapes` A list of boundary-conditioned shapes. + :attr:`logger` The performance logger, by default ``getLogger("qlbm")``. + ============================ ====================================================================== + + Example usage: + + .. plot:: + :include-source: + + from qlbm.components.lqlga import LQLGAReflectionOperator + from qlbm.lattice import LQLGALattice + + lattice = LQLGALattice( + { + "lattice": { + "dim": {"x": 7}, + "velocities": "D1Q3", + }, + "geometry": [{"shape": "cuboid", "x": [3, 5], "boundary": "bounceback"}], + }, + ) + reflection_operator = LQLGAReflectionOperator( + lattice, shapes=lattice.shapes["bounceback"] + ) + reflection_operator.draw("mpl") + + """ + + shapes: List[LQLGAShape] + """ + A list of shapes that require reflection at the boundaries. + """ + + def __init__( + self, + lattice: LQLGALattice, + shapes: List[Shape], + logger: Logger = getLogger("qlbm"), + ) -> None: + super().__init__(lattice, logger) + self.shapes = cast(List[LQLGAShape], shapes) + + self.logger.info(f"Creating circuit {str(self)}...") + circuit_creation_start_time = perf_counter_ns() + self.circuit = self.create_circuit() + self.logger.info( + f"Creating circuit {str(self)} took {perf_counter_ns() - circuit_creation_start_time} (ns)" + ) + + @override + def create_circuit(self) -> QuantumCircuit: + discretization = self.lattice.discretization + if discretization == LatticeDiscretization.D1Q2: + return self.__create_circuit_d1q2() + + elif discretization == LatticeDiscretization.D1Q3: + return self.__create_circuit_d1q3() + + raise CircuitException(f"Reflection Operator unsupported for {discretization}.") + + def __create_circuit_d1q2(self) -> QuantumCircuit: + circuit = self.lattice.circuit.copy() + + for shape in self.shapes: + for reflection_data in shape.get_lqlga_reflection_data_d1q2(): + circuit.swap( + self.lattice.velocity_index_tuple( + reflection_data.gridpoints[0], + reflection_data.velocity_indices_to_swap[0], + ), + self.lattice.velocity_index_tuple( + reflection_data.gridpoints[1], + reflection_data.velocity_indices_to_swap[1], + ), + ) + return circuit + + def __create_circuit_d1q3(self) -> QuantumCircuit: + circuit = self.lattice.circuit.copy() + + for shape in self.shapes: + for reflection_data in shape.get_lqlga_reflection_data_d1q3(): + circuit.swap( + self.lattice.velocity_index_tuple( + reflection_data.gridpoints[0], + reflection_data.velocity_indices_to_swap[0], + ), + self.lattice.velocity_index_tuple( + reflection_data.gridpoints[1], + reflection_data.velocity_indices_to_swap[1], + ), + ) + return circuit + + @override + def __str__(self) -> str: + return f"[PointWiseLQLGAReflectionOperator for lattice {self.lattice}, shapes {self.shapes}]" diff --git a/qlbm/components/lqlga/streaming.py b/qlbm/components/lqlga/streaming.py new file mode 100644 index 0000000..373b7a8 --- /dev/null +++ b/qlbm/components/lqlga/streaming.py @@ -0,0 +1,120 @@ +"""Streaming operators for the :class:`.SpaceTimeQLBM` algorithm :cite:`spacetime`.""" + +from logging import Logger, getLogger +from time import perf_counter_ns +from typing import List, Tuple + +from typing_extensions import override + +from qlbm.components.base import LQLGAOperator +from qlbm.lattice.lattices.lqlga_lattice import LQLGALattice + + +class LQLGAStreamingOperator(LQLGAOperator): + # TODO: Improve documentation + """ + Streaming operator for the :class:`.LQLGA` algorithm. + + Streaming is implemented by a series of swap gates as described in :cite:`spacetime`. + The number of gates scales linearly with size of the grid, + while the depth scales logarithmically. + + Example usage: + + .. plot:: + :include-source: + + from qlbm.components.lqlga import LQLGAStreamingOperator + from qlbm.lattice import LQLGALattice + + lattice = LQLGALattice( + { + "lattice": { + "dim": {"x": 4}, + "velocities": "D1Q3", + }, + "geometry": [], + }, + ) + streaming_operator = LQLGAStreamingOperator(lattice) + streaming_operator.draw("mpl") + """ + + def __init__( + self, + lattice: LQLGALattice, + logger: Logger = getLogger("qlbm"), + ) -> None: + super().__init__(lattice, logger) + self.lattice = lattice + + self.logger.info(f"Creating circuit {str(self)}...") + circuit_creation_start_time = perf_counter_ns() + self.circuit = self.create_circuit() + self.logger.info( + f"Creating circuit {str(self)} took {perf_counter_ns() - circuit_creation_start_time} (ns)" + ) + + @override + def create_circuit(self): + circuit = self.lattice.circuit.copy() + # ! TODO Generalize in 2 and 3D + num_gps = self.lattice.num_gridpoints[0] + 1 + + for direction, velocity_qubit_of_line in enumerate( + self.lattice.get_velocity_qubits_of_line(0) + ): + gridpoints_to_swap = self.logarithmic_depth_streaming_line_swaps( + num_gps, negative_direction=bool(direction) + ) + for layer in gridpoints_to_swap: + for i, j in layer: + circuit.swap( + self.lattice.velocity_index_flat(i, velocity_qubit_of_line), + self.lattice.velocity_index_flat(j, velocity_qubit_of_line), + ) + + return circuit + + def logarithmic_depth_streaming_line_swaps( + self, num_gridpoints: int, negative_direction: bool + ) -> List[List[Tuple[int, int]]]: + """ + + Implements the logarithmic depth streaming line permuation as described in Section 4 of :cite:`spacetime`. + + Parameters + ---------- + num_gridpoints : int + The number of gridpoints in the streaming line. + negative_direction : bool + Whether streaming occurs in the negative direction (i.e., from high to low indices). + + Returns + ------- + List[List[Tuple[int, int]]] + A list of layers, where each layer is a list of tuples representing pairs of gridpoints to swap. + """ + if num_gridpoints < 2: + return [] + + layers: List[List[Tuple[int, int]]] = [] + stride = 1 + + while stride < num_gridpoints: + layer: List[Tuple[int, int]] = [] + for i in range(0, num_gridpoints, 2 * stride): + if i + stride < num_gridpoints: + layer.append( + (i, i + stride) + if not negative_direction + else (num_gridpoints - 1 - i, num_gridpoints - 1 - i - stride) + ) + layers.append(layer) + stride *= 2 + + return layers + + @override + def __str__(self): + return f"[LQLGAStreamingOperator on lattice={self.lattice}]" diff --git a/qlbm/components/spacetime/__init__.py b/qlbm/components/spacetime/__init__.py index f743f4f..7aca5cc 100644 --- a/qlbm/components/spacetime/__init__.py +++ b/qlbm/components/spacetime/__init__.py @@ -1,6 +1,6 @@ """Modular qlbm quantum circuit components for the CQLBM algorithm :cite:p:`spacetime`.""" -from .collision import SpaceTimeCollisionOperator +from .collision.d2q4_old import SpaceTimeD2Q4CollisionOperator from .initial.pointwise import PointWiseSpaceTimeInitialConditions from .initial.volumetric import VolumetricSpaceTimeInitialConditions from .measurement import SpaceTimeGridVelocityMeasurement @@ -9,7 +9,7 @@ from .streaming import SpaceTimeStreamingOperator __all__ = [ - "SpaceTimeCollisionOperator", + "SpaceTimeD2Q4CollisionOperator", "PointWiseSpaceTimeInitialConditions", "VolumetricSpaceTimeInitialConditions", "SpaceTimeStreamingOperator", diff --git a/qlbm/components/spacetime/collision/__init__.py b/qlbm/components/spacetime/collision/__init__.py new file mode 100644 index 0000000..db8a47d --- /dev/null +++ b/qlbm/components/spacetime/collision/__init__.py @@ -0,0 +1,15 @@ +"""Classes implementing the logic of the collision operator in the Space-Time Data encoding.""" + +from ....lattice.eqc.eqc import EquivalenceClass +from ....lattice.eqc.eqc_generator import ( + EquivalenceClassGenerator, +) +from .d2q4_old import SpaceTimeD2Q4CollisionOperator +from .eqc_collision import GenericSpaceTimeCollisionOperator + +__all__ = [ + "SpaceTimeD2Q4CollisionOperator", + "EquivalenceClass", + "EquivalenceClassGenerator", + "GenericSpaceTimeCollisionOperator", +] diff --git a/qlbm/components/spacetime/collision.py b/qlbm/components/spacetime/collision/d2q4_old.py similarity index 93% rename from qlbm/components/spacetime/collision.py rename to qlbm/components/spacetime/collision/d2q4_old.py index 5b768be..b70b6a2 100644 --- a/qlbm/components/spacetime/collision.py +++ b/qlbm/components/spacetime/collision/d2q4_old.py @@ -1,4 +1,4 @@ -"""Collision operators for the :class:`.SpaceTimeQLBM` algorithm :cite:`spacetime`.""" +"""Collision operator for the :math:`D_2Q_4` discretization :class:`.SpaceTimeQLBM` algorithm as described in :cite:`spacetime`.""" from logging import Logger, getLogger from math import pi @@ -6,14 +6,14 @@ from qiskit import QuantumCircuit from qiskit.circuit import Gate -from qiskit.circuit.library import MCMT, RYGate +from qiskit.circuit.library import MCMTGate, RYGate from typing_extensions import override from qlbm.components.base import SpaceTimeOperator from qlbm.lattice.lattices.spacetime_lattice import SpaceTimeLattice -class SpaceTimeCollisionOperator(SpaceTimeOperator): +class SpaceTimeD2Q4CollisionOperator(SpaceTimeOperator): r"""An operator that performs collision part of the :class:`.SpaceTimeQLBM` algorithm. Collision is a local operation that is performed simultaneously on all velocity qubits corresponding to a grid location. @@ -46,20 +46,20 @@ class SpaceTimeCollisionOperator(SpaceTimeOperator): .. plot:: :include-source: - from qlbm.components.spacetime import SpaceTimeCollisionOperator + from qlbm.components.spacetime.collision.d2q4_old import SpaceTimeD2Q4CollisionOperator from qlbm.lattice import SpaceTimeLattice # Build an example lattice lattice = SpaceTimeLattice( num_timesteps=1, lattice_data={ - "lattice": {"dim": {"x": 4, "y": 8}, "velocities": {"x": 2, "y": 2}}, + "lattice": {"dim": {"x": 4, "y": 8}, "velocities": "D2Q4"}, "geometry": [], }, ) # Draw the collision operator for 1 time step - SpaceTimeCollisionOperator(lattice=lattice, timestep=1).draw("mpl") + SpaceTimeD2Q4CollisionOperator(lattice=lattice, timestep=1).draw("mpl") """ def __init__( @@ -87,7 +87,7 @@ def create_circuit(self) -> QuantumCircuit: collision_circuit = self.local_collision_circuit(reset_state=False) collision_circuit.compose( - MCMT( + MCMTGate( self.gate_to_apply, self.lattice.properties.get_num_velocities_per_point() - 1, 1, @@ -158,5 +158,4 @@ def local_collision_circuit(self, reset_state: bool) -> QuantumCircuit: @override def __str__(self) -> str: - # TODO: Implement return "Space Time Collision Operator" diff --git a/qlbm/components/spacetime/collision/eqc_collision.py b/qlbm/components/spacetime/collision/eqc_collision.py new file mode 100644 index 0000000..5ac8cd4 --- /dev/null +++ b/qlbm/components/spacetime/collision/eqc_collision.py @@ -0,0 +1,70 @@ +"""Collision operators for the :class:`.SpaceTimeQLBM` algorithm :cite:`spacetime`.""" + +from logging import Logger, getLogger +from time import perf_counter_ns + +from qiskit import QuantumCircuit +from typing_extensions import override + +from qlbm.components.base import SpaceTimeOperator +from qlbm.components.common.cbse_collision.cbse_collision import EQCCollisionOperator +from qlbm.lattice.lattices.spacetime_lattice import SpaceTimeLattice + + +class GenericSpaceTimeCollisionOperator(SpaceTimeOperator): + """ + A generic space-time collision operator that can be used to apply any gate to the velocities of a grid location. + + ========================= ====================================================================== + Attribute Summary + ========================= ====================================================================== + :attr:`lattice` The :class:`.SpaceTimeLattice` based on which the properties of the operator are inferred. + :attr:`gate` The gate to apply to the velocities. + :attr:`logger` The performance logger, by default ``getLogger("qlbm")``. + ========================= ====================================================================== + """ + + def __init__( + self, + lattice: SpaceTimeLattice, + timestep: int, + logger: Logger = getLogger("qlbm"), + ) -> None: + super().__init__(lattice, logger) + self.lattice = lattice + self.timestep = timestep + + self.logger.info(f"Creating circuit {str(self)}...") + circuit_creation_start_time = perf_counter_ns() + self.circuit = self.create_circuit() + self.logger.info( + f"Creating circuit {str(self)} took {perf_counter_ns() - circuit_creation_start_time} (ns)" + ) + + @override + def create_circuit(self) -> QuantumCircuit: + local_collision_circuit = EQCCollisionOperator( + self.lattice.properties.get_discretization() + ).circuit + circuit = self.lattice.circuit.copy() + + for velocity_qubit_indices in range( + self.lattice.properties.get_num_grid_qubits(), + self.lattice.properties.get_num_grid_qubits() + + self.lattice.properties.get_num_velocity_qubits(self.timestep), + self.lattice.properties.get_num_velocities_per_point(), + ): + circuit.compose( + local_collision_circuit, + inplace=True, + qubits=range( + velocity_qubit_indices, + velocity_qubit_indices + + self.lattice.properties.get_num_velocities_per_point(), + ), + ) + return circuit + + @override + def __str__(self) -> str: + return f"[GenericSpaceTimeCollisionOperator for discretization {self.lattice.properties.get_discretization()}]" diff --git a/qlbm/components/spacetime/initial/pointwise.py b/qlbm/components/spacetime/initial/pointwise.py index 923b8c4..e23451a 100644 --- a/qlbm/components/spacetime/initial/pointwise.py +++ b/qlbm/components/spacetime/initial/pointwise.py @@ -4,7 +4,7 @@ from typing import List, Tuple from qiskit import QuantumCircuit -from qiskit.circuit.library import MCMT, XGate +from qiskit.circuit.library import MCXGate from typing_extensions import override from qlbm.components.base import LBMPrimitive @@ -55,7 +55,7 @@ class PointWiseSpaceTimeInitialConditions(LBMPrimitive): lattice = SpaceTimeLattice( num_timesteps=1, lattice_data={ - "lattice": {"dim": {"x": 4, "y": 8}, "velocities": {"x": 2, "y": 2}}, + "lattice": {"dim": {"x": 4, "y": 8}, "velocities": "D2Q4"}, "geometry": [], }, ) @@ -93,26 +93,20 @@ def create_circuit(self) -> QuantumCircuit: if self.lattice.is_inside_an_obstacle(grid_point_data[0]): continue circuit.compose(self.set_grid_value(grid_point_data[0]), inplace=True) - circuit.compose( - MCMT( - XGate(), - num_ctrl_qubits=self.lattice.properties.get_num_grid_qubits(), - num_target_qubits=sum( - grid_point_data[1] - ), # The sum is equal to the number of velocities set to true - ), - qubits=list( - self.lattice.grid_index() - + flatten( - [ - self.lattice.velocity_index(0, c) - for c, is_velocity_enabled in enumerate(grid_point_data[1]) - if is_velocity_enabled - ] - ) - ), - inplace=True, - ) + for target_qubit in flatten( + [ + self.lattice.velocity_index(0, c) + for c, is_velocity_enabled in enumerate(grid_point_data[1]) + if is_velocity_enabled + ] + ): + circuit.compose( + MCXGate( + num_ctrl_qubits=self.lattice.properties.get_num_grid_qubits(), + ), + qubits=list(self.lattice.grid_index() + [target_qubit]), + inplace=True, + ) circuit.compose(self.set_grid_value(grid_point_data[0]), inplace=True) # Set the velocity state for neighbors in increasing velocity @@ -203,26 +197,20 @@ def set_neighbor_velocity( circuit.compose( self.set_grid_value(absolute_neighbor_coordinates), inplace=True ) - circuit.compose( - MCMT( - XGate(), - num_ctrl_qubits=self.lattice.properties.get_num_grid_qubits(), - num_target_qubits=sum( - velocity_values - ), # The sum is equal to the number of velocities set to true - ), - inplace=True, - qubits=list( - self.lattice.grid_index() - + flatten( - [ - self.lattice.velocity_index(neighbor.neighbor_index, c) - for c, is_velocity_enabled in enumerate(velocity_values) - if is_velocity_enabled - ] - ) - ), - ) + for target_qubit in flatten( + [ + self.lattice.velocity_index(neighbor.neighbor_index, c) + for c, is_velocity_enabled in enumerate(velocity_values) + if is_velocity_enabled + ] + ): + circuit.compose( + MCXGate( + num_ctrl_qubits=self.lattice.properties.get_num_grid_qubits(), + ), + qubits=list(self.lattice.grid_index() + [target_qubit]), + inplace=True, + ) circuit.compose( self.set_grid_value(absolute_neighbor_coordinates), inplace=True ) diff --git a/qlbm/components/spacetime/initial/volumetric.py b/qlbm/components/spacetime/initial/volumetric.py index f24b52b..9d28bc8 100644 --- a/qlbm/components/spacetime/initial/volumetric.py +++ b/qlbm/components/spacetime/initial/volumetric.py @@ -5,7 +5,7 @@ from typing import List, Tuple, cast from qiskit import QuantumCircuit -from qiskit.circuit.library import MCMT, XGate +from qiskit.circuit.library import MCMTGate, XGate from typing_extensions import override from qlbm.components.base import LBMPrimitive @@ -109,7 +109,7 @@ def create_circuit(self) -> QuantumCircuit: [any(pvb[1]) for pvb in periodic_volume_bounds] ): circuit.compose( - MCMT( + MCMTGate( XGate(), num_ctrl_qubits=len(control_qubit_sequence), num_target_qubits=sum( diff --git a/qlbm/components/spacetime/measurement.py b/qlbm/components/spacetime/measurement.py index cffbbe9..1095179 100644 --- a/qlbm/components/spacetime/measurement.py +++ b/qlbm/components/spacetime/measurement.py @@ -4,7 +4,7 @@ from typing import Tuple from qiskit import ClassicalRegister -from qiskit.circuit.library import MCMT, XGate +from qiskit.circuit.library import MCMTGate, XGate from typing_extensions import override from qlbm.components.base import SpaceTimeOperator @@ -39,7 +39,7 @@ class SpaceTimeGridVelocityMeasurement(SpaceTimeOperator): lattice = SpaceTimeLattice( num_timesteps=1, lattice_data={ - "lattice": {"dim": {"x": 4, "y": 8}, "velocities": {"x": 2, "y": 2}}, + "lattice": {"dim": {"x": 4, "y": 8}, "velocities": "D2Q4"}, "geometry": [], }, ) @@ -79,8 +79,7 @@ def create_circuit(self): @override def __str__(self) -> str: - # TODO: Implement - return "Space Gird Measurement" + return f"[SpaceTimeGridVelocityMeasurement for lattice {self.lattice}]" class SpaceTimePointWiseMassMeasurement(SpaceTimeOperator): @@ -132,7 +131,7 @@ def create_circuit(self): target_qubits = self.lattice.ancilla_mass_index() circuit.compose( - MCMT( + MCMTGate( XGate(), len(control_qubits), len(target_qubits), @@ -150,5 +149,4 @@ def create_circuit(self): @override def __str__(self) -> str: - # TODO: Implement - return "SpaceTimePointWiseMassMeasurement" + return f"[SpaceTimePointWiseMassMeasurement for lattice {self.lattice}, gridpoint {self.gridpoint}, velocity {self.velocity_index_to_measure}]" diff --git a/qlbm/components/spacetime/reflection/pointwise.py b/qlbm/components/spacetime/reflection/pointwise.py index 0be4024..8dad761 100644 --- a/qlbm/components/spacetime/reflection/pointwise.py +++ b/qlbm/components/spacetime/reflection/pointwise.py @@ -9,7 +9,7 @@ from qlbm.components.base import SpaceTimeOperator from qlbm.components.common.primitives import MCSwap -from qlbm.lattice.geometry.shapes.block import Block +from qlbm.lattice.geometry.shapes.base import Shape, SpaceTimeShape from qlbm.lattice.lattices.spacetime_lattice import SpaceTimeLattice from qlbm.lattice.spacetime.properties_base import LatticeDiscretization from qlbm.tools.exceptions import CircuitException @@ -25,7 +25,7 @@ class PointWiseSpaceTimeReflectionOperator(SpaceTimeOperator): ============================ ====================================================================== :attr:`lattice` The :class:`.SpaceTimeLattice` based on which the properties of the operator are inferred. :attr:`timestep` The timestep for to which to perform reflection. - :attr:`blocks` A list of :class:`.Block` objects for which to generate the BB boundary condition circuits. + :attr:`shapes` A list of :class:`.Shape` objects for which to generate the BB boundary condition circuits. :attr:`filter_inside_blocks` A ``bool`` that, when enabled, disregards operations that would have happened inside blocks. :attr:`logger` The performance logger, by default ``getLogger("qlbm")``. ============================ ====================================================================== @@ -36,7 +36,7 @@ def __init__( self, lattice: SpaceTimeLattice, timestep: int, - blocks: List[Block], + shapes: List[Shape], filter_inside_blocks: bool = True, logger: Logger = getLogger("qlbm"), ) -> None: @@ -48,7 +48,7 @@ def __init__( f"Invalid time step {timestep}, select a value between 1 and {lattice.num_timesteps}" ) - self.blocks = blocks + self.shapes = cast(List[SpaceTimeShape], shapes) self.filter_inside_blocks = filter_inside_blocks self.logger.info(f"Creating circuit {str(self)}...") @@ -71,7 +71,7 @@ def create_circuit(self) -> QuantumCircuit: def __create_circuit_d1q2(self) -> QuantumCircuit: circuit = self.lattice.circuit.copy() - for block in self.blocks: + for block in self.shapes: for reflection_data in block.get_spacetime_reflection_data_d1q2( self.lattice.properties, self.timestep ): @@ -116,7 +116,7 @@ def __create_circuit_d1q2(self) -> QuantumCircuit: def __create_circuit_d2q4(self) -> QuantumCircuit: circuit = self.lattice.circuit.copy() - for block in self.blocks: + for block in self.shapes: reflection_data_points = block.get_spacetime_reflection_data_d2q4( self.lattice.properties, self.timestep ) @@ -169,5 +169,4 @@ def __create_circuit_d2q4(self) -> QuantumCircuit: @override def __str__(self) -> str: - # TODO: Implement - return "Space Time Reflection Operator" + return f"[PointWiseSpaceTimeReflectionOperator for lattice {self.lattice}, blocks {self.shapes}, filtered {self.filter_inside_blocks}]" diff --git a/qlbm/components/spacetime/spacetime.py b/qlbm/components/spacetime/spacetime.py index 508d8c8..88771fd 100644 --- a/qlbm/components/spacetime/spacetime.py +++ b/qlbm/components/spacetime/spacetime.py @@ -8,8 +8,9 @@ from qlbm.components.base import LBMAlgorithm from qlbm.components.spacetime.reflection import PointWiseSpaceTimeReflectionOperator from qlbm.lattice.lattices.spacetime_lattice import SpaceTimeLattice +from qlbm.tools.exceptions import LatticeException -from .collision import SpaceTimeCollisionOperator +from .collision.d2q4_old import SpaceTimeD2Q4CollisionOperator from .streaming import SpaceTimeStreamingOperator @@ -23,7 +24,7 @@ class SpaceTimeQLBM(LBMAlgorithm): The algorithm is composed of two main steps, the implementation of which (in general) varies per individual time step: #. Streaming performed by the :class:`.SpaceTimeStreamingOperator` moves the particles on the grid by means of swap gates over velocity qubits. - #. Collision performed by the :class:`.SpaceTimeCollisionOperator` does not move particles on the grid, but locally alters the velocity qubits at each grid point, if applicable. + #. Collision performed by the :class:`.GenericSpaceTimeCollisionOperator` does not move particles on the grid, but locally alters the velocity qubits at each grid point, if applicable. ========================= ====================================================================== Attribute Summary @@ -44,7 +45,7 @@ class SpaceTimeQLBM(LBMAlgorithm): lattice = SpaceTimeLattice( num_timesteps=1, lattice_data={ - "lattice": {"dim": {"x": 4, "y": 8}, "velocities": {"x": 2, "y": 2}}, + "lattice": {"dim": {"x": 4, "y": 8}, "velocities": "D2Q4"}, "geometry": [], }, ) @@ -72,6 +73,12 @@ def create_circuit(self) -> QuantumCircuit: circuit = self.lattice.circuit.copy() for timestep in range(self.lattice.num_timesteps, 0, -1): + # Warn the user if there are any shapes that are NOT bounceback boundary conditions. + if self.lattice.shapes["specular"]: + raise LatticeException( + "Currently, the Space-Time QLBM algorithm only supports bounceback boundary conditions." + ) + circuit.compose( SpaceTimeStreamingOperator(self.lattice, timestep, self.logger).circuit, inplace=True, @@ -81,7 +88,7 @@ def create_circuit(self) -> QuantumCircuit: PointWiseSpaceTimeReflectionOperator( self.lattice, timestep, - self.lattice.blocks["bounceback"], + self.lattice.shapes["bounceback"], self.filter_inside_blocks, self.logger, ).circuit, @@ -91,7 +98,7 @@ def create_circuit(self) -> QuantumCircuit: # There is no collision in 1D if self.lattice.num_dims > 1: circuit.compose( - SpaceTimeCollisionOperator( + SpaceTimeD2Q4CollisionOperator( self.lattice, timestep, logger=self.logger ).circuit, inplace=True, @@ -101,4 +108,4 @@ def create_circuit(self) -> QuantumCircuit: @override def __str__(self) -> str: - return f"[Space Time QLBM on lattice={self.lattice}]" + return f"[SpaceTimes QLBM on lattice={self.lattice}]" diff --git a/qlbm/components/spacetime/streaming.py b/qlbm/components/spacetime/streaming.py index b4be9a5..88fc302 100644 --- a/qlbm/components/spacetime/streaming.py +++ b/qlbm/components/spacetime/streaming.py @@ -44,7 +44,7 @@ class SpaceTimeStreamingOperator(SpaceTimeOperator): lattice = SpaceTimeLattice( num_timesteps=1, lattice_data={ - "lattice": {"dim": {"x": 4, "y": 8}, "velocities": {"x": 2, "y": 2}}, + "lattice": {"dim": {"x": 4, "y": 8}, "velocities": "D2Q4"}, "geometry": [], }, ) @@ -165,5 +165,4 @@ def stream_lines( @override def __str__(self) -> str: - # TODO: Implement - return "Space Time Streaming Operator" + return f"[SpaceTimeStreamingOperator for lattice {self.lattice}]" diff --git a/qlbm/infra/__init__.py b/qlbm/infra/__init__.py index 132b189..f7f9839 100644 --- a/qlbm/infra/__init__.py +++ b/qlbm/infra/__init__.py @@ -7,7 +7,7 @@ from .compiler import CircuitCompiler from .result import CollisionlessResult, SpaceTimeResult -from .runner import CircuitRunner, QiskitRunner, QulacsRunner, SimulationConfig +from .runner import CircuitRunner, QiskitRunner, SimulationConfig __all__ = [ "CircuitCompiler", @@ -16,6 +16,5 @@ "CollisionlessResult", "SpaceTimeResult", "QiskitRunner", - "QulacsRunner", "SimulationConfig", ] diff --git a/qlbm/infra/compiler.py b/qlbm/infra/compiler.py index c280489..313610a 100644 --- a/qlbm/infra/compiler.py +++ b/qlbm/infra/compiler.py @@ -12,7 +12,6 @@ from qiskit import QuantumCircuit as QiskitQC from qiskit.compiler import transpile from qiskit_aer.backends.aerbackend import AerBackend -from qiskit_qulacs import QulacsProvider from qulacs import QuantumCircuit as QulacsQC from qlbm.components.base import QuantumComponent @@ -51,7 +50,7 @@ class CircuitCompiler: lattice = SpaceTimeLattice( num_timesteps=1, lattice_data={ - "lattice": {"dim": {"x": 4, "y": 8}, "velocities": {"x": 2, "y": 2}}, + "lattice": {"dim": {"x": 4, "y": 8}, "velocities": "D2Q4"}, "geometry": [], }, ) @@ -196,7 +195,10 @@ def compile( self.logger.warn( "Provided backend is ignored. The qiskit-qulacs backend is always used when compiling to Qulacs from Qiskit." ) - backend = QulacsProvider().get_backend("qulacs_simulator") + raise CompilerException( + "The qiksit-qulacs extension does not currently support qiskit>=2.0." + ) + # backend = QulacsProvider().get_backend("qulacs_simulator") if self.compiler_target == "QISKIT": if not isinstance(backend, AerBackend): diff --git a/qlbm/infra/reinitialize/__init__.py b/qlbm/infra/reinitialize/__init__.py index a834c09..1c06f5a 100644 --- a/qlbm/infra/reinitialize/__init__.py +++ b/qlbm/infra/reinitialize/__init__.py @@ -1,7 +1,7 @@ """Reinitialize objects for transitioning between time steps of the QLBM algorithm.""" from .base import Reinitializer -from .collisionless_reinitializer import CollisionlessReinitializer +from .identity_reinitializer import IdentityReinitializer from .spacetime_reinitializer import SpaceTimeReinitializer -__all__ = ["Reinitializer", "CollisionlessReinitializer", "SpaceTimeReinitializer"] +__all__ = ["Reinitializer", "IdentityReinitializer", "SpaceTimeReinitializer"] diff --git a/qlbm/infra/reinitialize/collisionless_reinitializer.py b/qlbm/infra/reinitialize/identity_reinitializer.py similarity index 82% rename from qlbm/infra/reinitialize/collisionless_reinitializer.py rename to qlbm/infra/reinitialize/identity_reinitializer.py index dfb4086..f38c434 100644 --- a/qlbm/infra/reinitialize/collisionless_reinitializer.py +++ b/qlbm/infra/reinitialize/identity_reinitializer.py @@ -1,4 +1,4 @@ -""":class:`.CQLBM`-specific implementation of the :class:`.Reinitializer`.""" +"""Identity reinitializer used for the :class:`.CQLBM` and :class:`.LQLGA` algorithms.""" from logging import Logger, getLogger @@ -11,15 +11,15 @@ from typing_extensions import override from qlbm.infra.compiler import CircuitCompiler -from qlbm.lattice import CollisionlessLattice +from qlbm.infra.reinitialize.base import Reinitializer +from qlbm.lattice.lattices.base import Lattice -from .base import Reinitializer - -class CollisionlessReinitializer(Reinitializer): +class IdentityReinitializer(Reinitializer): r""" - :class:`.CQLBM`-specific implementation of the :class:`.Reinitializer`. + Implementation of the :class:`.Reinitializer` that passes along the statevector to the following time step. + Useful for the :class:`.CQLBM` and :class:`.LQLGA` algorithms. Compatible with both :class:`.QiskitRunner`\ s and :class:`.QulacsRunner`\ s. To generate a new set of initial conditions for the CQLBM algorithm, the reinitializer simply returns the quantum state computed @@ -31,19 +31,19 @@ class CollisionlessReinitializer(Reinitializer): =========================== ====================================================================== Attribute Summary =========================== ====================================================================== - :attr:`lattice` The :class:`.CollisionlessLattice` of the simulated system. + :attr:`lattice` The :class:`.Lattice` of the simulated system. :attr:`compiler` The compiler that converts the novel initial conditions circuits. :attr:`logger` The performance logger, by default ``getLogger("qlbm")`` =========================== ====================================================================== """ - lattice: CollisionlessLattice + lattice: Lattice statevector: Statevector counts: Counts def __init__( self, - lattice: CollisionlessLattice, + lattice: Lattice, compiler: CircuitCompiler, logger: Logger = getLogger("qlbm"), ): @@ -87,4 +87,4 @@ def reinitialize( @override def requires_statevector(self) -> bool: - return True + return True \ No newline at end of file diff --git a/qlbm/infra/result/__init__.py b/qlbm/infra/result/__init__.py index bc012b2..78cfce1 100644 --- a/qlbm/infra/result/__init__.py +++ b/qlbm/infra/result/__init__.py @@ -2,10 +2,7 @@ from .base import QBMResult from .collisionless_result import CollisionlessResult +from .lqlga_result import LQLGAResult from .spacetime_result import SpaceTimeResult -__all__ = [ - "QBMResult", - "CollisionlessResult", - "SpaceTimeResult", -] +__all__ = ["QBMResult", "CollisionlessResult", "SpaceTimeResult", "LQLGAResult"] diff --git a/qlbm/infra/result/base.py b/qlbm/infra/result/base.py index d4b258b..ff3b39d 100644 --- a/qlbm/infra/result/base.py +++ b/qlbm/infra/result/base.py @@ -62,9 +62,9 @@ def visualize_geometry(self): Creates ``stl`` files for each block in the lattice. Output files are formatted as ``output_dir/paraview_dir/cube_.stl``. - The output is created through the :class:`.Block`'s :meth:`.Block.stl_mesh` method. + The output is created through the :class:`.Shape`'s :meth:`.Shape.stl_mesh` method. """ - for c, shape in enumerate(flatten(self.lattice.blocks.values())): + for c, shape in enumerate(flatten(self.lattice.shapes.values())): shape.stl_mesh().save(f"{self.paraview_dir}/{shape.name()}_{c}.stl") def save_timestep_array( diff --git a/qlbm/infra/result/collisionless_result.py b/qlbm/infra/result/collisionless_result.py index 070546c..ea9b88c 100644 --- a/qlbm/infra/result/collisionless_result.py +++ b/qlbm/infra/result/collisionless_result.py @@ -25,15 +25,21 @@ class CollisionlessResult(QBMResult): =========================== ====================================================================== :attr:`lattice` The :class:`.CollisionlessLattice` of the simulated system. :attr:`directory` The directory to which the results outputs data to. - :attr:`paraview_dir` The subdirectory under ``directory`` which stores the Paraview files. :attr:`output_file_name` The root name for files containing time step artifacts, by default "step". =========================== ====================================================================== """ num_steps: int + """The time step to which this result corresponds.""" + directory: str + """The output directory for the results.""" + output_file_name: str + """The name of the file to output the artifacts to.""" + lattice: CollisionlessLattice + """The lattice the result corresponds to.""" def __init__( self, diff --git a/qlbm/infra/result/lqlga_result.py b/qlbm/infra/result/lqlga_result.py new file mode 100644 index 0000000..5d45326 --- /dev/null +++ b/qlbm/infra/result/lqlga_result.py @@ -0,0 +1,123 @@ +""":class:`.LQLGA`-specific implementation of the :class:`.QBMResult`.""" + +import re +from os import listdir +from typing import Dict + +import numpy as np +import vtk +from typing_extensions import override +from vtkmodules.util import numpy_support + +from qlbm.lattice.lattices.lqlga_lattice import LQLGALattice +from qlbm.lattice.spacetime.properties_base import LatticeDiscretizationProperties + +from .base import QBMResult + + +class LQLGAResult(QBMResult): + """ + :class:`.LQLGA`-specific implementation of the :class:`.QBMResult`. + + Processes counts sampled from :class:`.LQLGAGridVelocityMeasurement` primitives. + + =========================== ====================================================================== + Attribute Summary + =========================== ====================================================================== + :attr:`lattice` The :class:`.LQLGA` of the simulated system. + :attr:`directory` The directory to which the results outputs data to. + :attr:`output_file_name` The root name for files containing time step artifacts, by default "step". + =========================== ====================================================================== + """ + + num_steps: int + """The time step to which this result corresponds.""" + + directory: str + """The output directory for the results.""" + + output_file_name: str + """The name of the file to output the artifacts to.""" + + lattice: LQLGALattice + """The lattice the result corresponds to.""" + + def __init__( + self, + lattice: LQLGALattice, + directory: str, + output_file_name: str = "step", + ) -> None: + super().__init__(lattice, directory, output_file_name) + + @override + def save_timestep_counts( + self, + counts: Dict[str, float], + timestep: int, + create_vis: bool = True, + save_array: bool = False, + ): + total_counts = sum(counts.values()) + if self.lattice.num_dims == 1: + # The second dimension is a dirty rendering trick for VTK and Paraview + count_history = np.zeros((self.lattice.num_gridpoints[0] + 1, 2)) + channel_masses = LatticeDiscretizationProperties.get_channel_masses( + self.lattice.discretization + ) + for count in counts: + count_inverse = count[::-1] + num_vel = self.lattice.num_velocities_per_point + for gp in range( + self.lattice.num_base_qubits + // self.lattice.num_velocities_per_point + ): + mass = np.dot( + np.array( + list( + map( + lambda x: float(x), + count_inverse[gp * num_vel : (gp + 1) * num_vel], + ) + ) + ), + channel_masses, + ) + pops = counts[count] * mass / total_counts + count_history[gp][0] += pops + count_history[gp][1] += pops + self.save_timestep_array( + np.transpose(count_history), + timestep, + create_vis=create_vis, + save_counts_array=save_array, + ) + + @override + def visualize_all_numpy_data(self): + # Filter the algorithm output files + r = re.compile("[a-zA-Z0-9]+_[0-9]+.csv") + for data_file_name in filter(r.match, listdir(self.directory)): + data = np.genfromtxt( + f"{self.directory}/{data_file_name}", + dtype=None, + delimiter=",", + autostrip=True, + ) + vtk_data = numpy_support.numpy_to_vtk( + num_array=data.flatten, deep=True, array_type=vtk.VTK_FLOAT + ) + img = vtk.vtkImageData() + img.SetDimensions( + self.lattice.num_gridpoints[0] + 1, + self.lattice.num_gridpoints[1] + 1 if self.lattice.num_dims > 1 else 1, + self.lattice.num_gridpoints[2] + 1 if self.lattice.num_dims > 2 else 1, + ) + img.GetPointData().SetScalars(vtk_data) + + writer = vtk.vtkXMLImageDataWriter() + writer.SetFileName( + f"{self.paraview_dir}/{data_file_name.split('.')[0]}.vti" + ) + writer.SetInputData(img) + writer.Write() diff --git a/qlbm/infra/result/spacetime_result.py b/qlbm/infra/result/spacetime_result.py index bed4691..60a3506 100644 --- a/qlbm/infra/result/spacetime_result.py +++ b/qlbm/infra/result/spacetime_result.py @@ -25,15 +25,21 @@ class SpaceTimeResult(QBMResult): =========================== ====================================================================== :attr:`lattice` The :class:`.SpaceTimeLattice` of the simulated system. :attr:`directory` The directory to which the results outputs data to. - :attr:`paraview_dir` The subdirectory under ``directory`` which stores the Paraview files. :attr:`output_file_name` The root name for files containing time step artifacts, by default "step". =========================== ====================================================================== """ num_steps: int + """The time step to which this result corresponds.""" + directory: str + """The output directory for the results.""" + output_file_name: str + """The name of the file to output the artifacts to.""" + lattice: SpaceTimeLattice + """The lattice the result corresponds to.""" def __init__( self, diff --git a/qlbm/infra/runner/base.py b/qlbm/infra/runner/base.py index acb84d2..0de08ca 100644 --- a/qlbm/infra/runner/base.py +++ b/qlbm/infra/runner/base.py @@ -10,12 +10,18 @@ from qiskit_aer import AerSimulator from qlbm.infra.reinitialize import ( - CollisionlessReinitializer, Reinitializer, SpaceTimeReinitializer, ) -from qlbm.infra.result import CollisionlessResult, QBMResult, SpaceTimeResult +from qlbm.infra.reinitialize.identity_reinitializer import IdentityReinitializer +from qlbm.infra.result import ( + CollisionlessResult, + LQLGAResult, + QBMResult, + SpaceTimeResult, +) from qlbm.lattice import CollisionlessLattice, Lattice +from qlbm.lattice.lattices.lqlga_lattice import LQLGALattice from qlbm.lattice.lattices.spacetime_lattice import SpaceTimeLattice from qlbm.tools.exceptions import CircuitException, ResultsException @@ -127,6 +133,10 @@ def new_result(self, output_directory: str, output_file_name: str) -> QBMResult: return SpaceTimeResult( cast(SpaceTimeLattice, self.lattice), output_directory, output_file_name ) + elif isinstance(self.lattice, LQLGALattice): + return LQLGAResult( + cast(LQLGALattice, self.lattice), output_directory, output_file_name + ) else: raise ResultsException(f"Unsupported lattice: {self.lattice}.") @@ -144,9 +154,11 @@ def new_reinitializer(self) -> Reinitializer: ResultsException If the underlying algorithm does not support reinitialization. """ - if isinstance(self.lattice, CollisionlessLattice): - return CollisionlessReinitializer( - cast(CollisionlessLattice, self.lattice), + if isinstance(self.lattice, CollisionlessLattice) or isinstance( + self.lattice, LQLGALattice + ): + return IdentityReinitializer( + self.lattice, self.config.get_execution_compiler(), self.logger, ) diff --git a/qlbm/infra/runner/qiskit_runner.py b/qlbm/infra/runner/qiskit_runner.py index 2e643b9..ec12b8b 100644 --- a/qlbm/infra/runner/qiskit_runner.py +++ b/qlbm/infra/runner/qiskit_runner.py @@ -4,6 +4,7 @@ from time import perf_counter_ns from qiskit import QuantumCircuit as QiskitQC +from qiskit import transpile from qiskit.circuit.library import Initialize from qiskit.quantum_info import Statevector from typing_extensions import override @@ -62,6 +63,7 @@ def run( output_directory: str, output_file_name: str = "step", statevector_snapshots: bool = False, + recompile_each_step: bool = False, # ! Document ) -> QBMResult: if (self.sampling_backend is None) and self.config.statevector_sampling: raise ExecutionException( @@ -89,6 +91,7 @@ def run( num_shots, initial_conditions, simulation_result, + recompile_each_step, ) if statevector_snapshots else self._run_time_loop( @@ -111,6 +114,7 @@ def _run_snapshot_time_loop( num_shots: int, initial_conditions: QiskitQC, simulation_result: QBMResult, + recompile_each_step: bool = False, # ! Document ) -> QBMResult: for step in range(num_steps + 1): step_start_time = perf_counter_ns() @@ -124,7 +128,15 @@ def _run_snapshot_time_loop( ) # The first step consists of just initial conditions if step > 0: - circuit.compose(self.config.algorithm, inplace=True) + if recompile_each_step: + circuit.compose(self.config.algorithm_copy, inplace=True) + circuit = transpile( + circuit, + backend=self.execution_backend, + optimization_level=0, + ) + else: + circuit.compose(self.config.algorithm, inplace=True) if ( self.reinitializer.requires_statevector() diff --git a/qlbm/infra/runner/simulation_config.py b/qlbm/infra/runner/simulation_config.py index 6df6a91..00d4dde 100644 --- a/qlbm/infra/runner/simulation_config.py +++ b/qlbm/infra/runner/simulation_config.py @@ -230,6 +230,11 @@ def __init__( # Circuits self.initial_conditions = initial_conditions self.algorithm = algorithm + self.algorithm_copy = ( + algorithm.circuit.copy() + if isinstance(algorithm, QuantumComponent) + else algorithm.copy() + ) self.postprocessing = postprocessing self.measurement = measurement diff --git a/qlbm/lattice/__init__.py b/qlbm/lattice/__init__.py index 368b435..b709602 100644 --- a/qlbm/lattice/__init__.py +++ b/qlbm/lattice/__init__.py @@ -9,16 +9,28 @@ from .geometry.shapes.block import ( Block, ) +from .geometry.shapes.circle import ( + Circle, +) from .lattices import CollisionlessLattice, Lattice +from .lattices.lqlga_lattice import LQLGALattice from .lattices.spacetime_lattice import SpaceTimeLattice +from .spacetime.properties_base import ( + LatticeDiscretization, + LatticeDiscretizationProperties, +) __all__ = [ "Lattice", "CollisionlessLattice", "SpaceTimeLattice", + "LQLGALattice", "DimensionalReflectionData", "ReflectionWall", "ReflectionPoint", "ReflectionResetEdge", "Block", + "Circle", + "LatticeDiscretization", + "LatticeDiscretizationProperties", ] diff --git a/qlbm/lattice/eqc/__init__.py b/qlbm/lattice/eqc/__init__.py new file mode 100644 index 0000000..fc8ac86 --- /dev/null +++ b/qlbm/lattice/eqc/__init__.py @@ -0,0 +1,9 @@ +"""Utilities for the equivalence class logic handling in LGA-based algorithms.""" + +from .eqc import EquivalenceClass +from .eqc_generator import EquivalenceClassGenerator + +__all__ = [ + "EquivalenceClass", + "EquivalenceClassGenerator", +] diff --git a/qlbm/lattice/eqc/eqc.py b/qlbm/lattice/eqc/eqc.py new file mode 100644 index 0000000..84c8d85 --- /dev/null +++ b/qlbm/lattice/eqc/eqc.py @@ -0,0 +1,163 @@ +"""Implementation of the equivalence class abstraction described in Section 4 of :cite:p:`spacetime2`.""" + +from typing import List, Set, Tuple, override + +import numpy as np + +from qlbm.lattice.spacetime.properties_base import ( + LatticeDiscretization, + LatticeDiscretizationProperties, +) +from qlbm.tools.exceptions import LatticeException + + +class EquivalenceClass: + """ + Class representing LGA equivalence classes. + + In ``qlbm``, an equivalence class is a set of velocity configurations that share the same mass and momentum. + For a more in depth explanation, consult Section 4 of :cite:p:`spacetime2`. + + .. list-table:: Constructor Attributes + :widths: 25 50 + :header-rows: 1 + + * - Attribute + - Description + * - :attr:`discretization` + - The :class:`.LatticeDiscretization` that the equivalence class belongs to. + * - :attr:`velocity_configurations` + - The :class:`Set[Tuple[bool, ...]]` that contains the velocity configurations of the equivalence class. Configurations are stored as ``q``-tuples where an entry is ``True`` if the velocity channel is occupied and ``False`` otherwise. + + .. list-table:: Class Attributes + :widths: 25 50 + :header-rows: 1 + + * - Attribute + - Description + * - :attr:`mass` + - The total mass of the equivalence class, which is the sum of all occupied velocity channels. + * - :attr:`momentum` + - The total momentum of the equivalence class, which is the vector sum of all occupied velocity channels multiplid by their :class:`.LatticeDiscretizationProperties` velocity contribution. + """ + + discretization: LatticeDiscretization + """The discretization of the equivalence class.""" + + velocity_configurations: Set[Tuple[bool, ...]] + """The velocity configurations of the equivalence class, represented as a set of tuples of boolean corresponding to velocity channel occupancies that all have the same mass and momentum.""" + + mass: int + """The total mass of the equivalence class, which is the sum of all occupied velocity channels.""" + + momentum: np.typing.NDArray + """The total momentum of the equivalence class, which is the vector sum of all occupied velocity channels multiplied by their :class:`.LatticeDiscretizationProperties` velocity contribution.""" + + def __init__( + self, + discretization: LatticeDiscretization, + velocity_configurations: Set[Tuple[bool, ...]], + ): + if len(velocity_configurations) < 2: + raise LatticeException( + f"Equivalence class must have at least two configurations. Provided velocity profiles are {velocity_configurations}." + ) + num_velocities = LatticeDiscretizationProperties.get_num_velocities( + discretization + ) + if not all(len(v) == num_velocities for v in velocity_configurations): + raise LatticeException( + f"All velocity configurations must have length {num_velocities}. Provided configurations are {velocity_configurations}." + ) + + self.discretization = discretization + self.velocity_configurations = velocity_configurations + channel_masses = LatticeDiscretizationProperties.get_channel_masses( + discretization + ) + self.mass = np.dot(channel_masses, list(velocity_configurations)[0]) + + velocity_vectors = LatticeDiscretizationProperties.get_velocity_vectors( + discretization + ) + + if not all( + (np.dot(velocity_cfg, channel_masses) == self.mass) + for velocity_cfg in velocity_configurations + ): + raise LatticeException("Velocity configurations have different masses.") + self.momentum = sum( + [ + velocity_vectors[i] * list(velocity_configurations)[0][i] + for i in range(len(list(velocity_configurations)[0])) + ] + ) + if not all( + ( + np.array_equal( + self.momentum, + np.sum( + [ + velocity_vectors[i] * v[i] + for i in range(len(velocity_vectors)) + ], + axis=0, + ), + ) + ) + for v in velocity_configurations + ): + raise LatticeException("Velocity configurations have different momenta.") + + def size(self) -> int: + """ + The number of velocity configurations in the equivalence class. + + Returns + ------- + int + The number of velocity configurations in the equivalence class. + """ + return len(self.velocity_configurations) + + def id(self) -> Tuple[int, List[float]]: + """ + The identifier of the equivalence class. + + For a given discretization, an equivalence class can be uniquely identified by the common mass and momentum of its velocity configurations. + + Returns + ------- + Tuple[int, List[float]] + The mass and momentum of the equivalence class. + """ + return (self.mass, self.momentum.tolist()) + + def get_bitstrings(self) -> List[str]: + """ + Returns the velocity configurations as bitstrings. + + Returns + ------- + List[str] + The velocity configurations of the equivalence class as bitstrings. + """ + return [ + "".join(["1" if x else "0" for x in cfg]) + for cfg in self.velocity_configurations + ] + + @override + def __eq__(self, value): + if not isinstance(value, EquivalenceClass): + return False + return ( + self.discretization == value.discretization + and self.velocity_configurations == value.velocity_configurations + and self.mass == value.mass + and np.array_equal(self.momentum, value.momentum) + ) + + @override + def __hash__(self): + return hash((self.discretization, tuple(self.velocity_configurations))) diff --git a/qlbm/lattice/eqc/eqc_generator.py b/qlbm/lattice/eqc/eqc_generator.py new file mode 100644 index 0000000..d9df93c --- /dev/null +++ b/qlbm/lattice/eqc/eqc_generator.py @@ -0,0 +1,83 @@ +"""Generator class for equivalence classes in LGA-based algorithms.""" + +from itertools import product +from typing import Dict, Set + +import numpy as np + +from qlbm.lattice.eqc.eqc import EquivalenceClass +from qlbm.lattice.spacetime.properties_base import ( + LatticeDiscretization, + LatticeDiscretizationProperties, +) + + +class EquivalenceClassGenerator: + """ + A class that generates equivalence classes for a given lattice discretization. + + .. list-table:: Constructor Attributes + :widths: 25 50 + :header-rows: 1 + + * - Attribute + - Description + * - :attr:`discretization` + - The :class:`.LatticeDiscretization` that the equivalence class belongs to. + + Example usage: + + .. code-block:: python + :linenos: + + from qlbm.lattice import LatticeDiscretization + from qlbm.lattice.eqc import EquivalenceClassGenerator + + # Generate some equivalence classes + eqcs = EquivalenceClassGenerator( + LatticeDiscretization.D3Q6 + ).generate_equivalence_classes() + + print(eqcs.pop().get_bitstrings()) + """ + + discretization: LatticeDiscretization + + def __init__(self, discretization): + self.discretization = discretization + + def generate_equivalence_classes(self) -> Set[EquivalenceClass]: + """ + Generates equivalence classes for the given lattice discretization. + + Returns + ------- + Set[EquivalenceClass] + All equivalence classes of the discretization. + """ + equivalence_classes: Dict = {} + for state in product( + [0, 1], + repeat=LatticeDiscretizationProperties.get_num_velocities( + self.discretization + ), + ): + velocity_vectors = LatticeDiscretizationProperties.get_velocity_vectors( + self.discretization + ) + state = np.array(state) # type: ignore + mass = state @ LatticeDiscretizationProperties.get_channel_masses( + self.discretization + ) + momentum = np.sum(state[:, None] * velocity_vectors, axis=0) # type: ignore + + key = (mass, tuple(momentum)) + if key not in equivalence_classes: + equivalence_classes[key] = [] + equivalence_classes[key].append(state) + + return { + EquivalenceClass(self.discretization, set(tuple(cfg.tolist()) for cfg in v)) + for _, v in equivalence_classes.items() + if len(v) > 1 + } diff --git a/qlbm/lattice/geometry/__init__.py b/qlbm/lattice/geometry/__init__.py index 54d9e72..d78094e 100644 --- a/qlbm/lattice/geometry/__init__.py +++ b/qlbm/lattice/geometry/__init__.py @@ -2,6 +2,8 @@ from .encodings import ( DimensionalReflectionData, + LQLGAPointwiseReflectionData, + LQLGAReflectionData, ReflectionPoint, ReflectionResetEdge, ReflectionWall, @@ -21,4 +23,6 @@ "SpaceTimeVolumetricReflectionData", "Block", "Circle", + "LQLGAPointwiseReflectionData", + "LQLGAReflectionData", ] diff --git a/qlbm/lattice/geometry/encodings/__init__.py b/qlbm/lattice/geometry/encodings/__init__.py index 293aa9d..462e32c 100644 --- a/qlbm/lattice/geometry/encodings/__init__.py +++ b/qlbm/lattice/geometry/encodings/__init__.py @@ -6,6 +6,7 @@ ReflectionResetEdge, ReflectionWall, ) +from .lqlga import LQLGAPointwiseReflectionData, LQLGAReflectionData from .spacetime import ( SpaceTimeDiagonalReflectionData, SpaceTimePWReflectionData, @@ -20,4 +21,6 @@ "SpaceTimePWReflectionData", "SpaceTimeVolumetricReflectionData", "SpaceTimeDiagonalReflectionData", + "LQLGAPointwiseReflectionData", + "LQLGAReflectionData", ] diff --git a/qlbm/lattice/geometry/encodings/lqlga.py b/qlbm/lattice/geometry/encodings/lqlga.py new file mode 100644 index 0000000..fb5ecdb --- /dev/null +++ b/qlbm/lattice/geometry/encodings/lqlga.py @@ -0,0 +1,33 @@ +"""Geometrical data encodings specific to the :class:`.LQLGA` algorithm.""" + +from abc import ABC +from typing import Tuple + + +class LQLGAReflectionData(ABC): + """Base class for all space-time reflectiopn data structures.""" + + pass + + +class LQLGAPointwiseReflectionData(LQLGAReflectionData): + """ + Data structure representing pointwise space-time reflection information for LQLGA geometries. + + Attributes + ---------- + gridpoints (Tuple[Tuple[int, ...], Tuple[int, ...]]): + A tuple containing two grid points (as tuples of integers) that are + related by a reflection operation. + velocity_indices_to_swap (Tuple[int, int]): + A tuple of two velocity indices that should be swapped during the + reflection operation. + """ + + def __init__( + self, + gridpoints: Tuple[Tuple[int, ...], Tuple[int, ...]], + velocity_indices_to_swap: Tuple[int, int], + ) -> None: + self.gridpoints = gridpoints + self.velocity_indices_to_swap = velocity_indices_to_swap diff --git a/qlbm/lattice/geometry/encodings/spacetime.py b/qlbm/lattice/geometry/encodings/spacetime.py index e437053..56cd294 100644 --- a/qlbm/lattice/geometry/encodings/spacetime.py +++ b/qlbm/lattice/geometry/encodings/spacetime.py @@ -7,7 +7,7 @@ class SpaceTimeReflectionData(ABC): - """Base class for all space-time reflectiopn data structures.""" + """Base class for all space-time reflection data structures.""" pass diff --git a/qlbm/lattice/geometry/shapes/base.py b/qlbm/lattice/geometry/shapes/base.py index 94b4369..0dd18f9 100644 --- a/qlbm/lattice/geometry/shapes/base.py +++ b/qlbm/lattice/geometry/shapes/base.py @@ -5,6 +5,7 @@ from stl import mesh +from qlbm.lattice.geometry.encodings.lqlga import LQLGAPointwiseReflectionData from qlbm.lattice.geometry.encodings.spacetime import ( SpaceTimePWReflectionData, SpaceTimeVolumetricReflectionData, @@ -65,6 +66,76 @@ def to_dict(self): pass +class LQLGAShape(Shape): + """Base class for all shapes compatible with the :class:`.LQLGA` algorithm.""" + + def __init__(self, num_grid_qubits: List[int], boundary_condition: str): + super().__init__(num_grid_qubits, boundary_condition) + + def get_lqlga_reflection_data_1d_from_points( + self, + gridpoints: List[Tuple[int, ...]], + before_reflection_velocity_index: int, + after_reflection_velocity_index: int, + reflection_increment_from_boundary: Tuple[int, ...], + max_grid_size: int, + ) -> List[LQLGAPointwiseReflectionData]: + """Calculate LQLGA reflection data for D1Q2 lattice from a list of points. + + Parameters + ---------- + gridpoints : List[Tuple[int, ...]] + List of grid points where reflections occur, inside the solid domain. + before_reflection_velocity_index : int + Velocity index before reflection. + after_reflection_velocity_index : int + Velocity index after reflection. + reflection_increment_from_boundary : Tuple[int, ...] + Increment vector from boundary point to reflection destination. + max_grid_size: int + The size of the grid on which reflection is performed, to account for periodicity. + + Returns + ------- + List[LQLGAPointwiseReflectionData] + List of reflection data for each point. + """ + return [ + LQLGAPointwiseReflectionData( + tuple( # type: ignore + [ + gridpoint, + tuple( + (a + b) % max_grid_size + for a, b in zip( + gridpoint, + reflection_increment_from_boundary, + ) + ), + ], + ), + tuple( # type: ignore + [before_reflection_velocity_index, after_reflection_velocity_index], + ), + ) + for gridpoint in gridpoints + ] + + @abstractmethod + def get_lqlga_reflection_data_d1q2( + self, + ) -> List[LQLGAPointwiseReflectionData]: + """Calculate space-time reflection data for :math:`D_1Q_2` :class:`.LQLGA`.""" + pass + + @abstractmethod + def get_lqlga_reflection_data_d1q3( + self, + ) -> List[LQLGAPointwiseReflectionData]: + """Calculate space-time reflection data for :math:`D_1Q_2` :class:`.LQLGA`.""" + pass + + class SpaceTimeShape(Shape): """Base class for all shapes compatible with the :class:`.SpaceTimeQLBM` algorithm.""" diff --git a/qlbm/lattice/geometry/shapes/block.py b/qlbm/lattice/geometry/shapes/block.py index e860ddc..3704c1c 100644 --- a/qlbm/lattice/geometry/shapes/block.py +++ b/qlbm/lattice/geometry/shapes/block.py @@ -18,12 +18,12 @@ SpaceTimePWReflectionData, SpaceTimeVolumetricReflectionData, ) -from qlbm.lattice.geometry.shapes.base import SpaceTimeShape +from qlbm.lattice.geometry.shapes.base import LQLGAShape, SpaceTimeShape from qlbm.lattice.spacetime.properties_base import SpaceTimeLatticeBuilder from qlbm.tools.utils import bit_value, dimension_letter, flatten, get_qubits_to_invert -class Block(SpaceTimeShape): +class Block(SpaceTimeShape, LQLGAShape): r""" Contains information required for the generation of boundary conditions for an axis-parallel cuboid obstacle. @@ -707,6 +707,30 @@ def get_d2q4_surfaces(self) -> List[List[List[Tuple[int, ...]]]]: return surfaces + @override + def get_lqlga_reflection_data_d1q2(self): + return self.get_lqlga_reflection_data_1d_from_points( + [tuple([self.bounds[0][0]])], + 0, + 1, + tuple([-1]), + 2 ** self.num_grid_qubits[0], + ) + self.get_lqlga_reflection_data_1d_from_points( + [tuple([self.bounds[0][1]])], 1, 0, tuple([1]), 2 ** self.num_grid_qubits[0] + ) + + @override + def get_lqlga_reflection_data_d1q3(self): + return self.get_lqlga_reflection_data_1d_from_points( + [tuple([self.bounds[0][0]])], + 1, + 2, + tuple([-1]), + 2 ** self.num_grid_qubits[0], + ) + self.get_lqlga_reflection_data_1d_from_points( + [tuple([self.bounds[0][1]])], 2, 1, tuple([1]), 2 ** self.num_grid_qubits[0] + ) + @override def contains_gridpoint(self, gridpoint: Tuple[int, ...]) -> bool: return all( diff --git a/qlbm/lattice/geometry/shapes/circle.py b/qlbm/lattice/geometry/shapes/circle.py index c173476..8979b34 100644 --- a/qlbm/lattice/geometry/shapes/circle.py +++ b/qlbm/lattice/geometry/shapes/circle.py @@ -433,9 +433,6 @@ def get_spacetime_reflection_data_d2q4( expanded_diag_gridpoints: List[Tuple[int, ...]] = ( Circle.expand_diagonal_segments([diag_segment]) ) - - if (4, 1) in expanded_diag_gridpoints: - print("ok!") for reflection_dim in [0, 1]: other_dim = 1 - reflection_dim streaming_line_velocities = [0, 2] if reflection_dim == 0 else [1, 3] diff --git a/qlbm/lattice/lattices/__init__.py b/qlbm/lattice/lattices/__init__.py index 9e74f3d..7a92929 100644 --- a/qlbm/lattice/lattices/__init__.py +++ b/qlbm/lattice/lattices/__init__.py @@ -2,6 +2,7 @@ from .base import Lattice from .collisionless_lattice import CollisionlessLattice +from .lqlga_lattice import LQLGALattice from .spacetime_lattice import SpaceTimeLattice -__all__ = ["Lattice", "CollisionlessLattice", "SpaceTimeLattice"] +__all__ = ["Lattice", "CollisionlessLattice", "SpaceTimeLattice", "LQLGALattice"] diff --git a/qlbm/lattice/lattices/base.py b/qlbm/lattice/lattices/base.py index fcc843b..6ae96bf 100644 --- a/qlbm/lattice/lattices/base.py +++ b/qlbm/lattice/lattices/base.py @@ -10,8 +10,12 @@ from qlbm.lattice.geometry.shapes.base import Shape from qlbm.lattice.geometry.shapes.block import Block from qlbm.lattice.geometry.shapes.circle import Circle +from qlbm.lattice.spacetime.properties_base import ( + LatticeDiscretization, + LatticeDiscretizationProperties, +) from qlbm.tools.exceptions import LatticeException -from qlbm.tools.utils import dimension_letter, flatten, is_two_pow +from qlbm.tools.utils import dimension_letter, flatten class Lattice(ABC): @@ -31,22 +35,17 @@ class Lattice(ABC): Attribute Summary =========================== ====================================================================== :attr:`num_dims` The number of dimensions of the lattice. - :attr:`num_gridpoints` A ``List[int]`` of the number of gridpoints of the lattice in each dimension. - **Important**\ : for easier compatibility with binary arithmetic, the number of gridpoints - specified in the input dicitionary is one larger than the one held in the ``Lattice``. - That is, for a ``16x64`` lattice, the ``num_gridpoints`` attribute will have the value ``[15, 63]``. - :attr:`num_grid_qubits` The total number of qubits required to encode the lattice grid. + :attr:`num_gridpoints` The number of gridpoints in each dimension. + :attr:`num_grid_qubits` The total number of qubits required to encode the grid. + :attr:`num_velocities` The number of discrete velocities in each dimension. :attr:`num_velocity_qubits` The total number of qubits required to encode the velocity discretization of the lattice. - :attr:`num_ancilla_qubits` The total number of ancilla (non-velocity, non-grid) qubits required for the quantum circuit to simulate this lattice. + :attr:`num_ancilla_qubits` The total number of ancllary qubits required for the quantum circuit to simulate this lattice. :attr:`num_total_qubits` The total number of qubits required for the quantum circuit to simulate the lattice. - This is the sum of the number of grid, velocity, and ancilla qubits. - :attr:`registers` A ``Tuple[qiskit.QuantumRegister, ...]`` that holds registers responsible for specific operations of the QLBM algorithm. - :attr:`circuit` An empty ``qiskit.QuantumCircuit`` with labeled registers that quantum components use as a base. - Each quantum component that is parameterized by a ``Lattice`` makes a copy of this quantum circuit - to which it appends its designated logic. - :attr:`blocks` A ``Dict[str, List[Block]]`` that contains all of the :class:`.Block`\ s encoding the solid geometry of the lattice. - The key of the dictionary is the specific kind of boundary condition of the obstacle (i.e., ``"bounceback"`` or ``"specular"``). - :attr:`logger` The performance logger, by default ``getLogger("qlbm")``. + :attr:`registers` The qubit registers of the quantum algorithm. + :attr:`circuit` The blueprint quantum circuit for all components of the algorithm. + :attr:`discretization` The discretization of the lattice, as an enum value of :class:`.LatticeDiscretization`. + :attr:`shapes` A list of the solid geometry objects. + :attr:`logger` The performance logger. =========================== ====================================================================== A lattice can be constructed from from either an input file or a Python dictionary. @@ -81,13 +80,81 @@ class Lattice(ABC): """ num_dims: int + """ + The number of dimensions of the system. + """ num_gridpoints: List[int] + """ + The number of gridpoints of the lattice in each dimension. + + .. warning:: + + For easier compatibility with binary arithmetic, the number of gridpoints + specified in the input dictionary is one larger than the one held in the :class:`.Lattice` s. + That is, for a :math:`16 \\times 64` lattice, the ``num_gridpoints`` attribute will have the value ``[15, 63]``. + """ + num_velocities: List[int] + """ + The number of velocities in each dimension. This will be refactored in the future to support :math:`D_dQ_q` discretizations. + + .. warning:: + + For easier compatibility with binary arithmetic, the number of velocities + specified in the input dictionary is one larger than the one held in the :class:`.Lattice` s. + If the numbers of discrete velocities are :math:`4` and :math:`2`, the ``num_velocities`` attribute will have the value ``[3, 1]``. + """ + + num_grid_qubits: int + """ + The number of qubits required to encode the grid. + """ + + num_velocity_qubits: int + """ + The number of qubits required to encode the velocity discretization of the lattice. + """ + + num_ancilla_qubits: int + """ + The number of ancilla (non-velocity, non-grid) qubits required for the quantum circuit to simulate this lattice. + """ + num_total_qubits: int - registers: Tuple[QuantumRegister, ...] - logger: Logger + """ + The total number of qubits required for the quantum circuit to simulate the lattice. This is the sum of the number of grid, velocity, and ancilla qubits. + """ + + velocity_register: Tuple[QuantumRegister, ...] + """ + A tuple that holds registers responsible for specific operations of the QLBM algorithm. + """ circuit: QuantumCircuit - blocks: Dict[str, List[Block]] + """ + An empty ``qiskit.QuantumCircuit`` with labeled registers that quantum components use as a base. + Each quantum component that is parameterized by a :class:`.Lattice` makes a copy of this quantum circuit + to which it appends its designated logic. + """ + + shapes: Dict[str, List[Shape]] + """ + Contains all of the :class:`.Shape`s encoding the solid geometry of the lattice. The key of the dictionary is the specific kind of boundary condition of the obstacle (i.e., ``"bounceback"`` or ``"specular"``). + """ + + logger: Logger + """ + The performance logger, by default ``getLogger("qlbm")``. + """ + + register: Tuple[List[QuantumRegister], ...] + """ + A tuple of lists of :class:`qiskit.QuantumRegister` s that are used to store the quantum information of the lattice. + """ + + discretization: LatticeDiscretization + """ + The discretization of the lattice, as an enum value of :class:`.LatticeDiscretization`. + """ def __init__( self, @@ -100,7 +167,7 @@ def __init__( def parse_input_data( self, lattice_data: str | Dict, # type: ignore - ) -> Tuple[List[int], List[int], Dict[str, List[Block]]]: + ) -> Tuple[List[int], List[int], Dict[str, List[Shape]], LatticeDiscretization]: r""" Parses the lattice input data, provided in either a file path or a dictionary. @@ -112,11 +179,12 @@ def parse_input_data( Returns ------- - Tuple[List[int], List[int], Dict[str, List[Block]]] + Tuple[List[int], List[int], Dict[str, List[Shape]], LatticeDiscretization] A tuple containing (i) a list of the number of gridpoints per dimension, (ii) a list of the number of velicities per dimension, - and (iii) a dictionary containing the solid :class:`.Block`\ s. + (iii) a dictionary containing the solid :class:`.Shape`\ s, + and (iv) the discretization enum of the lattice. The key of the dictionary is the specific kind of boundary condition of the obstacle (i.e., ``"bounceback"`` or ``"specular"``). @@ -137,7 +205,7 @@ def parse_input_data( def __parse_input_dict( self, input_dict: Dict, # type: ignore - ) -> Tuple[List[int], List[int], Dict[str, List[Shape]]]: + ) -> Tuple[List[int], List[int], Dict[str, List[Shape]], LatticeDiscretization]: r""" Parses the lattice input data, provided as a dictionary. @@ -148,7 +216,7 @@ def __parse_input_dict( Returns ------- - Tuple[List[int], List[int], Dict[str, List[Shape]]] + Tuple[List[int], List[int], Dict[str, List[Shape]], LatticeDiscretization] A tuple containing (i) a list of the number of gridpoints per dimension, (ii) a list of the number of velicities per dimension, @@ -175,11 +243,6 @@ def __parse_input_dict( 'Lattice configuration missing "velocities" properties.' ) - if len(lattice_dict["dim"]) != len(lattice_dict["velocities"]): # type: ignore - raise LatticeException( - "Lattice configuration dimensionality is inconsistent." - ) - num_dimensions = len(lattice_dict["dim"]) # type: ignore if num_dimensions not in [1, 2, 3]: @@ -187,36 +250,60 @@ def __parse_input_dict( f"Only 1, 2, and 3-dimensional lattices are supported. Provided lattice has {len(lattice_dict['dim'])} dimensions." # type: ignore ) - # Check whether the number of grid points and velocities is compatible - for dim in range(num_dimensions): - dim_index = dimension_letter(dim) + grid_list: List[int] = [ + # -1 because the bit_length() would "overshoot" for powers of 2 + lattice_dict["dim"][dimension_letter(dim)] - 1 + for dim in range(num_dimensions) + ] + + discretization: LatticeDiscretization = LatticeDiscretization.CFLDISCRETIZATION + velocity_list: List[int] = [] - if not is_two_pow(lattice_dict["dim"][dim_index]): # type: ignore + # Check if velocities is a string (DdQq format) or dict + if isinstance(lattice_dict["velocities"], str): # type: ignore + # Parse DdQq format (e.g., "D2Q4" means 2 dimensions, 4 velocities total) + velocity_spec = lattice_dict["velocities"].upper() # type: ignore + if not velocity_spec.startswith("D") or "Q" not in velocity_spec: raise LatticeException( - f"Lattice {dim_index}-dimension has a number of grid points that is not divisible by 2." + f"Invalid velocity specification format: {lattice_dict['velocities']}. Expected format like 'd2q4'." ) - if not is_two_pow(lattice_dict["velocities"][dim_index]): # type: ignore + try: + parts = velocity_spec[1:].split("Q") + spec_dims = int(parts[0]) + total_velocities = int(parts[1]) + velocity_list = [] + + if spec_dims != num_dimensions: + raise LatticeException( + f"Velocity specification dimensions ({spec_dims}) do not match lattice dimensions ({num_dimensions})." + ) + + discretization = LatticeDiscretizationProperties.get_discretization( + num_dimensions, total_velocities + ) + + except (ValueError, IndexError): raise LatticeException( - f"Lattice {dim_index}-dimension has a number of velocities that is not divisible by 2." + f"Invalid velocity specification format: {lattice_dict['velocities']}. Expected format like 'DQ'." ) - # The lattice properties are ok - grid_list: List[int] = [ - # -1 because the bit_length() would "overshoot" for powers of 2 - lattice_dict["dim"][dimension_letter(dim)] - 1 - for dim in range(num_dimensions) - ] - velocity_list: List[int] = [ - # -1 because the bit_length() would "overshoot" for powers of 2 - lattice_dict["velocities"][dimension_letter(dim)] - 1 - for dim in range(num_dimensions) - ] + else: + if len(lattice_dict["dim"]) != len(lattice_dict["velocities"]): # type: ignore + raise LatticeException( + "Lattice configuration dimensionality is inconsistent." + ) + + velocity_list = [ + # -1 because the bit_length() would "overshoot" for powers of 2 + lattice_dict["velocities"][dimension_letter(dim)] - 1 + for dim in range(num_dimensions) + ] parsed_obstacles: Dict[str, List[Shape]] = {"specular": [], "bounceback": []} if "geometry" not in input_dict: - return grid_list, velocity_list, parsed_obstacles + return grid_list, velocity_list, parsed_obstacles, discretization geometry_list: List[Dict[str, List[int]]] = input_dict["geometry"] # type: ignore @@ -315,7 +402,7 @@ def __parse_input_dict( ) ) - return grid_list, velocity_list, parsed_obstacles + return grid_list, velocity_list, parsed_obstacles, discretization def to_json(self) -> str: """ @@ -335,14 +422,18 @@ def to_json(self) -> str: "velocities": { dimension_letter(dim): self.num_velocities[dim] + 1 for dim in range(self.num_dims) - }, + } + if self.discretization == LatticeDiscretization.CFLDISCRETIZATION + else LatticeDiscretizationProperties.string_representation[ + self.discretization + ], # type: ignore }, } lattice_dict["geometry"] = flatten( [ - [block.to_dict() for block in self.blocks[boundary_type]] - for boundary_type in self.blocks + [block.to_dict() for block in self.shapes[boundary_type]] + for boundary_type in self.shapes ] ) diff --git a/qlbm/lattice/lattices/collisionless_lattice.py b/qlbm/lattice/lattices/collisionless_lattice.py index 0c04974..d02724e 100644 --- a/qlbm/lattice/lattices/collisionless_lattice.py +++ b/qlbm/lattice/lattices/collisionless_lattice.py @@ -6,9 +6,9 @@ from qiskit import QuantumCircuit, QuantumRegister from typing_extensions import override -from qlbm.lattice.geometry.shapes.block import Block +from qlbm.lattice.geometry.shapes.base import Shape from qlbm.tools.exceptions import LatticeException -from qlbm.tools.utils import dimension_letter, flatten +from qlbm.tools.utils import dimension_letter, flatten, is_two_pow from .base import Lattice @@ -23,7 +23,7 @@ class CollisionlessLattice(Lattice): :attr:`num_dims` The number of dimensions of the lattice. :attr:`num_gridpoints` A ``List[int]`` of the number of gridpoints of the lattice in each dimension. **Important**\ : for easier compatibility with binary arithmetic, the number of gridpoints - specified in the input dicitionary is one larger than the one held in the ``Lattice``. + specified in the input dictionary is one larger than the one held in the ``Lattice``. That is, for a ``16x64`` lattice, the ``num_gridpoints`` attribute will have the value ``[15, 63]``. :attr:`num_grid_qubits` The total number of qubits required to encode the lattice grid. :attr:`num_velocity_qubits` The total number of qubits required to encode the velocity discretization of the lattice. @@ -34,7 +34,7 @@ class CollisionlessLattice(Lattice): :attr:`circuit` An empty ``qiskit.QuantumCircuit`` with labeled registers that quantum components use as a base. Each quantum component that is parameterized by a ``Lattice`` makes a copy of this quantum circuit to which it appends its designated logic. - :attr:`blocks` A ``Dict[str, List[Block]]`` that contains all of the :class:`.Block`\ s encoding the solid geometry of the lattice. + :attr:`shapes` A ``Dict[str, List[Shape]]`` that contains all of the :class:`.Shape`\ s encoding the solid geometry of the lattice. The key of the dictionary is the specific kind of boundary condition of the obstacle (i.e., ``"bounceback"`` or ``"specular"``). :attr:`logger` The performance logger, by default ``getLogger("qlbm")``. =========================== ====================================================================== @@ -155,13 +155,26 @@ def __init__( logger: Logger = getLogger("qlbm"), ) -> None: super().__init__(lattice_data, logger) - dimensions, velocities, blocks = self.parse_input_data(lattice_data) # type: ignore + dimensions, velocities, shapes, self.discretization = self.parse_input_data( + lattice_data + ) # type: ignore self.num_dims = len(dimensions) self.num_gridpoints = dimensions + + for dim in range(self.num_dims): + if not is_two_pow(self.num_gridpoints[dim] + 1): # type: ignore + raise LatticeException( + f"Lattice has a number of grid points that is not divisible by 2 in dimension {dimension_letter(dim)}." + ) + if not is_two_pow(velocities[dim] + 1): # type: ignore + raise LatticeException( + f"Lattice has a number of velocities that is not divisible by 2 in dimension {dimension_letter(dim)}." + ) + self.num_velocities = velocities - self.blocks: Dict[str, List[Block]] = blocks - self.block_list: List[Block] = flatten(list(blocks.values())) + self.shapes: Dict[str, List[Shape]] = shapes # type: ignore + self.shape_list: List[Shape] = flatten(list(shapes.values())) self.num_comparator_qubits = 2 * (self.num_dims - 1) self.num_obstacle_qubits = self.__num_obstacle_qubits() self.num_ancilla_qubits = ( @@ -495,8 +508,8 @@ def get_registers(self) -> Tuple[List[QuantumRegister], ...]: def __num_obstacle_qubits(self) -> int: all_obstacle_bounceback: bool = len( - [b for b in self.block_list if b.boundary_condition == "bounceback"] - ) == len(self.block_list) + [b for b in self.shape_list if b.boundary_condition == "bounceback"] + ) == len(self.shape_list) if all_obstacle_bounceback: # A single qubit suffices to determine # Whether particles have streamed inside the object @@ -508,13 +521,13 @@ def __num_obstacle_qubits(self) -> int: @override def __str__(self) -> str: - return f"[Lattice with {self.num_gridpoints} gps, {self.num_velocities} vels, and {str(self.blocks)} blocks with {self.num_total_qubits} qubits]" + return f"[Lattice with {self.num_gridpoints} gps, {self.num_velocities} vels, and {str(self.shapes)} shapes with {self.num_total_qubits} qubits]" @override def logger_name(self) -> str: gp_string = "" for c, gp in enumerate(self.num_gridpoints): - gp_string += f"{gp+1}" + gp_string += f"{gp + 1}" if c < len(self.num_gridpoints) - 1: gp_string += "x" - return f"{self.num_dims}d-{gp_string}-{len(self.block_list)}-obstacle" + return f"{self.num_dims}d-{gp_string}-{len(self.shape_list)}-obstacle" diff --git a/qlbm/lattice/lattices/lqlga_lattice.py b/qlbm/lattice/lattices/lqlga_lattice.py new file mode 100644 index 0000000..41743c2 --- /dev/null +++ b/qlbm/lattice/lattices/lqlga_lattice.py @@ -0,0 +1,284 @@ +"""Implementation of the :class:`.Lattice` base specific to the 2D and 3D :class:`.LQLGA` algorithm.""" + +from itertools import product +from math import prod +from typing import Dict, List, Tuple, cast, override + +from qiskit import QuantumCircuit, QuantumRegister + +from qlbm.lattice.geometry.shapes.base import Shape +from qlbm.lattice.lattices.base import Lattice +from qlbm.lattice.spacetime.properties_base import ( + LatticeDiscretization, + LatticeDiscretizationProperties, +) +from qlbm.tools.exceptions import LatticeException + + +class LQLGALattice(Lattice): + r""" + Lattice class for the :class:`.LQLGA` algorithm. + + ================================= ======================================================================================== + Attribute Summary + ================================= ======================================================================================== + :attr:`num_gridpoints` The number of gridpoints in each dimension of the lattice. + :attr:`num_velocities` The number of discrete velocities in each dimension of the lattice. + :attr:`num_dims` The number of dimensions of the lattice. + :attr:`discretization` The discretization of the lattice. + :attr:`num_velocities_per_point` The number of discrete velocities per gridpoint. + :attr:`num_base_qubits` The number of qubits required to represent the lattice without velocities. + :attr:`num_total_qubits` The total number of qubits required to represent the lattice, including velocities. + :attr:`registers` The list of quantum registers for the lattice, one for each gridpoint. + :attr:`circuit` The quantum circuit representing the lattice, initialized with the registers. + ================================ ======================================================================================== + + The registers encoded in the lattice and their accessors are given below. + For the size of each register, + and :math:`d` is the total number of dimensions: 2 or 3. + + .. list-table:: Register allocation + :widths: 25 25 25 50 + :header-rows: 1 + + * - Register + - Size + - Access Method + - Description + * - :attr:`velocity_register` + - :math:`N_g \cdot q` + - :meth:`velocity_index_flat` and `:meth:`velocity_index_tuple` + - :math:`N_g` registers sized according to the number of discrete velocities of the lattice. + """ + + discretization: LatticeDiscretization + """The discretization of the lattice, one of :class:`.LatticeDiscretization`.""" + + num_gridpoints: List[int] + """The number of gridpoints in each dimension of the lattice. + **Important** : for easier compatibility with binary arithmetic, the number of gridpoints + specified in the input dictionary is one larger than the one held in the ``Lattice``.""" + + shapes: Dict[str, List[Shape]] + """The shapes of the lattice, which are used to define the geometry of the lattice. + The key consists of the type of the shape and the name of the shape, e.g. "bounceback" or "specular". + """ + + num_base_qubits: int + """The number of qubits required to represent the lattice.""" + + velocity_register: QuantumRegister + """The quantum register representing the velocities of the lattice.""" + + def __init__(self, lattice_data, logger=...): + super().__init__(lattice_data, logger) + + self.num_gridpoints, self.num_velocities, self.shapes, self.discretization = ( + self.parse_input_data(lattice_data) + ) # type: ignore + self.num_dims = len(self.num_gridpoints) + self.num_velocities_per_point = ( + LatticeDiscretizationProperties.get_num_velocities(self.discretization) + ) + + self.num_base_qubits = ( + prod(map(lambda x: x + 1, self.num_gridpoints)) + * self.num_velocities_per_point + ) + + self.num_total_qubits = self.num_base_qubits + + self.registers = self.get_registers() + + self.velocity_register = self.registers + + self.circuit = QuantumCircuit(*self.registers) + + @override + def get_registers(self) -> Tuple[List[QuantumRegister], ...]: + velocity_registers = [ + QuantumRegister( + LatticeDiscretizationProperties.get_num_velocities(self.discretization), + name=rf"v^{{{gp_tuple}}}", + ) + for gp_tuple in list( + product(*[range(ng + 1) for ng in self.num_gridpoints]) + ) + ] + + return velocity_registers # type: ignore + + def gridpoint_index_tuple(self, gridpoint: Tuple[int, ...]) -> int: + """ + Get the lexicographic index of a gridpoint in the lattice. + + Parameters + ---------- + gridpoint : Tuple[int, ...] + The gridpoint formatted as a tuple of indices for each dimension. + + Returns + ------- + int + The lexicographic index of the gridpoint in the lattice. + + Raises + ------ + LatticeException + If the gridpoint index is out of bounds for the lattice. + """ + if len(gridpoint) != self.num_dims: + raise LatticeException( + f"Gridpoint {gridpoint} has incorrect number of dimensions. Expected {self.num_dims}, got {len(gridpoint)}." + ) + if any( + gp < 0 or gp >= ng + 1 for gp, ng in zip(gridpoint, self.num_gridpoints) + ): + raise LatticeException( + f"Gridpoint {gridpoint} is out of bounds for the lattice with gridpoints {self.num_gridpoints}." + ) + + flat_index = 0 + multiplier = 1 + num_dims = len(self.num_gridpoints) + + for dim in reversed(range(num_dims)): + flat_index += gridpoint[dim] * multiplier + multiplier *= self.num_gridpoints[dim] + 1 + + return flat_index + + def gridpoint_index_flat(self, gridpoint: int) -> Tuple[int, ...]: + """ + Get the tuple representation of a gridpoint index in the lattice. + + Parameters + ---------- + gridpoint : int + The lexicographic index of the gridpoint in the lattice. + + Returns + ------- + Tuple[int, ...] + The tuple representation of the gridpoint index in the lattice. + + Raises + ------ + LatticeException + If the gridpoint index is out of bounds for the lattice. + """ + if ( + gridpoint < 0 + or gridpoint >= self.num_total_qubits // self.num_velocities_per_point + ): + raise LatticeException( + f"Gridpoint {gridpoint} is out of bounds for the lattice with {self.num_total_qubits // self.num_velocities_per_point} gridpoints." + ) + indices = [] + for size in reversed(self.num_gridpoints): + indices.append(gridpoint % (size + 1)) + gridpoint //= size + 1 + return cast(Tuple[int, ...], tuple(reversed(indices))) + + def velocity_index_flat(self, gridpoint: int, velocity: int) -> int: + """ + Get the index of a qubit representing a particular velocity channel of a given gridpoint. + + Parameters + ---------- + gridpoint : int + The lexicographic index of the gridpoint in the lattice. + velocity : int + The index of the velocity channel (0-indexed). + + Returns + ------- + int + The index of the qubit representing the velocity channel at the specified gridpoint. + + Raises + ------ + LatticeException + If the gridpoint or velocity indices are out of bounds for the lattice. + """ + if velocity < 0 or velocity >= self.num_velocities_per_point: + raise LatticeException( + f"Velocity {velocity} is out of bounds for the lattice with {self.num_velocities_per_point} velocities per point." + ) + + if ( + gridpoint < 0 + or gridpoint >= self.num_total_qubits // self.num_velocities_per_point + ): + raise LatticeException( + f"Gridpoint {gridpoint} is out of bounds for the lattice with {self.num_total_qubits // self.num_velocities_per_point} gridpoints." + ) + return gridpoint * self.num_velocities_per_point + velocity + + def velocity_index_tuple(self, gridpoint: Tuple[int, ...], velocity: int) -> int: + """ + Get the index of a qubit representing a particular velocity channel of a given gridpoint. + + Parameters + ---------- + gridpoint : Tuple[int, ...] + The gridpoint formatted as a tuple of indices for each dimension. + velocity : int + The index of the velocity channel (0-indexed). + + Returns + ------- + int + The index of the qubit representing the velocity channel at the specified gridpoint. + + Raises + ------ + LatticeException + If the gridpoint or velocity indices are out of bounds for the lattice. + """ + if velocity < 0 or velocity >= self.num_velocities_per_point: + raise LatticeException( + f"Velocity {velocity} is out of bounds for the lattice with {self.num_velocities_per_point} velocities per point." + ) + return ( + self.gridpoint_index_tuple(gridpoint) * self.num_velocities_per_point + + velocity + ) + + def get_velocity_qubits_of_line(self, line_index: int) -> Tuple[int, int]: + r""" + Returns the velocity qubits of the positive and negative directions of a streaming line. + + Assumes that the lattice follows a :math:`D_{d}Q_{q}` discretization, where if :math:`q` is even, there are + :math:`\lceil \frac{q}{2} \rceil` streaming lines. This is to be generalized in the future. + + Parameters + ---------- + line_index : int + The index of the line to get the velocity qubits for. Counted from 0 according to the regular discretization taxonomy. + + Returns + ------- + List[int] + The list of velocity qubits for the specified line. + """ + if line_index < 0 or line_index > self.num_velocities_per_point // 2: + raise LatticeException( + f"Streaming Line index {line_index} is out of bounds for the lattice with {self.num_velocities_per_point // 2} lines." + ) + + return ( + (self.num_velocities_per_point % 2) + line_index, + (self.num_velocities_per_point % 2) + + line_index + + self.num_velocities_per_point // 2, + ) + + @override + def logger_name(self) -> str: + gp_string = "" + for c, gp in enumerate(self.num_gridpoints): + gp_string += f"{gp + 1}" + if c < len(self.num_gridpoints) - 1: + gp_string += "x" + return f"lqlga-d{self.num_dims}-q{self.num_velocities_per_point}-{gp_string}" diff --git a/qlbm/lattice/lattices/spacetime_lattice.py b/qlbm/lattice/lattices/spacetime_lattice.py index 8c8d052..e7099e2 100644 --- a/qlbm/lattice/lattices/spacetime_lattice.py +++ b/qlbm/lattice/lattices/spacetime_lattice.py @@ -9,9 +9,13 @@ from qlbm.lattice.lattices.base import Lattice from qlbm.lattice.spacetime.d1q2 import D1Q2SpaceTimeLatticeBuilder from qlbm.lattice.spacetime.d2q4 import D2Q4SpaceTimeLatticeBuilder -from qlbm.lattice.spacetime.properties_base import SpaceTimeLatticeBuilder +from qlbm.lattice.spacetime.d3q6 import D3Q6SpaceTimeLatticeBuilder +from qlbm.lattice.spacetime.properties_base import ( + LatticeDiscretization, + SpaceTimeLatticeBuilder, +) from qlbm.tools.exceptions import LatticeException -from qlbm.tools.utils import flatten +from qlbm.tools.utils import flatten, is_two_pow class SpaceTimeLattice(Lattice): @@ -84,10 +88,7 @@ class SpaceTimeLattice(Lattice): "x": 16, "y": 16 }, - "velocities": { - "x": 2, - "y": 2 - } + "velocities": "D2Q4" }, "geometry": [] } @@ -102,7 +103,7 @@ class SpaceTimeLattice(Lattice): SpaceTimeLattice( num_timesteps=1, lattice_data={ - "lattice": {"dim": {"x": 4, "y": 8}, "velocities": {"x": 2, "y": 2}}, + "lattice": {"dim": {"x": 4, "y": 8}, "velocities": "D2Q4"}, "geometry": [], } ).circuit.draw("mpl") @@ -122,10 +123,17 @@ def __init__( self.include_measurement_qubit = include_measurement_qubit self.use_volumetric_ops = use_volumetric_ops - self.num_gridpoints, self.num_velocities, self.blocks = self.parse_input_data( - lattice_data + self.num_gridpoints, self.num_velocities, self.shapes, self.discretization = ( + self.parse_input_data(lattice_data) ) # type: ignore self.num_dims = len(self.num_gridpoints) + + for dim in range(self.num_dims): + if not is_two_pow(self.num_gridpoints[dim] + 1): # type: ignore + raise LatticeException( + f"Lattice has a number of grid points that is not divisible by 2 in dimension {dim}." + ) + self.num_timesteps = num_timesteps self.properties: SpaceTimeLatticeBuilder = self.__get_builder() @@ -151,34 +159,35 @@ def __init__( logger.info(self.__str__()) def __get_builder(self) -> SpaceTimeLatticeBuilder: - if self.num_dims == 1: - if self.num_velocities[0] == 1: - return D1Q2SpaceTimeLatticeBuilder( - self.num_timesteps, - self.num_gridpoints, - include_measurement_qubit=self.include_measurement_qubit, - use_volumetric_ops=self.use_volumetric_ops, - logger=self.logger, - ) - raise LatticeException( - f"Unsupported number of velocities for 1D: {self.num_velocities[0] + 1}. Only D1Q2 is supported at the moment." + if self.discretization == LatticeDiscretization.D1Q2: + return D1Q2SpaceTimeLatticeBuilder( + self.num_timesteps, + self.num_gridpoints, + include_measurement_qubit=self.include_measurement_qubit, + use_volumetric_ops=self.use_volumetric_ops, + logger=self.logger, ) - if self.num_dims == 2: - if self.num_velocities[0] == 1 and self.num_velocities[1] == 1: - return D2Q4SpaceTimeLatticeBuilder( - self.num_timesteps, - self.num_gridpoints, - include_measurement_qubit=self.include_measurement_qubit, - use_volumetric_ops=self.use_volumetric_ops, - logger=self.logger, - ) - raise LatticeException( - f"Unsupported number of velocities for 2D: {(self.num_velocities[0] + 1, self.num_velocities[1] + 1)}. Only D2Q4 is supported at the moment." + if self.discretization == LatticeDiscretization.D2Q4: + return D2Q4SpaceTimeLatticeBuilder( + self.num_timesteps, + self.num_gridpoints, + include_measurement_qubit=self.include_measurement_qubit, + use_volumetric_ops=self.use_volumetric_ops, + logger=self.logger, + ) + + if self.discretization == LatticeDiscretization.D3Q6: + return D3Q6SpaceTimeLatticeBuilder( + self.num_timesteps, + self.num_gridpoints, + include_measurement_qubit=self.include_measurement_qubit, + use_volumetric_ops=self.use_volumetric_ops, + logger=self.logger, ) raise LatticeException( - "Only 1D and 2D discretizations are currently available." + f"{self.discretization} not currently implemented for the Space-Time method. Only D1Q2 and D2Q4 are fully at the moment. D3Q6 only supports collision." ) def grid_index(self, dim: int | None = None) -> List[int]: @@ -400,7 +409,7 @@ def is_inside_an_obstacle(self, gridpoint: Tuple[int, ...]) -> bool: """ return any( block.contains_gridpoint(gridpoint) - for block in flatten(self.blocks.values()) + for block in flatten(self.shapes.values()) ) @override diff --git a/qlbm/lattice/spacetime/d1q2.py b/qlbm/lattice/spacetime/d1q2.py index 5244be4..546da08 100644 --- a/qlbm/lattice/spacetime/d1q2.py +++ b/qlbm/lattice/spacetime/d1q2.py @@ -1,4 +1,4 @@ -""":math:`D_2Q_4` STQBM builder.""" +""":math:`D_1Q_2` STQBM builder.""" from logging import Logger, getLogger from typing import Dict, List, Tuple, cast @@ -16,7 +16,7 @@ class D1Q2SpaceTimeLatticeBuilder(SpaceTimeLatticeBuilder): - """:math:`D_2Q_4` STQBM builder.""" + """:math:`D_1Q_2` STQBM builder.""" # Points to the right of the origin are marked as 0 # And points to the left are marked as 1 diff --git a/qlbm/lattice/spacetime/d2q4.py b/qlbm/lattice/spacetime/d2q4.py index 6a16103..6643c50 100644 --- a/qlbm/lattice/spacetime/d2q4.py +++ b/qlbm/lattice/spacetime/d2q4.py @@ -1,6 +1,5 @@ """:math:`D_2Q_4` STQBM builder.""" - from itertools import product from logging import Logger, getLogger from typing import Dict, List, Tuple, cast @@ -367,7 +366,7 @@ def coordinates_to_quadrant(self, distance: Tuple[int, int]) -> int: def get_reflected_index_of_velocity(self, velocity_index: int) -> int: if velocity_index not in list(range(4)): raise LatticeException( - f"D1Q2 discretization only supports velocities with indices {list(range(4))}. Index {velocity_index} unsupported." + f"D2Q4 discretization only supports velocities with indices {list(range(4))}. Index {velocity_index} unsupported." ) return self.velocity_reflection[velocity_index] @@ -376,7 +375,7 @@ def get_reflected_index_of_velocity(self, velocity_index: int) -> int: def get_reflection_increments(self, velocity_index: int) -> Tuple[int, ...]: if velocity_index not in list(range(4)): raise LatticeException( - f"D1Q2 discretization only supports velocities with indices {list(range(4))}. Index {velocity_index} unsupported." + f"D2Q4 discretization only supports velocities with indices {list(range(4))}. Index {velocity_index} unsupported." ) return self.extreme_point_classes[velocity_index][1] diff --git a/qlbm/lattice/spacetime/d3q6.py b/qlbm/lattice/spacetime/d3q6.py new file mode 100644 index 0000000..c6baf84 --- /dev/null +++ b/qlbm/lattice/spacetime/d3q6.py @@ -0,0 +1,158 @@ +""":math:`D_3Q_6` STQBM builder.""" + +from itertools import product +from logging import Logger, getLogger +from typing import Dict, List, Tuple + +from qiskit import QuantumRegister +from typing_extensions import override + +from qlbm.lattice.spacetime.properties_base import ( + LatticeDiscretization, + SpaceTimeLatticeBuilder, +) +from qlbm.tools.exceptions import LatticeException +from qlbm.tools.utils import dimension_letter + + +class D3Q6SpaceTimeLatticeBuilder(SpaceTimeLatticeBuilder): + """ + :math:`D_3Q_6` lattice builder. + + WIP. + """ + + velocity_reflection: Dict[int, int] = {0: 3, 1: 4, 2: 5, 3: 0, 4: 1, 5: 2} + + def __init__( + self, + num_timesteps: int, + num_gridpoints: List[int], + include_measurement_qubit: bool = False, + use_volumetric_ops: bool = False, + logger: Logger = getLogger("qlbm"), + ) -> None: + super().__init__( + num_timesteps, + include_measurement_qubit=include_measurement_qubit, + use_volumetric_ops=use_volumetric_ops, + logger=logger, + ) + self.num_gridpoints = num_gridpoints + + @override + def get_discretization(self) -> LatticeDiscretization: + return LatticeDiscretization.D3Q6 + + @override + def get_num_velocities_per_point(self) -> int: + return 6 + + @override + def get_num_ancilla_qubits(self) -> int: + return (6 if self.use_volumetric_ops else 0) + ( + 1 if self.include_measurement_qubit else 0 + ) + + @override + def get_num_grid_qubits(self) -> int: + return sum( + num_gridpoints_in_dim.bit_length() + for num_gridpoints_in_dim in self.num_gridpoints + ) + + @override + def get_num_previous_grid_qubits(self, dim: int) -> int: + # ! TODO add exception + return sum(self.num_gridpoints[i].bit_length() for i in range(dim)) + + @override + def get_num_velocity_qubits(self, num_timesteps: int | None = None) -> int: + total_gridpoints = ( + (sum(self.num_gridpoints) + len(self.num_gridpoints)) + * (sum(self.num_gridpoints) + len(self.num_gridpoints)) + * (sum(self.num_gridpoints) + len(self.num_gridpoints)) + ) + velocities_per_gp = self.get_num_velocities_per_point() + + if num_timesteps is None: + num_timesteps = self.num_timesteps + + return min( + total_gridpoints * velocities_per_gp, + int( + 8 * velocities_per_gp * velocities_per_gp * velocities_per_gp + + 12 * velocities_per_gp * velocities_per_gp + + 16 * velocities_per_gp + + velocities_per_gp + ), + ) + + @override + def get_registers(self) -> Tuple[List[QuantumRegister], ...]: + # Grid qubits + grid_registers = [ + QuantumRegister(gp.bit_length(), name=f"g_{dimension_letter(c)}") + for c, gp in enumerate(self.num_gridpoints) + ] + + # Velocity qubits + velocity_registers = [ + QuantumRegister( + self.get_num_velocity_qubits( + self.num_timesteps, + ), # The number of velocity qubits required at time t + name="v", + ) + ] + + ancilla_measurement_register = ( + [QuantumRegister(1, "a_m")] if self.include_measurement_qubit else [] + ) + + ancilla_comparator_registers = ( + [ + QuantumRegister(1, f"a_{bound}{dimension_letter(dim)}") + for dim, bound in product([0, 1], ["l", "u"]) + ] + if self.use_volumetric_ops + else [] + ) + + # Ancilla qubits + ancilla_registers = ancilla_measurement_register + ancilla_comparator_registers + + return (grid_registers, velocity_registers, ancilla_registers) + + @override + def get_index_of_neighbor(self, distance: Tuple[int, ...]) -> int: + raise LatticeException("Not implemented") + + @override + def get_streaming_lines( + self, dimension: int, direction: bool, timestep: int | None = None + ) -> List[List[int]]: + raise LatticeException("Not implemented") + + @override + def get_neighbor_indices(self): + # ! TODO!!! + return [], [] + + @override + def get_reflected_index_of_velocity(self, velocity_index: int) -> int: + if velocity_index not in list(range(4)): + raise LatticeException( + f"D3Q6 discretization only supports velocities with indices 0-5. Index {velocity_index} unsupported." + ) + + return self.velocity_reflection[velocity_index] + + @override + def get_reflection_increments(self, velocity_index: int) -> Tuple[int, ...]: + if velocity_index not in list(range(4)): + raise LatticeException( + f"D3Q6 discretization only supports velocities with indices 0-5. Index {velocity_index} unsupported." + ) + + raise LatticeException("Not implemented") diff --git a/qlbm/lattice/spacetime/properties_base.py b/qlbm/lattice/spacetime/properties_base.py index 1bec944..0d7d29a 100644 --- a/qlbm/lattice/spacetime/properties_base.py +++ b/qlbm/lattice/spacetime/properties_base.py @@ -10,6 +10,7 @@ from logging import Logger, getLogger from typing import Dict, List, Tuple, cast +import numpy as np from qiskit import QuantumRegister @@ -20,8 +21,188 @@ class LatticeDiscretization(Enum): The only supported discretizations currently are D1Q2 and D2Q4. """ + CFLDISCRETIZATION = (0,) D1Q2 = (1,) D2Q4 = (2,) + D3Q6 = (3,) + D1Q3 = (4,) + + +class LatticeDiscretizationProperties: + """ + Class containing properties of the lattice discretization in the :math:`D_dQ_q` taxonomy. + + .. list-table:: Attributes + :widths: 25 50 + :header-rows: 1 + + * - Attribute + - Description + * - :attr:`num_velocities` + - The number of velocities in the discretization. Stored as a ``Dict[LatticeDiscretization, int]``. + * - :attr:`velocity_vectors` + - The velocity profile each of the :math:`q` velocity channels. Each vector is :math:`d`-dimensional and each entry represents the velocity component of the channel in a particular dimension. Stored as a ``Dict[LatticeDiscretization, numpy.ndarray]``. + """ + + velocity_vectors: Dict[LatticeDiscretization, np.ndarray] = { + LatticeDiscretization.D1Q2: np.array([[1], [-1]]), + LatticeDiscretization.D2Q4: np.array([[1, 0], [0, 1], [-1, 0], [0, -1]]), + LatticeDiscretization.D3Q6: np.array( + [[1, 0, 0], [0, 1, 0], [0, 0, 1], [-1, 0, 0], [0, -1, 0], [0, 0, -1]] + ), + LatticeDiscretization.D1Q3: np.array([[0], [1], [-1]]), + } + + channel_masses: Dict[LatticeDiscretization, np.ndarray] = { + LatticeDiscretization.D1Q2: np.ones(2), + LatticeDiscretization.D2Q4: np.ones(4), + LatticeDiscretization.D3Q6: np.ones(6), + LatticeDiscretization.D1Q3: np.concatenate((np.array([2]), np.ones(2))), + } + + num_velocities: Dict[LatticeDiscretization, int] = { + LatticeDiscretization.D1Q2: 2, + LatticeDiscretization.D2Q4: 4, + LatticeDiscretization.D3Q6: 6, + LatticeDiscretization.D1Q3: 3, + } + + num_dimensions: Dict[LatticeDiscretization, int] = { + LatticeDiscretization.D1Q2: 1, + LatticeDiscretization.D2Q4: 2, + LatticeDiscretization.D3Q6: 3, + LatticeDiscretization.D1Q3: 1, + } + + string_representation: Dict[LatticeDiscretization, str] = { + LatticeDiscretization.D1Q2: "D1Q2", + LatticeDiscretization.D2Q4: "D2Q4", + LatticeDiscretization.D3Q6: "D3Q6", + LatticeDiscretization.D1Q3: "D1Q3", + } + + @staticmethod + def get_velocity_vectors( + discretization: LatticeDiscretization, + ) -> np.ndarray: + """ + Get the velocity vectors for a given discretization. + + Parameters + ---------- + discretization : LatticeDiscretization + The discretization for which to get the velocity vectors. + + Returns + ------- + np.ndarray + The velocity vectors corresponding to the discretization. + """ + return LatticeDiscretizationProperties.velocity_vectors[discretization] + + @staticmethod + def get_num_velocities( + discretization: LatticeDiscretization, + ) -> int: + """ + Get the number of velocities for a given discretization. + + Parameters + ---------- + discretization : LatticeDiscretization + The discretization for which to get the number of velocities. + + Returns + ------- + int + The number of velocities corresponding to the discretization. + """ + if discretization not in LatticeDiscretizationProperties.num_velocities: + raise ValueError(f"Discretization {discretization} is not supported.") + + return LatticeDiscretizationProperties.num_velocities[discretization] + + @staticmethod + def get_channel_masses( + discretization: LatticeDiscretization, + ) -> np.ndarray: + """ + Get the channel masses for a given discretization. + + Parameters + ---------- + discretization : LatticeDiscretization + The discretization for which to get the channel masses. + + Returns + ------- + np.ndarray + The channel masses corresponding to the discretization. + """ + return LatticeDiscretizationProperties.channel_masses[discretization] + + @staticmethod + def get_discretizations_of_dimensionality( + num_dimensions: int, + ) -> List[LatticeDiscretization]: + """ + Get all discretizations with a given number of dimensions. + + Parameters + ---------- + num_dimensions : int + The number of dimensions to filter by. + + Returns + ------- + List[LatticeDiscretization] + A list of discretizations with the specified number of dimensions. + """ + return [ + d + for d in LatticeDiscretization + if d != LatticeDiscretization.CFLDISCRETIZATION + and LatticeDiscretizationProperties.num_dimensions[d] == num_dimensions + ] + + @staticmethod + def get_discretization( + num_dimensions: int, num_velocities: int + ) -> LatticeDiscretization: + """ + Get the discretization for a given number of dimensions and velocities. + + Parameters + ---------- + num_dimensions : int + The number of dimensions. + num_velocities : int + The number of velocities. + + Returns + ------- + LatticeDiscretization + The discretization corresponding to the given parameters. + """ + compatible_discretizations = ( + LatticeDiscretizationProperties.get_discretizations_of_dimensionality( + num_dimensions + ) + ) + + compatible_discretizations = [ + discretization + for discretization in compatible_discretizations + if LatticeDiscretizationProperties.get_num_velocities(discretization) + == num_velocities + ] + + if not compatible_discretizations: + raise ValueError( + f"No discretization found for {num_dimensions} dimensions and {num_velocities} velocities." + ) + + return compatible_discretizations[0] class VonNeumannNeighborType(Enum): diff --git a/test/unit/blocks_3d_test.py b/test/unit/blocks_3d_test.py index babb025..d3353eb 100644 --- a/test/unit/blocks_3d_test.py +++ b/test/unit/blocks_3d_test.py @@ -38,8 +38,6 @@ def test_3d_mesh(simple_3d_block): assert all(a == b for a, b in zip(vertex_ull, mesh.vectors[1][1])) assert all(a == b for a, b in zip(vertex_ulu, mesh.vectors[1][2])) - # TODO: complete the remaining 5 faces - def test_3d_inner_x_wall_points(simple_3d_block): assert len(simple_3d_block.walls_inside) == 3 diff --git a/test/unit/collision/__init__.py b/test/unit/collision/__init__.py new file mode 100644 index 0000000..8b13789 --- /dev/null +++ b/test/unit/collision/__init__.py @@ -0,0 +1 @@ + diff --git a/test/unit/collision/eqc_collision_test.py b/test/unit/collision/eqc_collision_test.py new file mode 100644 index 0000000..2e09a7c --- /dev/null +++ b/test/unit/collision/eqc_collision_test.py @@ -0,0 +1,180 @@ +from typing import List + +import pytest +from qiskit import QuantumCircuit, transpile +from qiskit_aer import AerSimulator + +from qlbm.components.spacetime.collision.eqc_collision import ( + EQCCollisionOperator, +) +from qlbm.lattice.eqc.eqc_generator import ( + EquivalenceClassGenerator, +) +from qlbm.lattice.lattices.spacetime_lattice import SpaceTimeLattice +from qlbm.lattice.spacetime.properties_base import LatticeDiscretization +from qlbm.tools.utils import flatten + + +def int_to_bool_list(num, num_bits): + return bit_string_to_bool_list(format(num, f"0{num_bits}b")) + + +def bit_string_to_bool_list(bitstring): + return [x == "1" for x in bitstring] + + +@pytest.fixture +def d2q4_equivalence_class_bitstrings() -> List[List[str]]: + return [ + eqc.get_bitstrings() + for eqc in EquivalenceClassGenerator( + LatticeDiscretization.D2Q4 + ).generate_equivalence_classes() + ] + + +@pytest.fixture +def d3q6_equivalence_class_bitstrings() -> List[List[str]]: + return sorted( + [ + eqc.get_bitstrings() + for eqc in EquivalenceClassGenerator( + LatticeDiscretization.D3Q6 + ).generate_equivalence_classes() + ], + key=lambda x: x[0], + ) # Sort by first element + + +def verify_simulation_outcome( + input_state: str, + expected_outcomes: List[str], + collision_circuit: QuantumCircuit, + num_shots=128, + verify_negative_cases: bool = True, +): + init_circuit = QuantumCircuit(collision_circuit.num_qubits) + # for c, b in enumerate(bit_string_to_bool_list(input_state[::-1])): + # if b: + # init_circuit.x(collision_circuit.num_qubits - c - 1) + for c, b in enumerate(bit_string_to_bool_list(input_state)): + if b: + init_circuit.x(c) + qc_final = init_circuit.compose(collision_circuit) + qc_final.measure_all() + sim = AerSimulator() + counts = sim.run(transpile(qc_final, sim), shots=num_shots).result().get_counts() + + assert all(b[::-1] in counts for b in expected_outcomes), ( + f"{expected_outcomes} not covered by counts {counts}" + ) + + if verify_negative_cases: + all_biststrings = [ + format(i, f"0{collision_circuit.num_qubits}b") + for i in range(2**collision_circuit.num_qubits) + ] + assert all( + b[::-1] not in counts for b in all_biststrings if b not in expected_outcomes + ) + + +@pytest.mark.parametrize( + "equivalence_class_index", + list( + range( + len( + EquivalenceClassGenerator( + LatticeDiscretization.D2Q4 + ).generate_equivalence_classes() + ) + ) + ), +) +def test_d2q4_collision_positive_cases( + d2q4_equivalence_class_bitstrings, equivalence_class_index +): + lattice = SpaceTimeLattice( + 1, + { + "lattice": {"dim": {"x": 4, "y": 4}, "velocities": "D2Q4"}, + "geometry": [], + }, + ) + + local_circuit = EQCCollisionOperator( + lattice.properties.get_discretization() + ).circuit + assert local_circuit.num_qubits == 4 + eqc = d2q4_equivalence_class_bitstrings[equivalence_class_index] + for velocity_cfg in eqc: + verify_simulation_outcome( + velocity_cfg, eqc, local_circuit, verify_negative_cases=True + ) + + +def test_d2q4_collision_negative_cases(d2q4_equivalence_class_bitstrings): + lattice = SpaceTimeLattice( + 1, + { + "lattice": {"dim": {"x": 4, "y": 4}, "velocities": "D2Q4"}, + "geometry": [], + }, + ) + + local_circuit = EQCCollisionOperator( + lattice.properties.get_discretization() + ).circuit + assert local_circuit.num_qubits == 4 + + for b in [ + format(i, f"0{local_circuit.num_qubits}b") + for i in range(2**local_circuit.num_qubits) + if format(i, f"0{local_circuit.num_qubits}b") + not in flatten(d2q4_equivalence_class_bitstrings) + ]: + verify_simulation_outcome( + b, + [b], + local_circuit, + verify_negative_cases=True, + ) + + +@pytest.mark.parametrize( + "equivalence_class_index", + list( + range( + len( + EquivalenceClassGenerator( + LatticeDiscretization.D3Q6 + ).generate_equivalence_classes() + ) + ) + ), +) +def test_d3q6_collision_positive_cases( + d3q6_equivalence_class_bitstrings, equivalence_class_index +): + lattice = SpaceTimeLattice( + 1, + { + "lattice": { + "dim": {"x": 2, "y": 2, "z": 2}, + "velocities": "D3Q6", + }, + "geometry": [], + }, + ) + + local_circuit = EQCCollisionOperator( + lattice.properties.get_discretization() + ).circuit + assert local_circuit.num_qubits == 6 + eqc = d3q6_equivalence_class_bitstrings[equivalence_class_index] + for velocity_cfg in eqc: + verify_simulation_outcome( + velocity_cfg, eqc, local_circuit, verify_negative_cases=True + ) + + # [['010011', '100101'], ['010110', '001101'], ['011011', '110110', '101101'], ['100100', '001001', '010010'], ['100110', '001011'], ['101100', '011010'], ['110010', '101001'], ['110100', '011001']] diff --git a/test/unit/collision/eqc_generator_test.py b/test/unit/collision/eqc_generator_test.py new file mode 100644 index 0000000..66934e3 --- /dev/null +++ b/test/unit/collision/eqc_generator_test.py @@ -0,0 +1,101 @@ +import numpy as np + +from qlbm.lattice.eqc.eqc import EquivalenceClass +from qlbm.lattice.eqc.eqc_generator import ( + EquivalenceClassGenerator, +) +from qlbm.lattice.spacetime.properties_base import LatticeDiscretization + + +def test_eqc_generator_d1q2(): + generator = EquivalenceClassGenerator(LatticeDiscretization.D1Q2) + eqcs = generator.generate_equivalence_classes() + assert len(eqcs) == 0 + +def test_eqc_generator_d1q3(): + generator = EquivalenceClassGenerator(LatticeDiscretization.D1Q3) + eqcs = generator.generate_equivalence_classes() + assert len(eqcs) == 1 + +def test_eqc_generator_d2q4(): + generator = EquivalenceClassGenerator(LatticeDiscretization.D2Q4) + eqcs = generator.generate_equivalence_classes() + eqcs_expected = set( + [ + EquivalenceClass( + LatticeDiscretization.D2Q4, + {(True, False, True, False), (False, True, False, True)}, + ) + ] + ) + assert len(eqcs) == 1 + assert eqcs_expected == eqcs + + +def test_eqc_generator_d3q6(): + generator = EquivalenceClassGenerator(LatticeDiscretization.D3Q6) + eqcs = generator.generate_equivalence_classes() + eqcs_expected = set( + [ + EquivalenceClass( + LatticeDiscretization.D3Q6, + { + (True, False, False, True, False, False), # 100100 + (False, True, False, False, True, False), # 010010 + (False, False, True, False, False, True), # 001001 + }, + ), + EquivalenceClass( + LatticeDiscretization.D3Q6, + { + (True, True, False, True, True, False), # 110110 + (True, False, True, True, False, True), # 101101 + (False, True, True, False, True, True), # 011011 + }, + ), + EquivalenceClass( + LatticeDiscretization.D3Q6, # 3 + { + (True, True, False, False, True, False), # 110010 + (True, False, True, False, False, True), # 101001 + }, + ), + EquivalenceClass( + LatticeDiscretization.D3Q6, # 4 + { + (False, True, False, True, True, False), # 010110 + (False, False, True, True, False, True), # 001101 + }, + ), + EquivalenceClass( + LatticeDiscretization.D3Q6, # 5 + { + (True, True, False, True, False, False), # 110100 + (False, True, True, False, False, True), # 011001 + }, + ), + EquivalenceClass( + LatticeDiscretization.D3Q6, # 6 + { + (True, False, False, True, True, False), # 100110 + (False, False, True, False, True, True), # 001011 + }, + ), + EquivalenceClass( + LatticeDiscretization.D3Q6, # 7 + { + (True, False, True, True, False, False), # 101100 + (False, True, True, False, True, False), # 011010 + }, + ), + EquivalenceClass( + LatticeDiscretization.D3Q6, # 8 + { + (False, True, False, False, True, True), # 010011 + (True, False, False, True, False, True), # 100101 + }, + ), + ] + ) + assert len(eqcs) == 8 + assert all(eqc in eqcs_expected for eqc in eqcs) diff --git a/test/unit/collision/eqc_permutation_test.py b/test/unit/collision/eqc_permutation_test.py new file mode 100644 index 0000000..e69de29 diff --git a/test/unit/collision/eqc_velocity_discretization_test.py b/test/unit/collision/eqc_velocity_discretization_test.py new file mode 100644 index 0000000..4f78d9e --- /dev/null +++ b/test/unit/collision/eqc_velocity_discretization_test.py @@ -0,0 +1,51 @@ +from typing import Set, Tuple + +import pytest + +from qlbm.lattice.eqc.eqc import EquivalenceClass +from qlbm.lattice.eqc.eqc_generator import ( + EquivalenceClassGenerator, +) +from qlbm.lattice.spacetime.properties_base import LatticeDiscretization +from qlbm.tools.exceptions import LatticeException + + +def test_equivalence_class_bad_number_of_velocity_configurations(): + velocity_profile: Set[Tuple[bool, ...]] = {(True, True, True, True)} + with pytest.raises(LatticeException) as excinfo: + EquivalenceClass(LatticeDiscretization.D2Q4, velocity_profile) + assert ( + f"Equivalence class must have at least two configurations. Provided velocity profiles are {velocity_profile}." + == str(excinfo.value) + ) + + +def test_equivalence_class_bad_number_velocity_configurations(): + velocity_profile: Set[Tuple[bool, ...]] = {(True, True, True), (True, True, True, False)} + with pytest.raises(LatticeException) as excinfo: + EquivalenceClass(LatticeDiscretization.D2Q4, velocity_profile) + assert ( + f"All velocity configurations must have length 4. Provided configurations are {velocity_profile}." + == str(excinfo.value) + ) + + +def test_equivalence_class_bad_mass(): + velocity_profile: Set[Tuple[bool, ...]] = {(True, False, True, False), (True, True, True, False)} + with pytest.raises(LatticeException) as excinfo: + EquivalenceClass(LatticeDiscretization.D2Q4, velocity_profile) + assert "Velocity configurations have different masses." == str(excinfo.value) + + +def test_equivalence_class_bad_momentum(): + velocity_profile: Set[Tuple[bool, ...]] = {(True, False, True, False), (True, True, False, False)} + with pytest.raises(LatticeException) as excinfo: + EquivalenceClass(LatticeDiscretization.D2Q4, velocity_profile) + assert "Velocity configurations have different momenta." == str(excinfo.value) + + +def test_equivalence_class_ok(): + velocity_profile: Set[Tuple[bool, ...]] = {(True, False, True, False), (False, True, False, True)} + eqc = EquivalenceClass(LatticeDiscretization.D2Q4, velocity_profile) + assert eqc.mass == 2 + assert eqc.momentum.tolist() == [0, 0] diff --git a/test/unit/collisionles_lattice_test.py b/test/unit/collisionles_lattice_test.py index 7e92a67..f69fa66 100644 --- a/test/unit/collisionles_lattice_test.py +++ b/test/unit/collisionles_lattice_test.py @@ -183,7 +183,7 @@ def test_lattice_exception_mismatched_bad_dimensions(): ) assert ( - "Lattice y-dimension has a number of grid points that is not divisible by 2." + "Lattice has a number of grid points that is not divisible by 2 in dimension y." == str(excinfo.value) ) @@ -200,7 +200,7 @@ def test_lattice_exception_mismatched_bad_velocities(): ) assert ( - "Lattice y-dimension has a number of velocities that is not divisible by 2." + "Lattice has a number of velocities that is not divisible by 2 in dimension y." == str(excinfo.value) ) diff --git a/test/unit/incrementer_test.py b/test/unit/incrementer_test.py index 03c7980..5d012ff 100644 --- a/test/unit/incrementer_test.py +++ b/test/unit/incrementer_test.py @@ -17,10 +17,10 @@ def lattice_symmetric_small_2d(): def test_3d_incrementer_size(lattice_asymmetric_medium_3d): inc: ControlledIncrementer = ControlledIncrementer(lattice_asymmetric_medium_3d) - assert inc.circuit.size() == 32 + assert inc.circuit.size() == 78 def test_2d_incrementer_size(lattice_symmetric_small_2d): inc: ControlledIncrementer = ControlledIncrementer(lattice_symmetric_small_2d) - assert inc.circuit.size() == 24 + assert inc.circuit.size() == 68 diff --git a/test/unit/lqlga/__init__.py b/test/unit/lqlga/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/test/unit/lqlga/circuits/__init__.py b/test/unit/lqlga/circuits/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/test/unit/lqlga/circuits/conftest.py b/test/unit/lqlga/circuits/conftest.py new file mode 100644 index 0000000..e94a2e1 --- /dev/null +++ b/test/unit/lqlga/circuits/conftest.py @@ -0,0 +1,39 @@ +import pytest + +from qlbm.lattice.lattices.lqlga_lattice import LQLGALattice + + +@pytest.fixture +def lattice_d2q4_256_8() -> LQLGALattice: + return LQLGALattice( + { + "lattice": { + "dim": {"x": 256, "y": 8}, + "velocities": "D2Q4", + }, + }, + ) + + +@pytest.fixture +def lattice_d1q2_256() -> LQLGALattice: + return LQLGALattice( + { + "lattice": { + "dim": {"x": 256}, + "velocities": "D1Q2", + }, + }, + ) + +@pytest.fixture +def lattice_d1q2_8() -> LQLGALattice: + return LQLGALattice( + { + "lattice": { + "dim": {"x": 8}, + "velocities": "D1Q2", + }, + }, + ) + diff --git a/test/unit/lqlga/circuits/streaming_test.py b/test/unit/lqlga/circuits/streaming_test.py new file mode 100644 index 0000000..444d474 --- /dev/null +++ b/test/unit/lqlga/circuits/streaming_test.py @@ -0,0 +1,50 @@ +from typing import List, Tuple + +import pytest +from qlbm.components.lqlga.streaming import LQLGAStreamingOperator + + +def verify_streaming_line_correctness( + size: int, swaps: List[List[Tuple[int, int]]] +) -> bool: + configuration = list(range(size)) + for layer in swaps: + for i, j in layer: + configuration[i], configuration[j] = ( + configuration[j], + configuration[i], + ) + + return configuration == [size - 1] + list(range(size - 1)) + + +@pytest.mark.parametrize("size", list(map(lambda x: 2**x, list(range(10))))) +def test_streaming_line_creation_power_of_2(lattice_d1q2_8, size): + operator = LQLGAStreamingOperator(lattice_d1q2_8) + swaps = operator.logarithmic_depth_streaming_line_swaps(size, False) + + assert verify_streaming_line_correctness(size, swaps) + + +@pytest.mark.parametrize("size", list(map(lambda x: 2**x + 1, list(range(1, 10))))) +def test_streaming_line_creation_off_power_of_2_positive(lattice_d1q2_8, size): + operator = LQLGAStreamingOperator(lattice_d1q2_8) + swaps = operator.logarithmic_depth_streaming_line_swaps(size, False) + + assert verify_streaming_line_correctness(size, swaps) + + +@pytest.mark.parametrize("size", list(map(lambda x: 2**x - 1, list(range(3, 10))))) +def test_streaming_line_creation__off_power_2_negative(lattice_d1q2_8, size): + operator = LQLGAStreamingOperator(lattice_d1q2_8) + swaps = operator.logarithmic_depth_streaming_line_swaps(size, False) + + assert verify_streaming_line_correctness(size, swaps) + + +@pytest.mark.parametrize("size", list(range(7, 16))) +def test_streaming_line_creation_intermediate(lattice_d1q2_8, size): + operator = LQLGAStreamingOperator(lattice_d1q2_8) + swaps = operator.logarithmic_depth_streaming_line_swaps(size, False) + + assert verify_streaming_line_correctness(size, swaps) diff --git a/test/unit/lqlga/conftest.py b/test/unit/lqlga/conftest.py new file mode 100644 index 0000000..0bb459e --- /dev/null +++ b/test/unit/lqlga/conftest.py @@ -0,0 +1,27 @@ +import pytest + +from qlbm.lattice.lattices.lqlga_lattice import LQLGALattice + + +@pytest.fixture +def lattice_d2q4_256_8() -> LQLGALattice: + return LQLGALattice( + { + "lattice": { + "dim": {"x": 256, "y": 8}, + "velocities": "D2Q4", + }, + }, + ) + + +@pytest.fixture +def lattice_d1q2_256() -> LQLGALattice: + return LQLGALattice( + { + "lattice": { + "dim": {"x": 256}, + "velocities": "D1Q2", + }, + }, + ) diff --git a/test/unit/lqlga/lqlga_lattice_test.py b/test/unit/lqlga/lqlga_lattice_test.py new file mode 100644 index 0000000..15daded --- /dev/null +++ b/test/unit/lqlga/lqlga_lattice_test.py @@ -0,0 +1,51 @@ +import pytest + +from qlbm.lattice.lattices.lqlga_lattice import LQLGALattice +from qlbm.lattice.spacetime.properties_base import LatticeDiscretizationProperties + + +def test_lqlga_lattice_num_registers_d1q2(lattice_d1q2_256): + assert len(lattice_d1q2_256.registers) == 256 + assert all(reg.size == 2 for reg in lattice_d1q2_256.registers) + assert lattice_d1q2_256.circuit.num_qubits == 512 + + +def test_lqlga_lattice_num_registers_d2q4(lattice_d2q4_256_8): + assert len(lattice_d2q4_256_8.registers) == 256 * 8 + assert all(reg.size == 4 for reg in lattice_d2q4_256_8.registers) + assert lattice_d2q4_256_8.circuit.num_qubits == 256 * 8 * 4 + + +def test_lqlga_grid_index_mapping_edge(): + lattice = LQLGALattice( + { + "lattice": { + "dim": {"x": 64, "y": 8}, + "velocities": "D2Q4", + }, + }, + ) + + assert lattice.gridpoint_index_flat(0) == (0, 0) + assert lattice.gridpoint_index_flat(1) == (0, 1) + assert lattice.gridpoint_index_flat(7) == (0, 7) + assert lattice.gridpoint_index_flat(8) == (1, 0) + assert lattice.gridpoint_index_flat(9) == (1, 1) + assert lattice.gridpoint_index_flat(511) == (63, 7) + + +@pytest.mark.parametrize( + "lattice_name", + ["lattice_d2q4_256_8", "lattice_d1q2_256"], +) +def test_lqlga_grid_index_mapping_general(lattice_name, request): + lattice = request.getfixturevalue(lattice_name) + assert all( + lattice.gridpoint_index_tuple(lattice.gridpoint_index_flat(gp)) == gp + for gp in range( + lattice.num_total_qubits + // LatticeDiscretizationProperties.get_num_velocities( + lattice.discretization + ) + ) + ) diff --git a/test/unit/spacetime/circuits/mcswap_test.py b/test/unit/spacetime/circuits/mcswap_test.py index bed58b7..27fdd2e 100644 --- a/test/unit/spacetime/circuits/mcswap_test.py +++ b/test/unit/spacetime/circuits/mcswap_test.py @@ -1,4 +1,4 @@ -from qiskit.circuit.library import MCMT, XGate +from qiskit.circuit.library import MCMTGate, XGate from qlbm.components.common.primitives import MCSwap from qlbm.lattice.lattices.spacetime_lattice import SpaceTimeLattice @@ -10,7 +10,7 @@ def test_mcswap_13ctrl(): { "lattice": { "dim": {"x": 16, "y": 16}, - "velocities": {"x": 2, "y": 2}, + "velocities": "D2Q4", }, "geometry": [], }, @@ -21,13 +21,13 @@ def test_mcswap_13ctrl(): expected_circuit.cx(6, 5) expected_circuit.compose( - MCMT(XGate(), 4, 1), + MCMTGate(XGate(), 4, 1), qubits=[0, 2, 3, 5, 6], inplace=True, ) expected_circuit.cx(6, 5) - assert mcswap.circuit == expected_circuit + assert mcswap.circuit.decompose() == expected_circuit.decompose() def test_mcswap_grid_ctrl(): @@ -36,7 +36,7 @@ def test_mcswap_grid_ctrl(): { "lattice": { "dim": {"x": 16, "y": 16}, - "velocities": {"x": 2, "y": 2}, + "velocities": "D2Q4", }, "geometry": [], }, @@ -51,10 +51,10 @@ def test_mcswap_grid_ctrl(): expected_circuit.cx(11, 9) expected_circuit.compose( - MCMT(XGate(), 9, 1), + MCMTGate(XGate(), 9, 1), qubits=list(range(8)) + [9, 11], inplace=True, ) expected_circuit.cx(11, 9) - assert mcswap.circuit == expected_circuit + assert mcswap.circuit.decompose() == expected_circuit.decompose() diff --git a/test/unit/spacetime/conftest.py b/test/unit/spacetime/conftest.py index c098e93..52a7c34 100644 --- a/test/unit/spacetime/conftest.py +++ b/test/unit/spacetime/conftest.py @@ -13,7 +13,7 @@ def dummy_1d_lattice() -> SpaceTimeLattice: { "lattice": { "dim": {"x": 256}, - "velocities": {"x": 2}, + "velocities": "D1Q2", }, }, ) @@ -26,7 +26,7 @@ def lattice_1d_16_1_obstacle_1_timestep() -> SpaceTimeLattice: { "lattice": { "dim": {"x": 16}, - "velocities": {"x": 2}, + "velocities": "D1Q2" }, "geometry": [ {"shape": "cuboid", "x": [4, 6], "boundary": "bounceback"}, @@ -42,7 +42,7 @@ def lattice_1d_16_1_obstacle_2_timesteps() -> SpaceTimeLattice: { "lattice": { "dim": {"x": 16}, - "velocities": {"x": 2}, + "velocities": "D1Q2" }, "geometry": [ {"shape": "cuboid", "x": [4, 6], "boundary": "bounceback"}, @@ -58,7 +58,7 @@ def lattice_1d_16_1_obstacle_5_timesteps() -> SpaceTimeLattice: { "lattice": { "dim": {"x": 16}, - "velocities": {"x": 2}, + "velocities": "D1Q2" }, }, ) @@ -71,7 +71,7 @@ def volumetric_lattice_1d_16_1_obstacle_1_timestep() -> SpaceTimeLattice: { "lattice": { "dim": {"x": 16}, - "velocities": {"x": 2}, + "velocities": "D1Q2" }, }, use_volumetric_ops=True, @@ -86,7 +86,7 @@ def dummy_lattice() -> SpaceTimeLattice: { "lattice": { "dim": {"x": 32, "y": 32}, - "velocities": {"x": 2, "y": 2}, + "velocities": "D2Q4" }, }, ) @@ -99,7 +99,7 @@ def lattice_2d_16x16_1_obstacle_1_timestep() -> SpaceTimeLattice: { "lattice": { "dim": {"x": 16, "y": 16}, - "velocities": {"x": 2, "y": 2}, + "velocities": "D2Q4" }, "geometry": [ { @@ -120,7 +120,7 @@ def lattice_2d_16x16_1_obstacle_2_timesteps() -> SpaceTimeLattice: { "lattice": { "dim": {"x": 16, "y": 16}, - "velocities": {"x": 2, "y": 2}, + "velocities": "D2Q4" }, "geometry": [ { @@ -141,7 +141,7 @@ def lattice_2d_16x16_1_obstacle_5_timesteps() -> SpaceTimeLattice: { "lattice": { "dim": {"x": 16, "y": 16}, - "velocities": {"x": 2, "y": 2}, + "velocities": "D2Q4" }, "geometry": [ { @@ -162,7 +162,7 @@ def lattice_1d_16_1_obstacle_5_timesteps_1_obstacle() -> SpaceTimeLattice: { "lattice": { "dim": {"x": 16}, - "velocities": {"x": 2}, + "velocities": "D1Q2" }, "geometry": [ {"shape": "cuboid", "x": [5, 11], "boundary": "bounceback"}, @@ -178,7 +178,7 @@ def lattice_1d_16_1_obstacle_5_timesteps_large_obstacle() -> SpaceTimeLattice: { "lattice": { "dim": {"x": 16}, - "velocities": {"x": 2}, + "velocities": "D1Q2" }, "geometry": [ {"shape": "cuboid", "x": [2, 14], "boundary": "bounceback"}, diff --git a/test/unit/spacetime/d1q2_lattice_test.py b/test/unit/spacetime/d1q2_lattice_test.py index 99bd5e3..3a64214 100644 --- a/test/unit/spacetime/d1q2_lattice_test.py +++ b/test/unit/spacetime/d1q2_lattice_test.py @@ -239,7 +239,7 @@ def test_bad_lattice_specification_velocities(): ) assert ( - "Unsupported number of velocities for 1D: 4. Only D1Q2 is supported at the moment." + "LatticeDiscretization.CFLDISCRETIZATION not currently implemented for the Space-Time method. Only D1Q2 and D2Q4 are fully at the moment. D3Q6 only supports collision." == str(excinfo_measure.value) ) diff --git a/test/unit/spacetime/d2q4_lattice_test.py b/test/unit/spacetime/d2q4_lattice_test.py index 37448b5..c685bea 100644 --- a/test/unit/spacetime/d2q4_lattice_test.py +++ b/test/unit/spacetime/d2q4_lattice_test.py @@ -15,7 +15,7 @@ def lattice_2d_16x16_1_obstacle_1_timestep() -> SpaceTimeLattice: { "lattice": { "dim": {"x": 16, "y": 16}, - "velocities": {"x": 2, "y": 2}, + "velocities": "D2Q4", }, "geometry": [ {"shape": "cuboid", "x": [4, 6], "y": [3, 12], "boundary": "specular"}, @@ -31,7 +31,7 @@ def lattice_2d_16x16_1_obstacle_2_timesteps() -> SpaceTimeLattice: { "lattice": { "dim": {"x": 16, "y": 16}, - "velocities": {"x": 2, "y": 2}, + "velocities": "D2Q4", }, "geometry": [ {"shape": "cuboid", "x": [4, 6], "y": [3, 12], "boundary": "specular"}, @@ -47,7 +47,7 @@ def volumetric_lattice_2d_16x16_1_obstacle_1_timestep() -> SpaceTimeLattice: { "lattice": { "dim": {"x": 16, "y": 16}, - "velocities": {"x": 2, "y": 2}, + "velocities": "D2Q4", }, }, use_volumetric_ops=True, @@ -207,7 +207,7 @@ def test_adaptable_qubit_register_indexing_measure_off_volume_off(): { "lattice": { "dim": {"x": 16, "y": 16}, - "velocities": {"x": 2, "y": 2}, + "velocities": "D2Q4", }, }, ) @@ -237,7 +237,7 @@ def test_adaptable_qubit_register_indexing_measure_on_volume_off(): { "lattice": { "dim": {"x": 16, "y": 16}, - "velocities": {"x": 2, "y": 2}, + "velocities": "D2Q4", }, }, include_measurement_qubit=True, @@ -262,7 +262,7 @@ def test_adaptable_qubit_register_indexing_measure_on_volume_on(): { "lattice": { "dim": {"x": 16, "y": 16}, - "velocities": {"x": 2, "y": 2}, + "velocities": "D2Q4", }, }, include_measurement_qubit=True, @@ -1123,25 +1123,26 @@ def test_bad_lattice_specification_velocities(): ) assert ( - "Unsupported number of velocities for 2D: (4, 2). Only D2Q4 is supported at the moment." + "LatticeDiscretization.CFLDISCRETIZATION not currently implemented for the Space-Time method. Only D1Q2 and D2Q4 are fully at the moment. D3Q6 only supports collision." == str(excinfo_measure.value) ) -def test_bad_lattice_specification_velocities_3d(): +def test_bad_lattice_specification_velocities_4d(): with pytest.raises(LatticeException) as excinfo_measure: SpaceTimeLattice( 2, { "lattice": { - "dim": {"x": 16, "y": 16, "z": 16}, - "velocities": {"x": 2, "y": 2, "z": 2}, + "dim": {"x": 16, "y": 16, "z": 16, "a": 32}, + "velocities": {"x": 2, "y": 2, "z": 2, "a": 2}, }, }, ) - assert "Only 1D and 2D discretizations are currently available." == str( - excinfo_measure.value + assert ( + "Only 1, 2, and 3-dimensional lattices are supported. Provided lattice has 4 dimensions." + == str(excinfo_measure.value) ) @@ -1151,7 +1152,7 @@ def test_2d_lattice_grid_register(): { "lattice": { "dim": {"x": 256, "y": 256}, - "velocities": {"x": 2, "y": 2}, + "velocities": "D2Q4", }, }, )