Simulation

The simulate module provides tools for generating synthetic M/EEG data with known ground truth source activations. This is essential for validating and comparing inverse methods under controlled conditions.

Overview

Simulation in invertmeeg supports:

• Configurable source patterns: Single dipoles, extended patches, or multiple simultaneous sources
• Noise models: White noise, colored noise, and realistic sensor noise
• Covariance structures: Correlated sources for testing robustness
• Batch generation: Efficient generation of many samples for training or benchmarking

Quick Start

from invert.simulate import SimulationConfig, SimulationGenerator

# Create a simulation configuration
config = SimulationConfig(
    n_sources=2,
    snr=5.0,
    source_extent=10.0,  # mm
)

# Generate simulations
generator = SimulationGenerator(forward, config)
evoked, stc_true = generator.generate()

API Reference

SimulationConfig

Bases: BaseModel

Configuration for EEG simulation generator.

Attributes: batch_size: Number of simulations per batch batch_repetitions: Number of times to repeat each batch n_sources: Number of active sources (int or tuple for range) n_orders: Smoothness order(s) for spatial patterns amplitude_range: Min/max amplitude for source activity n_timepoints: Number of time samples per simulation snr_range: Signal-to-noise ratio range in dB n_timecourses: Number of pre-generated timecourses beta_range: Power-law exponent range for 1/f noise add_forward_error: Whether to add perturbations to leadfield forward_error: Magnitude of forward model error inter_source_correlation: Correlation between sources (float or range) diffusion_smoothing: Whether to use diffusion-based smoothing diffusion_parameter: Smoothing strength (alpha parameter) correlation_mode: Spatial noise correlation pattern noise_color_coeff: Spatial noise correlation strength noise_temporal_beta: 1/f^beta temporal coloring of sensor noise noise_rank_deficiency: Number of projected-out sensor dimensions apply_sensor_projector: Apply projector to both signal and noise noise_low_rank_dim: Rank for low-rank spatial noise model return_noise_cov: Include per-sample noise covariance matrices in metadata estimate_noise_cov: Estimate covariance from baseline noise samples noise_cov_n_baseline: Number of baseline samples used for covariance estimate noise_cov_shrinkage: Shrinkage factor for estimated covariance random_seed: Random seed for reproducibility normalize_leadfield: Whether to normalize leadfield columns verbose: Verbosity level simulation_mode: Simulation mode ('patches' or 'mixture') background_beta: 1/f^beta exponent for smooth background background_mixture_alpha: Mixing coefficient alpha (higher = more background) source_spatial_model: Spatial source model for patch generation source_extent: Number of dipoles per source patch patch_smoothness_sigma: Gaussian smoothness (graph-hop units) within a patch patch_rank: Temporal rank per patch (1=single latent, 2=two latent components)

Source code in invert/simulate/config.py

class SimulationConfig(BaseModel):
    """Configuration for EEG simulation generator.

    Attributes:
        batch_size: Number of simulations per batch
        batch_repetitions: Number of times to repeat each batch
        n_sources: Number of active sources (int or tuple for range)
        n_orders: Smoothness order(s) for spatial patterns
        amplitude_range: Min/max amplitude for source activity
        n_timepoints: Number of time samples per simulation
        snr_range: Signal-to-noise ratio range in dB
        n_timecourses: Number of pre-generated timecourses
        beta_range: Power-law exponent range for 1/f noise
        add_forward_error: Whether to add perturbations to leadfield
        forward_error: Magnitude of forward model error
        inter_source_correlation: Correlation between sources (float or range)
        diffusion_smoothing: Whether to use diffusion-based smoothing
        diffusion_parameter: Smoothing strength (alpha parameter)
        correlation_mode: Spatial noise correlation pattern
        noise_color_coeff: Spatial noise correlation strength
        noise_temporal_beta: 1/f^beta temporal coloring of sensor noise
        noise_rank_deficiency: Number of projected-out sensor dimensions
        apply_sensor_projector: Apply projector to both signal and noise
        noise_low_rank_dim: Rank for low-rank spatial noise model
        return_noise_cov: Include per-sample noise covariance matrices in metadata
        estimate_noise_cov: Estimate covariance from baseline noise samples
        noise_cov_n_baseline: Number of baseline samples used for covariance estimate
        noise_cov_shrinkage: Shrinkage factor for estimated covariance
        random_seed: Random seed for reproducibility
        normalize_leadfield: Whether to normalize leadfield columns
        verbose: Verbosity level
        simulation_mode: Simulation mode ('patches' or 'mixture')
        background_beta: 1/f^beta exponent for smooth background
        background_mixture_alpha: Mixing coefficient alpha (higher = more background)
        source_spatial_model: Spatial source model for patch generation
        source_extent: Number of dipoles per source patch
        patch_smoothness_sigma: Gaussian smoothness (graph-hop units) within a patch
        patch_rank: Temporal rank per patch (1=single latent, 2=two latent components)
    """

    batch_size: int = Field(default=1284, ge=1)
    batch_repetitions: int = Field(default=1, ge=1)
    n_sources: int | tuple[int, int] = Field(
        default=(1, 5), description="Single value or (min, max) tuple"
    )
    n_orders: int | tuple[int, int] = Field(
        default=(0, 3), description="Smoothness order or (min, max) tuple"
    )
    amplitude_range: tuple[float, float] = Field(
        default=(0.5, 1.0), description="Source amplitude range"
    )
    n_timepoints: int = Field(default=20, ge=1)
    snr_range: tuple[float, float] = Field(
        default=(-5.0, 5.0), description="SNR range in dB"
    )
    n_timecourses: int = Field(default=5000, ge=1)
    beta_range: tuple[float, float] = Field(default=(0.0, 3.0))
    add_forward_error: bool = Field(default=False)
    forward_error: float = Field(default=0.1, ge=0.0)
    inter_source_correlation: float | tuple[float, float] = Field(default=(0.25, 0.75))
    diffusion_smoothing: bool = Field(default=True)
    diffusion_parameter: float = Field(default=0.1, ge=0.0)
    correlation_mode: (
        Literal["auto", "cholesky", "banded", "diagonal", "low_rank"] | None
    ) = Field(default=None)
    noise_color_coeff: float | tuple[float, float] = Field(default=(0.25, 0.75))
    noise_temporal_beta: float | tuple[float, float] = Field(
        default=0.0,
        description="Power-law exponent for sensor noise (0=white).",
    )
    noise_rank_deficiency: int | tuple[int, int] = Field(
        default=0,
        description="Projected-out sensor dimensions (rank loss).",
    )
    apply_sensor_projector: bool = Field(
        default=True,
        description="Apply simulated projector to signal and sensor noise.",
    )
    noise_low_rank_dim: int | tuple[int, int] = Field(
        default=4,
        description="Latent dimension for low-rank spatial noise.",
    )
    return_noise_cov: bool = Field(
        default=False,
        description="Store per-sample noise covariance matrices in metadata.",
    )
    estimate_noise_cov: bool = Field(
        default=True,
        description="Estimate noise covariance from baseline noise samples.",
    )
    noise_cov_n_baseline: int = Field(
        default=200,
        ge=8,
        description="Baseline time samples used for covariance estimation.",
    )
    noise_cov_shrinkage: float = Field(
        default=0.05,
        ge=0.0,
        le=1.0,
        description="Shrinkage factor for covariance estimation.",
    )
    random_seed: int | None = Field(default=None)
    normalize_leadfield: bool = Field(default=False)
    verbose: int = Field(default=0, ge=0)

    # Mixture mode parameters (simple background + patches)
    simulation_mode: Literal["patches", "mixture"] = Field(
        default="patches",
        description="Simulation mode: 'patches' (sparse only) or 'mixture' (smooth background + patches)",
    )
    background_beta: float | tuple[float, float] = Field(
        default=(0.5, 2.0),
        description="1/f^beta exponent for smooth background temporal dynamics",
    )
    background_mixture_alpha: float | tuple[float, float] = Field(
        default=(0.7, 0.9),
        description="Mixing coefficient alpha: y = alpha*y_background + (1-alpha)*y_patches",
    )
    source_spatial_model: Literal["diffusion_basis", "contiguous_gaussian"] = Field(
        default="diffusion_basis",
        description="Spatial source model: legacy diffusion basis or contiguous Gaussian patches.",
    )
    source_extent: int | tuple[int, int] = Field(
        default=(1, 60),
        description="Source extent as number of dipoles (single value or min/max tuple).",
    )
    patch_smoothness_sigma: float | tuple[float, float] = Field(
        default=(0.8, 2.0),
        description="Gaussian smoothness in graph-hop units for contiguous patches.",
    )
    patch_rank: int | tuple[int, int] = Field(
        default=(1, 2),
        description="Temporal rank per patch (1 or 2).",
    )

    model_config = ConfigDict(frozen=False)

SimulationGenerator

Class-based EEG simulation generator with precomputed components.

This generator creates realistic EEG simulations by: 1. Generating spatially smooth source patterns 2. Assigning colored noise timecourses to sources 3. Projecting through the leadfield matrix 4. Adding spatially/temporally colored sensor noise

The class precomputes spatial smoothing operators and timecourses during initialization for faster batch generation.

Source code in invert/simulate/simulate.py

class SimulationGenerator:
    """Class-based EEG simulation generator with precomputed components.

    This generator creates realistic EEG simulations by:
    1. Generating spatially smooth source patterns
    2. Assigning colored noise timecourses to sources
    3. Projecting through the leadfield matrix
    4. Adding spatially/temporally colored sensor noise

    The class precomputes spatial smoothing operators and timecourses
    during initialization for faster batch generation.
    """

    def __init__(self, fwd, config: SimulationConfig | None = None, **kwargs):
        """Initialize the simulation generator.

        Parameters:
            fwd: MNE forward solution object
            config: SimulationConfig instance (optional)
            **kwargs: Configuration parameters (used if config is None)
        """
        # Initialize configuration
        if config is None:
            self.config = SimulationConfig(**kwargs)
        else:
            self.config = config

        # Store forward solution components
        self.fwd = fwd
        self.channel_types = np.array(fwd["info"].get_channel_types())
        self.leadfield_original = deepcopy(fwd["sol"]["data"])
        self.leadfield = deepcopy(self.leadfield_original)
        self.n_chans, self.n_dipoles = self.leadfield.shape

        # Initialize random number generator
        self.rng = np.random.default_rng(self.config.random_seed)

        # Parse n_sources parameter
        if isinstance(self.config.n_sources, (int, float)):
            n_sources_val = np.clip(self.config.n_sources, a_min=1, a_max=np.inf)
            self.min_sources, self.max_sources = int(n_sources_val), int(n_sources_val)
        else:
            self.min_sources, self.max_sources = self.config.n_sources

        # Precompute spatial smoothing operator
        self._precompute_spatial_operators()

        # Precompute timecourses
        self._precompute_timecourses()

        # Setup correlation sampling functions
        self._setup_correlation_samplers()

        # Normalize leadfield if requested
        if self.config.normalize_leadfield:
            self.leadfield /= np.linalg.norm(self.leadfield, axis=0)

    def _precompute_spatial_operators(self):
        """Precompute graph Laplacian and smoothing operators."""
        self.adjacency = build_adjacency(self.fwd, verbose=self.config.verbose)
        adjacency_csr = self.adjacency.tocsr()
        self.neighbors = [
            adjacency_csr.indices[
                adjacency_csr.indptr[i] : adjacency_csr.indptr[i + 1]
            ].astype(int)
            for i in range(self.n_dipoles)
        ]

        # Parse n_orders parameter
        if isinstance(self.config.n_orders, (tuple, list)):
            self.min_order, self.max_order = self.config.n_orders
            if self.min_order == self.max_order:
                self.max_order += 1
        else:
            self.min_order = 0
            self.max_order = self.config.n_orders

        self.sources, self.sources_dense, self.gradient = build_spatial_basis(
            self.adjacency,
            self.n_dipoles,
            self.min_order,
            self.max_order,
            diffusion_smoothing=self.config.diffusion_smoothing,
            diffusion_parameter=self.config.diffusion_parameter,
        )
        self.n_candidates = self.sources.shape[0]

    def _precompute_timecourses(self):
        """Precompute colored noise timecourses."""
        betas = self.rng.uniform(*self.config.beta_range, self.config.n_timecourses)

        time_courses = powerlaw_noise(
            betas,
            self.config.n_timepoints,
            n_signals=self.config.n_timecourses,
            rng=self.rng,
        )

        # Normalize to max(abs()) == 1
        self.time_courses = (time_courses.T / abs(time_courses).max(axis=1)).T

    def _setup_correlation_samplers(self):
        """Setup sampling functions for simulation hyperparameters."""

        def _sampler(value, n=1, dtype=float):
            if isinstance(value, (tuple, list)):
                lo, hi = value
                if np.issubdtype(np.dtype(dtype), np.integer):
                    return self.rng.integers(int(lo), int(hi) + 1, size=n)
                return self.rng.uniform(lo, hi, n)
            return np.full(n, value, dtype=dtype)

        isc = self.config.inter_source_correlation
        self.get_inter_source_correlation = lambda n=1: _sampler(isc, n=n, dtype=float)

        ncc = self.config.noise_color_coeff
        self.get_noise_color_coeff = lambda n=1: _sampler(ncc, n=n, dtype=float)

        ntb = self.config.noise_temporal_beta
        self.get_noise_temporal_beta = lambda n=1: _sampler(ntb, n=n, dtype=float)

        nrd = self.config.noise_rank_deficiency
        self.get_noise_rank_deficiency = lambda n=1: _sampler(nrd, n=n, dtype=int)

        nlr = self.config.noise_low_rank_dim
        self.get_noise_low_rank_dim = lambda n=1: _sampler(nlr, n=n, dtype=int)

        se = self.config.source_extent
        self.get_source_extent = lambda n=1: _sampler(se, n=n, dtype=int)

        pss = self.config.patch_smoothness_sigma
        self.get_patch_smoothness_sigma = lambda n=1: _sampler(pss, n=n, dtype=float)

        pr = self.config.patch_rank
        self.get_patch_rank = lambda n=1: _sampler(pr, n=n, dtype=int)

    def _generate_smooth_background(self, batch_size):
        """Generate smooth background activity with 1/f^beta temporal dynamics.

        Uses vectorized FFT-based colored noise generation instead of
        per-dipole Python loops.

        Parameters:
            batch_size: Number of simulations to generate

        Returns:
            y_background: [batch_size, n_dipoles, n_timepoints] background activity
        """
        # Sample beta parameters for background
        if isinstance(self.config.background_beta, tuple):
            betas = self.rng.uniform(*self.config.background_beta, batch_size)
        else:
            betas = np.full(batch_size, self.config.background_beta)

        y_background_all = np.empty(
            (batch_size, self.n_dipoles, self.config.n_timepoints)
        )

        for b_idx, beta in enumerate(betas):
            # Vectorized: generate all dipole timecourses at once
            background_timecourses = powerlaw_noise(
                beta, self.config.n_timepoints, n_signals=self.n_dipoles, rng=self.rng
            )  # [n_dipoles, n_timepoints]

            # Apply spatial smoothing using gradient operator
            background_smooth = self.gradient @ background_timecourses

            # Normalize
            max_val = np.max(np.abs(background_smooth))
            if max_val > 0:
                background_smooth = background_smooth / max_val
            y_background_all[b_idx] = background_smooth

        return y_background_all

    def _setup_leadfield(self):
        """Get the leadfield matrix, optionally with forward model error."""
        if self.config.add_forward_error:
            return add_error(
                self.leadfield_original,
                self.config.forward_error,
                self.gradient,
                self.rng,
            )
        return self.leadfield

    def _generate_diffusion_basis_patches(self, batch_size):
        """Generate patch activity from the legacy diffusion basis."""
        n_sources_batch = self.rng.integers(
            self.min_sources, self.max_sources + 1, batch_size
        )

        # Select source locations
        selection = [
            self.rng.integers(0, self.n_candidates, n) for n in n_sources_batch
        ]

        # Sample amplitudes and timecourses
        amplitude_values = [
            self.rng.uniform(*self.config.amplitude_range, n) for n in n_sources_batch
        ]
        timecourse_choices = [
            self.rng.choice(self.config.n_timecourses, n) for n in n_sources_batch
        ]
        amplitudes = [self.time_courses[choice].T for choice in timecourse_choices]

        # Apply inter-source correlations
        inter_source_correlations = self.get_inter_source_correlation(n=batch_size)
        noise_color_coeffs = self.get_noise_color_coeff(n=batch_size)

        source_covariances = [
            get_cov(n, isc)
            for n, isc in zip(n_sources_batch, inter_source_correlations, strict=False)
        ]
        amplitudes = [
            amp @ np.diag(amplitude_values[i]) @ cov
            for i, (amp, cov) in enumerate(
                zip(amplitudes, source_covariances, strict=False)
            )
        ]

        y_patches = np.stack(
            [
                (amplitudes[i] @ self.sources_dense[selection[i]]).T
                / len(amplitudes[i])
                for i in range(batch_size)
            ],
            axis=0,
        )

        source_meta = {"source_vertices": [None] * batch_size, "source_ranks": None}

        return (
            y_patches,
            n_sources_batch,
            selection,
            amplitude_values,
            inter_source_correlations,
            noise_color_coeffs,
            source_meta,
        )

    def _sample_connected_vertices(self, target_size):
        """Sample one contiguous region from the source adjacency graph."""
        target_size = int(np.clip(target_size, 1, self.n_dipoles))
        if target_size == 1:
            return np.array([int(self.rng.integers(0, self.n_dipoles))], dtype=int)

        best_region = None
        for _ in range(12):
            seed = int(self.rng.integers(0, self.n_dipoles))
            region = [seed]
            region_set = {seed}
            frontier = {int(v) for v in self.neighbors[seed]}

            while len(region) < target_size and frontier:
                candidate = int(self.rng.choice(np.fromiter(frontier, dtype=int)))
                frontier.remove(candidate)
                if candidate in region_set:
                    continue
                region.append(candidate)
                region_set.add(candidate)
                for nb in self.neighbors[candidate]:
                    nb = int(nb)
                    if nb not in region_set:
                        frontier.add(nb)

            if len(region) == target_size:
                return np.asarray(region, dtype=int)

            if best_region is None or len(region) > len(best_region):
                best_region = region

        return np.asarray(best_region, dtype=int)

    def _bfs_distances_within_region(self, center, region_vertices):
        """Compute hop distances from center restricted to region vertices."""
        region_set = {int(v) for v in region_vertices}
        dists = {int(center): 0}
        queue = [int(center)]
        q_idx = 0
        while q_idx < len(queue):
            node = queue[q_idx]
            q_idx += 1
            for nb in self.neighbors[node]:
                nb = int(nb)
                if nb not in region_set or nb in dists:
                    continue
                dists[nb] = dists[node] + 1
                queue.append(nb)

        return np.array([float(dists.get(int(v), np.inf)) for v in region_vertices])

    def _make_contiguous_gaussian_patch_components(self, region_vertices, rank, sigma):
        """Create rank-1/2 smooth components inside one contiguous source region."""
        rank = int(np.clip(rank, 1, 2))
        sigma = float(max(sigma, 1e-3))
        region_vertices = np.asarray(region_vertices, dtype=int)
        if region_vertices.size == 0:
            return []

        components = []
        used_centers = set()
        for _ in range(rank):
            if len(used_centers) < region_vertices.size:
                available = [v for v in region_vertices if int(v) not in used_centers]
                center = int(self.rng.choice(available))
            else:
                center = int(self.rng.choice(region_vertices))
            used_centers.add(center)

            hop_dists = self._bfs_distances_within_region(center, region_vertices)
            weights = np.exp(-(hop_dists**2) / (2.0 * sigma**2))
            if np.max(weights) > 0:
                weights = weights / np.max(weights)

            spatial_map = np.zeros(self.n_dipoles, dtype=float)
            spatial_map[region_vertices] = weights
            components.append(spatial_map)

        return components

    def _generate_contiguous_gaussian_patches(self, batch_size):
        """Generate contiguous, irregular source patches with Gaussian smoothness.

        Inter-source correlation is controlled at the source level and patch-internal
        activity can be rank-1 or rank-2.
        """
        n_sources_batch = self.rng.integers(
            self.min_sources, self.max_sources + 1, batch_size
        )
        inter_source_correlations = self.get_inter_source_correlation(n=batch_size)
        noise_color_coeffs = self.get_noise_color_coeff(n=batch_size)

        y_batch = np.zeros((batch_size, self.n_dipoles, self.config.n_timepoints))
        selection = []
        amplitude_values = []
        source_vertices_meta = []
        source_ranks_meta = []

        for i, n_sources in enumerate(n_sources_batch):
            extents = self.get_source_extent(n=n_sources).astype(int)
            extents = np.clip(extents, 1, self.n_dipoles)

            sigmas = self.get_patch_smoothness_sigma(n=n_sources).astype(float)
            sigmas = np.clip(sigmas, 1e-3, None)

            patch_ranks = self.get_patch_rank(n=n_sources).astype(int)
            patch_ranks = np.clip(patch_ranks, 1, 2)

            amps = self.rng.uniform(*self.config.amplitude_range, n_sources)
            amplitude_values.append(amps)

            cov = get_cov(int(n_sources), float(inter_source_correlations[i]))
            base_choices_1 = self.rng.choice(self.config.n_timecourses, n_sources)
            tc_rank1 = (self.time_courses[base_choices_1].T @ cov).T

            base_choices_2 = self.rng.choice(self.config.n_timecourses, n_sources)
            tc_rank2 = (self.time_courses[base_choices_2].T @ cov).T

            centers_i = []
            source_vertices_i = []
            source_ranks_i = []

            for s in range(n_sources):
                region_vertices = self._sample_connected_vertices(int(extents[s]))
                components = self._make_contiguous_gaussian_patch_components(
                    region_vertices=region_vertices,
                    rank=int(patch_ranks[s]),
                    sigma=float(sigmas[s]),
                )

                if len(components) == 0:
                    continue

                patch_ts = components[0][:, None] * tc_rank1[s][None, :]
                if int(patch_ranks[s]) >= 2 and len(components) > 1:
                    patch_ts += components[1][:, None] * tc_rank2[s][None, :]

                y_batch[i] += (amps[s] / int(patch_ranks[s])) * patch_ts
                center_vertex = int(
                    region_vertices[np.argmax(components[0][region_vertices])]
                )
                centers_i.append(center_vertex)
                source_vertices_i.append(region_vertices)
                source_ranks_i.append(int(patch_ranks[s]))

            selection.append(np.asarray(centers_i, dtype=int))
            source_vertices_meta.append(source_vertices_i)
            source_ranks_meta.append(np.asarray(source_ranks_i, dtype=int))

        source_meta = {
            "source_vertices": source_vertices_meta,
            "source_ranks": source_ranks_meta,
        }

        return (
            y_batch,
            n_sources_batch,
            selection,
            amplitude_values,
            inter_source_correlations,
            noise_color_coeffs,
            source_meta,
        )

    def _generate_patches(self, batch_size):
        """Generate patch-based source activity.

        Returns:
            y_patches: [batch_size, n_dipoles, n_timepoints] patch activity
            selection: list of source index arrays
            amplitude_values: list of amplitude arrays
            inter_source_correlations: array of correlation values
            noise_color_coeffs: array of noise color coefficients
        """
        if self.config.source_spatial_model == "contiguous_gaussian":
            return self._generate_contiguous_gaussian_patches(batch_size)

        return self._generate_diffusion_basis_patches(batch_size)

    def _generate_background(self, batch_size, y_patches):
        """Mix background activity with patches if in mixture mode.

        Returns:
            y: [batch_size, n_dipoles, n_timepoints] combined activity
            alphas: mixing coefficients or None
        """
        if self.config.simulation_mode == "mixture":
            y_background = self._generate_smooth_background(batch_size)

            if isinstance(self.config.background_mixture_alpha, tuple):
                alphas = self.rng.uniform(
                    *self.config.background_mixture_alpha, batch_size
                )
            else:
                alphas = np.full(batch_size, self.config.background_mixture_alpha)

            alphas_bc = alphas[:, np.newaxis, np.newaxis]
            y = alphas_bc * y_background + (1 - alphas_bc) * y_patches
        else:
            y = y_patches
            alphas = None

        return y, alphas

    def _apply_noise(self, x, batch_size, noise_color_coeffs, modes_batch):
        """Apply realistic sensor noise with optional projection/rank loss.

        Returns
        -------
        x_noisy : np.ndarray
            Noisy EEG data [batch_size, n_channels, n_timepoints].
        noise_meta : dict
            Per-sample noise metadata (SNR, projector, covariance statistics).
        """
        snr_levels = self.rng.uniform(
            low=self.config.snr_range[0],
            high=self.config.snr_range[1],
            size=batch_size,
        )
        temporal_betas = self.get_noise_temporal_beta(n=batch_size)
        rank_losses = self.get_noise_rank_deficiency(n=batch_size)
        low_rank_dims = self.get_noise_low_rank_dim(n=batch_size)

        x_noisy_list = []
        realized_snrs = []
        noise_scales = []
        projectors = []
        projector_ranks = []
        covariance_model = []
        covariance_true = []
        covariance_est = []
        covariance_est_scaled = []
        covariance_rank_true = []
        covariance_rank_est = []

        for (
            xx,
            snr_level,
            corr_mode,
            noise_color_level,
            temporal_beta,
            rank_loss,
            low_rank_dim,
        ) in zip(
            x,
            snr_levels,
            modes_batch,
            noise_color_coeffs,
            temporal_betas,
            rank_losses,
            low_rank_dims,
            strict=False,
        ):
            P, _ = make_rank_projector(
                self.n_chans, rank_deficiency=int(rank_loss), rng=self.rng
            )
            if self.config.apply_sensor_projector:
                signal = P @ xx
            else:
                signal = xx

            cov_spatial = make_sensor_noise_covariance(
                self.n_chans,
                mode=corr_mode,
                noise_color_coeff=float(noise_color_level),
                rng=self.rng,
                low_rank_dim=int(low_rank_dim),
            )

            noise = sample_sensor_noise(
                cov_spatial,
                self.config.n_timepoints,
                rng=self.rng,
                temporal_beta=float(temporal_beta),
            )
            if self.config.apply_sensor_projector:
                noise = P @ noise

            noise_scaled, _scale, realized_snr = scale_noise_to_snr(
                signal, noise, float(snr_level)
            )
            x_noisy = signal + noise_scaled

            # Theoretical covariance of the realized noise (up to sampling error).
            if self.config.apply_sensor_projector:
                cov_model = (_scale**2) * (P @ cov_spatial @ P.T)
            else:
                cov_model = (_scale**2) * cov_spatial

            # True covariance from realized noise in this sample.
            cov_true = empirical_covariance(noise_scaled, shrinkage=0.0, eps=0.0)

            # Estimated covariance from independent baseline noise.
            cov_est = None
            cov_est_s = None
            if self.config.estimate_noise_cov:
                baseline = sample_sensor_noise(
                    cov_spatial,
                    int(self.config.noise_cov_n_baseline),
                    rng=self.rng,
                    temporal_beta=float(temporal_beta),
                )
                if self.config.apply_sensor_projector:
                    baseline = P @ baseline
                cov_est = empirical_covariance(
                    baseline, shrinkage=float(self.config.noise_cov_shrinkage)
                )
                # Match the baseline estimate to the scaled noise actually added to the data.
                cov_est_s = (_scale**2) * cov_est

            x_noisy_list.append(x_noisy)
            realized_snrs.append(realized_snr)
            noise_scales.append(_scale)
            projectors.append(P)
            projector_ranks.append(int(np.linalg.matrix_rank(P)))
            covariance_model.append(cov_model)
            covariance_true.append(cov_true)
            covariance_rank_true.append(int(np.linalg.matrix_rank(cov_true)))
            covariance_est.append(cov_est)
            covariance_est_scaled.append(cov_est_s)
            covariance_rank_est.append(
                int(np.linalg.matrix_rank(cov_est)) if cov_est is not None else -1
            )

        x_noisy = np.stack(x_noisy_list, axis=0)
        noise_meta = {
            "snr_target": snr_levels,
            "snr_realized": np.asarray(realized_snrs, dtype=float),
            "noise_scale": np.asarray(noise_scales, dtype=float),
            "projector": projectors,
            "projector_rank": np.asarray(projector_ranks, dtype=int),
            "noise_cov_model": covariance_model,
            "noise_cov_true": covariance_true,
            "noise_cov_est": covariance_est,
            "noise_cov_est_scaled": covariance_est_scaled,
            "noise_cov_rank_true": np.asarray(covariance_rank_true, dtype=int),
            "noise_cov_rank_est": np.asarray(covariance_rank_est, dtype=int),
            "noise_temporal_beta": np.asarray(temporal_betas, dtype=float),
            "noise_rank_deficiency": np.asarray(rank_losses, dtype=int),
            "noise_low_rank_dim": np.asarray(low_rank_dims, dtype=int),
        }
        return x_noisy, noise_meta

    def _build_metadata(
        self,
        batch_size,
        n_sources_batch,
        amplitude_values,
        inter_source_correlations,
        noise_color_coeffs,
        modes_batch,
        selection,
        source_meta,
        alphas,
        noise_meta,
    ):
        """Build simulation metadata DataFrame."""
        info_dict = {
            "n_sources": n_sources_batch,
            "amplitudes": amplitude_values,
            "snr": noise_meta["snr_target"],
            "snr_realized": noise_meta["snr_realized"],
            "noise_scale": noise_meta["noise_scale"],
            "inter_source_correlations": inter_source_correlations,
            "n_orders": [[self.min_order, self.max_order]] * batch_size,
            "diffusion_parameter": [self.config.diffusion_parameter] * batch_size,
            "n_timepoints": [self.config.n_timepoints] * batch_size,
            "n_timecourses": [self.config.n_timecourses] * batch_size,
            "correlation_mode": modes_batch,
            "noise_color_coeff": noise_color_coeffs,
            "noise_temporal_beta": noise_meta["noise_temporal_beta"],
            "noise_rank_deficiency": noise_meta["noise_rank_deficiency"],
            "projector_rank": noise_meta["projector_rank"],
            "noise_cov_rank_true": noise_meta["noise_cov_rank_true"],
            "add_forward_error": [self.config.add_forward_error] * batch_size,
            "forward_error": [self.config.forward_error] * batch_size,
            "centers": selection,
            "simulation_mode": [self.config.simulation_mode] * batch_size,
            "source_spatial_model": [self.config.source_spatial_model] * batch_size,
        }

        if source_meta.get("source_vertices", None) is not None:
            info_dict["source_vertices"] = source_meta["source_vertices"]
        if source_meta.get("source_ranks", None) is not None:
            info_dict["source_ranks"] = source_meta["source_ranks"]

        if self.config.estimate_noise_cov:
            info_dict["noise_cov_rank_est"] = noise_meta["noise_cov_rank_est"]

        if self.config.return_noise_cov:
            info_dict["projector"] = noise_meta["projector"]
            info_dict["noise_cov_model"] = noise_meta["noise_cov_model"]
            info_dict["noise_cov_true"] = noise_meta["noise_cov_true"]
            if self.config.estimate_noise_cov:
                info_dict["noise_cov_est"] = noise_meta["noise_cov_est"]
                info_dict["noise_cov_est_scaled"] = noise_meta["noise_cov_est_scaled"]

        if self.config.simulation_mode == "mixture":
            info_dict.update(
                {
                    "background_beta": [self.config.background_beta] * batch_size,
                    "background_mixture_alpha": alphas,
                }
            )

        return pd.DataFrame(info_dict)

    def generate(self):
        """Generate batches of simulations.

        Yields:
            tuple: (x, y, info) where:
                - x: EEG data [batch_size, n_channels, n_timepoints]
                - y: Source activity [batch_size, n_dipoles, n_timepoints] (scaled)
                - info: DataFrame with simulation metadata
        """
        # Setup correlation modes
        if (
            isinstance(self.config.correlation_mode, str)
            and self.config.correlation_mode.lower() == "auto"
        ):
            correlation_modes = ["cholesky", "banded", "diagonal", "low_rank", None]
            modes_batch = self.rng.choice(correlation_modes, self.config.batch_size)
        else:
            modes_batch = [self.config.correlation_mode] * self.config.batch_size

        while True:
            leadfield = self._setup_leadfield()

            (
                y_patches,
                n_sources_batch,
                selection,
                amplitude_values,
                inter_source_correlations,
                noise_color_coeffs,
                source_meta,
            ) = self._generate_patches(self.config.batch_size)

            y, alphas = self._generate_background(self.config.batch_size, y_patches)

            # Vectorized leadfield projection
            x = np.einsum("cd,bdt->bct", leadfield, y)

            x, noise_meta = self._apply_noise(
                x, self.config.batch_size, noise_color_coeffs, modes_batch
            )

            info = self._build_metadata(
                self.config.batch_size,
                n_sources_batch,
                amplitude_values,
                inter_source_correlations,
                noise_color_coeffs,
                modes_batch,
                selection,
                source_meta,
                alphas,
                noise_meta,
            )

            output = (x, y, info)

            for _ in range(self.config.batch_repetitions):
                yield output

__init__

__init__(
    fwd, config: SimulationConfig | None = None, **kwargs
)

Initialize the simulation generator.

Parameters: fwd: MNE forward solution object config: SimulationConfig instance (optional) **kwargs: Configuration parameters (used if config is None)

Source code in invert/simulate/simulate.py

def __init__(self, fwd, config: SimulationConfig | None = None, **kwargs):
    """Initialize the simulation generator.

    Parameters:
        fwd: MNE forward solution object
        config: SimulationConfig instance (optional)
        **kwargs: Configuration parameters (used if config is None)
    """
    # Initialize configuration
    if config is None:
        self.config = SimulationConfig(**kwargs)
    else:
        self.config = config

    # Store forward solution components
    self.fwd = fwd
    self.channel_types = np.array(fwd["info"].get_channel_types())
    self.leadfield_original = deepcopy(fwd["sol"]["data"])
    self.leadfield = deepcopy(self.leadfield_original)
    self.n_chans, self.n_dipoles = self.leadfield.shape

    # Initialize random number generator
    self.rng = np.random.default_rng(self.config.random_seed)

    # Parse n_sources parameter
    if isinstance(self.config.n_sources, (int, float)):
        n_sources_val = np.clip(self.config.n_sources, a_min=1, a_max=np.inf)
        self.min_sources, self.max_sources = int(n_sources_val), int(n_sources_val)
    else:
        self.min_sources, self.max_sources = self.config.n_sources

    # Precompute spatial smoothing operator
    self._precompute_spatial_operators()

    # Precompute timecourses
    self._precompute_timecourses()

    # Setup correlation sampling functions
    self._setup_correlation_samplers()

    # Normalize leadfield if requested
    if self.config.normalize_leadfield:
        self.leadfield /= np.linalg.norm(self.leadfield, axis=0)

generate

generate()

Generate batches of simulations.

Yields: tuple: (x, y, info) where: - x: EEG data [batch_size, n_channels, n_timepoints] - y: Source activity [batch_size, n_dipoles, n_timepoints] (scaled) - info: DataFrame with simulation metadata

Source code in invert/simulate/simulate.py

def generate(self):
    """Generate batches of simulations.

    Yields:
        tuple: (x, y, info) where:
            - x: EEG data [batch_size, n_channels, n_timepoints]
            - y: Source activity [batch_size, n_dipoles, n_timepoints] (scaled)
            - info: DataFrame with simulation metadata
    """
    # Setup correlation modes
    if (
        isinstance(self.config.correlation_mode, str)
        and self.config.correlation_mode.lower() == "auto"
    ):
        correlation_modes = ["cholesky", "banded", "diagonal", "low_rank", None]
        modes_batch = self.rng.choice(correlation_modes, self.config.batch_size)
    else:
        modes_batch = [self.config.correlation_mode] * self.config.batch_size

    while True:
        leadfield = self._setup_leadfield()

        (
            y_patches,
            n_sources_batch,
            selection,
            amplitude_values,
            inter_source_correlations,
            noise_color_coeffs,
            source_meta,
        ) = self._generate_patches(self.config.batch_size)

        y, alphas = self._generate_background(self.config.batch_size, y_patches)

        # Vectorized leadfield projection
        x = np.einsum("cd,bdt->bct", leadfield, y)

        x, noise_meta = self._apply_noise(
            x, self.config.batch_size, noise_color_coeffs, modes_batch
        )

        info = self._build_metadata(
            self.config.batch_size,
            n_sources_batch,
            amplitude_values,
            inter_source_correlations,
            noise_color_coeffs,
            modes_batch,
            selection,
            source_meta,
            alphas,
            noise_meta,
        )

        output = (x, y, info)

        for _ in range(self.config.batch_repetitions):
            yield output

compute_covariance

compute_covariance(Y, cov_type='basic')

Compute the covariance matrix of the data.

Parameters:

Name	Type	Description	Default
Y	ndarray	The data matrix.	required
cov_type	str	The type of covariance matrix to compute. Options are 'basic' and 'SSM'. Default is 'basic'.	'basic'

Return

C : numpy.ndarray The covariance matrix.

Source code in invert/simulate/covariance.py

def compute_covariance(Y, cov_type="basic"):
    """Compute the covariance matrix of the data.

    Parameters
    ----------
    Y : numpy.ndarray
        The data matrix.
    cov_type : str
        The type of covariance matrix to compute. Options are 'basic' and 'SSM'. Default is 'basic'.

    Return
    ------
    C : numpy.ndarray
        The covariance matrix.
    """
    if cov_type == "basic":
        C = Y @ Y.T
    elif cov_type == "SSM":
        n_time = Y.shape[1]
        M_Y = Y.T @ Y
        YY = M_Y + 0.001 * (50 / n_time) * np.trace(M_Y) * np.eye(n_time)
        P_Y = (Y @ np.linalg.inv(YY)) @ Y.T
        C = P_Y.T @ P_Y
    elif cov_type == "riemann":
        msg = "Riemannian covariance is not yet implemented as a standalone function."
        raise NotImplementedError(msg)
    else:
        msg = "Covariance type not recognized. Use 'basic', 'SSM' or provide a custom covariance matrix."
        raise ValueError(msg)

    return C

gen_correlated_sources

gen_correlated_sources(corr_coeff, T, Q)

Generate Q correlated sources with a specified correlation coefficient. The sources are generated as sinusoids with random frequencies and phases.

Parameters:

Name	Type	Description	Default
corr_coeff	float	The correlation coefficient between the sources.	required
T	int	The number of time points in the sources.	required
Q	int	The number of sources to generate.	required

Returns:

Name	Type	Description
Y	ndarray	The generated sources.

Source code in invert/simulate/covariance.py

def gen_correlated_sources(corr_coeff, T, Q):
    """Generate Q correlated sources with a specified correlation coefficient.
    The sources are generated as sinusoids with random frequencies and phases.

    Parameters
    ----------
    corr_coeff : float
        The correlation coefficient between the sources.
    T : int
        The number of time points in the sources.
    Q : int
        The number of sources to generate.

    Returns
    -------
    Y : numpy.ndarray
        The generated sources.
    """
    Cov = np.ones((Q, Q)) * corr_coeff + np.diag(
        np.ones(Q) * (1 - corr_coeff)
    )  # required covariance matrix
    freq = np.random.randint(10, 31, Q)  # random frequencies between 10Hz to 30Hz

    phases = 2 * np.pi * np.random.rand(Q)  # random phases
    t = np.linspace(10 * np.pi / T, 10 * np.pi, T)
    Signals = np.sqrt(2) * np.cos(
        2 * np.pi * freq[:, None] * t + phases[:, None]
    )  # the basic signals

    if corr_coeff < 1:
        A = np.linalg.cholesky(Cov).T  # cholesky Decomposition
        Y = A @ Signals
    else:  # Coherent Sources
        Y = np.tile(Signals[0, :], (Q, 1))

    return Y

get_cov

get_cov(n, corr_coef)

Generate a covariance matrix that is symmetric along the diagonal that correlates sources to a specified degree.

Source code in invert/simulate/covariance.py

def get_cov(n, corr_coef):
    """Generate a covariance matrix that is symmetric along the
    diagonal that correlates sources to a specified degree."""
    if corr_coef < 1:
        cov = np.ones((n, n)) * corr_coef + np.eye(n) * (1 - corr_coef)
        cov = np.linalg.cholesky(cov)
    else:
        # Make all signals be exactly the first one (perfectly coherent)
        cov = np.zeros((n, n))
        cov[:, 0] = 1
    return cov.T

add_white_noise

add_white_noise(
    X_clean,
    snr,
    rng,
    channel_types,
    noise_color_coeff=0.5,
    correlation_mode=None,
)

Parameters:

Name	Type	Description	Default
X_clean	ndarray	The clean EEG data.	required
snr	float	The signal to noise ratio in dB.	required
correlation_mode	None / str	None implies no correlation between the noise in different channels. 'banded' : Colored banded noise, where channels closer to each other will be more correlated. 'diagonal' : Some channels have varying degrees of noise. 'cholesky' : A set correlation coefficient between each pair of channels	None
noise_color_coeff	float	The magnitude of spatial coloring of the noise (not the magnitude of noise overall!).	0.5

Source code in invert/simulate/noise.py

def add_white_noise(
    X_clean, snr, rng, channel_types, noise_color_coeff=0.5, correlation_mode=None
):
    """
    Parameters
    ----------
    X_clean : numpy.ndarray
        The clean EEG data.
    snr : float
        The signal to noise ratio in dB.
    correlation_mode : None/str
        None implies no correlation between the noise in different channels.
        'banded' : Colored banded noise, where channels closer to each other will be more correlated.
        'diagonal' : Some channels have varying degrees of noise.
        'cholesky' : A set correlation coefficient between each pair of channels
    noise_color_coeff : float
        The magnitude of spatial coloring of the noise (not the magnitude of noise overall!).
    """
    n_chans, n_time = X_clean.shape
    X_noise = rng.standard_normal((n_chans, n_time))
    snr_linear = 10 ** (snr / 10)

    if isinstance(channel_types, list):
        channel_types = np.array(channel_types)
    # Ensure the channel_types array is correct length
    assert len(channel_types) == n_chans, (
        "Length of channel_types must match the number of channels in X_clean"
    )

    unique_types = np.unique(channel_types)
    X_full = np.zeros_like(X_clean)

    for ch_type in unique_types:
        type_indices = np.where(channel_types == ch_type)[0]
        X_clean_type = X_clean[type_indices, :]
        X_noise_type = X_noise[type_indices, :]
        if isinstance(noise_color_coeff, str) and isinstance(
            correlation_mode, np.ndarray
        ):
            # Real Noise Covariance
            X_noise_type = (
                np.linalg.cholesky(correlation_mode[type_indices][:, type_indices])
                @ X_noise_type
            )
        elif correlation_mode == "cholesky":
            covariance_matrix = np.full(
                (len(type_indices), len(type_indices)), noise_color_coeff
            )
            np.fill_diagonal(covariance_matrix, 1)  # Set diagonal to 1 for variance

            # Cholesky decomposition
            X_noise_type = np.linalg.cholesky(covariance_matrix) @ X_noise_type
        elif correlation_mode == "banded":
            num_sensors = X_noise_type.shape[0]
            Y = np.zeros_like(X_noise_type)
            for i in range(num_sensors):
                Y[i, :] = X_noise_type[i, :]
                for j in range(num_sensors):
                    if abs(i - j) % num_sensors == 1:
                        Y[i, :] += (noise_color_coeff / np.sqrt(2)) * X_noise_type[j, :]
            X_noise_type = Y
        elif correlation_mode == "diagonal":
            X_noise_type[1::3, :] *= 1 - noise_color_coeff
            X_noise_type[2::3, :] *= 1 + noise_color_coeff
        elif correlation_mode is None:
            pass
        else:
            msg = f"correlation_mode can be either None, cholesky, banded or diagonal, but was {correlation_mode}"
            raise AttributeError(msg)

        rms_noise = rms(X_noise_type)
        rms_signal = rms(X_clean_type)
        scaler = rms_signal / (snr_linear * rms_noise)

        X_full[type_indices] = X_clean_type + X_noise_type * scaler

    return X_full

empirical_covariance

empirical_covariance(noise, shrinkage=0.0, eps=1e-12)

Estimate covariance from noise samples with optional Ledoit-style shrinkage.

Source code in invert/simulate/noise.py

def empirical_covariance(noise, shrinkage=0.0, eps=1e-12):
    """Estimate covariance from noise samples with optional Ledoit-style shrinkage."""
    noise = np.asarray(noise, dtype=float)
    n_chans, n_time = noise.shape
    demeaned = noise - noise.mean(axis=1, keepdims=True)
    cov = (demeaned @ demeaned.T) / max(n_time - 1, 1)
    cov = 0.5 * (cov + cov.T)

    gamma = float(np.clip(shrinkage, 0.0, 1.0))
    if gamma > 0:
        target = np.trace(cov) / n_chans
        cov = (1.0 - gamma) * cov + gamma * target * np.eye(n_chans)

    if eps > 0:
        cov = cov + eps * np.eye(n_chans)
    return cov

make_rank_projector

make_rank_projector(n_chans, rank_deficiency=0, rng=None)

Create an orthogonal projector P = I - U U^T with controlled rank loss.

Source code in invert/simulate/noise.py

def make_rank_projector(n_chans, rank_deficiency=0, rng=None):
    """Create an orthogonal projector P = I - U U^T with controlled rank loss."""
    if rng is None:
        rng = np.random.default_rng()

    k = int(np.clip(rank_deficiency, 0, n_chans - 1))
    if k == 0:
        return np.eye(n_chans), np.zeros((n_chans, 0))

    basis = rng.standard_normal((n_chans, k))
    q, _ = np.linalg.qr(basis)
    U = q[:, :k]
    P = np.eye(n_chans) - U @ U.T
    P = 0.5 * (P + P.T)
    return P, U

make_sensor_noise_covariance

make_sensor_noise_covariance(
    n_chans,
    mode=None,
    noise_color_coeff=0.5,
    rng=None,
    low_rank_dim=4,
    eps=1e-12,
)

Create a PSD spatial covariance matrix for sensor noise.

Parameters:

Name	Type	Description	Default
n_chans	int	Number of channels.	required
mode	None | str	One of {None, "cholesky", "banded", "diagonal", "low_rank"}.	None
noise_color_coeff	float	Strength of spatial correlation/coloring.	0.5
rng	Generator | None	Random number generator.	None
low_rank_dim	int	Latent rank used when mode="low_rank".	4
eps	float	Diagonal jitter for numerical safety.	1e-12

Source code in invert/simulate/noise.py

def make_sensor_noise_covariance(
    n_chans,
    mode=None,
    noise_color_coeff=0.5,
    rng=None,
    low_rank_dim=4,
    eps=1e-12,
):
    """Create a PSD spatial covariance matrix for sensor noise.

    Parameters
    ----------
    n_chans : int
        Number of channels.
    mode : None | str
        One of {None, "cholesky", "banded", "diagonal", "low_rank"}.
    noise_color_coeff : float
        Strength of spatial correlation/coloring.
    rng : numpy.random.Generator | None
        Random number generator.
    low_rank_dim : int
        Latent rank used when ``mode="low_rank"``.
    eps : float
        Diagonal jitter for numerical safety.
    """
    if rng is None:
        rng = np.random.default_rng()

    coeff = float(np.clip(noise_color_coeff, 0.0, 0.999))

    if mode is None:
        cov = np.eye(n_chans)
    elif mode == "cholesky":
        cov = np.full((n_chans, n_chans), coeff)
        np.fill_diagonal(cov, 1.0)
    elif mode == "banded":
        idx = np.arange(n_chans)
        dist = np.abs(idx[:, None] - idx[None, :])
        cov = coeff**dist
        np.fill_diagonal(cov, 1.0)
    elif mode == "diagonal":
        variances = np.ones(n_chans)
        variances[1::3] *= max(1.0 - coeff, eps)
        variances[2::3] *= 1.0 + coeff
        cov = np.diag(variances)
    elif mode == "low_rank":
        latent_rank = int(np.clip(low_rank_dim, 1, max(1, n_chans - 1)))
        basis = rng.standard_normal((n_chans, latent_rank))
        low_rank = (basis @ basis.T) / max(latent_rank, 1)
        low_rank = low_rank / max(np.trace(low_rank) / n_chans, eps)
        cov = (1.0 - coeff) * np.eye(n_chans) + coeff * low_rank
    else:
        msg = (
            "mode must be one of None, 'cholesky', 'banded', 'diagonal', "
            f"'low_rank', but got {mode!r}"
        )
        raise ValueError(msg)

    cov = 0.5 * (cov + cov.T)
    # Normalize to unit average variance to keep coefficient effects interpretable.
    cov = cov / max(np.trace(cov) / n_chans, eps)
    cov = cov + eps * np.eye(n_chans)
    return cov

powerlaw_noise

powerlaw_noise(beta, n_timepoints, n_signals=1, rng=None)

Generate 1/f^beta colored noise via FFT spectral shaping.

Parameters:

Name	Type	Description	Default
beta	float or array - like	Power-law exponent(s). 0=white, 1=pink, 2=brown. If array, must have length n_signals.	required
n_timepoints	int	Number of time samples.	required
n_signals	int	Number of independent signals to generate.	1
rng	Generator or None	Random number generator.	None

Returns:

Name	Type	Description
signals	(ndarray, shape(n_signals, n_timepoints))	Colored noise signals.

Source code in invert/simulate/noise.py

def powerlaw_noise(beta, n_timepoints, n_signals=1, rng=None):
    """Generate 1/f^beta colored noise via FFT spectral shaping.

    Parameters
    ----------
    beta : float or array-like
        Power-law exponent(s). 0=white, 1=pink, 2=brown.
        If array, must have length n_signals.
    n_timepoints : int
        Number of time samples.
    n_signals : int
        Number of independent signals to generate.
    rng : numpy.random.Generator or None
        Random number generator.

    Returns
    -------
    signals : ndarray, shape (n_signals, n_timepoints)
        Colored noise signals.
    """
    if rng is None:
        rng = np.random.default_rng()

    beta = np.atleast_1d(np.asarray(beta, dtype=float))
    if beta.shape[0] == 1:
        beta = np.broadcast_to(beta, (n_signals,))

    # Generate white noise and FFT
    white = rng.standard_normal((n_signals, n_timepoints))
    spectrum = np.fft.rfft(white)

    # Build frequency-domain filter: 1/f^(beta/2) (power spectrum goes as 1/f^beta)
    freqs = np.fft.rfftfreq(n_timepoints)
    # Skip DC (index 0) to avoid division by zero
    freq_filter = np.ones((n_signals, len(freqs)))
    freq_filter[:, 1:] = freqs[1:][np.newaxis, :] ** (-beta[:, np.newaxis] / 2.0)

    spectrum *= freq_filter
    signals = np.fft.irfft(spectrum, n=n_timepoints)

    return signals

sample_sensor_noise

sample_sensor_noise(
    covariance,
    n_timepoints,
    rng=None,
    temporal_beta=0.0,
    eps=1e-12,
)

Sample sensor noise with desired spatial covariance and temporal color.

Source code in invert/simulate/noise.py

def sample_sensor_noise(
    covariance,
    n_timepoints,
    rng=None,
    temporal_beta=0.0,
    eps=1e-12,
):
    """Sample sensor noise with desired spatial covariance and temporal color."""
    if rng is None:
        rng = np.random.default_rng()

    covariance = np.asarray(covariance, dtype=float)
    n_chans = covariance.shape[0]

    if np.isscalar(temporal_beta):
        beta = np.full(n_chans, float(temporal_beta))
    else:
        beta = np.asarray(temporal_beta, dtype=float)
        if beta.shape[0] != n_chans:
            msg = f"temporal_beta length {beta.shape[0]} must match n_chans={n_chans}"
            raise ValueError(msg)

    if np.allclose(beta, 0.0):
        base = rng.standard_normal((n_chans, n_timepoints))
    else:
        base = powerlaw_noise(beta, n_timepoints, n_signals=n_chans, rng=rng)

    # Normalize each channel before spatial mixing.
    std = np.std(base, axis=1, keepdims=True)
    base = base / np.maximum(std, eps)

    # Numerically safe covariance factorization.
    covariance = 0.5 * (covariance + covariance.T)
    evals, evecs = np.linalg.eigh(covariance)
    evals = np.maximum(evals, eps)
    sqrt_cov = evecs @ np.diag(np.sqrt(evals)) @ evecs.T
    noise = sqrt_cov @ base
    return noise

scale_noise_to_snr

scale_noise_to_snr(signal, noise, snr_db, eps=1e-12)

Scale noise to match target SNR (in dB) for a given signal.

Source code in invert/simulate/noise.py

def scale_noise_to_snr(signal, noise, snr_db, eps=1e-12):
    """Scale noise to match target SNR (in dB) for a given signal."""
    signal_power = float(np.mean(np.asarray(signal, dtype=float) ** 2))
    noise_power = float(np.mean(np.asarray(noise, dtype=float) ** 2))
    snr_linear = float(10 ** (float(snr_db) / 10.0))

    scale = np.sqrt(signal_power / max(snr_linear * noise_power, eps))
    noise_scaled = noise * scale
    realized_snr = 10.0 * np.log10(
        signal_power / max(float(np.mean(noise_scaled**2)), eps)
    )
    return noise_scaled, scale, realized_snr

build_adjacency

build_adjacency(forward, verbose=0)

Build sparse adjacency matrix from an MNE forward solution.

Parameters:

Name	Type	Description	Default
forward	dict	MNE forward solution object.	required
verbose	int	Verbosity level.	0

Returns:

Name	Type	Description
adjacency	csr_matrix	Sparse adjacency matrix.

Source code in invert/simulate/spatial.py

def build_adjacency(forward, verbose=0):
    """Build sparse adjacency matrix from an MNE forward solution.

    Parameters
    ----------
    forward : dict
        MNE forward solution object.
    verbose : int
        Verbosity level.

    Returns
    -------
    adjacency : csr_matrix
        Sparse adjacency matrix.
    """
    return build_source_adjacency(forward["src"], verbose=verbose)

build_spatial_basis

build_spatial_basis(
    adjacency,
    n_dipoles,
    min_order,
    max_order,
    diffusion_smoothing=True,
    diffusion_parameter=0.1,
)

Build multi-order spatial basis from adjacency matrix.

Parameters:

Name	Type	Description	Default
adjacency	csr_matrix	Sparse adjacency matrix.	required
n_dipoles	int	Number of dipoles (source space size).	required
min_order	int	Minimum smoothing order to include.	required
max_order	int	Maximum smoothing order (exclusive).	required
diffusion_smoothing	bool	Whether to use diffusion smoothing (True) or absolute Laplacian (False).	True
diffusion_parameter	float	Diffusion parameter alpha for smoothing.	0.1

Returns:

Name	Type	Description
sources	csr_matrix	Stacked spatial basis (sparse).
sources_dense	ndarray	Dense version for fast indexing.
gradient	csr_matrix	The gradient/smoothing operator.

Source code in invert/simulate/spatial.py

def build_spatial_basis(
    adjacency,
    n_dipoles,
    min_order,
    max_order,
    diffusion_smoothing=True,
    diffusion_parameter=0.1,
):
    """Build multi-order spatial basis from adjacency matrix.

    Parameters
    ----------
    adjacency : csr_matrix
        Sparse adjacency matrix.
    n_dipoles : int
        Number of dipoles (source space size).
    min_order : int
        Minimum smoothing order to include.
    max_order : int
        Maximum smoothing order (exclusive).
    diffusion_smoothing : bool
        Whether to use diffusion smoothing (True) or absolute Laplacian (False).
    diffusion_parameter : float
        Diffusion parameter alpha for smoothing.

    Returns
    -------
    sources : csr_matrix
        Stacked spatial basis (sparse).
    sources_dense : ndarray
        Dense version for fast indexing.
    gradient : csr_matrix
        The gradient/smoothing operator.
    """
    if diffusion_smoothing:
        gradient = np.identity(n_dipoles) - diffusion_parameter * laplacian(adjacency)
    else:
        gradient = abs(laplacian(adjacency))

    gradient = csr_matrix(gradient)

    # Build multi-order source basis using iterative sparse multiplication
    last_block = csr_matrix(np.identity(n_dipoles))
    blocks = [last_block]

    for _i in range(1, max_order):
        last_block = last_block @ gradient

        # Normalize each column (with guard against zero max)
        row_maxes = last_block.max(axis=0).toarray().flatten()
        row_maxes[row_maxes == 0] = 1.0
        last_block = last_block / row_maxes[np.newaxis]

        blocks.append(last_block)

    sources = vstack(blocks[min_order:])
    sources_dense = sources.toarray()

    return sources, sources_dense, gradient

SimulationConfig

invert.simulate.SimulationConfig

Bases: BaseModel

Configuration for EEG simulation generator.

Attributes: batch_size: Number of simulations per batch batch_repetitions: Number of times to repeat each batch n_sources: Number of active sources (int or tuple for range) n_orders: Smoothness order(s) for spatial patterns amplitude_range: Min/max amplitude for source activity n_timepoints: Number of time samples per simulation snr_range: Signal-to-noise ratio range in dB n_timecourses: Number of pre-generated timecourses beta_range: Power-law exponent range for 1/f noise add_forward_error: Whether to add perturbations to leadfield forward_error: Magnitude of forward model error inter_source_correlation: Correlation between sources (float or range) diffusion_smoothing: Whether to use diffusion-based smoothing diffusion_parameter: Smoothing strength (alpha parameter) correlation_mode: Spatial noise correlation pattern noise_color_coeff: Spatial noise correlation strength noise_temporal_beta: 1/f^beta temporal coloring of sensor noise noise_rank_deficiency: Number of projected-out sensor dimensions apply_sensor_projector: Apply projector to both signal and noise noise_low_rank_dim: Rank for low-rank spatial noise model return_noise_cov: Include per-sample noise covariance matrices in metadata estimate_noise_cov: Estimate covariance from baseline noise samples noise_cov_n_baseline: Number of baseline samples used for covariance estimate noise_cov_shrinkage: Shrinkage factor for estimated covariance random_seed: Random seed for reproducibility normalize_leadfield: Whether to normalize leadfield columns verbose: Verbosity level simulation_mode: Simulation mode ('patches' or 'mixture') background_beta: 1/f^beta exponent for smooth background background_mixture_alpha: Mixing coefficient alpha (higher = more background) source_spatial_model: Spatial source model for patch generation source_extent: Number of dipoles per source patch patch_smoothness_sigma: Gaussian smoothness (graph-hop units) within a patch patch_rank: Temporal rank per patch (1=single latent, 2=two latent components)

Source code in invert/simulate/config.py

class SimulationConfig(BaseModel):
    """Configuration for EEG simulation generator.

    Attributes:
        batch_size: Number of simulations per batch
        batch_repetitions: Number of times to repeat each batch
        n_sources: Number of active sources (int or tuple for range)
        n_orders: Smoothness order(s) for spatial patterns
        amplitude_range: Min/max amplitude for source activity
        n_timepoints: Number of time samples per simulation
        snr_range: Signal-to-noise ratio range in dB
        n_timecourses: Number of pre-generated timecourses
        beta_range: Power-law exponent range for 1/f noise
        add_forward_error: Whether to add perturbations to leadfield
        forward_error: Magnitude of forward model error
        inter_source_correlation: Correlation between sources (float or range)
        diffusion_smoothing: Whether to use diffusion-based smoothing
        diffusion_parameter: Smoothing strength (alpha parameter)
        correlation_mode: Spatial noise correlation pattern
        noise_color_coeff: Spatial noise correlation strength
        noise_temporal_beta: 1/f^beta temporal coloring of sensor noise
        noise_rank_deficiency: Number of projected-out sensor dimensions
        apply_sensor_projector: Apply projector to both signal and noise
        noise_low_rank_dim: Rank for low-rank spatial noise model
        return_noise_cov: Include per-sample noise covariance matrices in metadata
        estimate_noise_cov: Estimate covariance from baseline noise samples
        noise_cov_n_baseline: Number of baseline samples used for covariance estimate
        noise_cov_shrinkage: Shrinkage factor for estimated covariance
        random_seed: Random seed for reproducibility
        normalize_leadfield: Whether to normalize leadfield columns
        verbose: Verbosity level
        simulation_mode: Simulation mode ('patches' or 'mixture')
        background_beta: 1/f^beta exponent for smooth background
        background_mixture_alpha: Mixing coefficient alpha (higher = more background)
        source_spatial_model: Spatial source model for patch generation
        source_extent: Number of dipoles per source patch
        patch_smoothness_sigma: Gaussian smoothness (graph-hop units) within a patch
        patch_rank: Temporal rank per patch (1=single latent, 2=two latent components)
    """

    batch_size: int = Field(default=1284, ge=1)
    batch_repetitions: int = Field(default=1, ge=1)
    n_sources: int | tuple[int, int] = Field(
        default=(1, 5), description="Single value or (min, max) tuple"
    )
    n_orders: int | tuple[int, int] = Field(
        default=(0, 3), description="Smoothness order or (min, max) tuple"
    )
    amplitude_range: tuple[float, float] = Field(
        default=(0.5, 1.0), description="Source amplitude range"
    )
    n_timepoints: int = Field(default=20, ge=1)
    snr_range: tuple[float, float] = Field(
        default=(-5.0, 5.0), description="SNR range in dB"
    )
    n_timecourses: int = Field(default=5000, ge=1)
    beta_range: tuple[float, float] = Field(default=(0.0, 3.0))
    add_forward_error: bool = Field(default=False)
    forward_error: float = Field(default=0.1, ge=0.0)
    inter_source_correlation: float | tuple[float, float] = Field(default=(0.25, 0.75))
    diffusion_smoothing: bool = Field(default=True)
    diffusion_parameter: float = Field(default=0.1, ge=0.0)
    correlation_mode: (
        Literal["auto", "cholesky", "banded", "diagonal", "low_rank"] | None
    ) = Field(default=None)
    noise_color_coeff: float | tuple[float, float] = Field(default=(0.25, 0.75))
    noise_temporal_beta: float | tuple[float, float] = Field(
        default=0.0,
        description="Power-law exponent for sensor noise (0=white).",
    )
    noise_rank_deficiency: int | tuple[int, int] = Field(
        default=0,
        description="Projected-out sensor dimensions (rank loss).",
    )
    apply_sensor_projector: bool = Field(
        default=True,
        description="Apply simulated projector to signal and sensor noise.",
    )
    noise_low_rank_dim: int | tuple[int, int] = Field(
        default=4,
        description="Latent dimension for low-rank spatial noise.",
    )
    return_noise_cov: bool = Field(
        default=False,
        description="Store per-sample noise covariance matrices in metadata.",
    )
    estimate_noise_cov: bool = Field(
        default=True,
        description="Estimate noise covariance from baseline noise samples.",
    )
    noise_cov_n_baseline: int = Field(
        default=200,
        ge=8,
        description="Baseline time samples used for covariance estimation.",
    )
    noise_cov_shrinkage: float = Field(
        default=0.05,
        ge=0.0,
        le=1.0,
        description="Shrinkage factor for covariance estimation.",
    )
    random_seed: int | None = Field(default=None)
    normalize_leadfield: bool = Field(default=False)
    verbose: int = Field(default=0, ge=0)

    # Mixture mode parameters (simple background + patches)
    simulation_mode: Literal["patches", "mixture"] = Field(
        default="patches",
        description="Simulation mode: 'patches' (sparse only) or 'mixture' (smooth background + patches)",
    )
    background_beta: float | tuple[float, float] = Field(
        default=(0.5, 2.0),
        description="1/f^beta exponent for smooth background temporal dynamics",
    )
    background_mixture_alpha: float | tuple[float, float] = Field(
        default=(0.7, 0.9),
        description="Mixing coefficient alpha: y = alpha*y_background + (1-alpha)*y_patches",
    )
    source_spatial_model: Literal["diffusion_basis", "contiguous_gaussian"] = Field(
        default="diffusion_basis",
        description="Spatial source model: legacy diffusion basis or contiguous Gaussian patches.",
    )
    source_extent: int | tuple[int, int] = Field(
        default=(1, 60),
        description="Source extent as number of dipoles (single value or min/max tuple).",
    )
    patch_smoothness_sigma: float | tuple[float, float] = Field(
        default=(0.8, 2.0),
        description="Gaussian smoothness in graph-hop units for contiguous patches.",
    )
    patch_rank: int | tuple[int, int] = Field(
        default=(1, 2),
        description="Temporal rank per patch (1 or 2).",
    )

    model_config = ConfigDict(frozen=False)

SimulationGenerator

invert.simulate.SimulationGenerator

Class-based EEG simulation generator with precomputed components.

This generator creates realistic EEG simulations by: 1. Generating spatially smooth source patterns 2. Assigning colored noise timecourses to sources 3. Projecting through the leadfield matrix 4. Adding spatially/temporally colored sensor noise

The class precomputes spatial smoothing operators and timecourses during initialization for faster batch generation.

Source code in invert/simulate/simulate.py

class SimulationGenerator:
    """Class-based EEG simulation generator with precomputed components.

    This generator creates realistic EEG simulations by:
    1. Generating spatially smooth source patterns
    2. Assigning colored noise timecourses to sources
    3. Projecting through the leadfield matrix
    4. Adding spatially/temporally colored sensor noise

    The class precomputes spatial smoothing operators and timecourses
    during initialization for faster batch generation.
    """

    def __init__(self, fwd, config: SimulationConfig | None = None, **kwargs):
        """Initialize the simulation generator.

        Parameters:
            fwd: MNE forward solution object
            config: SimulationConfig instance (optional)
            **kwargs: Configuration parameters (used if config is None)
        """
        # Initialize configuration
        if config is None:
            self.config = SimulationConfig(**kwargs)
        else:
            self.config = config

        # Store forward solution components
        self.fwd = fwd
        self.channel_types = np.array(fwd["info"].get_channel_types())
        self.leadfield_original = deepcopy(fwd["sol"]["data"])
        self.leadfield = deepcopy(self.leadfield_original)
        self.n_chans, self.n_dipoles = self.leadfield.shape

        # Initialize random number generator
        self.rng = np.random.default_rng(self.config.random_seed)

        # Parse n_sources parameter
        if isinstance(self.config.n_sources, (int, float)):
            n_sources_val = np.clip(self.config.n_sources, a_min=1, a_max=np.inf)
            self.min_sources, self.max_sources = int(n_sources_val), int(n_sources_val)
        else:
            self.min_sources, self.max_sources = self.config.n_sources

        # Precompute spatial smoothing operator
        self._precompute_spatial_operators()

        # Precompute timecourses
        self._precompute_timecourses()

        # Setup correlation sampling functions
        self._setup_correlation_samplers()

        # Normalize leadfield if requested
        if self.config.normalize_leadfield:
            self.leadfield /= np.linalg.norm(self.leadfield, axis=0)

    def _precompute_spatial_operators(self):
        """Precompute graph Laplacian and smoothing operators."""
        self.adjacency = build_adjacency(self.fwd, verbose=self.config.verbose)
        adjacency_csr = self.adjacency.tocsr()
        self.neighbors = [
            adjacency_csr.indices[
                adjacency_csr.indptr[i] : adjacency_csr.indptr[i + 1]
            ].astype(int)
            for i in range(self.n_dipoles)
        ]

        # Parse n_orders parameter
        if isinstance(self.config.n_orders, (tuple, list)):
            self.min_order, self.max_order = self.config.n_orders
            if self.min_order == self.max_order:
                self.max_order += 1
        else:
            self.min_order = 0
            self.max_order = self.config.n_orders

        self.sources, self.sources_dense, self.gradient = build_spatial_basis(
            self.adjacency,
            self.n_dipoles,
            self.min_order,
            self.max_order,
            diffusion_smoothing=self.config.diffusion_smoothing,
            diffusion_parameter=self.config.diffusion_parameter,
        )
        self.n_candidates = self.sources.shape[0]

    def _precompute_timecourses(self):
        """Precompute colored noise timecourses."""
        betas = self.rng.uniform(*self.config.beta_range, self.config.n_timecourses)

        time_courses = powerlaw_noise(
            betas,
            self.config.n_timepoints,
            n_signals=self.config.n_timecourses,
            rng=self.rng,
        )

        # Normalize to max(abs()) == 1
        self.time_courses = (time_courses.T / abs(time_courses).max(axis=1)).T

    def _setup_correlation_samplers(self):
        """Setup sampling functions for simulation hyperparameters."""

        def _sampler(value, n=1, dtype=float):
            if isinstance(value, (tuple, list)):
                lo, hi = value
                if np.issubdtype(np.dtype(dtype), np.integer):
                    return self.rng.integers(int(lo), int(hi) + 1, size=n)
                return self.rng.uniform(lo, hi, n)
            return np.full(n, value, dtype=dtype)

        isc = self.config.inter_source_correlation
        self.get_inter_source_correlation = lambda n=1: _sampler(isc, n=n, dtype=float)

        ncc = self.config.noise_color_coeff
        self.get_noise_color_coeff = lambda n=1: _sampler(ncc, n=n, dtype=float)

        ntb = self.config.noise_temporal_beta
        self.get_noise_temporal_beta = lambda n=1: _sampler(ntb, n=n, dtype=float)

        nrd = self.config.noise_rank_deficiency
        self.get_noise_rank_deficiency = lambda n=1: _sampler(nrd, n=n, dtype=int)

        nlr = self.config.noise_low_rank_dim
        self.get_noise_low_rank_dim = lambda n=1: _sampler(nlr, n=n, dtype=int)

        se = self.config.source_extent
        self.get_source_extent = lambda n=1: _sampler(se, n=n, dtype=int)

        pss = self.config.patch_smoothness_sigma
        self.get_patch_smoothness_sigma = lambda n=1: _sampler(pss, n=n, dtype=float)

        pr = self.config.patch_rank
        self.get_patch_rank = lambda n=1: _sampler(pr, n=n, dtype=int)

    def _generate_smooth_background(self, batch_size):
        """Generate smooth background activity with 1/f^beta temporal dynamics.

        Uses vectorized FFT-based colored noise generation instead of
        per-dipole Python loops.

        Parameters:
            batch_size: Number of simulations to generate

        Returns:
            y_background: [batch_size, n_dipoles, n_timepoints] background activity
        """
        # Sample beta parameters for background
        if isinstance(self.config.background_beta, tuple):
            betas = self.rng.uniform(*self.config.background_beta, batch_size)
        else:
            betas = np.full(batch_size, self.config.background_beta)

        y_background_all = np.empty(
            (batch_size, self.n_dipoles, self.config.n_timepoints)
        )

        for b_idx, beta in enumerate(betas):
            # Vectorized: generate all dipole timecourses at once
            background_timecourses = powerlaw_noise(
                beta, self.config.n_timepoints, n_signals=self.n_dipoles, rng=self.rng
            )  # [n_dipoles, n_timepoints]

            # Apply spatial smoothing using gradient operator
            background_smooth = self.gradient @ background_timecourses

            # Normalize
            max_val = np.max(np.abs(background_smooth))
            if max_val > 0:
                background_smooth = background_smooth / max_val
            y_background_all[b_idx] = background_smooth

        return y_background_all

    def _setup_leadfield(self):
        """Get the leadfield matrix, optionally with forward model error."""
        if self.config.add_forward_error:
            return add_error(
                self.leadfield_original,
                self.config.forward_error,
                self.gradient,
                self.rng,
            )
        return self.leadfield

    def _generate_diffusion_basis_patches(self, batch_size):
        """Generate patch activity from the legacy diffusion basis."""
        n_sources_batch = self.rng.integers(
            self.min_sources, self.max_sources + 1, batch_size
        )

        # Select source locations
        selection = [
            self.rng.integers(0, self.n_candidates, n) for n in n_sources_batch
        ]

        # Sample amplitudes and timecourses
        amplitude_values = [
            self.rng.uniform(*self.config.amplitude_range, n) for n in n_sources_batch
        ]
        timecourse_choices = [
            self.rng.choice(self.config.n_timecourses, n) for n in n_sources_batch
        ]
        amplitudes = [self.time_courses[choice].T for choice in timecourse_choices]

        # Apply inter-source correlations
        inter_source_correlations = self.get_inter_source_correlation(n=batch_size)
        noise_color_coeffs = self.get_noise_color_coeff(n=batch_size)

        source_covariances = [
            get_cov(n, isc)
            for n, isc in zip(n_sources_batch, inter_source_correlations, strict=False)
        ]
        amplitudes = [
            amp @ np.diag(amplitude_values[i]) @ cov
            for i, (amp, cov) in enumerate(
                zip(amplitudes, source_covariances, strict=False)
            )
        ]

        y_patches = np.stack(
            [
                (amplitudes[i] @ self.sources_dense[selection[i]]).T
                / len(amplitudes[i])
                for i in range(batch_size)
            ],
            axis=0,
        )

        source_meta = {"source_vertices": [None] * batch_size, "source_ranks": None}

        return (
            y_patches,
            n_sources_batch,
            selection,
            amplitude_values,
            inter_source_correlations,
            noise_color_coeffs,
            source_meta,
        )

    def _sample_connected_vertices(self, target_size):
        """Sample one contiguous region from the source adjacency graph."""
        target_size = int(np.clip(target_size, 1, self.n_dipoles))
        if target_size == 1:
            return np.array([int(self.rng.integers(0, self.n_dipoles))], dtype=int)

        best_region = None
        for _ in range(12):
            seed = int(self.rng.integers(0, self.n_dipoles))
            region = [seed]
            region_set = {seed}
            frontier = {int(v) for v in self.neighbors[seed]}

            while len(region) < target_size and frontier:
                candidate = int(self.rng.choice(np.fromiter(frontier, dtype=int)))
                frontier.remove(candidate)
                if candidate in region_set:
                    continue
                region.append(candidate)
                region_set.add(candidate)
                for nb in self.neighbors[candidate]:
                    nb = int(nb)
                    if nb not in region_set:
                        frontier.add(nb)

            if len(region) == target_size:
                return np.asarray(region, dtype=int)

            if best_region is None or len(region) > len(best_region):
                best_region = region

        return np.asarray(best_region, dtype=int)

    def _bfs_distances_within_region(self, center, region_vertices):
        """Compute hop distances from center restricted to region vertices."""
        region_set = {int(v) for v in region_vertices}
        dists = {int(center): 0}
        queue = [int(center)]
        q_idx = 0
        while q_idx < len(queue):
            node = queue[q_idx]
            q_idx += 1
            for nb in self.neighbors[node]:
                nb = int(nb)
                if nb not in region_set or nb in dists:
                    continue
                dists[nb] = dists[node] + 1
                queue.append(nb)

        return np.array([float(dists.get(int(v), np.inf)) for v in region_vertices])

    def _make_contiguous_gaussian_patch_components(self, region_vertices, rank, sigma):
        """Create rank-1/2 smooth components inside one contiguous source region."""
        rank = int(np.clip(rank, 1, 2))
        sigma = float(max(sigma, 1e-3))
        region_vertices = np.asarray(region_vertices, dtype=int)
        if region_vertices.size == 0:
            return []

        components = []
        used_centers = set()
        for _ in range(rank):
            if len(used_centers) < region_vertices.size:
                available = [v for v in region_vertices if int(v) not in used_centers]
                center = int(self.rng.choice(available))
            else:
                center = int(self.rng.choice(region_vertices))
            used_centers.add(center)

            hop_dists = self._bfs_distances_within_region(center, region_vertices)
            weights = np.exp(-(hop_dists**2) / (2.0 * sigma**2))
            if np.max(weights) > 0:
                weights = weights / np.max(weights)

            spatial_map = np.zeros(self.n_dipoles, dtype=float)
            spatial_map[region_vertices] = weights
            components.append(spatial_map)

        return components

    def _generate_contiguous_gaussian_patches(self, batch_size):
        """Generate contiguous, irregular source patches with Gaussian smoothness.

        Inter-source correlation is controlled at the source level and patch-internal
        activity can be rank-1 or rank-2.
        """
        n_sources_batch = self.rng.integers(
            self.min_sources, self.max_sources + 1, batch_size
        )
        inter_source_correlations = self.get_inter_source_correlation(n=batch_size)
        noise_color_coeffs = self.get_noise_color_coeff(n=batch_size)

        y_batch = np.zeros((batch_size, self.n_dipoles, self.config.n_timepoints))
        selection = []
        amplitude_values = []
        source_vertices_meta = []
        source_ranks_meta = []

        for i, n_sources in enumerate(n_sources_batch):
            extents = self.get_source_extent(n=n_sources).astype(int)
            extents = np.clip(extents, 1, self.n_dipoles)

            sigmas = self.get_patch_smoothness_sigma(n=n_sources).astype(float)
            sigmas = np.clip(sigmas, 1e-3, None)

            patch_ranks = self.get_patch_rank(n=n_sources).astype(int)
            patch_ranks = np.clip(patch_ranks, 1, 2)

            amps = self.rng.uniform(*self.config.amplitude_range, n_sources)
            amplitude_values.append(amps)

            cov = get_cov(int(n_sources), float(inter_source_correlations[i]))
            base_choices_1 = self.rng.choice(self.config.n_timecourses, n_sources)
            tc_rank1 = (self.time_courses[base_choices_1].T @ cov).T

            base_choices_2 = self.rng.choice(self.config.n_timecourses, n_sources)
            tc_rank2 = (self.time_courses[base_choices_2].T @ cov).T

            centers_i = []
            source_vertices_i = []
            source_ranks_i = []

            for s in range(n_sources):
                region_vertices = self._sample_connected_vertices(int(extents[s]))
                components = self._make_contiguous_gaussian_patch_components(
                    region_vertices=region_vertices,
                    rank=int(patch_ranks[s]),
                    sigma=float(sigmas[s]),
                )

                if len(components) == 0:
                    continue

                patch_ts = components[0][:, None] * tc_rank1[s][None, :]
                if int(patch_ranks[s]) >= 2 and len(components) > 1:
                    patch_ts += components[1][:, None] * tc_rank2[s][None, :]

                y_batch[i] += (amps[s] / int(patch_ranks[s])) * patch_ts
                center_vertex = int(
                    region_vertices[np.argmax(components[0][region_vertices])]
                )
                centers_i.append(center_vertex)
                source_vertices_i.append(region_vertices)
                source_ranks_i.append(int(patch_ranks[s]))

            selection.append(np.asarray(centers_i, dtype=int))
            source_vertices_meta.append(source_vertices_i)
            source_ranks_meta.append(np.asarray(source_ranks_i, dtype=int))

        source_meta = {
            "source_vertices": source_vertices_meta,
            "source_ranks": source_ranks_meta,
        }

        return (
            y_batch,
            n_sources_batch,
            selection,
            amplitude_values,
            inter_source_correlations,
            noise_color_coeffs,
            source_meta,
        )

    def _generate_patches(self, batch_size):
        """Generate patch-based source activity.

        Returns:
            y_patches: [batch_size, n_dipoles, n_timepoints] patch activity
            selection: list of source index arrays
            amplitude_values: list of amplitude arrays
            inter_source_correlations: array of correlation values
            noise_color_coeffs: array of noise color coefficients
        """
        if self.config.source_spatial_model == "contiguous_gaussian":
            return self._generate_contiguous_gaussian_patches(batch_size)

        return self._generate_diffusion_basis_patches(batch_size)

    def _generate_background(self, batch_size, y_patches):
        """Mix background activity with patches if in mixture mode.

        Returns:
            y: [batch_size, n_dipoles, n_timepoints] combined activity
            alphas: mixing coefficients or None
        """
        if self.config.simulation_mode == "mixture":
            y_background = self._generate_smooth_background(batch_size)

            if isinstance(self.config.background_mixture_alpha, tuple):
                alphas = self.rng.uniform(
                    *self.config.background_mixture_alpha, batch_size
                )
            else:
                alphas = np.full(batch_size, self.config.background_mixture_alpha)

            alphas_bc = alphas[:, np.newaxis, np.newaxis]
            y = alphas_bc * y_background + (1 - alphas_bc) * y_patches
        else:
            y = y_patches
            alphas = None

        return y, alphas

    def _apply_noise(self, x, batch_size, noise_color_coeffs, modes_batch):
        """Apply realistic sensor noise with optional projection/rank loss.

        Returns
        -------
        x_noisy : np.ndarray
            Noisy EEG data [batch_size, n_channels, n_timepoints].
        noise_meta : dict
            Per-sample noise metadata (SNR, projector, covariance statistics).
        """
        snr_levels = self.rng.uniform(
            low=self.config.snr_range[0],
            high=self.config.snr_range[1],
            size=batch_size,
        )
        temporal_betas = self.get_noise_temporal_beta(n=batch_size)
        rank_losses = self.get_noise_rank_deficiency(n=batch_size)
        low_rank_dims = self.get_noise_low_rank_dim(n=batch_size)

        x_noisy_list = []
        realized_snrs = []
        noise_scales = []
        projectors = []
        projector_ranks = []
        covariance_model = []
        covariance_true = []
        covariance_est = []
        covariance_est_scaled = []
        covariance_rank_true = []
        covariance_rank_est = []

        for (
            xx,
            snr_level,
            corr_mode,
            noise_color_level,
            temporal_beta,
            rank_loss,
            low_rank_dim,
        ) in zip(
            x,
            snr_levels,
            modes_batch,
            noise_color_coeffs,
            temporal_betas,
            rank_losses,
            low_rank_dims,
            strict=False,
        ):
            P, _ = make_rank_projector(
                self.n_chans, rank_deficiency=int(rank_loss), rng=self.rng
            )
            if self.config.apply_sensor_projector:
                signal = P @ xx
            else:
                signal = xx

            cov_spatial = make_sensor_noise_covariance(
                self.n_chans,
                mode=corr_mode,
                noise_color_coeff=float(noise_color_level),
                rng=self.rng,
                low_rank_dim=int(low_rank_dim),
            )

            noise = sample_sensor_noise(
                cov_spatial,
                self.config.n_timepoints,
                rng=self.rng,
                temporal_beta=float(temporal_beta),
            )
            if self.config.apply_sensor_projector:
                noise = P @ noise

            noise_scaled, _scale, realized_snr = scale_noise_to_snr(
                signal, noise, float(snr_level)
            )
            x_noisy = signal + noise_scaled

            # Theoretical covariance of the realized noise (up to sampling error).
            if self.config.apply_sensor_projector:
                cov_model = (_scale**2) * (P @ cov_spatial @ P.T)
            else:
                cov_model = (_scale**2) * cov_spatial

            # True covariance from realized noise in this sample.
            cov_true = empirical_covariance(noise_scaled, shrinkage=0.0, eps=0.0)

            # Estimated covariance from independent baseline noise.
            cov_est = None
            cov_est_s = None
            if self.config.estimate_noise_cov:
                baseline = sample_sensor_noise(
                    cov_spatial,
                    int(self.config.noise_cov_n_baseline),
                    rng=self.rng,
                    temporal_beta=float(temporal_beta),
                )
                if self.config.apply_sensor_projector:
                    baseline = P @ baseline
                cov_est = empirical_covariance(
                    baseline, shrinkage=float(self.config.noise_cov_shrinkage)
                )
                # Match the baseline estimate to the scaled noise actually added to the data.
                cov_est_s = (_scale**2) * cov_est

            x_noisy_list.append(x_noisy)
            realized_snrs.append(realized_snr)
            noise_scales.append(_scale)
            projectors.append(P)
            projector_ranks.append(int(np.linalg.matrix_rank(P)))
            covariance_model.append(cov_model)
            covariance_true.append(cov_true)
            covariance_rank_true.append(int(np.linalg.matrix_rank(cov_true)))
            covariance_est.append(cov_est)
            covariance_est_scaled.append(cov_est_s)
            covariance_rank_est.append(
                int(np.linalg.matrix_rank(cov_est)) if cov_est is not None else -1
            )

        x_noisy = np.stack(x_noisy_list, axis=0)
        noise_meta = {
            "snr_target": snr_levels,
            "snr_realized": np.asarray(realized_snrs, dtype=float),
            "noise_scale": np.asarray(noise_scales, dtype=float),
            "projector": projectors,
            "projector_rank": np.asarray(projector_ranks, dtype=int),
            "noise_cov_model": covariance_model,
            "noise_cov_true": covariance_true,
            "noise_cov_est": covariance_est,
            "noise_cov_est_scaled": covariance_est_scaled,
            "noise_cov_rank_true": np.asarray(covariance_rank_true, dtype=int),
            "noise_cov_rank_est": np.asarray(covariance_rank_est, dtype=int),
            "noise_temporal_beta": np.asarray(temporal_betas, dtype=float),
            "noise_rank_deficiency": np.asarray(rank_losses, dtype=int),
            "noise_low_rank_dim": np.asarray(low_rank_dims, dtype=int),
        }
        return x_noisy, noise_meta

    def _build_metadata(
        self,
        batch_size,
        n_sources_batch,
        amplitude_values,
        inter_source_correlations,
        noise_color_coeffs,
        modes_batch,
        selection,
        source_meta,
        alphas,
        noise_meta,
    ):
        """Build simulation metadata DataFrame."""
        info_dict = {
            "n_sources": n_sources_batch,
            "amplitudes": amplitude_values,
            "snr": noise_meta["snr_target"],
            "snr_realized": noise_meta["snr_realized"],
            "noise_scale": noise_meta["noise_scale"],
            "inter_source_correlations": inter_source_correlations,
            "n_orders": [[self.min_order, self.max_order]] * batch_size,
            "diffusion_parameter": [self.config.diffusion_parameter] * batch_size,
            "n_timepoints": [self.config.n_timepoints] * batch_size,
            "n_timecourses": [self.config.n_timecourses] * batch_size,
            "correlation_mode": modes_batch,
            "noise_color_coeff": noise_color_coeffs,
            "noise_temporal_beta": noise_meta["noise_temporal_beta"],
            "noise_rank_deficiency": noise_meta["noise_rank_deficiency"],
            "projector_rank": noise_meta["projector_rank"],
            "noise_cov_rank_true": noise_meta["noise_cov_rank_true"],
            "add_forward_error": [self.config.add_forward_error] * batch_size,
            "forward_error": [self.config.forward_error] * batch_size,
            "centers": selection,
            "simulation_mode": [self.config.simulation_mode] * batch_size,
            "source_spatial_model": [self.config.source_spatial_model] * batch_size,
        }

        if source_meta.get("source_vertices", None) is not None:
            info_dict["source_vertices"] = source_meta["source_vertices"]
        if source_meta.get("source_ranks", None) is not None:
            info_dict["source_ranks"] = source_meta["source_ranks"]

        if self.config.estimate_noise_cov:
            info_dict["noise_cov_rank_est"] = noise_meta["noise_cov_rank_est"]

        if self.config.return_noise_cov:
            info_dict["projector"] = noise_meta["projector"]
            info_dict["noise_cov_model"] = noise_meta["noise_cov_model"]
            info_dict["noise_cov_true"] = noise_meta["noise_cov_true"]
            if self.config.estimate_noise_cov:
                info_dict["noise_cov_est"] = noise_meta["noise_cov_est"]
                info_dict["noise_cov_est_scaled"] = noise_meta["noise_cov_est_scaled"]

        if self.config.simulation_mode == "mixture":
            info_dict.update(
                {
                    "background_beta": [self.config.background_beta] * batch_size,
                    "background_mixture_alpha": alphas,
                }
            )

        return pd.DataFrame(info_dict)

    def generate(self):
        """Generate batches of simulations.

        Yields:
            tuple: (x, y, info) where:
                - x: EEG data [batch_size, n_channels, n_timepoints]
                - y: Source activity [batch_size, n_dipoles, n_timepoints] (scaled)
                - info: DataFrame with simulation metadata
        """
        # Setup correlation modes
        if (
            isinstance(self.config.correlation_mode, str)
            and self.config.correlation_mode.lower() == "auto"
        ):
            correlation_modes = ["cholesky", "banded", "diagonal", "low_rank", None]
            modes_batch = self.rng.choice(correlation_modes, self.config.batch_size)
        else:
            modes_batch = [self.config.correlation_mode] * self.config.batch_size

        while True:
            leadfield = self._setup_leadfield()

            (
                y_patches,
                n_sources_batch,
                selection,
                amplitude_values,
                inter_source_correlations,
                noise_color_coeffs,
                source_meta,
            ) = self._generate_patches(self.config.batch_size)

            y, alphas = self._generate_background(self.config.batch_size, y_patches)

            # Vectorized leadfield projection
            x = np.einsum("cd,bdt->bct", leadfield, y)

            x, noise_meta = self._apply_noise(
                x, self.config.batch_size, noise_color_coeffs, modes_batch
            )

            info = self._build_metadata(
                self.config.batch_size,
                n_sources_batch,
                amplitude_values,
                inter_source_correlations,
                noise_color_coeffs,
                modes_batch,
                selection,
                source_meta,
                alphas,
                noise_meta,
            )

            output = (x, y, info)

            for _ in range(self.config.batch_repetitions):
                yield output

__init__

__init__(
    fwd, config: SimulationConfig | None = None, **kwargs
)

Initialize the simulation generator.

Parameters: fwd: MNE forward solution object config: SimulationConfig instance (optional) **kwargs: Configuration parameters (used if config is None)

Source code in invert/simulate/simulate.py

def __init__(self, fwd, config: SimulationConfig | None = None, **kwargs):
    """Initialize the simulation generator.

    Parameters:
        fwd: MNE forward solution object
        config: SimulationConfig instance (optional)
        **kwargs: Configuration parameters (used if config is None)
    """
    # Initialize configuration
    if config is None:
        self.config = SimulationConfig(**kwargs)
    else:
        self.config = config

    # Store forward solution components
    self.fwd = fwd
    self.channel_types = np.array(fwd["info"].get_channel_types())
    self.leadfield_original = deepcopy(fwd["sol"]["data"])
    self.leadfield = deepcopy(self.leadfield_original)
    self.n_chans, self.n_dipoles = self.leadfield.shape

    # Initialize random number generator
    self.rng = np.random.default_rng(self.config.random_seed)

    # Parse n_sources parameter
    if isinstance(self.config.n_sources, (int, float)):
        n_sources_val = np.clip(self.config.n_sources, a_min=1, a_max=np.inf)
        self.min_sources, self.max_sources = int(n_sources_val), int(n_sources_val)
    else:
        self.min_sources, self.max_sources = self.config.n_sources

    # Precompute spatial smoothing operator
    self._precompute_spatial_operators()

    # Precompute timecourses
    self._precompute_timecourses()

    # Setup correlation sampling functions
    self._setup_correlation_samplers()

    # Normalize leadfield if requested
    if self.config.normalize_leadfield:
        self.leadfield /= np.linalg.norm(self.leadfield, axis=0)

generate

generate()

Generate batches of simulations.

Yields: tuple: (x, y, info) where: - x: EEG data [batch_size, n_channels, n_timepoints] - y: Source activity [batch_size, n_dipoles, n_timepoints] (scaled) - info: DataFrame with simulation metadata

Source code in invert/simulate/simulate.py

def generate(self):
    """Generate batches of simulations.

    Yields:
        tuple: (x, y, info) where:
            - x: EEG data [batch_size, n_channels, n_timepoints]
            - y: Source activity [batch_size, n_dipoles, n_timepoints] (scaled)
            - info: DataFrame with simulation metadata
    """
    # Setup correlation modes
    if (
        isinstance(self.config.correlation_mode, str)
        and self.config.correlation_mode.lower() == "auto"
    ):
        correlation_modes = ["cholesky", "banded", "diagonal", "low_rank", None]
        modes_batch = self.rng.choice(correlation_modes, self.config.batch_size)
    else:
        modes_batch = [self.config.correlation_mode] * self.config.batch_size

    while True:
        leadfield = self._setup_leadfield()

        (
            y_patches,
            n_sources_batch,
            selection,
            amplitude_values,
            inter_source_correlations,
            noise_color_coeffs,
            source_meta,
        ) = self._generate_patches(self.config.batch_size)

        y, alphas = self._generate_background(self.config.batch_size, y_patches)

        # Vectorized leadfield projection
        x = np.einsum("cd,bdt->bct", leadfield, y)

        x, noise_meta = self._apply_noise(
            x, self.config.batch_size, noise_color_coeffs, modes_batch
        )

        info = self._build_metadata(
            self.config.batch_size,
            n_sources_batch,
            amplitude_values,
            inter_source_correlations,
            noise_color_coeffs,
            modes_batch,
            selection,
            source_meta,
            alphas,
            noise_meta,
        )

        output = (x, y, info)

        for _ in range(self.config.batch_repetitions):
            yield output

Noise Functions

invert.simulate.add_white_noise

add_white_noise(
    X_clean,
    snr,
    rng,
    channel_types,
    noise_color_coeff=0.5,
    correlation_mode=None,
)

Parameters:

Name	Type	Description	Default
X_clean	ndarray	The clean EEG data.	required
snr	float	The signal to noise ratio in dB.	required
correlation_mode	None / str	None implies no correlation between the noise in different channels. 'banded' : Colored banded noise, where channels closer to each other will be more correlated. 'diagonal' : Some channels have varying degrees of noise. 'cholesky' : A set correlation coefficient between each pair of channels	None
noise_color_coeff	float	The magnitude of spatial coloring of the noise (not the magnitude of noise overall!).	0.5

Source code in invert/simulate/noise.py

def add_white_noise(
    X_clean, snr, rng, channel_types, noise_color_coeff=0.5, correlation_mode=None
):
    """
    Parameters
    ----------
    X_clean : numpy.ndarray
        The clean EEG data.
    snr : float
        The signal to noise ratio in dB.
    correlation_mode : None/str
        None implies no correlation between the noise in different channels.
        'banded' : Colored banded noise, where channels closer to each other will be more correlated.
        'diagonal' : Some channels have varying degrees of noise.
        'cholesky' : A set correlation coefficient between each pair of channels
    noise_color_coeff : float
        The magnitude of spatial coloring of the noise (not the magnitude of noise overall!).
    """
    n_chans, n_time = X_clean.shape
    X_noise = rng.standard_normal((n_chans, n_time))
    snr_linear = 10 ** (snr / 10)

    if isinstance(channel_types, list):
        channel_types = np.array(channel_types)
    # Ensure the channel_types array is correct length
    assert len(channel_types) == n_chans, (
        "Length of channel_types must match the number of channels in X_clean"
    )

    unique_types = np.unique(channel_types)
    X_full = np.zeros_like(X_clean)

    for ch_type in unique_types:
        type_indices = np.where(channel_types == ch_type)[0]
        X_clean_type = X_clean[type_indices, :]
        X_noise_type = X_noise[type_indices, :]
        if isinstance(noise_color_coeff, str) and isinstance(
            correlation_mode, np.ndarray
        ):
            # Real Noise Covariance
            X_noise_type = (
                np.linalg.cholesky(correlation_mode[type_indices][:, type_indices])
                @ X_noise_type
            )
        elif correlation_mode == "cholesky":
            covariance_matrix = np.full(
                (len(type_indices), len(type_indices)), noise_color_coeff
            )
            np.fill_diagonal(covariance_matrix, 1)  # Set diagonal to 1 for variance

            # Cholesky decomposition
            X_noise_type = np.linalg.cholesky(covariance_matrix) @ X_noise_type
        elif correlation_mode == "banded":
            num_sensors = X_noise_type.shape[0]
            Y = np.zeros_like(X_noise_type)
            for i in range(num_sensors):
                Y[i, :] = X_noise_type[i, :]
                for j in range(num_sensors):
                    if abs(i - j) % num_sensors == 1:
                        Y[i, :] += (noise_color_coeff / np.sqrt(2)) * X_noise_type[j, :]
            X_noise_type = Y
        elif correlation_mode == "diagonal":
            X_noise_type[1::3, :] *= 1 - noise_color_coeff
            X_noise_type[2::3, :] *= 1 + noise_color_coeff
        elif correlation_mode is None:
            pass
        else:
            msg = f"correlation_mode can be either None, cholesky, banded or diagonal, but was {correlation_mode}"
            raise AttributeError(msg)

        rms_noise = rms(X_noise_type)
        rms_signal = rms(X_clean_type)
        scaler = rms_signal / (snr_linear * rms_noise)

        X_full[type_indices] = X_clean_type + X_noise_type * scaler

    return X_full

invert.simulate.powerlaw_noise

powerlaw_noise(beta, n_timepoints, n_signals=1, rng=None)

Generate 1/f^beta colored noise via FFT spectral shaping.

Parameters:

Name	Type	Description	Default
beta	float or array - like	Power-law exponent(s). 0=white, 1=pink, 2=brown. If array, must have length n_signals.	required
n_timepoints	int	Number of time samples.	required
n_signals	int	Number of independent signals to generate.	1
rng	Generator or None	Random number generator.	None

Returns:

Name	Type	Description
signals	(ndarray, shape(n_signals, n_timepoints))	Colored noise signals.

Source code in invert/simulate/noise.py

def powerlaw_noise(beta, n_timepoints, n_signals=1, rng=None):
    """Generate 1/f^beta colored noise via FFT spectral shaping.

    Parameters
    ----------
    beta : float or array-like
        Power-law exponent(s). 0=white, 1=pink, 2=brown.
        If array, must have length n_signals.
    n_timepoints : int
        Number of time samples.
    n_signals : int
        Number of independent signals to generate.
    rng : numpy.random.Generator or None
        Random number generator.

    Returns
    -------
    signals : ndarray, shape (n_signals, n_timepoints)
        Colored noise signals.
    """
    if rng is None:
        rng = np.random.default_rng()

    beta = np.atleast_1d(np.asarray(beta, dtype=float))
    if beta.shape[0] == 1:
        beta = np.broadcast_to(beta, (n_signals,))

    # Generate white noise and FFT
    white = rng.standard_normal((n_signals, n_timepoints))
    spectrum = np.fft.rfft(white)

    # Build frequency-domain filter: 1/f^(beta/2) (power spectrum goes as 1/f^beta)
    freqs = np.fft.rfftfreq(n_timepoints)
    # Skip DC (index 0) to avoid division by zero
    freq_filter = np.ones((n_signals, len(freqs)))
    freq_filter[:, 1:] = freqs[1:][np.newaxis, :] ** (-beta[:, np.newaxis] / 2.0)

    spectrum *= freq_filter
    signals = np.fft.irfft(spectrum, n=n_timepoints)

    return signals

Covariance Utilities

invert.simulate.compute_covariance

compute_covariance(Y, cov_type='basic')

Compute the covariance matrix of the data.

Parameters:

Name	Type	Description	Default
Y	ndarray	The data matrix.	required
cov_type	str	The type of covariance matrix to compute. Options are 'basic' and 'SSM'. Default is 'basic'.	'basic'

Return

C : numpy.ndarray The covariance matrix.

Source code in invert/simulate/covariance.py

def compute_covariance(Y, cov_type="basic"):
    """Compute the covariance matrix of the data.

    Parameters
    ----------
    Y : numpy.ndarray
        The data matrix.
    cov_type : str
        The type of covariance matrix to compute. Options are 'basic' and 'SSM'. Default is 'basic'.

    Return
    ------
    C : numpy.ndarray
        The covariance matrix.
    """
    if cov_type == "basic":
        C = Y @ Y.T
    elif cov_type == "SSM":
        n_time = Y.shape[1]
        M_Y = Y.T @ Y
        YY = M_Y + 0.001 * (50 / n_time) * np.trace(M_Y) * np.eye(n_time)
        P_Y = (Y @ np.linalg.inv(YY)) @ Y.T
        C = P_Y.T @ P_Y
    elif cov_type == "riemann":
        msg = "Riemannian covariance is not yet implemented as a standalone function."
        raise NotImplementedError(msg)
    else:
        msg = "Covariance type not recognized. Use 'basic', 'SSM' or provide a custom covariance matrix."
        raise ValueError(msg)

    return C

invert.simulate.gen_correlated_sources

gen_correlated_sources(corr_coeff, T, Q)

Generate Q correlated sources with a specified correlation coefficient. The sources are generated as sinusoids with random frequencies and phases.

Parameters:

Name	Type	Description	Default
corr_coeff	float	The correlation coefficient between the sources.	required
T	int	The number of time points in the sources.	required
Q	int	The number of sources to generate.	required

Returns:

Name	Type	Description
Y	ndarray	The generated sources.

Source code in invert/simulate/covariance.py

def gen_correlated_sources(corr_coeff, T, Q):
    """Generate Q correlated sources with a specified correlation coefficient.
    The sources are generated as sinusoids with random frequencies and phases.

    Parameters
    ----------
    corr_coeff : float
        The correlation coefficient between the sources.
    T : int
        The number of time points in the sources.
    Q : int
        The number of sources to generate.

    Returns
    -------
    Y : numpy.ndarray
        The generated sources.
    """
    Cov = np.ones((Q, Q)) * corr_coeff + np.diag(
        np.ones(Q) * (1 - corr_coeff)
    )  # required covariance matrix
    freq = np.random.randint(10, 31, Q)  # random frequencies between 10Hz to 30Hz

    phases = 2 * np.pi * np.random.rand(Q)  # random phases
    t = np.linspace(10 * np.pi / T, 10 * np.pi, T)
    Signals = np.sqrt(2) * np.cos(
        2 * np.pi * freq[:, None] * t + phases[:, None]
    )  # the basic signals

    if corr_coeff < 1:
        A = np.linalg.cholesky(Cov).T  # cholesky Decomposition
        Y = A @ Signals
    else:  # Coherent Sources
        Y = np.tile(Signals[0, :], (Q, 1))

    return Y