Simulation & Monte Carlo API

simulation

╔══════════════════════════════════════════════════════════════════════════════╗ ║ MHRQI - Monte Carlo Simulation Backend ║ ║ Efficient shot-based sampling for hierarchical quantum image processing ║ ║ ║ ║ Author: Keno S. Jose ║ ║ License: Apache 2.0 ║ ╚══════════════════════════════════════════════════════════════════════════════╝

HierarchicalMeasurementAggregator

Aggregates Monte Carlo measurements across hierarchical levels, computing statistics for binning and reconstruction.

Source code in mhrqi/core/simulation.py

class HierarchicalMeasurementAggregator:
    """
    Aggregates Monte Carlo measurements across hierarchical levels,
    computing statistics for binning and reconstruction.
    """

    def __init__(self, hierarchical_coord_matrix: list, bit_depth: int = 8):
        """
        Initialize aggregator.

        Args:
            hierarchical_coord_matrix: List of HCV vectors (one per pixel).
            bit_depth: Bits per intensity value.
        """
        self.hc_matrix = hierarchical_coord_matrix
        self.bit_depth = bit_depth
        self.bins = defaultdict(list)

    def aggregate(self, measurement_counts: Dict[str, int]) -> None:
        """
        Aggregate measurement outcomes into hierarchical bins.

        Args:
            measurement_counts: Dictionary from Monte Carlo sampling.
        """
        for outcome_str, count in measurement_counts.items():
            # Parse outcome: position bits + intensity bits
            n_position_bits = len(self.hc_matrix[0]) if self.hc_matrix else 0

            if len(outcome_str) >= n_position_bits + self.bit_depth:
                pos_bits = outcome_str[:n_position_bits]
                intensity_bits = outcome_str[n_position_bits:n_position_bits + self.bit_depth]

                # Store multiple copies (one per shot)
                for _ in range(count):
                    self.bins[pos_bits].append(intensity_bits)

    def get_statistics(self, pos_bits: str) -> Dict[str, float]:
        """
        Compute statistics for a hierarchical position.

        Args:
            pos_bits: Position bits (binary string).

        Returns:
            Dictionary with mean, std, median, mode.
        """
        if pos_bits not in self.bins or len(self.bins[pos_bits]) == 0:
            return {'mean': 0.0, 'std': 0.0, 'median': 0.0, 'mode': '0'}

        # Convert intensity bits to integer values
        intensity_values = np.array([
            int(bits, 2) for bits in self.bins[pos_bits]
        ])

        # Compute statistics
        from scipy import stats
        mode_result = stats.mode(intensity_values, keepdims=True)

        return {
            'mean': float(np.mean(intensity_values)),
            'std': float(np.std(intensity_values)),
            'median': float(np.median(intensity_values)),
            'mode': int(mode_result.mode[0]),
            'samples': len(self.bins[pos_bits]),
        }

    def confidence_weighted_value(self, pos_bits: str, context_value: float = None) -> float:
        """
        Compute confidence-weighted reconstruction value.

        Args:
            pos_bits: Position bits.
            context_value: Value from hierarchical context (sibling average, etc.).

        Returns:
            Blended reconstruction value.
        """
        stats = self.get_statistics(pos_bits)

        if stats['samples'] == 0:
            return context_value if context_value is not None else 0.0

        # Confidence based on measurement variance
        # Lower variance → higher confidence in measurement
        variance = stats['std'] ** 2
        # Normalize by bit scale
        max_variance = (2 ** self.bit_depth) ** 2
        confidence = np.exp(-variance / max_variance)

        measurement = stats['mean'] / (2 ** self.bit_depth - 1)  # Normalize to [0, 1]

        if context_value is not None:
            return confidence * measurement + (1 - confidence) * context_value

        return measurement

__init__(hierarchical_coord_matrix, bit_depth=8)

Initialize aggregator.

Parameters:

Name	Type	Description	Default
hierarchical_coord_matrix	list	List of HCV vectors (one per pixel).	required
bit_depth	int	Bits per intensity value.	8

Source code in mhrqi/core/simulation.py

def __init__(self, hierarchical_coord_matrix: list, bit_depth: int = 8):
    """
    Initialize aggregator.

    Args:
        hierarchical_coord_matrix: List of HCV vectors (one per pixel).
        bit_depth: Bits per intensity value.
    """
    self.hc_matrix = hierarchical_coord_matrix
    self.bit_depth = bit_depth
    self.bins = defaultdict(list)

aggregate(measurement_counts)

Aggregate measurement outcomes into hierarchical bins.

Parameters:

Name	Type	Description	Default
measurement_counts	Dict[str, int]	Dictionary from Monte Carlo sampling.	required

Source code in mhrqi/core/simulation.py

def aggregate(self, measurement_counts: Dict[str, int]) -> None:
    """
    Aggregate measurement outcomes into hierarchical bins.

    Args:
        measurement_counts: Dictionary from Monte Carlo sampling.
    """
    for outcome_str, count in measurement_counts.items():
        # Parse outcome: position bits + intensity bits
        n_position_bits = len(self.hc_matrix[0]) if self.hc_matrix else 0

        if len(outcome_str) >= n_position_bits + self.bit_depth:
            pos_bits = outcome_str[:n_position_bits]
            intensity_bits = outcome_str[n_position_bits:n_position_bits + self.bit_depth]

            # Store multiple copies (one per shot)
            for _ in range(count):
                self.bins[pos_bits].append(intensity_bits)

confidence_weighted_value(pos_bits, context_value=None)

Compute confidence-weighted reconstruction value.

Parameters:

Name	Type	Description	Default
pos_bits	str	Position bits.	required
context_value	float	Value from hierarchical context (sibling average, etc.).	None

Returns:

Type	Description
float	Blended reconstruction value.

Source code in mhrqi/core/simulation.py

def confidence_weighted_value(self, pos_bits: str, context_value: float = None) -> float:
    """
    Compute confidence-weighted reconstruction value.

    Args:
        pos_bits: Position bits.
        context_value: Value from hierarchical context (sibling average, etc.).

    Returns:
        Blended reconstruction value.
    """
    stats = self.get_statistics(pos_bits)

    if stats['samples'] == 0:
        return context_value if context_value is not None else 0.0

    # Confidence based on measurement variance
    # Lower variance → higher confidence in measurement
    variance = stats['std'] ** 2
    # Normalize by bit scale
    max_variance = (2 ** self.bit_depth) ** 2
    confidence = np.exp(-variance / max_variance)

    measurement = stats['mean'] / (2 ** self.bit_depth - 1)  # Normalize to [0, 1]

    if context_value is not None:
        return confidence * measurement + (1 - confidence) * context_value

    return measurement

get_statistics(pos_bits)

Compute statistics for a hierarchical position.

Parameters:

Name	Type	Description	Default
pos_bits	str	Position bits (binary string).	required

Returns:

Type	Description
Dict[str, float]	Dictionary with mean, std, median, mode.

Source code in mhrqi/core/simulation.py

def get_statistics(self, pos_bits: str) -> Dict[str, float]:
    """
    Compute statistics for a hierarchical position.

    Args:
        pos_bits: Position bits (binary string).

    Returns:
        Dictionary with mean, std, median, mode.
    """
    if pos_bits not in self.bins or len(self.bins[pos_bits]) == 0:
        return {'mean': 0.0, 'std': 0.0, 'median': 0.0, 'mode': '0'}

    # Convert intensity bits to integer values
    intensity_values = np.array([
        int(bits, 2) for bits in self.bins[pos_bits]
    ])

    # Compute statistics
    from scipy import stats
    mode_result = stats.mode(intensity_values, keepdims=True)

    return {
        'mean': float(np.mean(intensity_values)),
        'std': float(np.std(intensity_values)),
        'median': float(np.median(intensity_values)),
        'mode': int(mode_result.mode[0]),
        'samples': len(self.bins[pos_bits]),
    }

MonteCarloSimulator

Efficient Monte Carlo sampler for MHRQI quantum simulations.

Provides configurable shot-based sampling with GPU acceleration support, statistical aggregation, and reproducibility via seed management.

Source code in mhrqi/core/simulation.py

class MonteCarloSimulator:
    """
    Efficient Monte Carlo sampler for MHRQI quantum simulations.

    Provides configurable shot-based sampling with GPU acceleration support,
    statistical aggregation, and reproducibility via seed management.
    """

    def __init__(self, seed: Optional[int] = None, use_gpu: bool = True):
        """
        Initialize Monte Carlo simulator.

        Args:
            seed: Random seed for reproducibility. If None, uses system entropy.
            use_gpu: Attempt GPU acceleration via cuStateVec if available.
        """
        self.seed = seed
        self.use_gpu = use_gpu
        self.rng = np.random.RandomState(seed)
        self._gpu_available = False

        # Check GPU availability
        if use_gpu:
            try:
                import pycuda.driver as cuda
                cuda.init()
                self._gpu_available = True
            except (ImportError, RuntimeError):
                self._gpu_available = False
                if use_gpu:
                    warnings.warn("GPU not available; falling back to CPU simulation.", stacklevel=2)

    @property
    def gpu_available(self) -> bool:
        """Check if GPU acceleration is available."""
        return self._gpu_available

    def sample_statevector(
        self,
        statevector: np.ndarray,
        shots: int,
        measured_qubits: Optional[list] = None,
    ) -> Dict[str, int]:
        """
        Sample measurement outcomes from statevector via Monte Carlo method.

        Args:
            statevector: Full quantum state vector (shape: 2^n_qubits).
            shots: Number of measurement repetitions.
            measured_qubits: Indices of qubits to measure. If None, measures all.

        Returns:
            Dictionary mapping measurement outcomes (binary strings) to counts.
        """
        if measured_qubits is None:
            n_qubits = int(np.log2(len(statevector)))
            measured_qubits = list(range(n_qubits))

        # Compute probabilities
        probabilities = np.abs(statevector) ** 2

        # Monte Carlo sampling
        outcomes = self.rng.choice(
            len(probabilities),
            size=shots,
            p=probabilities / probabilities.sum()
        )

        # Convert indices to binary strings and aggregate
        counts = defaultdict(int)
        for outcome_idx in outcomes:
            outcome_binary = format(outcome_idx, f'0{len(measured_qubits)}b')
            counts[outcome_binary] += 1

        return dict(counts)

    def bootstrap_statistics(
        self,
        samples: np.ndarray,
        n_bootstraps: int = 1000,
        confidence: float = 0.95,
    ) -> Dict[str, float]:
        """
        Compute bootstrap confidence intervals for measurement statistics.

        Args:
            samples: Measured sample values.
            n_bootstraps: Number of bootstrap resamples.
            confidence: Confidence level (e.g., 0.95 for 95% CI).

        Returns:
            Dictionary with mean, std, CI_lower, CI_upper.
        """
        bootstrap_means = []
        for _ in range(n_bootstraps):
            resample = self.rng.choice(samples, size=len(samples), replace=True)
            bootstrap_means.append(np.mean(resample))

        bootstrap_means = np.array(bootstrap_means)
        alpha = (1 - confidence) / 2

        return {
            'mean': np.mean(samples),
            'std': np.std(samples),
            'ci_lower': np.percentile(bootstrap_means, alpha * 100),
            'ci_upper': np.percentile(bootstrap_means, (1 - alpha) * 100),
            'bootstrap_samples': bootstrap_means,
        }

    def estimate_convergence(
        self,
        statevector: np.ndarray,
        max_shots: int = 10000,
        step_size: int = 100,
    ) -> Dict[int, float]:
        """
        Estimate Monte Carlo convergence by increasing shot count.

        Args:
            statevector: Full quantum state vector.
            max_shots: Maximum number of shots to simulate.
            step_size: Shots per convergence step.

        Returns:
            Dictionary mapping shot counts to estimated errors.
        """
        true_probs = np.abs(statevector) ** 2
        convergence = {}

        for shots in range(step_size, max_shots + 1, step_size):
            # Estimate probabilities via sampling
            outcomes = self.rng.choice(
                len(statevector),
                size=shots,
                p=true_probs
            )
            estimated_probs = np.bincount(outcomes, minlength=len(statevector)) / shots

            # Compute KL divergence (proxy for estimation error)
            kl_div = np.sum(
                true_probs[true_probs > 0] * 
                np.log(true_probs[true_probs > 0] / (estimated_probs[true_probs > 0] + 1e-10))
            )
            convergence[shots] = kl_div

        return convergence

    def adaptive_shot_allocation(
        self,
        statevector: np.ndarray,
        target_error: float = 0.01,
        max_shots: int = 100000,
    ) -> int:
        """
        Determine minimum shots needed to achieve target error via binary search.

        Args:
            statevector: Full quantum state vector.
            target_error: Target KL divergence.
            max_shots: Maximum shots to try.

        Returns:
            Recommended shot count.
        """
        convergence = self.estimate_convergence(statevector, max_shots // 100, max_shots // 100)

        for shots in sorted(convergence.keys()):
            if convergence[shots] <= target_error:
                return shots

        return max_shots

gpu_available property

Check if GPU acceleration is available.

__init__(seed=None, use_gpu=True)

Initialize Monte Carlo simulator.

Parameters:

Name	Type	Description	Default
seed	Optional[int]	Random seed for reproducibility. If None, uses system entropy.	None
use_gpu	bool	Attempt GPU acceleration via cuStateVec if available.	True

Source code in mhrqi/core/simulation.py

def __init__(self, seed: Optional[int] = None, use_gpu: bool = True):
    """
    Initialize Monte Carlo simulator.

    Args:
        seed: Random seed for reproducibility. If None, uses system entropy.
        use_gpu: Attempt GPU acceleration via cuStateVec if available.
    """
    self.seed = seed
    self.use_gpu = use_gpu
    self.rng = np.random.RandomState(seed)
    self._gpu_available = False

    # Check GPU availability
    if use_gpu:
        try:
            import pycuda.driver as cuda
            cuda.init()
            self._gpu_available = True
        except (ImportError, RuntimeError):
            self._gpu_available = False
            if use_gpu:
                warnings.warn("GPU not available; falling back to CPU simulation.", stacklevel=2)

adaptive_shot_allocation(statevector, target_error=0.01, max_shots=100000)

Determine minimum shots needed to achieve target error via binary search.

Parameters:

Name	Type	Description	Default
statevector	ndarray	Full quantum state vector.	required
target_error	float	Target KL divergence.	0.01
max_shots	int	Maximum shots to try.	100000

Returns:

Type	Description
int	Recommended shot count.

Source code in mhrqi/core/simulation.py

def adaptive_shot_allocation(
    self,
    statevector: np.ndarray,
    target_error: float = 0.01,
    max_shots: int = 100000,
) -> int:
    """
    Determine minimum shots needed to achieve target error via binary search.

    Args:
        statevector: Full quantum state vector.
        target_error: Target KL divergence.
        max_shots: Maximum shots to try.

    Returns:
        Recommended shot count.
    """
    convergence = self.estimate_convergence(statevector, max_shots // 100, max_shots // 100)

    for shots in sorted(convergence.keys()):
        if convergence[shots] <= target_error:
            return shots

    return max_shots

bootstrap_statistics(samples, n_bootstraps=1000, confidence=0.95)

Compute bootstrap confidence intervals for measurement statistics.

Parameters:

Name	Type	Description	Default
samples	ndarray	Measured sample values.	required
n_bootstraps	int	Number of bootstrap resamples.	1000
confidence	float	Confidence level (e.g., 0.95 for 95% CI).	0.95

Returns:

Type	Description
Dict[str, float]	Dictionary with mean, std, CI_lower, CI_upper.

Source code in mhrqi/core/simulation.py

def bootstrap_statistics(
    self,
    samples: np.ndarray,
    n_bootstraps: int = 1000,
    confidence: float = 0.95,
) -> Dict[str, float]:
    """
    Compute bootstrap confidence intervals for measurement statistics.

    Args:
        samples: Measured sample values.
        n_bootstraps: Number of bootstrap resamples.
        confidence: Confidence level (e.g., 0.95 for 95% CI).

    Returns:
        Dictionary with mean, std, CI_lower, CI_upper.
    """
    bootstrap_means = []
    for _ in range(n_bootstraps):
        resample = self.rng.choice(samples, size=len(samples), replace=True)
        bootstrap_means.append(np.mean(resample))

    bootstrap_means = np.array(bootstrap_means)
    alpha = (1 - confidence) / 2

    return {
        'mean': np.mean(samples),
        'std': np.std(samples),
        'ci_lower': np.percentile(bootstrap_means, alpha * 100),
        'ci_upper': np.percentile(bootstrap_means, (1 - alpha) * 100),
        'bootstrap_samples': bootstrap_means,
    }

estimate_convergence(statevector, max_shots=10000, step_size=100)

Estimate Monte Carlo convergence by increasing shot count.

Parameters:

Name	Type	Description	Default
statevector	ndarray	Full quantum state vector.	required
max_shots	int	Maximum number of shots to simulate.	10000
step_size	int	Shots per convergence step.	100

Returns:

Type	Description
Dict[int, float]	Dictionary mapping shot counts to estimated errors.

Source code in mhrqi/core/simulation.py

def estimate_convergence(
    self,
    statevector: np.ndarray,
    max_shots: int = 10000,
    step_size: int = 100,
) -> Dict[int, float]:
    """
    Estimate Monte Carlo convergence by increasing shot count.

    Args:
        statevector: Full quantum state vector.
        max_shots: Maximum number of shots to simulate.
        step_size: Shots per convergence step.

    Returns:
        Dictionary mapping shot counts to estimated errors.
    """
    true_probs = np.abs(statevector) ** 2
    convergence = {}

    for shots in range(step_size, max_shots + 1, step_size):
        # Estimate probabilities via sampling
        outcomes = self.rng.choice(
            len(statevector),
            size=shots,
            p=true_probs
        )
        estimated_probs = np.bincount(outcomes, minlength=len(statevector)) / shots

        # Compute KL divergence (proxy for estimation error)
        kl_div = np.sum(
            true_probs[true_probs > 0] * 
            np.log(true_probs[true_probs > 0] / (estimated_probs[true_probs > 0] + 1e-10))
        )
        convergence[shots] = kl_div

    return convergence

sample_statevector(statevector, shots, measured_qubits=None)

Sample measurement outcomes from statevector via Monte Carlo method.

Parameters:

Name	Type	Description	Default
statevector	ndarray	Full quantum state vector (shape: 2^n_qubits).	required
shots	int	Number of measurement repetitions.	required
measured_qubits	Optional[list]	Indices of qubits to measure. If None, measures all.	None

Returns:

Type	Description
Dict[str, int]	Dictionary mapping measurement outcomes (binary strings) to counts.

Source code in mhrqi/core/simulation.py

def sample_statevector(
    self,
    statevector: np.ndarray,
    shots: int,
    measured_qubits: Optional[list] = None,
) -> Dict[str, int]:
    """
    Sample measurement outcomes from statevector via Monte Carlo method.

    Args:
        statevector: Full quantum state vector (shape: 2^n_qubits).
        shots: Number of measurement repetitions.
        measured_qubits: Indices of qubits to measure. If None, measures all.

    Returns:
        Dictionary mapping measurement outcomes (binary strings) to counts.
    """
    if measured_qubits is None:
        n_qubits = int(np.log2(len(statevector)))
        measured_qubits = list(range(n_qubits))

    # Compute probabilities
    probabilities = np.abs(statevector) ** 2

    # Monte Carlo sampling
    outcomes = self.rng.choice(
        len(probabilities),
        size=shots,
        p=probabilities / probabilities.sum()
    )

    # Convert indices to binary strings and aggregate
    counts = defaultdict(int)
    for outcome_idx in outcomes:
        outcome_binary = format(outcome_idx, f'0{len(measured_qubits)}b')
        counts[outcome_binary] += 1

    return dict(counts)

configure_backend(use_gpu=True, seed=None)

Factory function to configure optimal simulation backend.

Parameters:

Name	Type	Description	Default
use_gpu	bool	Attempt GPU acceleration.	True
seed	Optional[int]	Random seed.	None

Returns:

Type	Description
MonteCarloSimulator	Configured MonteCarloSimulator instance.

Source code in mhrqi/core/simulation.py

def configure_backend(use_gpu: bool = True, seed: Optional[int] = None) -> MonteCarloSimulator:
    """
    Factory function to configure optimal simulation backend.

    Args:
        use_gpu: Attempt GPU acceleration.
        seed: Random seed.

    Returns:
        Configured MonteCarloSimulator instance.
    """
    return MonteCarloSimulator(seed=seed, use_gpu=use_gpu)

monte_carlo

╔══════════════════════════════════════════════════════════════════════════════╗ ║ MHRQI - Monte Carlo Utilities ║ ║ Helper functions for Monte Carlo sampling and statistical analysis ║ ║ ║ ║ Author: Keno S. Jose ║ ║ License: Apache 2.0 ║ ╚══════════════════════════════════════════════════════════════════════════════╝

aggregate_multi_run_results(results, method='average')

Aggregate multiple Monte Carlo runs for improved statistics.

Parameters:

Name	Type	Description	Default
results	List[Dict[int, int]]	List of count dictionaries from separate runs.	required
method	str	'average' (mean estimates) or 'pooled' (combine all samples).	'average'

Returns:

Type	Description
Dict[int, float]	Aggregated probability estimates.

Source code in mhrqi/utils/monte_carlo.py

def aggregate_multi_run_results(
    results: List[Dict[int, int]],
    method: str = 'average',
) -> Dict[int, float]:
    """
    Aggregate multiple Monte Carlo runs for improved statistics.

    Args:
        results: List of count dictionaries from separate runs.
        method: 'average' (mean estimates) or 'pooled' (combine all samples).

    Returns:
        Aggregated probability estimates.
    """
    if method == 'pooled':
        # Combine all shots as single distribution
        combined = {}
        total_shots = 0
        for result in results:
            for outcome, count in result.items():
                combined[outcome] = combined.get(outcome, 0) + count
                total_shots += count

        for outcome in combined:
            combined[outcome] /= total_shots

        return combined

    elif method == 'average':
        # Average probabilities across runs
        averaged = {}
        for result in results:
            total = sum(result.values())
            for outcome, count in result.items():
                if outcome not in averaged:
                    averaged[outcome] = []
                averaged[outcome].append(count / total)

        final = {}
        for outcome, probs in averaged.items():
            final[outcome] = np.mean(probs)

        return final

    else:
        raise ValueError(f"Unknown aggregation method: {method}")

auto_correlate_length(samples, max_lag=None)

Estimate autocorrelation length to assess sampling efficiency.

Parameters:

Name	Type	Description	Default
samples	ndarray	Time series of measurements.	required
max_lag	int	Maximum lag to check. If None, uses 10% of series length.	None

Returns:

Type	Description
float	Autocorrelation length (ACL).

Source code in mhrqi/utils/monte_carlo.py

def auto_correlate_length(samples: np.ndarray, max_lag: int = None) -> float:
    """
    Estimate autocorrelation length to assess sampling efficiency.

    Args:
        samples: Time series of measurements.
        max_lag: Maximum lag to check. If None, uses 10% of series length.

    Returns:
        Autocorrelation length (ACL).
    """
    if max_lag is None:
        max_lag = len(samples) // 10

    mean = np.mean(samples)
    c0 = np.mean((samples - mean) ** 2)

    acl = 0.5
    for lag in range(1, max_lag):
        c_lag = np.mean((samples[:-lag] - mean) * (samples[lag:] - mean))
        rho_lag = c_lag / c0
        acl += rho_lag

        if rho_lag < 0.05:  # Stop when autocorrelation negligible
            break

    return acl

compute_statistical_error(samples, confidence=0.95)

Compute mean, standard error, and confidence interval.

Parameters:

Name	Type	Description	Default
samples	ndarray	Measurement samples.	required
confidence	float	Confidence level.	0.95

Returns:

Type	Description
Tuple[float, float, float]	Tuple of (mean, std_error, ci_half_width).

Source code in mhrqi/utils/monte_carlo.py

def compute_statistical_error(
    samples: np.ndarray,
    confidence: float = 0.95,
) -> Tuple[float, float, float]:
    """
    Compute mean, standard error, and confidence interval.

    Args:
        samples: Measurement samples.
        confidence: Confidence level.

    Returns:
        Tuple of (mean, std_error, ci_half_width).
    """
    from scipy import stats

    mean = np.mean(samples)
    se = np.std(samples) / np.sqrt(len(samples))
    ci_hw = stats.t.ppf((1 + confidence) / 2, len(samples) - 1) * se

    return mean, se, ci_hw

effective_sample_size(weights)

Compute effective sample size for importance-weighted samples.

ESS = (∑ w_i)² / ∑ w_i²

Parameters:

Name	Type	Description	Default
weights	ndarray	Importance weights.	required

Returns:

Type	Description
float	Effective sample size.

Source code in mhrqi/utils/monte_carlo.py

def effective_sample_size(weights: np.ndarray) -> float:
    """
    Compute effective sample size for importance-weighted samples.

    ESS = (∑ w_i)² / ∑ w_i²

    Args:
        weights: Importance weights.

    Returns:
        Effective sample size.
    """
    weights = np.asarray(weights)
    weights_normalized = weights / np.sum(weights)
    ess = np.sum(weights_normalized) ** 2 / np.sum(weights_normalized ** 2)
    return ess

estimate_required_shots(target_accuracy=0.95, target_error=0.05, num_outcomes=None)

Estimate Monte Carlo shots needed for target accuracy.

Uses Hoeffding inequality: P(|hat_p - p| > ε) ≤ 2 * exp(-2 * n * ε^2)

Parameters:

Name	Type	Description	Default
target_accuracy	float	Desired confidence level (e.g., 0.95).	0.95
target_error	float	Maximum error tolerance (e.g., 0.05).	0.05
num_outcomes	int	Number of possible measurement outcomes. If None, computes conservative estimate.	None

Returns:

Type	Description
int	Recommended number of shots.

Source code in mhrqi/utils/monte_carlo.py

def estimate_required_shots(
    target_accuracy: float = 0.95,
    target_error: float = 0.05,
    num_outcomes: int = None,
) -> int:
    """
    Estimate Monte Carlo shots needed for target accuracy.

    Uses Hoeffding inequality: P(|hat_p - p| > ε) ≤ 2 * exp(-2 * n * ε^2)

    Args:
        target_accuracy: Desired confidence level (e.g., 0.95).
        target_error: Maximum error tolerance (e.g., 0.05).
        num_outcomes: Number of possible measurement outcomes. If None,
                      computes conservative estimate.

    Returns:
        Recommended number of shots.
    """
    # Hoeffding bound: n ≥ ln(2/δ) / (2ε²)
    delta = 1 - target_accuracy
    shots = int(np.log(2 / delta) / (2 * target_error ** 2))
    return shots

importance_sampling(target_probs, proposal_probs, shots, rng=None)

Importance sampling to estimate target distribution from proposal.

Parameters:

Name	Type	Description	Default
target_probs	ndarray	Target probability distribution.	required
proposal_probs	ndarray	Proposal (sampling) distribution.	required
shots	int	Number of proposal samples.	required
rng	RandomState	Random state for reproducibility.	None

Returns:

Type	Description
Dict[int, float]	Dictionary of weighted estimates.

Source code in mhrqi/utils/monte_carlo.py

def importance_sampling(
    target_probs: np.ndarray,
    proposal_probs: np.ndarray,
    shots: int,
    rng: np.random.RandomState = None,
) -> Dict[int, float]:
    """
    Importance sampling to estimate target distribution from proposal.

    Args:
        target_probs: Target probability distribution.
        proposal_probs: Proposal (sampling) distribution.
        shots: Number of proposal samples.
        rng: Random state for reproducibility.

    Returns:
        Dictionary of weighted estimates.
    """
    if rng is None:
        rng = np.random.RandomState()

    # Sample from proposal
    outcomes = rng.choice(
        len(proposal_probs),
        size=shots,
        p=proposal_probs / proposal_probs.sum()
    )

    # Compute importance weights
    weights = {}
    for outcome in outcomes:
        if proposal_probs[outcome] > 0:
            w = target_probs[outcome] / proposal_probs[outcome]
            weights[outcome] = weights.get(outcome, 0.0) + w

    # Normalize
    total_weight = sum(weights.values())
    for outcome in weights:
        weights[outcome] /= total_weight

    return weights

shot_count_for_precision(std_dev, target_uncertainty, confidence=0.95)

Determine shots needed for desired measurement precision.

Uses standard error formula: SE = σ / √n

Parameters:

Name	Type	Description	Default
std_dev	float	Estimated standard deviation of measurements.	required
target_uncertainty	float	Desired standard error.	required
confidence	float	Confidence level (affects z-score).	0.95

Returns:

Type	Description
int	Recommended number of shots.

Source code in mhrqi/utils/monte_carlo.py

def shot_count_for_precision(
    std_dev: float,
    target_uncertainty: float,
    confidence: float = 0.95,
) -> int:
    """
    Determine shots needed for desired measurement precision.

    Uses standard error formula: SE = σ / √n

    Args:
        std_dev: Estimated standard deviation of measurements.
        target_uncertainty: Desired standard error.
        confidence: Confidence level (affects z-score).

    Returns:
        Recommended number of shots.
    """
    from scipy import stats

    z_score = stats.norm.ppf((1 + confidence) / 2)
    shots = int((z_score * std_dev / target_uncertainty) ** 2)
    return max(shots, 10)  # Minimum 10 shots

stratified_sampling(probabilities, shots, rng=None)

Stratified Monte Carlo sampling to reduce variance.

Parameters:

Name	Type	Description	Default
probabilities	ndarray	Probability distribution over outcomes.	required
shots	int	Number of shots.	required
rng	RandomState	Random state for reproducibility.	None

Returns:

Type	Description
Dict[int, int]	Dictionary mapping outcomes to counts.

Source code in mhrqi/utils/monte_carlo.py

def stratified_sampling(
    probabilities: np.ndarray,
    shots: int,
    rng: np.random.RandomState = None,
) -> Dict[int, int]:
    """
    Stratified Monte Carlo sampling to reduce variance.

    Args:
        probabilities: Probability distribution over outcomes.
        shots: Number of shots.
        rng: Random state for reproducibility.

    Returns:
        Dictionary mapping outcomes to counts.
    """
    if rng is None:
        rng = np.random.RandomState()

    n_outcomes = len(probabilities)
    shots_per_stratum = shots // n_outcomes
    remainder = shots % n_outcomes

    counts = {}
    for outcome_idx in range(n_outcomes):
        stratum_shots = shots_per_stratum
        if outcome_idx < remainder:
            stratum_shots += 1
        counts[outcome_idx] = stratum_shots

    # Resample within stratified buckets for variance reduction
    outcome_indices = rng.choice(
        n_outcomes,
        size=shots,
        p=probabilities / probabilities.sum()
    )

    final_counts = {}
    for idx in outcome_indices:
        final_counts[idx] = final_counts.get(idx, 0) + 1

    return final_counts

visualize_convergence(convergence_data, metric_name='KL Divergence')

Generate ASCII plot of Monte Carlo convergence.

Parameters:

Name	Type	Description	Default
convergence_data	Dict[int, float]	Dictionary mapping shot counts to error values.	required
metric_name	str	Name of metric being tracked.	'KL Divergence'

Returns:

Type	Description
str	ASCII plot string.

Source code in mhrqi/utils/monte_carlo.py

def visualize_convergence(
    convergence_data: Dict[int, float],
    metric_name: str = "KL Divergence",
) -> str:
    """
    Generate ASCII plot of Monte Carlo convergence.

    Args:
        convergence_data: Dictionary mapping shot counts to error values.
        metric_name: Name of metric being tracked.

    Returns:
        ASCII plot string.
    """
    shots = sorted(convergence_data.keys())
    errors = [convergence_data[s] for s in shots]

    max_error = max(errors)
    min_error = min(errors)

    plot = f"\n{metric_name} vs. Shot Count\n"
    plot += "=" * 50 + "\n"

    for shot_count, error in zip(shots, errors):
        normalized = (error - min_error) / (max_error - min_error + 1e-10)
        bar_length = int(normalized * 40)
        plot += f"{shot_count:6d}: {'█' * bar_length} {error:.4f}\n"

    plot += "=" * 50 + "\n"
    return plot