Pcvl 733 slos C++

Bench_backends.py

@@ -0,0 +1,148 @@
+import gc
+import time
+import perceval as pcvl
+import psutil
+
+from perceval.backends import BACKEND_LIST
+from perceval.backends._slos_v2 import SLOSV2Backend
+from perceval.backends._slos_v3 import SLOSV3Backend
+from perceval.utils.dist_metrics import tvd_dist
+from perceval.utils.postselect import PostSelect
+import exqalibur as xq
+
+from perceval.utils import get_logger
+get_logger().set_level(pcvl.logging.level.warn, pcvl.logging.channel.general)
+# get_logger().set_level(pcvl.logging.level.info, pcvl.logging.channel.general)
+
+# BACKEND_LIST[ "SLOS_V2_PS" ] = SLOSV2Backend
+# BACKEND_LIST[ "SLOS_V3_PS" ] = SLOSV3Backend
+
+def human_readable_size(size, decimal_places=2):
+    for unit in ['B', 'KiB', 'MiB', 'GiB', 'TiB', 'PiB']:
+        if size < 1024.0 or unit == 'PiB':
+            break
+        size /= 1024.0
+    return f"{size:.{decimal_places}f} {unit}"
+
+def do_compute(backend, method):
+    if method == 1:
+        return backend.all_prob_ampli()
+    if method == 2:
+        return backend.all_prob()
+    if method == 3:
+        return backend.prob_distribution()
+
+def set_mask(backend, mask_str, m, n):
+    if mask_str:
+        backend.set_mask((mask_str + "*" * m)[:m], n)
+    # if backend.name == "SLOS_V2_PS" or backend.name == "SLOS_V3_PS":
+    #     backend.set_post_select(ps)
+
+def bench(m, n, backends, mask_str, method):
+    u1 = pcvl.Matrix.random_unitary(m)
+    photons = [ 0 ] * m
+    photons[0] = n
+    bs = pcvl.BasicState(photons)
+    photons[0] = n - 1
+    photons[1] = 1
+    bs2 = pcvl.BasicState(photons)
+    circuit = pcvl.Circuit(m) // pcvl.components.Unitary(u1)
+
+    tvds = {}
+    result = []
+    refbsd = None
+    ref = "-"
+    for backend_name in backends:
+        # get_logger().warn(backend_name)
+        if backend_name == "SLAP" and n > 15:
+            print("\t-- ", end='', flush=True)
+            continue
+        if backend_name == "SLOS" and n > 18:
+            print("\t-- ", end='', flush=True)
+            continue
+        start = time.time()
+        backend = pcvl.backends.BackendFactory.get_backend(backend_name)
+        set_mask(backend, mask_str, m, n)
+        backend.set_circuit(circuit)
+        backend.set_input_state(bs)
+        bsd = do_compute(backend, method)
+        end = time.time()
+
+        backend.set_input_state(bs2)
+        bsd = do_compute(backend, method)
+        end2 = time.time()
+
+        # if backend_name == "SLOS_V2" or backend_name == "SLOS_V2_PS":
+        #     print(bsd)
+        # if refbsd:
+        #     tvds[backend_name] = tvd_dist(bsd, refbsd)
+        #     pass
+        # else:
+        #     refbsd = bsd
+        #     ref = backend_name
+        print(f"\t{end-start:.4f} {end2-end:.4f}", end='', flush=True)
+        result.append(end - start)
+    for k, v in tvds.items():
+        print(f"\n{k} TVD against {ref}: {v}", end='')
+
+    return result
+
+
+def bench_mem(m, n, backend_name, mask_str, method):
+    u1 = pcvl.Matrix.random_unitary(m)
+    photons = [ 0 ] * m
+    photons[0] = n
+    bs = pcvl.BasicState(photons)
+    circuit = pcvl.Circuit(m) // pcvl.components.Unitary(u1)
+
+    # if backend_name == "SLOS" and n > 12:
+    #     break
+    gc.collect()
+    mem_1 = psutil.Process().memory_info().rss
+    backend = pcvl.backends.BackendFactory.get_backend(backend_name)
+    set_mask(backend, mask_str, m, n)
+    backend.set_circuit(circuit)
+    backend.set_input_state(bs)
+    bsd = do_compute(backend, method)
+    mem_2 = psutil.Process().memory_info().rss
+    backend = None
+    bsd = None
+    gc.collect()
+    print(f"{backend_name:12}\t{human_readable_size(mem_2-mem_1)}")
+
+if __name__ == "__main__":
+    import argparse
+    parser = argparse.ArgumentParser(description='Benchmarks of strong simulation backends')
+    parser.add_argument('--modes', '-m',
+                        type=int, required=True, action='store',
+                        help='number of modes')
+    parser.add_argument('--nphotons', '-n',
+                        type=int, action='store', default=20,
+                        help='number of modes')
+    parser.add_argument('--minphotons', '-i',
+                        type=int, action='store', default=4,
+                        help='number of modes')
+    parser.add_argument('--memusage', '-u', action='store_true', default=False)
+    parser.add_argument('--backends', '-b',
+                        type=str, action='store', default='SLOS SLAP SLOS_CPP SLOS_V2 SLOS_V3',
+                        help='backend used')
+    parser.add_argument('--mask', '-k',
+                        type=str, action='store', default='',
+                        help='mask')
+    parser.add_argument('--compute', '-c',
+                        type=int, action='store', default='3',
+                        help='Computations: 0 -> all_prob_ampli / 1 -> all_prob / 2 -> prob_distribution')
+    args = parser.parse_args()
+
+    backend_list = args.backends.split(" ")
+    if args.memusage:
+        bench_mem(args.modes, args.nphotons, backend_list[0], args.mask, args.compute)
+    else:
+        print('\t', end = '', flush=False)
+        for b in backend_list:
+            print(f"{b:<16}", end = '', flush=False)
+        print('', flush=True)
+        for n in range(args.minphotons, args.nphotons + 1):
+            print(f"{n}", end='')
+            bench(args.modes, n, backend_list, args.mask, args.compute)
+            print("")

benchmarks-memory.sh

@@ -0,0 +1,46 @@
+
+#header#{{{#
+scriptPath="$( cd "$( dirname "${BASH_SOURCE[0]}" )" >/dev/null 2>&1 && pwd )"
+advancedPath="${scriptPath}/.."
+#}}}#
+
+##                                                                                                ##
+## arguments                                                                                      ##
+##                                                                                                ##
+
+modes=""
+nphotons="4"
+backends="SLOS SLAP SLOS_CPP SLOS_V2 SLOS_V3"
+mask=""
+compute_method=1
+
+# message strings
+USAGESTR="Usage: $0 [--modes/-m] 20 [--nphotons/-n] 10 [--backends/-b] \"SLOS_V2 SLOS_V3\" [--mask/-k] \"1010**\" [--compute/-c] 2"
+ERRORSTR="Failed to parse options... exiting."
+
+# option strings
+SHORTOPTIONS="m:n:i:b:k:c:"
+LONGOPTIONS="modes:,nphotons:,minphotons:,backends:,mask:,compute:"
+
+# read the options
+ARGS=$(getopt --options $SHORTOPTIONS --longoptions $LONGOPTIONS --name "$0" -- "$@") || exit 1
+eval "set -- $ARGS"
+
+# extract options and their arguments into variables
+while true ; do
+  case "$1" in
+    -h | --help ) echo "$USAGESTR" ; exit 0 ;;
+    -m | --modes ) modes="$2" ; shift 2 ;;
+    -n | --nphotons ) nphotons="$2" ; shift 2 ;;
+    -b | --backends ) backends="$2" ; shift 2 ;;
+    -k | --mask ) mask="$2" ; shift 2 ;;
+    -c | --compute ) compute="$2" ; shift 2 ;;
+    -- ) shift ; break ;;
+    *) echo "$ERRORSTR" >&2 ; exit 1 ;;
+  esac
+done
+
+for backend in $backends
+do
+    python ${scriptPath}/Bench_backends.py -u -m ${modes} -n ${nphotons} -b "${backend}" -c ${compute} -k "$mask"
+done

perceval/backends/__init__.py

@@ -34,6 +34,9 @@
 from ._naive_approx import NaiveApproxBackend
 from ._slos import SLOSBackend
 from ._slap import SLAPBackend
+from ._slos_v2 import SLOSV2Backend
+from ._slos_v3 import SLOSV3Backend
+from ._slos_cpp import SLOSCPPBackend
 
 
 BACKEND_LIST = {
@@ -42,7 +45,10 @@
     "Naive": NaiveBackend,
     "NaiveApprox": NaiveApproxBackend,
     "SLAP": SLAPBackend,
-    "SLOS": SLOSBackend
+    "SLOS": SLOSBackend,
+    "SLOS_CPP": SLOSCPPBackend,
+    "SLOS_V2": SLOSV2Backend,
+    "SLOS_V3": SLOSV3Backend
 }
 
 

perceval/backends/_slap.py

@@ -45,41 +45,41 @@ def __init__(self, mask=None):
 
     def set_circuit(self, circuit: ACircuit):
         super().set_circuit(circuit)  # Computes circuit unitary as _umat
-        self._slap.reset_feed_forward()
+        # self._slap.reset_feed_forward()
         self._slap.set_unitary(self._umat)
 
     def _init_mask(self):
         super()._init_mask()
-        if self._mask:
-            self._slap.set_mask(self._mask)
-        else:
-            self._slap.reset_mask()
+        self._slap.set_mask(self._mask)
 
     def prob_amplitude(self, output_state: FockState) -> complex:
         istate = self._input_state
-        all_pa = self._slap.all_prob_ampli(istate)
+        self._slap.set_input_state(self._input_state)
+        all_pa = self._slap.all_amplitudes()
         if self._mask:
             return all_pa[xq.FSArray(self._input_state.m, self._input_state.n, self._mask).find(output_state)]
         else:
             return all_pa[xq.FSArray(self._input_state.m, self._input_state.n).find(output_state)]
 
     def prob_distribution(self) -> BSDistribution:
-        return self._slap.prob_distribution(self._input_state)
+        self._slap.set_input_state(self._input_state)
+        return self._slap.distribution()
+
+    def all_prob_ampli(self):
+        self._slap.set_input_state(self._input_state)
+        return self._slap.all_amplitudes()
 
     @property
     def name(self) -> str:
         return "SLAP"
 
     def all_prob(self, input_state: FockState = None):
-        if input_state is not None:
-            self.set_input_state(input_state)
-        else:
-            input_state = self._input_state
-        return self._slap.all_prob(input_state)
+        self._slap.set_input_state(input_state or self._input_state)
+        return self._slap.all_probabilities()
 
     def evolve(self) -> StateVector:
-        istate = self._input_state
-        all_pa = self._slap.all_prob_ampli(istate)
+        self._slap.set_input_state(self._input_state)
+        all_pa = self._slap.all_amplitudes()
         res = StateVector()
         for output_state, pa in zip(self._get_iterator(self._input_state), all_pa):
             res += output_state * pa

perceval/backends/_slos.py

@@ -158,10 +158,10 @@ def _deploy(self, input_list: list[FockState]):
             n = input_state.n
             if n < len(self._fsms) and n not in self._fsas:
                 # we are missing the intermediate states - let us retrieve/load it back
-                current_fsa = xq.FSArray(m, n, self._mask) if self._mask else xq.FSArray(m, n)
+                current_fsa = xq.FSArray(m, n, self._mask) # if self._mask else xq.FSArray(m, n)
             for k in range(len(self._fsms), n + 1):
                 fsa_n_m1 = current_fsa
-                current_fsa = xq.FSArray(m, k, self._mask) if self._mask else xq.FSArray(m, k)
+                current_fsa = xq.FSArray(m, k, self._mask) # if self._mask else xq.FSArray(m, k)
                 self._mk_l.append(current_fsa.count())
                 self._fsms.append(xq.FSMap(current_fsa, fsa_n_m1, True))
             if n not in self._fsas:
@@ -192,6 +192,9 @@ def prob_amplitude(self, output_state: FockState) -> complex:
         result = self._state_mapping[self._input_state].coefs[output_idx, 0] * math.sqrt(output_state.prodnfact() / self._input_state.prodnfact())
         return result if self._symb else complex(result)
 
+    def all_prob_ampli(self):
+        return self._state_mapping[self._input_state].coefs.reshape(self._fsas[self._input_state.n].count())
+
     def prob_distribution(self) -> BSDistribution:
         istate = self._input_state
         c = self._state_mapping[istate].coefs.reshape(self._fsas[istate.n].count())

perceval/backends/_slos_cpp.py

@@ -0,0 +1,89 @@
+# MIT License
+#
+# Copyright (c) 2022 Quandela
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in all
+# copies or substantial portions of the Software.
+#
+# As a special exception, the copyright holders of exqalibur library give you
+# permission to combine exqalibur with code included in the standard release of
+# Perceval under the MIT license (or modified versions of such code). You may
+# copy and distribute such a combined system following the terms of the MIT
+# license for both exqalibur and Perceval. This exception for the usage of
+# exqalibur is limited to the python bindings used by Perceval.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+# SOFTWARE.
+import exqalibur as xq
+from perceval.utils import FockState, BSDistribution, StateVector
+from perceval.components import ACircuit
+from perceval.utils.states import BasicState
+
+from ._abstract_backends import AStrongSimulationBackend
+
+
+class SLOSCPPBackend(AStrongSimulationBackend):
+
+    def __init__(self, mask=None):
+        super().__init__()
+        self._slos = xq.SLOS()
+        self._fock_space = None
+        if mask is not None:
+            self.set_mask(mask)
+
+    def set_circuit(self, circuit: ACircuit):
+        super().set_circuit(circuit)  # Computes circuit unitary as _umat
+        self._slos.set_unitary(self._umat)
+
+    def set_input_state(self, input_state: BasicState):
+        super().set_input_state(input_state)
+        if self._fock_space is None or self._fock_space.m != input_state.m or self._fock_space.n != input_state.n:
+            self._fock_space = xq.FSArray(input_state.m, input_state.n, self._mask)
+
+    def _init_mask(self):
+        super()._init_mask()
+        self._slos.set_mask(self._mask)
+
+    def prob_amplitude(self, output_state: FockState) -> complex:
+        self._slos.set_input_state(self._input_state)
+        all_pa = self._slos.all_amplitudes()
+        return all_pa[self._fock_space.find(output_state)]
+
+    def prob_distribution(self) -> BSDistribution:
+        self._slos.set_input_state(self._input_state)
+
+        return self._slos.distribution()
+
+    def all_prob_ampli(self):
+        self._slos.set_input_state(self._input_state)
+        return self._slos.all_amplitudes()
+
+    @property
+    def name(self) -> str:
+        return "SLOS_CPP"
+
+    def all_prob(self, input_state: FockState = None):
+        self._slos.set_input_state(input_state or self._input_state)
+        return self._slos.all_probabilities()
+
+    def evolve(self) -> StateVector:
+        self._slos.set_input_state(self._input_state)
+        all_pa = self._slos.all_amplitudes()
+        res = StateVector()
+        for output_state, pa in zip(self._fock_space, all_pa):
+            # Utterly non-optimized. Mask management should be added in the computation
+            if self._mask is None or self._mask.match(output_state):
+                res += output_state * pa
+        return res