Source Reconstruction Examples¶

pyAVS provides source reconstruction via LCMV beamforming: forward modeling, per-session filter computation, ROI-based data extraction, and population code computation for encoding analyses. This page walks through the two real, runnable example scripts that demonstrate the full pipeline, plus the lower-level functions they’re built from.

Quickstart: Synthetic Data¶

examples/simple_source_reconstruction.py is the fastest way to see the pipeline run end to end. It fabricates synthetic MEG epochs and a synthetic forward model, so it needs no real AVS data and runs standalone:

#!/usr/bin/env python3
"""
Simple pyAVS Source Reconstruction Example (synthetic-data quickstart)

A minimal example showing how to:
1. Create synthetic MEG data
2. Perform source reconstruction
3. Save data using pyAVS I/O system

This is the quickstart to reach for when you don't have (or don't yet want to
configure) real AVS data - it runs standalone against synthetic MEG/forward-model
data. For a complete, real-data, config-driven workflow, see
compute_population_codes_example.py instead.
"""

import numpy as np
import pandas as pd
import mne

import pyavs
from pyavs.io.write import save_population_codes_h5, save_epochs
from pyavs.source.reconstruction import compute_beamformer_filters, apply_beamformer

The core reconstruction step is a two-call sequence – compute beamformer filters, then apply them:

# Compute beamformer filters
filters = compute_beamformer_filters(
    epochs=epochs,
    forward=forward,
    reg=0.05,
    weight_norm='unit-noise-gain',
    verbose=False
)

# Apply beamformer
source_data = apply_beamformer(
    epochs=epochs,
    filters=filters,
    verbose=False
)

print(f"Source reconstruction complete. Shape: {source_data.shape}")

See the full script (examples/simple_source_reconstruction.py) for the surrounding synthetic-data setup and how results are saved via pyavs.io.write.save_population_codes_h5().

Real-Data, Config-Driven Workflow¶

examples/compute_population_codes_example.py mirrors the full AVS-machine-room population-code pipeline against real subject data: MEG loading and ICA via AVSComposer, eye-tracking event epoching, per-session LCMV filter loading/computation, ROI extraction, and HDF5 storage. It uses the PyAVSConfig system (see Configuration System) to drive all analysis parameters rather than hardcoding them.

The source-level reconstruction step loads or computes per-session LCMV filters and applies them to the current session’s epochs:

def compute_source_population_codes(epochs, composer, subject_id, session_num, 
                                  source_rois, method, pick_ori, output_dir, event_type, data_path,
                                  tmin, tmax, filter_params, resample_to_hz, hemi, blocks=None):
    """Compute population codes for source-level data using per-session LCMV filters."""
    
    logger.info("Setting up source reconstruction with per-session LCMV filters...")
    
    try:
        # Import the new filter management system
        from pyavs.source.filters import load_or_compute_lcmv_filters, apply_lcmv_to_epochs
        
        # Load or compute per-session LCMV filters with full parameter set
        logger.info(f"Loading/computing LCMV filters for event type: {event_type}")
        
        # Get filter kwargs from config for consistent parameter usage
        from pyavs.config import get_config
        config = get_config()
        filter_kwargs = config.get_filter_kwargs()
        
        # Load or compute filters with config-derived parameters
        filters = load_or_compute_lcmv_filters(
            data_path=data_path,
            subject_id=subject_id,
            sessions=[session_num],  # Only need current session for application
            event_type=event_type,
            **filter_kwargs
        )
        
        # Apply beamformer filters to epochs
        logger.info(f"Applying LCMV beamformer for session {session_num}...")
        stcs = apply_lcmv_to_epochs(epochs, filters, session_num)
        
        # Extract ROI data
        population_codes = {}
        for roi in source_rois:
            if roi == "stc":
                # Full source space
                population_codes[roi] = np.array([stc.data for stc in stcs])
            else:
                # Specific ROI (would need label files)
                try:
                    label_fname = os.path.join(composer.subject_dir, "label", f"lh.L_{roi}_ROI.label")
                    label = mne.read_label(label_fname, subject=f"as{subject_id:02d}")
                    roi_data = []
                    for stc in stcs:
                        stc_roi = stc.in_label(label)
                        roi_data.append(stc_roi.data)
                    population_codes[roi] = np.array(roi_data)
                except:
                    logger.warning(f"Could not load ROI {roi}, creating mock data")
                    n_sources = 50  # Mock number of sources in ROI
                    population_codes[roi] = np.random.randn(len(epochs), n_sources, len(epochs.times))
        
        return population_codes
        
    except Exception as e:
        logger.error(f"Source reconstruction failed: {e}")
        raise

Note that this uses load_or_compute_lcmv_filters() (which wraps compute_cross_session_data_covariance() and compute_per_session_lcmv_filters() – see Cross-Session Beamformer Filters for computing those filters as a standalone step) and apply_lcmv_to_epochs(), not the lower-level apply_beamformer() shown above – the per-session-filter route is what the real pipeline uses when filters need to be shared/reused across a session, while apply_beamformer is the more direct call for a one-off epochs/forward/filters triple.

ROI-Based Extraction¶

Once you have source-space data (an array of shape (n_epochs, n_sources, n_times)), use pyavs.extract_roi_data() to average or select within named regions:

import pyavs

# Glasser atlas ROI names for one area category at a time
# (area is a single string: 'all', 'high_visual', 'early_visual', or 'intermediate_visual')
visual_rois = pyavs.get_glasser_roi_labels(area='early_visual')

roi_data = pyavs.extract_roi_data(
    source_data,          # np.ndarray, shape (n_epochs, n_sources, n_times)
    forward_model['src'],
    roi_labels=visual_rois,
    method='mean',         # average activity within each ROI
    verbose=True,
)

for roi_name, data in roi_data.items():
    print(f"{roi_name}: {data.shape}")  # (n_epochs, n_times)

Population Codes¶

pyavs.compute_population_codes() turns source-space data plus per-epoch metadata into condition-averaged population codes:

population_codes = pyavs.compute_population_codes(
    source_data,             # np.ndarray, shape (n_epochs, n_sources, n_times)
    events_metadata=epochs.metadata,
    conditions=['scene_id'],  # column(s) in metadata defining conditions
    time_window=(0.0, 0.3),
    times=epochs.times,
)