DynDEM: validation fixes and updated tests

Signed-off-by: Tracy Shannon <tshannon@nvidia.com>

Author: tlshannon

libs/qec/include/cudaq/qec/code_matrices.h

@@ -115,13 +115,69 @@ inline void check_num_cols(const sparse_binary_matrix &m, std::size_t n,
         ") does not match n_qubits (" + std::to_string(n) + ")");
 }
 
-/// Throw if per_qubit is non-empty with a length that does not match n.
+/// Throw if rates is non-empty with a length that does not match n.
+/// size_label names the dimension n counts, so a per-check vector does not
+/// report a per-qubit mismatch.
+inline void check_rate_vector_size(const std::vector<double> &rates,
+                                   std::size_t n, const char *label,
+                                   const char *size_label) {
+  if (!rates.empty() && rates.size() != n)
+    throw std::invalid_argument(std::string(label) + " has " +
+                                std::to_string(rates.size()) + " entries but " +
+                                size_label + "=" + std::to_string(n));
+}
+
+/// Throw if per_qubit is non-empty with a length that does not match n_qubits.
 inline void check_per_qubit_size(const std::vector<double> &per_qubit,
                                  std::size_t n, const char *label) {
-  if (!per_qubit.empty() && per_qubit.size() != n)
-    throw std::invalid_argument(std::string(label) + " has " +
-                                std::to_string(per_qubit.size()) +
-                                " entries but n_qubits=" + std::to_string(n));
+  check_rate_vector_size(per_qubit, n, label, "n_qubits");
+}
+
+/// Throw unless p is a probability. Written as a negated range test so that
+/// NaN, which compares false against everything, is rejected too.
+inline void check_probability(double p, const std::string &label) {
+  if (!(p >= 0.0 && p <= 1.0))
+    throw std::invalid_argument(label + " must be a probability in [0, 1]," +
+                                " got " + std::to_string(p));
+}
+
+/// Throw unless every scalar and per-element rate is a probability. Zero rates
+/// are legal and simply produce no DEM column; the point of this check is that
+/// a negative or NaN rate is silently inactive under that same rule, so a
+/// mistyped configuration would otherwise build a smaller DEM instead of
+/// failing.
+inline void validate_noise_rates(const css_noise_params &noise) {
+  check_probability(noise.px, "px");
+  check_probability(noise.py, "py");
+  check_probability(noise.pz, "pz");
+  check_probability(noise.pm, "pm");
+
+  const auto check_each = [](const std::vector<double> &rates,
+                             const char *label) {
+    for (std::size_t i = 0; i < rates.size(); ++i)
+      check_probability(rates[i],
+                        std::string(label) + "[" + std::to_string(i) + "]");
+  };
+  check_each(noise.px_per_qubit, "px_per_qubit");
+  check_each(noise.py_per_qubit, "py_per_qubit");
+  check_each(noise.pz_per_qubit, "pz_per_qubit");
+  check_each(noise.pm_per_check, "pm_per_check");
+}
+
+/// True when any noise rate is nonzero, i.e. the model asks for at least one
+/// fault mechanism. Call only after validate_noise_rates(), so that a nonzero
+/// rate here really is a positive probability.
+inline bool has_any_noise(const css_noise_params &noise) {
+  const auto any_nonzero = [](const std::vector<double> &rates) {
+    for (const double p : rates)
+      if (p != 0.0)
+        return true;
+    return false;
+  };
+  return noise.px != 0.0 || noise.py != 0.0 || noise.pz != 0.0 ||
+         noise.pm != 0.0 || any_nonzero(noise.px_per_qubit) ||
+         any_nonzero(noise.py_per_qubit) || any_nonzero(noise.pz_per_qubit) ||
+         any_nonzero(noise.pm_per_check);
 }
 
 /// Nested CSC for m padded to outer size n; new entries are empty vectors.

libs/qec/include/cudaq/qec/dem_construction.h

@@ -7,9 +7,13 @@
  ******************************************************************************/
 
 // This file declares dem_from_css_matrices(), which builds a T-round
-// code-capacity detector_error_model directly from CSS generator matrices
-// and a depolarizing noise model, without requiring a stabilizer circuit or
-// Stim round-trip. num_rounds defaults to 1 for the single-round case.
+// detector_error_model directly from CSS generator matrices and a noise
+// model, without requiring a stabilizer circuit or Stim round-trip.
+// num_rounds defaults to 1 for the single-round case.
+//
+// The model is code-capacity when only data-qubit rates are set, and
+// phenomenological when a measurement error rate (pm / pm_per_check) is set
+// as well, which adds one fault column per active check per round.
 
 #pragma once
 
@@ -18,8 +22,8 @@
 
 namespace cudaq::qec {
 
-/// @brief Build a T-round code-capacity detector_error_model from CSS
-/// generator matrices and a depolarizing noise model.
+/// @brief Build a T-round detector_error_model from CSS generator matrices
+/// and a depolarizing (optionally phenomenological) noise model.
 ///
 /// Models T independent syndrome measurement rounds. Detectors are syndrome
 /// *differences* between consecutive rounds: detector[r] fires when the
@@ -36,29 +40,35 @@ namespace cudaq::qec {
 /// @param code       CSS code matrices. All non-empty matrices must share the
 ///                   same num_cols() value (n_qubits). Default-constructed
 ///                   (zero-row, zero-column) matrices are treated as empty.
-/// @param noise      Depolarizing noise rates applied identically to every
-///                   round. Per-qubit vectors override the scalar rates;
-///                   qubits with effective rate 0 produce no column.
+/// @param noise      Noise rates applied identically to every round.
+///                   Per-element vectors override the scalar rates; elements
+///                   with effective rate 0 produce no column.
 /// @param num_rounds Number of syndrome measurement rounds T (default 1,
 ///                   must be >= 1).
 /// @return           detector_error_model with:
 ///   - detector_error_matrix: [T*d x e] where
 ///       d = hz.num_rows() + hx.num_rows(),
-///       e = T * (|active_X| + |active_Z| + |active_Y|).
+///       e = T * (|active_X| + |active_Z| + |active_Y| + |active_checks|).
+///     |active_checks| is the number of checks with a nonzero measurement
+///     error rate, and is 0 unless pm or pm_per_check is set.
 ///     Round r occupies rows r*d .. (r+1)*d-1. Within each round:
 ///       rows 0 .. hz.num_rows()-1 are Z-type detectors (X and Y faults);
 ///       rows hz.num_rows() .. d-1 are X-type detectors (Z and Y faults).
+///     A measurement-error column touches the one detector row of its own
+///     check in rounds r and r+1, whichever type that check is.
 ///   - observables_flips_matrix: [k x e],
 ///       k = lz.num_rows() + lx.num_rows().
 ///     Faults in any round flip the same observable rows (the logical
-///     measurement is taken once at the end of the experiment).
+///     measurement is taken once at the end of the experiment). Measurement
+///     errors flip no observable.
 ///   - error_rates: column layout is
 ///     [round 0 faults | ... | round T-1 faults]; within each round:
-///     [active X qubits | active Z qubits | active Y qubits], each in
-///     ascending qubit-index order.
+///     [active X qubits | active Z qubits | active Y qubits |
+///      active checks], each in ascending index order.
 /// @throws std::invalid_argument if num_rounds is 0, if two non-empty
-///         matrices have inconsistent num_cols() values, or if a per-qubit
-///         rate vector length does not equal n_qubits.
+///         matrices have inconsistent num_cols() values, if a per-qubit rate
+///         vector length does not equal n_qubits, if pm_per_check's length
+///         does not equal d, or if any rate is not a probability in [0, 1].
 detector_error_model dem_from_css_matrices(const css_code_matrices &code,
                                            const css_noise_params &noise,
                                            std::size_t num_rounds = 1);

libs/qec/include/cudaq/qec/extended_dem.h

@@ -87,6 +87,23 @@ struct extended_dem {
 
   /// Number of outgoing-seam rows, i.e. rows of out_syndrome.
   uint32_t num_out_seam_rows() const;
+
+  /// @brief Throw std::invalid_argument unless this chunk is internally
+  /// consistent, i.e. safe to stitch, close, or merge.
+  ///
+  /// num_faults() reports in_syndrome's width alone, so nothing else about a
+  /// chunk is self-describing: a block that disagrees with it, or a prior list
+  /// that does not have one entry per fault, would otherwise be read at the
+  /// wrong width and silently misalign columns (or walk past the end of a
+  /// nested column list). Checks that
+  ///   - every block is num_faults() columns wide, including zero-row blocks
+  ///     (use a width-n empty matrix, not a default-constructed one, when the
+  ///     chunk has faults),
+  ///   - fault_priors has one entry per fault, and
+  ///   - each tag vector has one entry per row of the seam it names.
+  ///
+  /// @param context Prefix for the error message, naming the caller.
+  void validate(const char *context) const;
 };
 
 /// @brief Build a one-round extended_dem from CSS matrices and noise.
@@ -141,9 +158,10 @@ struct dem_chunk_spec {
   /// True when nothing has been set, used to detect an omitted phase.
   bool is_empty() const;
 
-  /// @brief Check internal consistency: a positive fault count, one error rate
-  /// per fault with each in [0, 1], and index lists that are -1 terminated
-  /// with every index in [0, num_faults).
+  /// @brief Check internal consistency: a positive fault count that fits in
+  /// uint32_t (sparse matrix column index width), one error rate per fault with
+  /// each in [0, 1], and index lists that are -1 terminated with every index
+  /// in [0, num_faults).
   /// @param context Prefix for error messages, e.g. "dem_chunks.init".
   /// @throws std::invalid_argument on the first violation.
   void validate(const std::string &context) const;
@@ -260,12 +278,18 @@ extended_dem dem_stitch_all(const std::vector<extended_dem> &dem_chunks);
 /// into one. The merged prior is computed from the individual priors using
 /// one of these two rules:
 ///
-///   - or_combine (default): p_merged = 1 - prod_i(1 - p_i)
-///     Exact probability that at least one independent event fires.
-///     Use for physical fault mechanisms.
-///
-///   - sum_combine: p_merged = sum_i(p_i)
-///     Linear approximation valid when all p_i are small.
+///   - or_combine (default): p_merged = 1/2 * (1 - prod_i(1 - 2 p_i))
+///     Exact probability that an odd number of independent events fire, which
+///     is the net GF(2) effect of identical DEM columns (even counts cancel).
+///     Pairwise this is P(A xor B) = p + q - 2 p q, matching
+///     detector_error_model canonicalization. Prefer this for physical fault
+///     mechanisms.
+///
+///   - sum_combine: p_merged = min(1, sum_i(p_i))
+///     Linear approximation valid when all p_i are small. The sum of several
+///     larger priors can exceed 1, which is not a probability any decoder can
+///     use, so the result is clamped; prefer or_combine when the priors are
+///     not small.
 enum class prior_combine_mode { or_combine, sum_combine };
 
 /// @brief Merge fault columns with identical row support into single columns.
@@ -424,7 +448,10 @@ dem_chunks_to_o_sparse(const std::vector<extended_dem> &dem_chunks);
 ///
 /// Seams only have to contract pairwise, so a phase decomposition works: the
 /// first chunk may have no incoming seam (its interior carries round 0) and the
-/// last may have no outgoing one.
+/// last may have no outgoing one. As with dem_close(), the last chunk's
+/// out_syndrome is discarded: any detector that should appear in the closed
+/// DEM must already live in some chunk's in_syndrome or interior (for a
+/// dem_chunks_spec final phase, that means H_in_sparse / H_mid_sparse).
 ///
 /// @param dem_chunks Non-empty sequence of chunks in round order. Each chunk's
 ///               out_syndrome must match the next one's in_syndrome, and all
@@ -457,7 +484,15 @@ dem_chunks_to_pcm(const std::vector<extended_dem> &dem_chunks);
 ///   - detector_error_matrix: [in_syndrome stacked above interior]
 ///   - observables_flips_matrix: observables
 ///   - error_rates: fault_priors
-///   - out_syndrome rows are dropped (no final data measurement assumed).
+///
+/// out_syndrome is intentionally dropped. Closing models a terminated
+/// experiment: there is no later round for the outgoing seam to differ against,
+/// matching dem_from_css_matrices (final-round faults touch only the last
+/// detector band). Put any detector that must survive closing into in_syndrome
+/// or interior instead — for example a final data-readout boundary belongs in
+/// the last chunk's in_syndrome / interior (dem_chunks.final.H_in_sparse /
+/// H_mid_sparse), never only in out_syndrome. dem_chunks_spec::validate()
+/// already rejects a nonempty final.H_out_sparse for this reason.
 ///
 /// Invariant (up to canonicalization):
 ///   dem_close(dem_stitch_all(T one-round chunks))

libs/qec/include/cudaq/qec/realtime/decoding_config.h

@@ -91,6 +91,14 @@ struct decoder_config {
   /// requires knowing the round count up front; these phases describe one
   /// round each so the round count can be chosen (or grown) at run time. See
   /// cudaq::qec::dem_chunks_from_spec() for expansion to a chunk sequence.
+  ///
+  /// A configuration that also has a nonempty H_sparse is flat, and that
+  /// matrix is the one decoders are built from -- the phases are then only a
+  /// record of where it came from. Form selection keys off H_sparse.empty(),
+  /// so an omitted H_sparse and an explicit empty list both count as chunk
+  /// form. Nonempty H_sparse is exactly the state expand_dem_chunks() leaves
+  /// behind, which is what lets an expanded configuration round-trip through
+  /// YAML; it is not a way to override individual rounds.
   std::optional<dem_chunks_spec> dem_chunks;
   /// How many rounds to expand `dem_chunks` into. Required with dem_chunks and
   /// rejected without it. This is the round count the flat form would otherwise
@@ -190,8 +198,9 @@ class multi_decoder_config {
 /// therefore only has to understand the flat form.
 ///
 /// Does nothing to a configuration that is already flat (one whose `H_sparse`
-/// is set, or which carries no `dem_chunks` at all), so it is safe to call
-/// unconditionally.
+/// is nonempty, or which carries no `dem_chunks` at all), so it is safe to
+/// call unconditionally. An empty H_sparse with dem_chunks present is still
+/// treated as chunk form.
 ///
 /// @return The closed DEM the flat fields were derived from, so a caller that
 ///         also wants its per-fault priors does not have to expand a second

libs/qec/lib/dem_construction.cpp

@@ -29,6 +29,10 @@ detector_error_model dem_from_css_matrices(const css_code_matrices &code,
   if (num_rounds == 0)
     throw std::invalid_argument("num_rounds must be >= 1");
 
+  // Checked before the n == 0 early-out below so that a malformed rate is
+  // reported even when the code matrices describe no qubits to apply it to.
+  validate_noise_rates(noise);
+
   detector_error_model result;
   const std::size_t n = resolve_num_qubits(code);
 
@@ -62,7 +66,7 @@ detector_error_model dem_from_css_matrices(const css_code_matrices &code,
     const std::size_t d = nz + nx; // detector rows per round = total checks
 
     // pm_per_check must have length d (= nz + nx) when non-empty.
-    check_per_qubit_size(noise.pm_per_check, d, "pm_per_check");
+    check_rate_vector_size(noise.pm_per_check, d, "pm_per_check", "n_checks");
 
     const std::size_t n_detectors = num_rounds * d;
     const std::size_t n_observables = kz + kx;