AVSComposer Guide

AVSComposer is the recommended entry point for MEG + eye-tracking data fusion in pyAVS: it loads MEG blocks, applies ICA artifact removal and filtering, concatenates blocks per session, finds MEG trigger events, and aligns eye-tracking events (fixations, saccades, blinks, or whole-scene trials) to build epoched MEG data with rich per-epoch metadata.

Full Worked Example

examples/avs_composer_example.py runs the complete pipeline end to end – MEG loading, filtering, ICA, event finding, eye-tracking annotation and epoching for every event type, and a simple median-ERF visualization:

"""
Example demonstrating the use of AVS Composer for MEG-ET data fusion.

This example shows how to use the AVSComposer class to replicate the functionality
of the original AVS-machine-room composer workflow in the pyAVS package.

NEW FEATURE: ET Event Offset Timing
- Set onset_offset=True to use event end timing (time_in_trial + duration) 
  instead of event onset for MEG epoch creation
- This is a fringe feature for specialized analyses requiring offset timing

Author: P. Sulewski (phsulewski@gmail.com)
"""

import os
import numpy as np
import pyavs
from pyavs.utils.logging import get_logger, configure_logging

def main():
    """Run AVS Composer example."""
    
    # Configure logging for better output formatting
    configure_logging(level='DEBUG', console=True)
    logger = get_logger('examples.avs_composer')
    
    logger.info("=== pyAVS AVS Composer Example ===")
    
    # Configuration
    subject_id = 4
    session = 2
    data_path = pyavs.get_data_path()
    if data_path is None:
        logger.error("No data path configured. Run: pyavs configure --data-path /path/to/data")
        return

    # NEW: ET event offset parameter for fringe feature
    # When True, uses time_in_trial + duration as epoch timing instead of event onset
    onset_offset = False  # Set to True to use event offset timing instead of onset
    recording = "scene"  # Use "caption" for captioned data, "scene" for scene data
    
    # Initialize AVS Composer
    logger.info("\n1. Initializing AVS Composer...")

    composer = pyavs.AVSComposer(
        subject=subject_id,
        session_num=session,
        data_path=data_path,
        output_path=data_path,
        et_path=data_path,
        preprocessed=True,
        recompute_prepro=False,
        max_block=None,  
        min_block=1,
        verbose=True,
        interpolate_bad_channels=True,
        use_precomputed_ica=True,  # Enable ICA artifact removal with precomputed solutions
        apply_ica=False,  # Set to True to compute ICA on-the-fly instead
        l_freq=0.2,  # Low-pass frequency for filtering
        h_freq=200,  # High-pass frequency for filtering
        causal_filter=False,  # Use causal filtering for temporal order preservation
        resample_freq=500.0  # Target sampling frequency
    )
    logger.info(f"   AVS Composer initialized for subject {subject_id}, session {session}")
    logger.info(f"   Selected blocks: {composer.blocks_this_session}")

    
    # Load MEG data
    logger.info("\n2. Loading MEG data...")
    try:
        composer.load_meg_data(compute_missing_prepro=False)
        logger.info(f"   Loaded MEG data for blocks: {list(composer.raws_dict.keys())}")
        logger.info(f"   Empty room recordings available: {composer.empty_room_available}")
    except Exception as e:
        logger.error(f"   Error loading MEG data: {e}")
        return
    
    
    # Filter MEG data (note: filtering is handled by preprocess_meg_block when recompute_prepro=True)
    logger.info("\n3. MEG preprocessing and filtering...")
    if composer.recompute_prepro:
        logger.info("   Filtering handled automatically by preprocess_meg_block during data loading")
        logger.info(f"   Applied {composer.l_freq}-{composer.h_freq} Hz band-pass filter ({'causal' if composer.causal_filter else 'non-causal'})")
    else:
        try:
            composer.filter_meg_data(ignore_existing_filter=True)  # Uses instance variables as defaults
            logger.info(f"   Applied {composer.l_freq}-{composer.h_freq} Hz band-pass filter ({'causal' if composer.causal_filter else 'non-causal'})")
        except Exception as e:
            logger.error(f"   Error filtering MEG data: {e}")
            return
        
    # Apply ICA artifact removal to raw data blocks
    logger.info("\n4. Applying ICA artifact removal...")
    composer.apply_ica_to_blocks()
    logger.info("   ICA artifact removal completed for all blocks")

    # Concatenate MEG blocks
    logger.info("\n5. Concatenating MEG blocks...")
    try:
        composer.concatenate_raws_per_session()
        logger.info(f"   Concatenated {len(composer.raws_dict)} MEG blocks")
        logger.info(f"   Total channels: {composer.raws_concatenated.info['nchan']}")
        logger.info(f"   Total samples: {len(composer.raws_concatenated.times)}")
        logger.info(f"   Duration: {composer.raws_concatenated.times[-1]:.2f} seconds")
    except Exception as e:
        logger.error(f"   Error concatenating MEG blocks: {e}")
        return
    
    # Find MEG events
    logger.info("\n6. Finding MEG events...")
    try:
        composer.find_events_in_raw()
        logger.info(f"   Found {len(composer.meg_trigger_events)} MEG events")
        unique_event_ids = np.unique(composer.meg_trigger_events[:, 2])
        logger.info(f"   Unique event IDs: {unique_event_ids}")
    except Exception as e:
        logger.error(f"   Error finding MEG events: {e}")
        return
    
    # Get eye tracking annotations and create epochs
    logger.info("\n7. Processing eye tracking data...")
    
    # Process each event type separately (new pyAVS approach)
    event_types = ["fixation", "blink", "saccade", "scene"]
    epochs_results = {}
    
    for event_type in event_types:
        logger.info(f"\n   Processing {event_type} events...")

        # Get annotations for this event type
        # Use onset_offset parameter to control timing (onset vs offset)
        composer.get_et_annotations(
            event_type=event_type,
            recording=recording,
            exclude_last_fixation=True,
            add_cross_event_info=True,
            preprocessed=True,
            onset_offset="offset" if onset_offset else "onset"
        )
        logger.info(f"   Loaded {len(composer.et_events)} {event_type} events")
        logger.info(f"   Added {event_type} annotations to MEG data")
        logger.info(f"   Total annotations: {len(composer.raws_annotated.annotations)}")
        
        # Create epochs for this event type
        timing_mode = "event offset" if onset_offset else "event onset"
        logger.info(f"      Creating {event_type} epochs using {timing_mode} timing...")
        if event_type == "scene":
            tmax = 0.5
        else:
            tmax = 0.5# Longer time window for scene events
        composer.make_et_event_epochs(
            tmin=-0.2,
            tmax=tmax,
            event_type=event_type,
            recording=recording,
            get_metadata=True,
            baseline=None
        )
        
        # Store results
        epochs_results[event_type] = composer.et_epochs
        n_epochs = len(composer.et_epochs)
        logger.info(f"   Created {n_epochs} {event_type} epochs")
        
        # Show some metadata columns
        if hasattr(composer.et_epochs, 'metadata') and composer.et_epochs.metadata is not None:
            metadata_cols = list(composer.et_epochs.metadata.columns)[:5]
            logger.info(f"   Metadata columns (first 5): {metadata_cols}")
        
   
    
    # Create simple median ERF plots
    logger.info("\n8. Creating median ERF plots...")
    try:
        import matplotlib.pyplot as plt
        
        # Create figure with subplots for each event type
        fig, axes = plt.subplots(len(epochs_results),1, figsize=(10, 3 * len(epochs_results)))
        axes = axes.flatten() if len(epochs_results) > 1 else [axes]  # Ensure axes is always iterable
        if len(epochs_results) == 1:
            axes = [axes]
        
        for idx, (event_type, epochs) in enumerate(epochs_results.items()):
            # Calculate median ERF across all epochs
            # Use magnetometers for cleaner visualization
            mag_picks = epochs.copy().pick_types(meg='mag')
            if len(mag_picks) > 0:
                evoked_median = mag_picks.average()
                
                # Plot the median ERF
                evoked_median.plot(axes=axes[idx], show=False, time_unit='ms')
                axes[idx].set_title(f'{event_type.capitalize()} Median ERF\n({len(epochs)} epochs)')
                axes[idx].set_ylabel('Magnetic Field (fT)')
                axes[idx].grid(True, alpha=0.3)
                
                logger.info(f"   Created median ERF plot for {event_type} ({len(epochs)} epochs)")
            else:
                logger.warning(f"   No magnetometer data found for {event_type}")
        
        plt.tight_layout()
        plt.savefig(f'avs_composer_median_erf_subject_{subject_id}_session_{session}_recording_{recording}.png', 
                   dpi=150, bbox_inches='tight')
        plt.close()
        
        logger.info("   Saved median ERF plots to file")
        
    except Exception as e:
        logger.error(f"   Error creating ERF plots: {e}")
        
        
    # report the time in trial of the event types
    logger.info("\n   Reporting time in trial for each event type...")
    for event_type, epochs in epochs_results.items():
        if hasattr(epochs, 'metadata') and 'time_in_trial' in epochs.metadata.columns:
            time_in_trial = np.nanmean(epochs.metadata['time_in_trial'].values)
            logger.info(f"   Average time in trial for {event_type}: {time_in_trial:.2f} seconds")
        elif hasattr(epochs, 'metadata') and 'time_to_first_event' in epochs.metadata.columns:
            # Scene epochs are trial-level; there's no single 'time_in_trial',
            # but 'time_to_first_event' gives the onset of the first fixation/
            # saccade/blink in the trial.
            time_to_first_event = np.nanmean(epochs.metadata['time_to_first_event'].values)
            logger.info(f"   Average time to first event for {event_type}: {time_to_first_event:.2f} seconds")
        else:
            logger.warning(f"   No 'time_in_trial' metadata found for {event_type} epochs")
    
    # Get data summary
    logger.info("\n9. Data summary...")
    try:
        summary = composer.get_data_summary()
        logger.info(f"   Subject: {summary['subject']}")
        logger.info(f"   Session: {summary['session']}")
        logger.info(f"   Blocks loaded: {summary['blocks_loaded']}")
        logger.info(f"   MEG channels: {summary['meg_channels']}")
        logger.info(f"   MEG duration: {summary['meg_duration']:.2f} seconds")
        logger.info(f"   Eye events: {summary['eye_events']}")
        logger.info(f"   Epochs created: {summary['epochs_created']}")
        logger.info(f"   Annotations: {summary['annotations']}")
    except Exception as e:
        logger.error(f"   Error getting data summary: {e}")
        return
    
    logger.info("\n=== AVS Composer Example Complete ===")
    logger.info("This example demonstrated:")
    logger.info("- MEG data loading and preprocessing using pyAVS meg.py functions")
    logger.info("- Standalone ICA artifact removal applied to unconcatenated blocks")
    logger.info("- Eye tracking data integration with single event type processing")
    logger.info("- Trigger-based MEG-ET alignment")
    logger.info("- Epoch creation with metadata for multiple event types")
    if onset_offset:
        logger.info("- ET event offset timing: using time_in_trial + duration for epoch timing")
    logger.info("- Simple median ERF visualization for different ET event types")
    logger.info("- Replication of AVS-machine-room composer functionality in pyAVS")
    logger.info("- Modular preprocessing pipeline with separated ICA processing")



if __name__ == "__main__":
    # Run full example
    main()
    

Advanced Initialization

The composer accepts many optional parameters beyond the basics shown above – data paths, block range, ICA source (precomputed vs. computed on the fly), and filter/resample settings:

composer = pyavs.AVSComposer(
    subject=1,
    session_num=1,

    # Data paths (default to the globally configured data path if omitted)
    data_path="/path/to/avs/dataset",
    output_path="/path/to/outputs",
    et_path="/path/to/eye_tracking",

    # Processing options
    preprocessed=True,
    recompute_prepro=False,      # set True to recompute preprocessing
    max_block=10,                # process blocks 1-10
    min_block=1,

    # Bad channel handling
    interpolate_bad_channels=True,

    # ICA artifact removal
    apply_ica=False,             # compute ICA on the fly
    use_precomputed_ica=True,    # use existing ICA solutions
    ica_solutions_path="/path/to/ica",
    ica_exclusions_file="/path/to/exclusions.json",

    # Filtering
    l_freq=0.2,                  # high-pass frequency
    h_freq=200.0,                # low-pass frequency
    causal_filter=True,          # causal filtering preserves timing

    # Resampling
    resample_freq=500.0,

    # Misc
    n_jobs=4,
    random_state=42,
    verbose=True,
    write_output=True,
)

Working with Empty-Room Recordings

The composer detects and separates empty-room recordings automatically during load_meg_data():

composer = pyavs.AVSComposer(subject=1, session_num=1)
composer.load_meg_data()

if composer.empty_room_available:
    print("Empty room recordings found!")
    print(f"Empty room blocks: {list(composer.raws_dict_empty_room.keys())}")
    composer.concatenate_raws_per_session()
    print(f"Empty room duration: {composer.raws_concatenated_empty_room.times[-1]:.1f}s")
else:
    print("No empty room recordings available")

Data Summary

get_data_summary() returns a dict summarizing what’s been loaded/processed so far:

summary = composer.get_data_summary()
print(f"Subject: {summary['subject']}, Session: {summary['session']}")
print(f"Blocks loaded: {summary['blocks_loaded']}")
print(f"MEG channels: {summary['meg_channels']}, duration: {summary['meg_duration']:.1f}s")
print(f"Eye events: {summary['eye_events']}, epochs created: {summary['epochs_created']}")
print(f"Empty room available: {summary['empty_room_available']}")

Integration with Source Reconstruction

Composer epochs feed directly into source reconstruction:

composer = pyavs.AVSComposer(subject=1, session_num=1)
composer.load_meg_data()
composer.apply_ica_to_blocks()
composer.concatenate_raws_per_session()
composer.find_events_in_raw()
composer.get_et_annotations(event_type="fixation")
composer.make_et_event_epochs(tmin=-0.2, tmax=0.5, event_type="fixation")

epochs = composer.et_epochs

forward_model = pyavs.load_forward_model(subject_id=1, session=1)
source_data = pyavs.apply_source_reconstruction(
    epochs, forward_model, method='beamformer'
)
roi_data = pyavs.extract_roi_data(
    source_data, forward_model['src'],
    roi_labels=pyavs.get_glasser_roi_labels(area='early_visual'),
)

print(f"Source data: {source_data.shape} (epochs, sources, timepoints)")

See Source Reconstruction and Population Codes and Source Reconstruction Examples for the full pipeline, including population code computation and storage.

See Also