We are currently solving the following issue within our repository. Here is the issue text:
--- BEGIN ISSUE ---
Bug: incompatible matrix is accepted as input when creating local operator
I caught the following bug:

```python
>>> import netket as nk
/Users/vladimirvargas/.local/share/virtualenvs/physics-jfSzH1Jb/lib/python3.9/site-packages/jax/_src/lib/__init__.py:33: UserWarning: JAX on Mac ARM machines is experimental and minimally tested. Please see https://github.com/google/jax/issues/5501 in the event of problems.
  warnings.warn("JAX on Mac ARM machines is experimental and minimally tested. "
>>> hi = nk.hilbert.Fock(5, N=3)
>>> ham = nk.operator.boson.number(hi, 0) * nk.operator.boson.number(hi, 1)
>>> sparse_ham = ham.to_sparse()
>>> ham2 = nk.operator.LocalOperator(hi, sparse_ham, [0, 1])
>>> sparse_ham.shape
(216, 216)
>>> 6**3
216
```

You see, the Hilbert space has three parts. First I create `ham` acting on parts 0 and 1 only. However, when I invoke the sparse form of this Hamiltonian `ham`, I get a matrix for the whole Hilbert space, i.e. this is an operator that acts on the three parts. When I create another operator whose matrix is `sparse_ham` acting on the parts 0 and 1 of the Hilbert space, an error should be raised.

I guess this can be solved by checking that the product of the elements of the input matrix shape must equal the product of the local sizes onto which the matrix acts on in the Hilbert space.
Bug: incompatible matrix is accepted as input when creating local operator
I caught the following bug:

```python
>>> import netket as nk
/Users/vladimirvargas/.local/share/virtualenvs/physics-jfSzH1Jb/lib/python3.9/site-packages/jax/_src/lib/__init__.py:33: UserWarning: JAX on Mac ARM machines is experimental and minimally tested. Please see https://github.com/google/jax/issues/5501 in the event of problems.
  warnings.warn("JAX on Mac ARM machines is experimental and minimally tested. "
>>> hi = nk.hilbert.Fock(5, N=3)
>>> ham = nk.operator.boson.number(hi, 0) * nk.operator.boson.number(hi, 1)
>>> sparse_ham = ham.to_sparse()
>>> ham2 = nk.operator.LocalOperator(hi, sparse_ham, [0, 1])
>>> sparse_ham.shape
(216, 216)
>>> 6**3
216
```

You see, the Hilbert space has three parts. First I create `ham` acting on parts 0 and 1 only. However, when I invoke the sparse form of this Hamiltonian `ham`, I get a matrix for the whole Hilbert space, i.e. this is an operator that acts on the three parts. When I create another operator whose matrix is `sparse_ham` acting on the parts 0 and 1 of the Hilbert space, an error should be raised.

I guess this can be solved by checking that the product of the elements of the input matrix shape must equal the product of the local sizes onto which the matrix acts on in the Hilbert space.

--- END ISSUE ---

Below are some code segments, each from a relevant file. One or more of these files may contain bugs.
--- BEGIN FILES ---
Path: `netket/operator/_local_operator_helpers.py`
Content:
```
1 # Copyright 2022 The NetKet Authors - All rights reserved.
2 #
3 # Licensed under the Apache License, Version 2.0 (the "License");
4 # you may not use this file except in compliance with the License.
5 # You may obtain a copy of the License at
6 #
7 #    http://www.apache.org/licenses/LICENSE-2.0
8 #
9 # Unless required by applicable law or agreed to in writing, software
10 # distributed under the License is distributed on an "AS IS" BASIS,
11 # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12 # See the License for the specific language governing permissions and
13 # limitations under the License.
14 
15 from typing import Tuple, Union
16 
17 import functools
18 import numbers
19 
20 import numpy as np
21 
22 from scipy.sparse import spmatrix
23 
24 from netket.hilbert import AbstractHilbert, Fock
25 from netket.utils.types import DType, Array
26 
27 from ._abstract_operator import AbstractOperator
28 from ._discrete_operator import DiscreteOperator
29 
30 
31 def _dtype(obj: Union[numbers.Number, Array, AbstractOperator]) -> DType:
32     """
33     Returns the dtype of the input object
34     """
35     if isinstance(obj, numbers.Number):
36         return type(obj)
37     elif isinstance(obj, DiscreteOperator):
38         return obj.dtype
39     elif isinstance(obj, np.ndarray):
40         return obj.dtype
41     else:
42         raise TypeError(f"cannot deduce dtype of object type {type(obj)}: {obj}")
43 
44 
45 def cast_operator_matrix_dtype(matrix, dtype):
46     """
47     Changes the dtype of a matrix, without changing the structural type of the object.
48 
49     This makes sure that if you pass sparse arrays to a LocalOperator, they remain
50     sparse even if you change the dtype
51     """
52     # must copy
53     # return np.asarray(matrix, dtype=dtype)
54     return matrix.astype(dtype)
55 
56 
57 def _standardize_matrix_input_type(op):
58     """
59     Standardize the structural type of operators stored in LocalOperator.
60 
61     Eventually, we could also support spmatrices (but some work will be needed.)
62     """
63     if isinstance(op, list):
64         return np.asarray(op)
65     elif isinstance(op, spmatrix):
66         return op.todense()
67     else:
68         return op
69 
70 
71 def canonicalize_input(
72     hilbert: AbstractHilbert, operators, acting_on, constant, *, dtype=None
73 ):
74     """
75     Takes as inputs the inputs to the constructor of LocalOperator and canonicalizes
76     them by ensuring the following holds:
77      - acting_on is a list of list
78      - acting_on[i] are sorted
79      - operators is list of matrices
80      - all dtypes match
81 
82     Args:
83         hilbert: The hilbert space
84 
85     Returns:
86         List of operators, acting ons and dtypes.
87     """
88     # check if passing a single operator or a list of operators
89     if isinstance(acting_on, numbers.Number):
90         acting_on = [acting_on]
91 
92     is_nested = any(hasattr(i, "__len__") for i in acting_on)
93     if not is_nested:
94         operators = [operators]
95         acting_on = [acting_on]
96 
97     if all(len(aon) == 0 for aon in acting_on):
98         operators = []
99         acting_on = []
100     else:
101         if max(map(max, acting_on)) >= hilbert.size or min(map(min, acting_on)) < 0:
102             raise ValueError("An operator acts on an invalid set of sites.")
103 
104     acting_on = [tuple(aon) for aon in acting_on]
105     # operators = [np.asarray(operator) for operator in operators]
106     operators = [_standardize_matrix_input_type(op) for op in operators]
107 
108     # If we asked for a specific dtype, enforce it.
109     if dtype is None:
110         dtype = np.promote_types(np.float32, _dtype(constant))
111         dtype = functools.reduce(
112             lambda dt, op: np.promote_types(dt, op.dtype), operators, dtype
113         )
114     dtype = np.empty((), dtype=dtype).dtype
115 
116     canonicalized_operators = []
117     canonicalized_acting_on = []
118     for (operator, acting_on) in zip(operators, acting_on):
119         if operator.dtype is not dtype:
120             operator = cast_operator_matrix_dtype(operator, dtype=dtype)
121 
122         # re-sort the operator
123         operator, acting_on = _reorder_kronecker_product(hilbert, operator, acting_on)
124         canonicalized_operators.append(operator)
125         canonicalized_acting_on.append(acting_on)
126 
127     return canonicalized_operators, canonicalized_acting_on, dtype
128 
129 
130 # TODO: support sparse arrays without returning dense arrays
131 def _reorder_kronecker_product(hi, mat, acting_on) -> Tuple[Array, Tuple]:
132     """
133     Reorders the matrix resulting from a kronecker product of several
134     operators in such a way to sort acting_on.
135 
136     A conceptual example is the following:
137     if `mat = Â ⊗ B̂ ⊗ Ĉ` and `acting_on = [[2],[1],[3]`
138     you will get `result = B̂ ⊗ Â ⊗ Ĉ, [[1], [2], [3]].
139 
140     However, essentially, A,B,C represent some operators acting on
141     thei sub-space acting_on[1], [2] and [3] of the hilbert space.
142 
143     This function also handles any possible set of values in acting_on.
144 
145     The inner logic uses the Fock.all_states(), number_to_state and
146     state_to_number to perform the re-ordering.
147     """
148     acting_on_sorted = np.sort(acting_on)
149     if np.array_equal(acting_on_sorted, acting_on):
150         return mat, acting_on
151 
152     # could write custom binary <-> int logic instead of using Fock...
153     # Since i need to work with bit-strings (where instead of bits i
154     # have integers, in order to support arbitrary size spaces) this
155     # is exactly what hilbert.to_number() and viceversa do.
156 
157     # target ordering binary representation
158     hi_subspace = Fock(hi.shape[acting_on_sorted[0]] - 1)
159     for site in acting_on_sorted[1:]:
160         hi_subspace = hi_subspace * Fock(hi.shape[site] - 1)
161 
162     hi_unsorted_subspace = Fock(hi.shape[acting_on[0]] - 1)
163     for site in acting_on[1:]:
164         hi_unsorted_subspace = hi_unsorted_subspace * Fock(hi.shape[site] - 1)
165 
166     # find how to map target ordering back to unordered
167     acting_on_unsorted_ids = np.zeros(len(acting_on), dtype=np.intp)
168     for (i, site) in enumerate(acting_on):
169         acting_on_unsorted_ids[i] = np.argmax(site == acting_on_sorted)
170 
171     # now it is valid that
172     # acting_on_sorted == acting_on[acting_on_unsorted_ids]
173 
174     # generate n-bit strings in the target ordering
175     v = hi_subspace.all_states()
176 
177     # convert them to origin (unordered) ordering
178     v_unsorted = v[:, acting_on_unsorted_ids]
179     # convert the unordered bit-strings to numbers in the target space.
180     n_unsorted = hi_unsorted_subspace.states_to_numbers(v_unsorted)
181 
182     # reorder the matrix
183     mat_sorted = mat[n_unsorted, :][:, n_unsorted]
184 
185     return mat_sorted, tuple(acting_on_sorted)
186 
187 
188 # TODO: support sparse arrays without returning dense arrays
189 def _multiply_operators(
190     hilbert, support_A: Tuple, A: Array, support_B: Tuple, B: Array, *, dtype
191 ) -> Tuple[Tuple, Array]:
192     """
193     Returns the `Tuple[acting_on, Matrix]` representing the operator obtained by
194     multiplying the two input operators A and B.
195     """
196     support_A = np.asarray(support_A)
197     support_B = np.asarray(support_B)
198 
199     inters = np.intersect1d(support_A, support_B, return_indices=False)
200 
201     if support_A.size == support_B.size and np.array_equal(support_A, support_B):
202         return tuple(support_A), A @ B
203     elif inters.size == 0:
204         # disjoint supports
205         support = tuple(np.concatenate([support_A, support_B]))
206         operator = np.kron(A, B)
207         operator, support = _reorder_kronecker_product(hilbert, operator, support)
208         return tuple(support), operator
209     else:
210         _support_A = list(support_A)
211         _support_B = list(support_B)
212         _A = A.copy()
213         _B = B.copy()
214 
215         # expand _act to match _act_i
216         supp_B_min = min(support_B)
217         for site in support_A:
218             if site not in support_B:
219                 I = np.eye(hilbert.shape[site], dtype=dtype)
220                 if site < supp_B_min:
221                     _support_B = [site] + _support_B
222                     _B = np.kron(I, _B)
223                 else:  # site > actmax
224                     _support_B = _support_B + [site]
225                     _B = np.kron(_B, I)
226 
227         supp_A_min = min(support_A)
228         for site in support_B:
229             if site not in support_A:
230                 I = np.eye(hilbert.shape[site], dtype=dtype)
231                 if site < supp_A_min:
232                     _support_A = [site] + _support_A
233                     _A = np.kron(I, _A)
234                 else:  #  site > actmax
235                     _support_A = _support_A + [site]
236                     _A = np.kron(_A, I)
237 
238         # reorder
239         _A, _support_A = _reorder_kronecker_product(hilbert, _A, _support_A)
240         _B, _support_B = _reorder_kronecker_product(hilbert, _B, _support_B)
241 
242         if len(_support_A) == len(_support_B) and np.array_equal(
243             _support_A, _support_B
244         ):
245             # back to the case of non-interesecting with same support
246             return tuple(_support_A), _A @ _B
247         else:
248             raise ValueError("Something failed")
249
```
--- END FILES ---

Please first localize the bug based on the issue statement, and then generate a patch according to the `git diff` format fenced by three backticks.

Here is an example:

```diff
diff --git a/examples/server_async.py b/examples/server_async.py
--- a/examples/server_async.py
+++ b/examples/server_async.py
@@ -313,4 +313,4 @@
 
 
 if __name__ == "__main__":
-    asyncio.run(run_async_server("."), debug=True)
+    asyncio.run(run_async_server(), debug=True)
diff --git a/examples/server_sync.py b/examples/server_sync.py
--- a/examples/server_sync.py
+++ b/examples/server_sync.py
@@ -313,5 +313,5 @@
 
 
 if __name__ == "__main__":
-    server = run_sync_server(".")
+    server = run_sync_server()
     server.shutdown()

```

diff --git a/netket/operator/_local_operator_helpers.py b/netket/operator/_local_operator_helpers.py
--- a/netket/operator/_local_operator_helpers.py
+++ b/netket/operator/_local_operator_helpers.py
@@ -116,6 +116,8 @@
     canonicalized_operators = []
     canonicalized_acting_on = []
     for (operator, acting_on) in zip(operators, acting_on):
+        check_valid_opmatrix(hilbert, operator, acting_on)
+
         if operator.dtype is not dtype:
             operator = cast_operator_matrix_dtype(operator, dtype=dtype)
 
@@ -127,6 +129,18 @@
     return canonicalized_operators, canonicalized_acting_on, dtype
 
 
+def check_valid_opmatrix(hi, mat, acting_on):
+    """ """
+    expected_size = np.prod([hi.shape[aon] for aon in acting_on])
+
+    if mat.shape != (expected_size, expected_size):
+        raise ValueError(
+            f"The matrix of the sub-operator acting on sites {acting_on} "
+            f"must have shape {expected_size, expected_size}, "
+            f"but it has shape {mat.shape}."
+        )
+
+
 # TODO: support sparse arrays without returning dense arrays
 def _reorder_kronecker_product(hi, mat, acting_on) -> Tuple[Array, Tuple]:
     """

{"golden_diff": "diff --git a/netket/operator/_local_operator_helpers.py b/netket/operator/_local_operator_helpers.py\n--- a/netket/operator/_local_operator_helpers.py\n+++ b/netket/operator/_local_operator_helpers.py\n@@ -116,6 +116,8 @@\n     canonicalized_operators = []\n     canonicalized_acting_on = []\n     for (operator, acting_on) in zip(operators, acting_on):\n+        check_valid_opmatrix(hilbert, operator, acting_on)\n+\n         if operator.dtype is not dtype:\n             operator = cast_operator_matrix_dtype(operator, dtype=dtype)\n \n@@ -127,6 +129,18 @@\n     return canonicalized_operators, canonicalized_acting_on, dtype\n \n \n+def check_valid_opmatrix(hi, mat, acting_on):\n+    \"\"\" \"\"\"\n+    expected_size = np.prod([hi.shape[aon] for aon in acting_on])\n+\n+    if mat.shape != (expected_size, expected_size):\n+        raise ValueError(\n+            f\"The matrix of the sub-operator acting on sites {acting_on} \"\n+            f\"must have shape {expected_size, expected_size}, \"\n+            f\"but it has shape {mat.shape}.\"\n+        )\n+\n+\n # TODO: support sparse arrays without returning dense arrays\n def _reorder_kronecker_product(hi, mat, acting_on) -> Tuple[Array, Tuple]:\n     \"\"\"\n", "issue": "Bug: incompatible matrix is accepted as input when creating local operator\nI caught the following bug:\r\n\r\n```python\r\n>>> import netket as nk\r\n/Users/vladimirvargas/.local/share/virtualenvs/physics-jfSzH1Jb/lib/python3.9/site-packages/jax/_src/lib/__init__.py:33: UserWarning: JAX on Mac ARM machines is experimental and minimally tested. Please see https://github.com/google/jax/issues/5501 in the event of problems.\r\n  warnings.warn(\"JAX on Mac ARM machines is experimental and minimally tested. \"\r\n>>> hi = nk.hilbert.Fock(5, N=3)\r\n>>> ham = nk.operator.boson.number(hi, 0) * nk.operator.boson.number(hi, 1)\r\n>>> sparse_ham = ham.to_sparse()\r\n>>> ham2 = nk.operator.LocalOperator(hi, sparse_ham, [0, 1])\r\n>>> sparse_ham.shape\r\n(216, 216)\r\n>>> 6**3\r\n216\r\n```\r\n\r\nYou see, the Hilbert space has three parts. First I create `ham` acting on parts 0 and 1 only. However, when I invoke the sparse form of this Hamiltonian `ham`, I get a matrix for the whole Hilbert space, i.e. this is an operator that acts on the three parts. When I create another operator whose matrix is `sparse_ham` acting on the parts 0 and 1 of the Hilbert space, an error should be raised.\r\n\r\nI guess this can be solved by checking that the product of the elements of the input matrix shape must equal the product of the local sizes onto which the matrix acts on in the Hilbert space.\nBug: incompatible matrix is accepted as input when creating local operator\nI caught the following bug:\r\n\r\n```python\r\n>>> import netket as nk\r\n/Users/vladimirvargas/.local/share/virtualenvs/physics-jfSzH1Jb/lib/python3.9/site-packages/jax/_src/lib/__init__.py:33: UserWarning: JAX on Mac ARM machines is experimental and minimally tested. Please see https://github.com/google/jax/issues/5501 in the event of problems.\r\n  warnings.warn(\"JAX on Mac ARM machines is experimental and minimally tested. \"\r\n>>> hi = nk.hilbert.Fock(5, N=3)\r\n>>> ham = nk.operator.boson.number(hi, 0) * nk.operator.boson.number(hi, 1)\r\n>>> sparse_ham = ham.to_sparse()\r\n>>> ham2 = nk.operator.LocalOperator(hi, sparse_ham, [0, 1])\r\n>>> sparse_ham.shape\r\n(216, 216)\r\n>>> 6**3\r\n216\r\n```\r\n\r\nYou see, the Hilbert space has three parts. First I create `ham` acting on parts 0 and 1 only. However, when I invoke the sparse form of this Hamiltonian `ham`, I get a matrix for the whole Hilbert space, i.e. this is an operator that acts on the three parts. When I create another operator whose matrix is `sparse_ham` acting on the parts 0 and 1 of the Hilbert space, an error should be raised.\r\n\r\nI guess this can be solved by checking that the product of the elements of the input matrix shape must equal the product of the local sizes onto which the matrix acts on in the Hilbert space.\n", "before_files": [{"content": "# Copyright 2022 The NetKet Authors - All rights reserved.\n#\n# Licensed under the Apache License, Version 2.0 (the \"License\");\n# you may not use this file except in compliance with the License.\n# You may obtain a copy of the License at\n#\n#    http://www.apache.org/licenses/LICENSE-2.0\n#\n# Unless required by applicable law or agreed to in writing, software\n# distributed under the License is distributed on an \"AS IS\" BASIS,\n# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.\n# See the License for the specific language governing permissions and\n# limitations under the License.\n\nfrom typing import Tuple, Union\n\nimport functools\nimport numbers\n\nimport numpy as np\n\nfrom scipy.sparse import spmatrix\n\nfrom netket.hilbert import AbstractHilbert, Fock\nfrom netket.utils.types import DType, Array\n\nfrom ._abstract_operator import AbstractOperator\nfrom ._discrete_operator import DiscreteOperator\n\n\ndef _dtype(obj: Union[numbers.Number, Array, AbstractOperator]) -> DType:\n    \"\"\"\n    Returns the dtype of the input object\n    \"\"\"\n    if isinstance(obj, numbers.Number):\n        return type(obj)\n    elif isinstance(obj, DiscreteOperator):\n        return obj.dtype\n    elif isinstance(obj, np.ndarray):\n        return obj.dtype\n    else:\n        raise TypeError(f\"cannot deduce dtype of object type {type(obj)}: {obj}\")\n\n\ndef cast_operator_matrix_dtype(matrix, dtype):\n    \"\"\"\n    Changes the dtype of a matrix, without changing the structural type of the object.\n\n    This makes sure that if you pass sparse arrays to a LocalOperator, they remain\n    sparse even if you change the dtype\n    \"\"\"\n    # must copy\n    # return np.asarray(matrix, dtype=dtype)\n    return matrix.astype(dtype)\n\n\ndef _standardize_matrix_input_type(op):\n    \"\"\"\n    Standardize the structural type of operators stored in LocalOperator.\n\n    Eventually, we could also support spmatrices (but some work will be needed.)\n    \"\"\"\n    if isinstance(op, list):\n        return np.asarray(op)\n    elif isinstance(op, spmatrix):\n        return op.todense()\n    else:\n        return op\n\n\ndef canonicalize_input(\n    hilbert: AbstractHilbert, operators, acting_on, constant, *, dtype=None\n):\n    \"\"\"\n    Takes as inputs the inputs to the constructor of LocalOperator and canonicalizes\n    them by ensuring the following holds:\n     - acting_on is a list of list\n     - acting_on[i] are sorted\n     - operators is list of matrices\n     - all dtypes match\n\n    Args:\n        hilbert: The hilbert space\n\n    Returns:\n        List of operators, acting ons and dtypes.\n    \"\"\"\n    # check if passing a single operator or a list of operators\n    if isinstance(acting_on, numbers.Number):\n        acting_on = [acting_on]\n\n    is_nested = any(hasattr(i, \"__len__\") for i in acting_on)\n    if not is_nested:\n        operators = [operators]\n        acting_on = [acting_on]\n\n    if all(len(aon) == 0 for aon in acting_on):\n        operators = []\n        acting_on = []\n    else:\n        if max(map(max, acting_on)) >= hilbert.size or min(map(min, acting_on)) < 0:\n            raise ValueError(\"An operator acts on an invalid set of sites.\")\n\n    acting_on = [tuple(aon) for aon in acting_on]\n    # operators = [np.asarray(operator) for operator in operators]\n    operators = [_standardize_matrix_input_type(op) for op in operators]\n\n    # If we asked for a specific dtype, enforce it.\n    if dtype is None:\n        dtype = np.promote_types(np.float32, _dtype(constant))\n        dtype = functools.reduce(\n            lambda dt, op: np.promote_types(dt, op.dtype), operators, dtype\n        )\n    dtype = np.empty((), dtype=dtype).dtype\n\n    canonicalized_operators = []\n    canonicalized_acting_on = []\n    for (operator, acting_on) in zip(operators, acting_on):\n        if operator.dtype is not dtype:\n            operator = cast_operator_matrix_dtype(operator, dtype=dtype)\n\n        # re-sort the operator\n        operator, acting_on = _reorder_kronecker_product(hilbert, operator, acting_on)\n        canonicalized_operators.append(operator)\n        canonicalized_acting_on.append(acting_on)\n\n    return canonicalized_operators, canonicalized_acting_on, dtype\n\n\n# TODO: support sparse arrays without returning dense arrays\ndef _reorder_kronecker_product(hi, mat, acting_on) -> Tuple[Array, Tuple]:\n    \"\"\"\n    Reorders the matrix resulting from a kronecker product of several\n    operators in such a way to sort acting_on.\n\n    A conceptual example is the following:\n    if `mat = \u00c2 \u2297 B\u0302 \u2297 \u0108` and `acting_on = [[2],[1],[3]`\n    you will get `result = B\u0302 \u2297 \u00c2 \u2297 \u0108, [[1], [2], [3]].\n\n    However, essentially, A,B,C represent some operators acting on\n    thei sub-space acting_on[1], [2] and [3] of the hilbert space.\n\n    This function also handles any possible set of values in acting_on.\n\n    The inner logic uses the Fock.all_states(), number_to_state and\n    state_to_number to perform the re-ordering.\n    \"\"\"\n    acting_on_sorted = np.sort(acting_on)\n    if np.array_equal(acting_on_sorted, acting_on):\n        return mat, acting_on\n\n    # could write custom binary <-> int logic instead of using Fock...\n    # Since i need to work with bit-strings (where instead of bits i\n    # have integers, in order to support arbitrary size spaces) this\n    # is exactly what hilbert.to_number() and viceversa do.\n\n    # target ordering binary representation\n    hi_subspace = Fock(hi.shape[acting_on_sorted[0]] - 1)\n    for site in acting_on_sorted[1:]:\n        hi_subspace = hi_subspace * Fock(hi.shape[site] - 1)\n\n    hi_unsorted_subspace = Fock(hi.shape[acting_on[0]] - 1)\n    for site in acting_on[1:]:\n        hi_unsorted_subspace = hi_unsorted_subspace * Fock(hi.shape[site] - 1)\n\n    # find how to map target ordering back to unordered\n    acting_on_unsorted_ids = np.zeros(len(acting_on), dtype=np.intp)\n    for (i, site) in enumerate(acting_on):\n        acting_on_unsorted_ids[i] = np.argmax(site == acting_on_sorted)\n\n    # now it is valid that\n    # acting_on_sorted == acting_on[acting_on_unsorted_ids]\n\n    # generate n-bit strings in the target ordering\n    v = hi_subspace.all_states()\n\n    # convert them to origin (unordered) ordering\n    v_unsorted = v[:, acting_on_unsorted_ids]\n    # convert the unordered bit-strings to numbers in the target space.\n    n_unsorted = hi_unsorted_subspace.states_to_numbers(v_unsorted)\n\n    # reorder the matrix\n    mat_sorted = mat[n_unsorted, :][:, n_unsorted]\n\n    return mat_sorted, tuple(acting_on_sorted)\n\n\n# TODO: support sparse arrays without returning dense arrays\ndef _multiply_operators(\n    hilbert, support_A: Tuple, A: Array, support_B: Tuple, B: Array, *, dtype\n) -> Tuple[Tuple, Array]:\n    \"\"\"\n    Returns the `Tuple[acting_on, Matrix]` representing the operator obtained by\n    multiplying the two input operators A and B.\n    \"\"\"\n    support_A = np.asarray(support_A)\n    support_B = np.asarray(support_B)\n\n    inters = np.intersect1d(support_A, support_B, return_indices=False)\n\n    if support_A.size == support_B.size and np.array_equal(support_A, support_B):\n        return tuple(support_A), A @ B\n    elif inters.size == 0:\n        # disjoint supports\n        support = tuple(np.concatenate([support_A, support_B]))\n        operator = np.kron(A, B)\n        operator, support = _reorder_kronecker_product(hilbert, operator, support)\n        return tuple(support), operator\n    else:\n        _support_A = list(support_A)\n        _support_B = list(support_B)\n        _A = A.copy()\n        _B = B.copy()\n\n        # expand _act to match _act_i\n        supp_B_min = min(support_B)\n        for site in support_A:\n            if site not in support_B:\n                I = np.eye(hilbert.shape[site], dtype=dtype)\n                if site < supp_B_min:\n                    _support_B = [site] + _support_B\n                    _B = np.kron(I, _B)\n                else:  # site > actmax\n                    _support_B = _support_B + [site]\n                    _B = np.kron(_B, I)\n\n        supp_A_min = min(support_A)\n        for site in support_B:\n            if site not in support_A:\n                I = np.eye(hilbert.shape[site], dtype=dtype)\n                if site < supp_A_min:\n                    _support_A = [site] + _support_A\n                    _A = np.kron(I, _A)\n                else:  # \u00a0site > actmax\n                    _support_A = _support_A + [site]\n                    _A = np.kron(_A, I)\n\n        # reorder\n        _A, _support_A = _reorder_kronecker_product(hilbert, _A, _support_A)\n        _B, _support_B = _reorder_kronecker_product(hilbert, _B, _support_B)\n\n        if len(_support_A) == len(_support_B) and np.array_equal(\n            _support_A, _support_B\n        ):\n            # back to the case of non-interesecting with same support\n            return tuple(_support_A), _A @ _B\n        else:\n            raise ValueError(\"Something failed\")\n", "path": "netket/operator/_local_operator_helpers.py"}], "after_files": [{"content": "# Copyright 2022 The NetKet Authors - All rights reserved.\n#\n# Licensed under the Apache License, Version 2.0 (the \"License\");\n# you may not use this file except in compliance with the License.\n# You may obtain a copy of the License at\n#\n#    http://www.apache.org/licenses/LICENSE-2.0\n#\n# Unless required by applicable law or agreed to in writing, software\n# distributed under the License is distributed on an \"AS IS\" BASIS,\n# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.\n# See the License for the specific language governing permissions and\n# limitations under the License.\n\nfrom typing import Tuple, Union\n\nimport functools\nimport numbers\n\nimport numpy as np\n\nfrom scipy.sparse import spmatrix\n\nfrom netket.hilbert import AbstractHilbert, Fock\nfrom netket.utils.types import DType, Array\n\nfrom ._abstract_operator import AbstractOperator\nfrom ._discrete_operator import DiscreteOperator\n\n\ndef _dtype(obj: Union[numbers.Number, Array, AbstractOperator]) -> DType:\n    \"\"\"\n    Returns the dtype of the input object\n    \"\"\"\n    if isinstance(obj, numbers.Number):\n        return type(obj)\n    elif isinstance(obj, DiscreteOperator):\n        return obj.dtype\n    elif isinstance(obj, np.ndarray):\n        return obj.dtype\n    else:\n        raise TypeError(f\"cannot deduce dtype of object type {type(obj)}: {obj}\")\n\n\ndef cast_operator_matrix_dtype(matrix, dtype):\n    \"\"\"\n    Changes the dtype of a matrix, without changing the structural type of the object.\n\n    This makes sure that if you pass sparse arrays to a LocalOperator, they remain\n    sparse even if you change the dtype\n    \"\"\"\n    # must copy\n    # return np.asarray(matrix, dtype=dtype)\n    return matrix.astype(dtype)\n\n\ndef _standardize_matrix_input_type(op):\n    \"\"\"\n    Standardize the structural type of operators stored in LocalOperator.\n\n    Eventually, we could also support spmatrices (but some work will be needed.)\n    \"\"\"\n    if isinstance(op, list):\n        return np.asarray(op)\n    elif isinstance(op, spmatrix):\n        return op.todense()\n    else:\n        return op\n\n\ndef canonicalize_input(\n    hilbert: AbstractHilbert, operators, acting_on, constant, *, dtype=None\n):\n    \"\"\"\n    Takes as inputs the inputs to the constructor of LocalOperator and canonicalizes\n    them by ensuring the following holds:\n     - acting_on is a list of list\n     - acting_on[i] are sorted\n     - operators is list of matrices\n     - all dtypes match\n\n    Args:\n        hilbert: The hilbert space\n\n    Returns:\n        List of operators, acting ons and dtypes.\n    \"\"\"\n    # check if passing a single operator or a list of operators\n    if isinstance(acting_on, numbers.Number):\n        acting_on = [acting_on]\n\n    is_nested = any(hasattr(i, \"__len__\") for i in acting_on)\n    if not is_nested:\n        operators = [operators]\n        acting_on = [acting_on]\n\n    if all(len(aon) == 0 for aon in acting_on):\n        operators = []\n        acting_on = []\n    else:\n        if max(map(max, acting_on)) >= hilbert.size or min(map(min, acting_on)) < 0:\n            raise ValueError(\"An operator acts on an invalid set of sites.\")\n\n    acting_on = [tuple(aon) for aon in acting_on]\n    # operators = [np.asarray(operator) for operator in operators]\n    operators = [_standardize_matrix_input_type(op) for op in operators]\n\n    # If we asked for a specific dtype, enforce it.\n    if dtype is None:\n        dtype = np.promote_types(np.float32, _dtype(constant))\n        dtype = functools.reduce(\n            lambda dt, op: np.promote_types(dt, op.dtype), operators, dtype\n        )\n    dtype = np.empty((), dtype=dtype).dtype\n\n    canonicalized_operators = []\n    canonicalized_acting_on = []\n    for (operator, acting_on) in zip(operators, acting_on):\n        check_valid_opmatrix(hilbert, operator, acting_on)\n\n        if operator.dtype is not dtype:\n            operator = cast_operator_matrix_dtype(operator, dtype=dtype)\n\n        # re-sort the operator\n        operator, acting_on = _reorder_kronecker_product(hilbert, operator, acting_on)\n        canonicalized_operators.append(operator)\n        canonicalized_acting_on.append(acting_on)\n\n    return canonicalized_operators, canonicalized_acting_on, dtype\n\n\ndef check_valid_opmatrix(hi, mat, acting_on):\n    \"\"\" \"\"\"\n    expected_size = np.prod([hi.shape[aon] for aon in acting_on])\n\n    if mat.shape != (expected_size, expected_size):\n        raise ValueError(\n            f\"The matrix of the sub-operator acting on sites {acting_on} \"\n            f\"must have shape {expected_size, expected_size}, \"\n            f\"but it has shape {mat.shape}.\"\n        )\n\n\n# TODO: support sparse arrays without returning dense arrays\ndef _reorder_kronecker_product(hi, mat, acting_on) -> Tuple[Array, Tuple]:\n    \"\"\"\n    Reorders the matrix resulting from a kronecker product of several\n    operators in such a way to sort acting_on.\n\n    A conceptual example is the following:\n    if `mat = \u00c2 \u2297 B\u0302 \u2297 \u0108` and `acting_on = [[2],[1],[3]`\n    you will get `result = B\u0302 \u2297 \u00c2 \u2297 \u0108, [[1], [2], [3]].\n\n    However, essentially, A,B,C represent some operators acting on\n    thei sub-space acting_on[1], [2] and [3] of the hilbert space.\n\n    This function also handles any possible set of values in acting_on.\n\n    The inner logic uses the Fock.all_states(), number_to_state and\n    state_to_number to perform the re-ordering.\n    \"\"\"\n    acting_on_sorted = np.sort(acting_on)\n    if np.array_equal(acting_on_sorted, acting_on):\n        return mat, acting_on\n\n    # could write custom binary <-> int logic instead of using Fock...\n    # Since i need to work with bit-strings (where instead of bits i\n    # have integers, in order to support arbitrary size spaces) this\n    # is exactly what hilbert.to_number() and viceversa do.\n\n    # target ordering binary representation\n    hi_subspace = Fock(hi.shape[acting_on_sorted[0]] - 1)\n    for site in acting_on_sorted[1:]:\n        hi_subspace = hi_subspace * Fock(hi.shape[site] - 1)\n\n    hi_unsorted_subspace = Fock(hi.shape[acting_on[0]] - 1)\n    for site in acting_on[1:]:\n        hi_unsorted_subspace = hi_unsorted_subspace * Fock(hi.shape[site] - 1)\n\n    # find how to map target ordering back to unordered\n    acting_on_unsorted_ids = np.zeros(len(acting_on), dtype=np.intp)\n    for (i, site) in enumerate(acting_on):\n        acting_on_unsorted_ids[i] = np.argmax(site == acting_on_sorted)\n\n    # now it is valid that\n    # acting_on_sorted == acting_on[acting_on_unsorted_ids]\n\n    # generate n-bit strings in the target ordering\n    v = hi_subspace.all_states()\n\n    # convert them to origin (unordered) ordering\n    v_unsorted = v[:, acting_on_unsorted_ids]\n    # convert the unordered bit-strings to numbers in the target space.\n    n_unsorted = hi_unsorted_subspace.states_to_numbers(v_unsorted)\n\n    # reorder the matrix\n    mat_sorted = mat[n_unsorted, :][:, n_unsorted]\n\n    return mat_sorted, tuple(acting_on_sorted)\n\n\n# TODO: support sparse arrays without returning dense arrays\ndef _multiply_operators(\n    hilbert, support_A: Tuple, A: Array, support_B: Tuple, B: Array, *, dtype\n) -> Tuple[Tuple, Array]:\n    \"\"\"\n    Returns the `Tuple[acting_on, Matrix]` representing the operator obtained by\n    multiplying the two input operators A and B.\n    \"\"\"\n    support_A = np.asarray(support_A)\n    support_B = np.asarray(support_B)\n\n    inters = np.intersect1d(support_A, support_B, return_indices=False)\n\n    if support_A.size == support_B.size and np.array_equal(support_A, support_B):\n        return tuple(support_A), A @ B\n    elif inters.size == 0:\n        # disjoint supports\n        support = tuple(np.concatenate([support_A, support_B]))\n        operator = np.kron(A, B)\n        operator, support = _reorder_kronecker_product(hilbert, operator, support)\n        return tuple(support), operator\n    else:\n        _support_A = list(support_A)\n        _support_B = list(support_B)\n        _A = A.copy()\n        _B = B.copy()\n\n        # expand _act to match _act_i\n        supp_B_min = min(support_B)\n        for site in support_A:\n            if site not in support_B:\n                I = np.eye(hilbert.shape[site], dtype=dtype)\n                if site < supp_B_min:\n                    _support_B = [site] + _support_B\n                    _B = np.kron(I, _B)\n                else:  # site > actmax\n                    _support_B = _support_B + [site]\n                    _B = np.kron(_B, I)\n\n        supp_A_min = min(support_A)\n        for site in support_B:\n            if site not in support_A:\n                I = np.eye(hilbert.shape[site], dtype=dtype)\n                if site < supp_A_min:\n                    _support_A = [site] + _support_A\n                    _A = np.kron(I, _A)\n                else:  # \u00a0site > actmax\n                    _support_A = _support_A + [site]\n                    _A = np.kron(_A, I)\n\n        # reorder\n        _A, _support_A = _reorder_kronecker_product(hilbert, _A, _support_A)\n        _B, _support_B = _reorder_kronecker_product(hilbert, _B, _support_B)\n\n        if len(_support_A) == len(_support_B) and np.array_equal(\n            _support_A, _support_B\n        ):\n            # back to the case of non-interesecting with same support\n            return tuple(_support_A), _A @ _B\n        else:\n            raise ValueError(\"Something failed\")\n", "path": "netket/operator/_local_operator_helpers.py"}]}