Define max_step from minimum variations of samples, nsteps (#666)

* Define max_step from minimum variations of samples

* Determine nsteps

* Add false positive test

* Delete draw

* Fix final_state

pulser-simulation/pulser_simulation/simulation.py

@@ -32,7 +32,7 @@
 from pulser.devices._device_datacls import BaseDevice
 from pulser.register.base_register import BaseRegister
 from pulser.result import SampledResult
-from pulser.sampler.samples import SequenceSamples
+from pulser.sampler.samples import ChannelSamples, SequenceSamples
 from pulser.sequence._seq_drawer import draw_samples, draw_sequence
 from pulser_simulation.hamiltonian import Hamiltonian
 from pulser_simulation.qutip_result import QutipResult
@@ -492,19 +492,39 @@ def run(
 
                 .. _docs: https://bit.ly/3il9A2u
         """
-        if "max_step" not in options:
-            pulse_durations = [
-                slot.tf - slot.ti
-                for ch_sample in self.samples_obj.samples_list
-                for slot in ch_sample.slots
-                if not (
-                    np.all(np.isclose(ch_sample.amp[slot.ti : slot.tf], 0))
-                    and np.all(np.isclose(ch_sample.det[slot.ti : slot.tf], 0))
+
+        def get_min_variation(ch_sample: ChannelSamples) -> int:
+            end_point = ch_sample.duration - 1
+            min_variations: list[int] = []
+            for sample in (ch_sample.amp, ch_sample.det):
+                min_variations.append(
+                    int(
+                        np.min(
+                            np.diff(
+                                np.nonzero(np.diff(sample)),
+                                prepend=-1,
+                                append=end_point,
+                            )
+                        )
+                    )
                 )
-            ]
-            if pulse_durations:
-                options["max_step"] = 0.5 * min(pulse_durations) / 1000
 
+            return min(min_variations)
+
+        if "max_step" not in options:
+            options["max_step"] = (
+                min(
+                    [
+                        get_min_variation(ch_sample)
+                        for ch_sample in self.samples_obj.samples_list
+                    ]
+                )
+                / 1000
+            )
+        if "nsteps" not in options:
+            options["nsteps"] = max(
+                1000, self._tot_duration // options["max_step"]
+            )
         solv_ops = qutip.Options(**options)
 
         meas_errors: Optional[Mapping[str, float]] = None

tests/test_simresults.py

@@ -19,9 +19,13 @@
 import qutip
 from qutip.piqs import isdiagonal
 
-from pulser import Pulse, Register, Sequence
+from pulser import AnalogDevice, Pulse, Register, Sequence
 from pulser.devices import DigitalAnalogDevice, MockDevice
-from pulser.waveforms import BlackmanWaveform
+from pulser.waveforms import (
+    BlackmanWaveform,
+    CompositeWaveform,
+    ConstantWaveform,
+)
 from pulser_simulation import QutipEmulator, SimConfig
 from pulser_simulation.simresults import CoherentResults, NoisyResults
 
@@ -213,7 +217,7 @@ def test_get_state_float_time(results):
                 [0.76522907 + 0.0j],
                 [0.08339973 - 0.39374219j],
                 [0.08339973 - 0.39374219j],
-                [-0.27977623 - 0.1103308j],
+                [-0.27977172 - 0.11031832j],
             ]
         ),
     ).all()
@@ -386,3 +390,23 @@ def test_results_xy(reg, pi_pulse):
     # Check that measurement projectors are correct
     assert results_._meas_projector(0) == qutip.basis(2, 0).proj()
     assert results_._meas_projector(1) == qutip.basis(2, 1).proj()
+
+
+def test_false_positive():
+    """Breaks for pulser version < v0.18."""
+    seq = Sequence(Register.square(2, 5), AnalogDevice)
+    seq.declare_channel("ryd_glob", "rydberg_global")
+    seq.add(
+        Pulse.ConstantDetuning(
+            CompositeWaveform(
+                ConstantWaveform(2500, 0.0),
+                BlackmanWaveform(1000, np.pi),
+                ConstantWaveform(500, 0.0),
+            ),
+            0,
+            0,
+        ),
+        channel="ryd_glob",
+    )
+    sim = QutipEmulator.from_sequence(seq)
+    assert sim.run().get_final_state() != sim.initial_state