+ end
+ end
+ end
+
+ if isempty(files)
+ % Return a placeholder so the test class can still be instantiated
+ files = {'__no_golden_files__'};
+ end
+ end
+ end
+end
+
+
+function out = extractOutput(processed)
+% Extract key output fields for golden comparison
+out = struct();
+fields = {'HbO', 'HbR', 'HbTotal', 'HbDiff', 'CBSI', 'units', 'DPF_factor'};
+for i = 1:length(fields)
+ f = fields{i};
+ if isfield(processed, f)
+ out.(f) = processed.(f);
+ end
+end
+end
+
+
+function result = compareOutputs(expected, actual, tolerance)
+% Compare two output structs field by field
+result = struct('passed', true, 'failures', {{}});
+
+fields = fieldnames(expected);
+for i = 1:length(fields)
+ fname = fields{i};
+ if ~isfield(actual, fname)
+ result.passed = false;
+ result.failures{end+1} = sprintf('Missing field: %s', fname);
+ continue;
+ end
+
+ exp = expected.(fname);
+ act = actual.(fname);
+
+ if isnumeric(exp) && isnumeric(act)
+ if ~isequal(size(exp), size(act))
+ result.passed = false;
+ result.failures{end+1} = sprintf('%s: size mismatch [%s] vs [%s]', ...
+ fname, mat2str(size(exp)), mat2str(size(act)));
+ continue;
+ end
+ maxDiff = max(abs(exp(:) - act(:)), [], 'omitnan');
+ if maxDiff > tolerance
+ result.passed = false;
+ result.failures{end+1} = sprintf('%s: max diff %.2e exceeds tolerance %.2e', ...
+ fname, maxDiff, tolerance);
+ end
+ elseif ischar(exp) && ischar(act)
+ if ~strcmp(exp, act)
+ result.passed = false;
+ result.failures{end+1} = sprintf('%s: string mismatch', fname);
+ end
+ elseif ~isequal(exp, act)
+ result.passed = false;
+ result.failures{end+1} = sprintf('%s: values do not match', fname);
+ end
+end
+end
diff --git a/+pf2_base/+tests/+synthetic/generateFNIRS.m b/+pf2_base/+tests/+synthetic/generateFNIRS.m
index 2826f481..06403fed 100644
--- a/+pf2_base/+tests/+synthetic/generateFNIRS.m
+++ b/+pf2_base/+tests/+synthetic/generateFNIRS.m
@@ -67,8 +67,8 @@
% time - Time vector [T x 1] in seconds
% fs - Sampling frequency in Hz
% fchMask - Channel mask [1 x nChannels], all ones (all good)
-% markers - Event markers [M x 3] from hrfOnsets
-% Format: [time, value, duration]
+% markers - Event markers [M x 4] from hrfOnsets
+% Format: [time, value, duration, amplitude]
% info - Metadata struct with:
% .header.filename = 'synthetic'
% .probename = 'synthetic'
@@ -237,10 +237,11 @@
end
% Create markers from HRF onsets
- markers = zeros(length(hrfOnsets), 3);
+ markers = zeros(length(hrfOnsets), 4);
markers(:, 1) = hrfOnsets(:); % Time
markers(:, 2) = 1; % Marker value
markers(:, 3) = 0; % Duration (impulse)
+ markers(:, 4) = 1; % Amplitude (default)
end
%% Add cardiac artifact (heartbeat)
diff --git a/+pf2_base/+tests/+synthetic/generateHemoglobin.m b/+pf2_base/+tests/+synthetic/generateHemoglobin.m
index 74a0b6c8..3fb9ae3b 100644
--- a/+pf2_base/+tests/+synthetic/generateHemoglobin.m
+++ b/+pf2_base/+tests/+synthetic/generateHemoglobin.m
@@ -64,8 +64,8 @@
% .time - Time vector [T x 1] in seconds
% .fs - Sampling frequency [scalar] in Hz
% .fchMask - Channel mask [1 x nChannels] all ones (good)
-% .markers - Event markers [M x 3] from hrfOnsets if provided
-% Format: [time, value, duration]
+% .markers - Event markers [M x 4] from hrfOnsets if provided
+% Format: [time, value, duration, amplitude]
% .info - Metadata struct with generation parameters
% .HbO - Oxygenated hemoglobin [T x nChannels] in uM
% .HbR - Deoxygenated hemoglobin [T x nChannels] in uM
@@ -208,7 +208,7 @@
% Create markers from onsets
if ~isempty(validOnsets)
- markers = [validOnsets(:), ones(length(validOnsets), 1), zeros(length(validOnsets), 1)];
+ markers = [validOnsets(:), ones(length(validOnsets), 1), zeros(length(validOnsets), 1), ones(length(validOnsets), 1)];
end
end
diff --git a/+pf2_base/+tests/+unit/BlockDefinitionTest.m b/+pf2_base/+tests/+unit/BlockDefinitionTest.m
new file mode 100644
index 00000000..90ab64f9
--- /dev/null
+++ b/+pf2_base/+tests/+unit/BlockDefinitionTest.m
@@ -0,0 +1,560 @@
+classdef BlockDefinitionTest < matlab.unittest.TestCase
+ % BLOCKDEFINITIONTEST Unit tests for pf2.data.defineBlocks and pf2.data.extractBlocks
+ %
+ % Tests block definition from markers (three modes) and block extraction
+ % to cell arrays. Covers ConditionMap, InfoTable, InfoFields, duration
+ % filtering, baseline subtraction, time shifting, and info merging.
+ %
+ % Run all tests:
+ % results = runtests('pf2_base.tests.unit.BlockDefinitionTest');
+ %
+ % See also: pf2.data.defineBlocks, pf2.data.extractBlocks
+
+ properties
+ processedData % Processed fNIRS sample data
+ end
+
+ methods (TestClassSetup)
+ function loadSampleData(testCase)
+ raw = pf2.import.sampleData.fNIR2000();
+ testCase.processedData = processFNIRS2(raw, 'ShowGUI', false);
+ end
+ end
+
+ methods (Static)
+ function data = makeSyntheticData()
+ % Create minimal fNIRS struct with synthetic markers for testing
+ fs = 10;
+ T = 300; % 300 seconds
+ nSamples = T * fs;
+ nCh = 4;
+ t = (0:nSamples-1)' / fs;
+
+ data.time = t;
+ data.fs = fs;
+ data.HbO = randn(nSamples, nCh) * 0.01;
+ data.HbR = randn(nSamples, nCh) * 0.005;
+ data.HbDiff = data.HbO - data.HbR;
+ data.HbTotal = data.HbO + data.HbR;
+ data.CBSI = randn(nSamples, nCh) * 0.008;
+ data.fchMask = ones(1, nCh);
+ data.info = struct('SubjectID', 'S01', 'Task', 'Test');
+
+ % Markers: [time, code, duration, amplitude]
+ % 3 blocks of marker 49 at t=30, 90, 180 with durations 20, 25, 30
+ data.markers = [
+ 30, 49, 20, 1;
+ 90, 49, 25, 1;
+ 180, 49, 30, 1;
+ ];
+ end
+
+ function data = makePairedMarkerData()
+ % Create data with start/end marker pairs
+ fs = 10;
+ T = 300;
+ nSamples = T * fs;
+ nCh = 4;
+ t = (0:nSamples-1)' / fs;
+
+ data.time = t;
+ data.fs = fs;
+ data.HbO = randn(nSamples, nCh) * 0.01;
+ data.HbR = randn(nSamples, nCh) * 0.005;
+ data.HbDiff = data.HbO - data.HbR;
+ data.HbTotal = data.HbO + data.HbR;
+ data.CBSI = randn(nSamples, nCh) * 0.008;
+ data.fchMask = ones(1, nCh);
+ data.info = struct('SubjectID', 'S01');
+
+ % Start/end pairs: condition A (50->52), condition B (51->53)
+ data.markers = [
+ 20, 50, 0, 1; % Start A block 1
+ 45, 52, 0, 1; % End A block 1
+ 60, 51, 0, 1; % Start B block 1
+ 90, 53, 0, 1; % End B block 1
+ 120, 50, 0, 1; % Start A block 2
+ 155, 52, 0, 1; % End A block 2
+ 200, 51, 0, 1; % Start B block 2
+ 240, 53, 0, 1; % End B block 2
+ ];
+ end
+ end
+
+ %% defineBlocks - Marker + Fixed Duration
+ methods (Test)
+ function testMarkerFixedDuration(testCase)
+ % Marker code + fixed duration creates correct blocks
+ data = pf2_base.tests.unit.BlockDefinitionTest.makeSyntheticData();
+ blocks = pf2.data.defineBlocks(data, 'MarkerCode', 49, 'Duration', 15);
+
+ testCase.verifyLength(blocks, 3, 'Should find 3 blocks');
+ testCase.verifyEqual(blocks(1).startTime, 30, 'First block starts at marker time');
+ testCase.verifyEqual(blocks(1).endTime, 45, 'First block ends at start + duration');
+ testCase.verifyEqual(blocks(1).duration, 15, 'Duration should be 15s');
+ testCase.verifyEqual(blocks(1).markerCode, 49, 'Marker code preserved');
+ end
+
+ function testMarkerFixedDurationBlockNumbers(testCase)
+ % BlockNumber auto-assigned sequentially
+ data = pf2_base.tests.unit.BlockDefinitionTest.makeSyntheticData();
+ blocks = pf2.data.defineBlocks(data, 'MarkerCode', 49, 'Duration', 15);
+
+ for k = 1:length(blocks)
+ testCase.verifyEqual(blocks(k).info.BlockNumber, k, ...
+ sprintf('Block %d should have BlockNumber = %d', k, k));
+ end
+ end
+ end
+
+ %% defineBlocks - Marker Duration from Column 3
+ methods (Test)
+ function testMarkerUseDuration(testCase)
+ % UseDuration flag reads duration from markers column 3
+ data = pf2_base.tests.unit.BlockDefinitionTest.makeSyntheticData();
+ blocks = pf2.data.defineBlocks(data, 'MarkerCode', 49, 'UseDuration', true);
+
+ testCase.verifyLength(blocks, 3);
+ testCase.verifyEqual(blocks(1).duration, 20, 'Block 1 duration from col 3');
+ testCase.verifyEqual(blocks(2).duration, 25, 'Block 2 duration from col 3');
+ testCase.verifyEqual(blocks(3).duration, 30, 'Block 3 duration from col 3');
+ testCase.verifyEqual(blocks(1).endTime, 50, 'endTime = startTime + duration');
+ end
+ end
+
+ %% defineBlocks - Start/End Pairs
+ methods (Test)
+ function testStartEndPairs(testCase)
+ % Paired start/end markers produce correct windows
+ data = pf2_base.tests.unit.BlockDefinitionTest.makePairedMarkerData();
+ blocks = pf2.data.defineBlocks(data, ...
+ 'StartMarker', [50; 51], 'EndMarker', [52; 53]);
+
+ testCase.verifyLength(blocks, 4, 'Should find 4 blocks (2 per condition)');
+
+ % Sorted by time: A1(20-45), B1(60-90), A2(120-155), B2(200-240)
+ testCase.verifyEqual(blocks(1).startTime, 20);
+ testCase.verifyEqual(blocks(1).endTime, 45);
+ testCase.verifyEqual(blocks(1).duration, 25);
+ testCase.verifyEqual(blocks(2).startTime, 60);
+ testCase.verifyEqual(blocks(2).endTime, 90);
+ end
+ end
+
+ %% defineBlocks - ConditionMap
+ methods (Test)
+ function testConditionMap(testCase)
+ % ConditionMap assigns condition labels per marker code
+ data = pf2_base.tests.unit.BlockDefinitionTest.makePairedMarkerData();
+ blocks = pf2.data.defineBlocks(data, ...
+ 'StartMarker', [50; 51], 'EndMarker', [52; 53], ...
+ 'ConditionMap', {50, 'Natural'; 51, 'Synthetic'});
+
+ % Sorted by time: A1(50), B1(51), A2(50), B2(51)
+ testCase.verifyEqual(blocks(1).info.Condition, 'Natural');
+ testCase.verifyEqual(blocks(2).info.Condition, 'Synthetic');
+ testCase.verifyEqual(blocks(3).info.Condition, 'Natural');
+ testCase.verifyEqual(blocks(4).info.Condition, 'Synthetic');
+ end
+
+ function testConditionMapWithMarkerCode(testCase)
+ % ConditionMap works with MarkerCode mode using OR logic
+ data = pf2_base.tests.unit.BlockDefinitionTest.makePairedMarkerData();
+
+ % Use column vector for OR: find markers 50 or 51
+ blocks = pf2.data.defineBlocks(data, ...
+ 'MarkerCode', [50; 51], 'Duration', 10, ...
+ 'ConditionMap', {50, 'CondA'; 51, 'CondB'});
+
+ testCase.verifyLength(blocks, 4, 'Should find 4 markers (2x50 + 2x51)');
+ testCase.verifyEqual(blocks(1).info.Condition, 'CondA');
+ testCase.verifyEqual(blocks(2).info.Condition, 'CondB');
+ end
+ end
+
+ %% defineBlocks - InfoTable
+ methods (Test)
+ function testInfoTable(testCase)
+ % Per-block table columns become .info fields
+ data = pf2_base.tests.unit.BlockDefinitionTest.makeSyntheticData();
+ scores = table([85; 92; 78], {'Easy';'Hard';'Easy'}, ...
+ 'VariableNames', {'Score','Difficulty'});
+
+ blocks = pf2.data.defineBlocks(data, 'MarkerCode', 49, ...
+ 'Duration', 15, 'InfoTable', scores);
+
+ testCase.verifyEqual(blocks(1).info.Score, 85);
+ testCase.verifyEqual(blocks(2).info.Difficulty, 'Hard');
+ testCase.verifyEqual(blocks(3).info.Score, 78);
+ end
+ end
+
+ %% defineBlocks - InfoFields
+ methods (Test)
+ function testInfoFields(testCase)
+ % Constant InfoFields applied to all blocks
+ data = pf2_base.tests.unit.BlockDefinitionTest.makeSyntheticData();
+ blocks = pf2.data.defineBlocks(data, 'MarkerCode', 49, ...
+ 'Duration', 15, 'InfoFields', struct('Task', 'Stroop', 'Group', 'Control'));
+
+ for k = 1:length(blocks)
+ testCase.verifyEqual(blocks(k).info.Task, 'Stroop');
+ testCase.verifyEqual(blocks(k).info.Group, 'Control');
+ end
+ end
+ end
+
+ %% defineBlocks - MinDuration Filter
+ methods (Test)
+ function testMinDurationFilter(testCase)
+ % Blocks shorter than MinDuration are rejected
+ data = pf2_base.tests.unit.BlockDefinitionTest.makeSyntheticData();
+ % Durations are 20, 25, 30 from column 3
+ blocks = pf2.data.defineBlocks(data, 'MarkerCode', 49, ...
+ 'UseDuration', true, 'MinDuration', 22);
+
+ testCase.verifyLength(blocks, 2, ...
+ 'Should reject block with duration=20 (< MinDuration=22)');
+ testCase.verifyEqual(blocks(1).duration, 25);
+ testCase.verifyEqual(blocks(2).duration, 30);
+ % BlockNumbers should be renumbered after filtering
+ testCase.verifyEqual(blocks(1).info.BlockNumber, 1);
+ testCase.verifyEqual(blocks(2).info.BlockNumber, 2);
+ end
+
+ function testMaxDurationFilter(testCase)
+ % Blocks longer than MaxDuration are rejected
+ data = pf2_base.tests.unit.BlockDefinitionTest.makeSyntheticData();
+ blocks = pf2.data.defineBlocks(data, 'MarkerCode', 49, ...
+ 'UseDuration', true, 'MaxDuration', 27);
+
+ testCase.verifyLength(blocks, 2, ...
+ 'Should reject block with duration=30 (> MaxDuration=27)');
+ end
+ end
+
+ %% defineBlocks - Empty Markers
+ methods (Test)
+ function testEmptyMarkers(testCase)
+ % No matching markers returns empty struct array
+ data = pf2_base.tests.unit.BlockDefinitionTest.makeSyntheticData();
+ blocks = pf2.data.defineBlocks(data, 'MarkerCode', 999, 'Duration', 10);
+
+ testCase.verifyEmpty(blocks, 'Should return empty for no matching markers');
+ testCase.verifyTrue(isstruct(blocks), 'Empty result should still be struct');
+ end
+
+ function testEmptyMarkerArray(testCase)
+ % Data with empty markers returns empty
+ data = pf2_base.tests.unit.BlockDefinitionTest.makeSyntheticData();
+ data.markers = zeros(0, 4);
+ blocks = pf2.data.defineBlocks(data, 'MarkerCode', 49, 'Duration', 10);
+
+ testCase.verifyEmpty(blocks);
+ end
+ end
+
+ %% extractBlocks - Basic Extraction
+ methods (Test)
+ function testBasicExtraction(testCase)
+ % Extract blocks produces correct number of segments
+ data = pf2_base.tests.unit.BlockDefinitionTest.makeSyntheticData();
+ blocks = pf2.data.defineBlocks(data, 'MarkerCode', 49, 'Duration', 20);
+ segments = pf2.data.extractBlocks(data, blocks, 'SetT0', false);
+
+ testCase.verifyLength(segments, 3, 'Should extract 3 segments');
+
+ % Each segment should have time within block bounds
+ for k = 1:length(segments)
+ seg = segments{k};
+ testCase.verifyTrue(isfield(seg, 'HbO'), 'Segment should have HbO');
+ testCase.verifyTrue(isfield(seg, 'time'), 'Segment should have time');
+ testCase.verifyGreaterThanOrEqual(min(seg.time), blocks(k).startTime - 0.5, ...
+ 'Segment time should start near block start');
+ testCase.verifyLessThanOrEqual(max(seg.time), blocks(k).endTime + 0.5, ...
+ 'Segment time should end near block end');
+ end
+ end
+ end
+
+ %% extractBlocks - PreTime/PostTime
+ methods (Test)
+ function testPreTimePostTime(testCase)
+ % PreTime and PostTime extend extraction window
+ data = pf2_base.tests.unit.BlockDefinitionTest.makeSyntheticData();
+ blocks = pf2.data.defineBlocks(data, 'MarkerCode', 49, 'Duration', 20);
+
+ % Extract with 5s pre and 3s post, no T0 shift for easier verification
+ segments = pf2.data.extractBlocks(data, blocks, ...
+ 'PreTime', 5, 'PostTime', 3, 'SetT0', false);
+
+ seg = segments{1};
+ expectedStart = blocks(1).startTime - 5;
+ expectedEnd = blocks(1).endTime + 3;
+
+ testCase.verifyLessThanOrEqual(min(seg.time), expectedStart + 0.2, ...
+ 'Segment should start at block start - PreTime');
+ testCase.verifyGreaterThanOrEqual(max(seg.time), expectedEnd - 0.2, ...
+ 'Segment should end at block end + PostTime');
+ end
+ end
+
+ %% extractBlocks - BaselineWindow
+ methods (Test)
+ function testBaselineWindow(testCase)
+ % BaselineWindow applies baseline subtraction
+ data = pf2_base.tests.unit.BlockDefinitionTest.makeSyntheticData();
+ % Add a known offset to HbO so baseline subtraction is visible
+ data.HbO = data.HbO + 5;
+
+ blocks = pf2.data.defineBlocks(data, 'MarkerCode', 49, 'Duration', 20);
+
+ % Extract without baseline
+ segsNoBL = pf2.data.extractBlocks(data, blocks, 'SetT0', false);
+
+ % Extract with baseline: [-5, 0] relative to block start
+ % Need PreTime >= 5 to include baseline period in extraction
+ segsBL = pf2.data.extractBlocks(data, blocks, ...
+ 'PreTime', 5, 'BaselineWindow', [-5, 0], 'SetT0', false);
+
+ % Baseline-corrected data should have lower mean than uncorrected
+ meanNoBL = mean(segsNoBL{1}.HbO(:));
+ meanBL = mean(segsBL{1}.HbO(:));
+ testCase.verifyLessThan(abs(meanBL), abs(meanNoBL), ...
+ 'Baseline correction should reduce mean offset');
+ end
+ end
+
+ %% extractBlocks - SetT0
+ methods (Test)
+ function testSetT0(testCase)
+ % SetT0 shifts time so block start = 0
+ data = pf2_base.tests.unit.BlockDefinitionTest.makeSyntheticData();
+ blocks = pf2.data.defineBlocks(data, 'MarkerCode', 49, 'Duration', 20);
+
+ segments = pf2.data.extractBlocks(data, blocks, 'SetT0', true);
+
+ for k = 1:length(segments)
+ seg = segments{k};
+ testCase.verifyEqual(min(seg.time), 0, 'AbsTol', 0.2, ...
+ 'Time should start near 0 after SetT0');
+ end
+ end
+
+ function testSetT0WithPreTime(testCase)
+ % With PreTime, time starts at negative value
+ data = pf2_base.tests.unit.BlockDefinitionTest.makeSyntheticData();
+ blocks = pf2.data.defineBlocks(data, 'MarkerCode', 49, 'Duration', 20);
+
+ segments = pf2.data.extractBlocks(data, blocks, ...
+ 'PreTime', 5, 'SetT0', true);
+
+ seg = segments{1};
+ testCase.verifyLessThan(min(seg.time), 0, ...
+ 'With PreTime, time should start before 0');
+ testCase.verifyEqual(min(seg.time), -5, 'AbsTol', 0.2, ...
+ 'Time should start near -PreTime');
+ end
+ end
+
+ %% extractBlocks - Info Merging
+ methods (Test)
+ function testInfoMerging(testCase)
+ % Parent data.info merged with block.info
+ data = pf2_base.tests.unit.BlockDefinitionTest.makeSyntheticData();
+ blocks = pf2.data.defineBlocks(data, 'MarkerCode', 49, ...
+ 'Duration', 20, 'ConditionMap', {49, 'StroopTask'});
+
+ segments = pf2.data.extractBlocks(data, blocks);
+
+ seg = segments{1};
+ % Parent info fields preserved
+ testCase.verifyEqual(seg.info.SubjectID, 'S01', ...
+ 'Parent SubjectID should be copied');
+ testCase.verifyEqual(seg.info.Task, 'Test', ...
+ 'Parent Task field should be preserved');
+ % Block info fields overlaid
+ testCase.verifyTrue(isfield(seg.info, 'BlockNumber'), ...
+ 'BlockNumber should be present');
+ testCase.verifyEqual(seg.info.Condition, 'StroopTask', ...
+ 'Condition from ConditionMap should be present');
+ end
+
+ function testCopyInfoFalse(testCase)
+ % CopyInfo=false skips parent data.info
+ data = pf2_base.tests.unit.BlockDefinitionTest.makeSyntheticData();
+ blocks = pf2.data.defineBlocks(data, 'MarkerCode', 49, 'Duration', 20);
+
+ segments = pf2.data.extractBlocks(data, blocks, 'CopyInfo', false);
+
+ seg = segments{1};
+ testCase.verifyFalse(isfield(seg.info, 'SubjectID'), ...
+ 'Parent SubjectID should not be present when CopyInfo=false');
+ testCase.verifyTrue(isfield(seg.info, 'BlockNumber'), ...
+ 'Block-level fields should still be present');
+ end
+ end
+
+ %% extractBlocks - SkipInvalid
+ methods (Test)
+ function testSkipInvalid(testCase)
+ % Blocks outside data range are skipped
+ data = pf2_base.tests.unit.BlockDefinitionTest.makeSyntheticData();
+
+ % Create blocks manually, one out of range
+ blocks(1).startTime = 30;
+ blocks(1).endTime = 50;
+ blocks(1).duration = 20;
+ blocks(1).markerCode = 49;
+ blocks(1).markerIndex = 1;
+ blocks(1).info = struct('BlockNumber', 1);
+
+ blocks(2).startTime = 500; % Beyond data range (300s)
+ blocks(2).endTime = 520;
+ blocks(2).duration = 20;
+ blocks(2).markerCode = 49;
+ blocks(2).markerIndex = 2;
+ blocks(2).info = struct('BlockNumber', 2);
+
+ segments = pf2.data.extractBlocks(data, blocks, 'SkipInvalid', true);
+
+ testCase.verifyLength(segments, 1, ...
+ 'Out-of-range block should be skipped');
+ end
+ end
+
+ %% Integration: defineBlocks -> extractBlocks -> Experiment
+ methods (Test)
+ function testEndToEndPipeline(testCase)
+ % Full pipeline: define blocks, extract, verify cell array structure
+ data = pf2_base.tests.unit.BlockDefinitionTest.makeSyntheticData();
+
+ blocks = pf2.data.defineBlocks(data, 'MarkerCode', 49, ...
+ 'UseDuration', true, ...
+ 'ConditionMap', {49, 'StroopTask'}, ...
+ 'InfoFields', struct('Group', 'Control'));
+
+ segments = pf2.data.extractBlocks(data, blocks, ...
+ 'PreTime', 5, 'PostTime', 2, 'SetT0', true);
+
+ testCase.verifyLength(segments, 3, 'Should have 3 segments');
+
+ % Verify each segment is a valid fNIRS struct
+ for k = 1:length(segments)
+ seg = segments{k};
+ testCase.verifyTrue(isfield(seg, 'HbO'), 'Must have HbO');
+ testCase.verifyTrue(isfield(seg, 'time'), 'Must have time');
+ testCase.verifyTrue(isfield(seg, 'info'), 'Must have info');
+ testCase.verifyEqual(seg.info.Condition, 'StroopTask');
+ testCase.verifyEqual(seg.info.Group, 'Control');
+ testCase.verifyEqual(seg.info.SubjectID, 'S01');
+ testCase.verifyEqual(seg.info.BlockNumber, k);
+
+ % Time should start near -5 (PreTime) after SetT0
+ testCase.verifyEqual(min(seg.time), -5, 'AbsTol', 0.2);
+ end
+ end
+
+ function testWithRealSampleData(testCase)
+ % Verify with real processed sample data
+ data = testCase.processedData;
+
+ % Add synthetic markers to real data
+ timeVec = data.time;
+ minT = min(timeVec);
+ data.markers = [
+ minT + 30, 10, 20, 1;
+ minT + 80, 10, 20, 1;
+ ];
+
+ blocks = pf2.data.defineBlocks(data, 'MarkerCode', 10, 'Duration', 20);
+ segments = pf2.data.extractBlocks(data, blocks, 'SetT0', true);
+
+ testCase.verifyLength(segments, 2);
+ testCase.verifyTrue(isfield(segments{1}, 'HbO'));
+ testCase.verifyEqual(size(segments{1}.HbO, 2), size(data.HbO, 2), ...
+ 'Channel count should be preserved');
+ end
+ end
+
+ %% defineBlocks - Positional API
+ methods (Test)
+ function testPositionalCodesAndDuration(testCase)
+ % Simple positional syntax: defineBlocks(data, codes, duration)
+ data = pf2_base.tests.unit.BlockDefinitionTest.makeSyntheticData();
+ blocks = pf2.data.defineBlocks(data, [49], 15);
+
+ testCase.verifyLength(blocks, 3, 'Should find 3 blocks');
+ testCase.verifyEqual(blocks(1).duration, 15);
+ testCase.verifyEqual(blocks(1).startTime, 30);
+ end
+
+ function testPositionalMultipleCodes(testCase)
+ % Multiple codes as row vector: defineBlocks(data, [49, 50], 30)
+ data = pf2_base.tests.unit.BlockDefinitionTest.makePairedMarkerData();
+ blocks = pf2.data.defineBlocks(data, [50, 51], 10);
+
+ testCase.verifyLength(blocks, 4, 'Should find 4 blocks (2x50 + 2x51)');
+ end
+
+ function testPositionalAutoDuration(testCase)
+ % No duration given: auto-detect from marker column 3
+ data = pf2_base.tests.unit.BlockDefinitionTest.makeSyntheticData();
+ blocks = pf2.data.defineBlocks(data, 49);
+
+ testCase.verifyLength(blocks, 3);
+ testCase.verifyEqual(blocks(1).duration, 20, ...
+ 'Should auto-use duration from marker column 3');
+ testCase.verifyEqual(blocks(2).duration, 25);
+ testCase.verifyEqual(blocks(3).duration, 30);
+ end
+
+ function testPositionalWithNameValue(testCase)
+ % Positional codes + name-value params mixed
+ data = pf2_base.tests.unit.BlockDefinitionTest.makeSyntheticData();
+ blocks = pf2.data.defineBlocks(data, 49, 15, ...
+ 'ConditionMap', {49, 'Stroop'});
+
+ testCase.verifyLength(blocks, 3);
+ testCase.verifyEqual(blocks(1).info.Condition, 'Stroop');
+ testCase.verifyEqual(blocks(1).duration, 15);
+ end
+ end
+
+ %% defineBlocks - PrePad/PostPad
+ methods (Test)
+ function testPrePadPostPad(testCase)
+ % PrePad and PostPad extend block boundaries
+ data = pf2_base.tests.unit.BlockDefinitionTest.makeSyntheticData();
+ blocks = pf2.data.defineBlocks(data, 49, 20, ...
+ 'PrePad', 5, 'PostPad', 3);
+
+ testCase.verifyEqual(blocks(1).startTime, 25, ...
+ 'startTime should be 30 - 5 PrePad = 25');
+ testCase.verifyEqual(blocks(1).endTime, 53, ...
+ 'endTime should be 50 + 3 PostPad = 53');
+ testCase.verifyEqual(blocks(1).duration, 28, ...
+ 'duration = 20 + 5 + 3 = 28');
+ end
+ end
+
+ %% Error Handling
+ methods (Test)
+ function testNoModeError(testCase)
+ % Must specify either MarkerCode or StartMarker+EndMarker
+ data = pf2_base.tests.unit.BlockDefinitionTest.makeSyntheticData();
+ testCase.verifyError(...
+ @() pf2.data.defineBlocks(data, 'Duration', 10), ...
+ 'pf2:defineBlocks:noMode');
+ end
+
+ function testAmbiguousModeError(testCase)
+ % Cannot specify both MarkerCode and StartMarker
+ data = pf2_base.tests.unit.BlockDefinitionTest.makeSyntheticData();
+ testCase.verifyError(...
+ @() pf2.data.defineBlocks(data, 'MarkerCode', 49, ...
+ 'StartMarker', 50, 'EndMarker', 51), ...
+ 'pf2:defineBlocks:ambiguousMode');
+ end
+ end
+end
diff --git a/+pf2_base/+tests/+unit/ConnectivityTest.m b/+pf2_base/+tests/+unit/ConnectivityTest.m
new file mode 100644
index 00000000..6b534621
--- /dev/null
+++ b/+pf2_base/+tests/+unit/ConnectivityTest.m
@@ -0,0 +1,719 @@
+classdef ConnectivityTest < matlab.unittest.TestCase
+ % CONNECTIVITYTEST Unit tests for connectivity and hyperscanning modules
+ %
+ % Tests cover:
+ % - Coupling functions: pearson, spearman, xcorr, coherence, wcoherence
+ % - Connectivity matrix computation
+ % - Subject pairing for hyperscanning
+ % - Dyad and group computation
+ % - Permutation testing
+ % - Export to table
+ % - Plotting: plotWcoherence, plotWindowed, plotGroup
+ %
+ % Example:
+ % results = runtests('pf2_base.tests.unit.ConnectivityTest');
+ % disp(results);
+
+ properties
+ fs % Sampling frequency
+ T % Number of samples
+ nChannels % Number of channels
+ end
+
+ methods (TestClassSetup)
+ function setupParams(testCase)
+ testCase.fs = 10;
+ testCase.T = 1000;
+ testCase.nChannels = 8;
+ end
+ end
+
+
+ %% Coupling Functions
+ methods (Test)
+
+ function testPearsonCorrelatedSignals(testCase)
+ % Two correlated signals should have r near target
+ rng(42);
+ target_r = 0.7;
+ [x, y] = generateCorrelatedPair(testCase.T, target_r);
+ result = exploreFNIRS.coupling.pearson(x, y, testCase.fs);
+
+ testCase.verifyEqual(result.method, 'pearson');
+ testCase.verifyFalse(result.windowed);
+ testCase.verifyGreaterThan(result.value, 0.5);
+ testCase.verifyLessThan(result.pvalue, 0.01);
+ end
+
+ function testPearsonUncorrelatedSignals(testCase)
+ % Two independent signals should have r near zero
+ rng(42);
+ x = randn(testCase.T, 1);
+ y = randn(testCase.T, 1);
+ result = exploreFNIRS.coupling.pearson(x, y, testCase.fs);
+
+ testCase.verifyLessThan(abs(result.value), 0.15);
+ end
+
+ function testPearsonWindowed(testCase)
+ % Windowed: high corr in first half, low in second
+ rng(42);
+ T = testCase.T;
+ shared = randn(T, 1);
+ x = [shared(1:T/2) + randn(T/2, 1) * 0.3; randn(T/2, 1)];
+ y = [shared(1:T/2) + randn(T/2, 1) * 0.3; randn(T/2, 1)];
+
+ result = exploreFNIRS.coupling.pearson(x, y, testCase.fs, ...
+ 'WindowSize', 10); % 10-second windows
+
+ testCase.verifyTrue(result.windowed);
+ testCase.verifyGreaterThan(length(result.value), 1);
+ testCase.verifyTrue(isfield(result, 'windowTimes'));
+
+ % First quarter should have higher correlation than last quarter
+ nWin = length(result.value);
+ firstQ = mean(result.value(1:floor(nWin/4)), 'omitnan');
+ lastQ = mean(result.value(ceil(3*nWin/4):end), 'omitnan');
+ testCase.verifyGreaterThan(firstQ, lastQ);
+ end
+
+ function testSpearmanCorrelatedSignals(testCase)
+ rng(42);
+ [x, y] = generateCorrelatedPair(testCase.T, 0.7);
+ result = exploreFNIRS.coupling.spearman(x, y, testCase.fs);
+
+ testCase.verifyEqual(result.method, 'spearman');
+ testCase.verifyGreaterThan(result.value, 0.4);
+ testCase.verifyLessThan(result.pvalue, 0.01);
+ end
+
+ function testXcorrPeakLag(testCase)
+ % Signal y is a delayed version of x - should detect lag
+ rng(42);
+ lagSamples = 5; % 0.5 seconds at 10 Hz
+ T = testCase.T;
+ x = randn(T, 1);
+ y = [zeros(lagSamples, 1); x(1:end-lagSamples)];
+
+ result = exploreFNIRS.coupling.xcorr(x, y, testCase.fs, ...
+ 'MaxLag', 2);
+
+ testCase.verifyEqual(result.method, 'xcorr');
+ testCase.verifyGreaterThan(abs(result.value), 0.8);
+ testCase.verifyEqual(abs(result.lag), lagSamples / testCase.fs, ...
+ 'AbsTol', 1/testCase.fs);
+ end
+
+ function testCoherenceCorrelatedSignals(testCase)
+ % Low-frequency coherent signals
+ rng(42);
+ t = (0:testCase.T-1)' / testCase.fs;
+ shared = sin(2*pi*0.05*t); % 0.05 Hz shared oscillation
+ x = shared + randn(testCase.T, 1) * 0.3;
+ y = shared + randn(testCase.T, 1) * 0.3;
+
+ result = exploreFNIRS.coupling.coherence(x, y, testCase.fs, ...
+ 'FreqRange', [0.01, 0.1]);
+
+ testCase.verifyEqual(result.method, 'coherence');
+ testCase.verifyGreaterThan(result.value, 0.3);
+ testCase.verifyTrue(isfield(result, 'spectrum'));
+ testCase.verifyTrue(isfield(result, 'freqs'));
+ end
+
+ function testPearsonLengthMismatch(testCase)
+ x = randn(100, 1);
+ y = randn(50, 1);
+ testCase.verifyError(@() exploreFNIRS.coupling.pearson(x, y, 10), ...
+ 'exploreFNIRS:coupling:pearson');
+ end
+
+ end
+
+
+ %% Connectivity Matrix
+ methods (Test)
+
+ function testComputeMatrixBasic(testCase)
+ rng(42);
+ data = createSyntheticSubject(testCase, 'correlated');
+
+ result = exploreFNIRS.connectivity.computeMatrix(data, ...
+ 'Method', 'pearson', 'Biomarker', 'HbO');
+
+ nCh = length(result.channels);
+ testCase.verifyEqual(size(result.matrix), [nCh, nCh]);
+ testCase.verifyEqual(size(result.pmatrix), [nCh, nCh]);
+ testCase.verifyEqual(result.method, 'pearson');
+ testCase.verifyEqual(result.biomarker, 'HbO');
+
+ % Diagonal should be 1
+ testCase.verifyEqual(diag(result.matrix), ones(nCh, 1), 'AbsTol', 1e-10);
+
+ % Matrix should be symmetric
+ testCase.verifyEqual(result.matrix, result.matrix', 'AbsTol', 1e-10);
+ end
+
+ function testComputeMatrixCorrelation(testCase)
+ % Channels 1-2 should be correlated, channel 3 independent
+ rng(42);
+ T = testCase.T;
+ shared = randn(T, 1);
+ data.HbO = [shared + randn(T,1)*0.2, shared + randn(T,1)*0.2, randn(T,1)];
+ data.time = (0:T-1)' / testCase.fs;
+ data.fs = testCase.fs;
+ data.fchMask = [1, 1, 1];
+
+ result = exploreFNIRS.connectivity.computeMatrix(data, ...
+ 'Method', 'pearson', 'Biomarker', 'HbO');
+
+ % Ch1-Ch2 should be highly correlated
+ testCase.verifyGreaterThan(result.matrix(1, 2), 0.7);
+ % Ch1-Ch3 and Ch2-Ch3 should be weakly correlated
+ testCase.verifyLessThan(abs(result.matrix(1, 3)), 0.3);
+ testCase.verifyLessThan(abs(result.matrix(2, 3)), 0.3);
+ end
+
+ function testComputeMatrixTimeWindow(testCase)
+ rng(42);
+ data = createSyntheticSubject(testCase, 'random');
+ result = exploreFNIRS.connectivity.computeMatrix(data, ...
+ 'TimeWindow', [10, 50]);
+
+ expectedSamples = sum(data.time >= 10 & data.time <= 50);
+ testCase.verifyEqual(result.nSamples, expectedSamples);
+ end
+
+ function testComputeMatrixChannelSubset(testCase)
+ rng(42);
+ data = createSyntheticSubject(testCase, 'random');
+ result = exploreFNIRS.connectivity.computeMatrix(data, ...
+ 'Channels', [1, 3, 5]);
+
+ testCase.verifyEqual(result.channels, [1, 3, 5]);
+ testCase.verifyEqual(size(result.matrix), [3, 3]);
+ end
+
+ end
+
+
+ %% Hyperscanning - Pairing
+ methods (Test)
+
+ function testPairSubjectsByMetadata(testCase)
+ data = createDyadData(testCase, 3);
+ pairs = exploreFNIRS.hyperscanning.pairSubjects(data);
+
+ testCase.verifyEqual(length(pairs), 3);
+ for d = 1:3
+ testCase.verifyEqual(length(pairs(d).indices), 2);
+ testCase.verifyNotEmpty(pairs(d).dyadID);
+ end
+ end
+
+ function testPairSubjectsManual(testCase)
+ data = createDyadData(testCase, 2);
+ pairs = exploreFNIRS.hyperscanning.pairSubjects(data, ...
+ 'ManualPairs', {{1,2}, {3,4}});
+
+ testCase.verifyEqual(length(pairs), 2);
+ testCase.verifyEqual(pairs(1).indices, [1, 2]);
+ testCase.verifyEqual(pairs(2).indices, [3, 4]);
+ end
+
+ function testPairSubjectsMismatchWarning(testCase)
+ % Create data with one incomplete dyad
+ data = createDyadData(testCase, 2);
+ % Remove one member of dyad 2
+ data = data(1:3); % only 3 of 4 subjects
+ testCase.verifyWarning( ...
+ @() exploreFNIRS.hyperscanning.pairSubjects(data), ...
+ 'exploreFNIRS:hyperscanning:pairSubjects');
+ end
+
+ end
+
+
+ %% Hyperscanning - Dyad Computation
+ methods (Test)
+
+ function testComputeDyadIdentical(testCase)
+ % Identical signals should yield coupling ~1
+ rng(42);
+ data = createSyntheticSubject(testCase, 'random');
+ result = exploreFNIRS.hyperscanning.computeDyad(data, data, ...
+ 'Method', 'pearson', 'Biomarker', 'HbO', 'ChannelPairing', 'same');
+
+ testCase.verifyEqual(result.pairing, 'same');
+ testCase.verifyEqual(length(result.values), testCase.nChannels);
+ % All channels should have r ~ 1
+ testCase.verifyGreaterThan(min(result.values), 0.95);
+ end
+
+ function testComputeDyadCorrelated(testCase)
+ % Shared signal across subjects should yield positive coupling
+ rng(42);
+ T = testCase.T;
+ nCh = testCase.nChannels;
+ shared = randn(T, nCh) * 0.5;
+
+ dataA = createSyntheticSubject(testCase, 'random');
+ dataB = createSyntheticSubject(testCase, 'random');
+ dataA.HbO = shared + randn(T, nCh) * 0.5;
+ dataB.HbO = shared + randn(T, nCh) * 0.5;
+
+ result = exploreFNIRS.hyperscanning.computeDyad(dataA, dataB, ...
+ 'Method', 'pearson', 'ChannelPairing', 'same');
+
+ testCase.verifyGreaterThan(mean(result.values, 'omitnan'), 0.2);
+ end
+
+ function testComputeDyadAllPairing(testCase)
+ rng(42);
+ dataA = createSyntheticSubject(testCase, 'random');
+ dataB = createSyntheticSubject(testCase, 'random');
+
+ result = exploreFNIRS.hyperscanning.computeDyad(dataA, dataB, ...
+ 'ChannelPairing', 'all');
+
+ testCase.verifyEqual(result.pairing, 'all');
+ nCh = testCase.nChannels;
+ testCase.verifyEqual(size(result.values), [nCh, nCh]);
+ end
+
+ end
+
+
+ %% Hyperscanning - Group Computation
+ methods (Test)
+
+ function testComputeGroupBasic(testCase)
+ rng(42);
+ nDyads = 5;
+ data = createDyadData(testCase, nDyads);
+ pairs = exploreFNIRS.hyperscanning.pairSubjects(data);
+
+ result = exploreFNIRS.hyperscanning.computeGroup(data, pairs, ...
+ 'Method', 'pearson', 'Biomarker', 'HbO');
+
+ testCase.verifyTrue(isfield(result, 'Mean'));
+ testCase.verifyTrue(isfield(result, 'SD'));
+ testCase.verifyTrue(isfield(result, 'SEM'));
+ testCase.verifyTrue(isfield(result, 'N'));
+ testCase.verifyTrue(isfield(result, 'tstat'));
+ testCase.verifyTrue(isfield(result, 'pvalue'));
+ testCase.verifyEqual(length(result.dyads), nDyads);
+ end
+
+ function testComputeGroupDetectsSignal(testCase)
+ % Dyads with shared signal should have positive mean coupling
+ rng(42);
+ nDyads = 5;
+ data = createDyadDataWithSharedSignal(testCase, nDyads, 0.5);
+ pairs = exploreFNIRS.hyperscanning.pairSubjects(data);
+
+ result = exploreFNIRS.hyperscanning.computeGroup(data, pairs, ...
+ 'Method', 'pearson', 'Biomarker', 'HbO');
+
+ testCase.verifyGreaterThan(mean(result.Mean, 'omitnan'), 0.15);
+ end
+
+ end
+
+
+ %% Permutation Testing
+ methods (Test)
+
+ function testPermutationTestSignificant(testCase)
+ % Strong shared signal should survive permutation
+ rng(42);
+ nDyads = 5;
+ data = createDyadDataWithSharedSignal(testCase, nDyads, 0.7);
+ pairs = exploreFNIRS.hyperscanning.pairSubjects(data);
+
+ result = exploreFNIRS.hyperscanning.permutationTest(data, pairs, ...
+ 'Permutations', 100, 'PThreshold', 0.05, ...
+ 'Method', 'pearson', 'Biomarker', 'HbO');
+
+ testCase.verifyTrue(isfield(result, 'pvalue'));
+ testCase.verifyTrue(isfield(result, 'significant'));
+ testCase.verifyTrue(isfield(result, 'nullDist'));
+ testCase.verifyTrue(isfield(result, 'zScore'));
+
+ % At least some channels should be significant
+ testCase.verifyTrue(any(result.significant(:)));
+ end
+
+ function testPermutationTestNull(testCase)
+ % Independent signals should not be significant
+ rng(42);
+ nDyads = 5;
+ data = createDyadData(testCase, nDyads); % no shared signal
+ pairs = exploreFNIRS.hyperscanning.pairSubjects(data);
+
+ result = exploreFNIRS.hyperscanning.permutationTest(data, pairs, ...
+ 'Permutations', 50, 'PThreshold', 0.05, ...
+ 'Method', 'pearson', 'Biomarker', 'HbO');
+
+ % Most channels should NOT be significant
+ testCase.verifyLessThan(sum(result.significant(:)) / numel(result.significant), 0.3);
+ end
+
+ end
+
+
+ %% Export
+ methods (Test)
+
+ function testConnectivityMatrixExport(testCase)
+ rng(42);
+ data = createSyntheticSubject(testCase, 'random');
+ result = exploreFNIRS.connectivity.computeMatrix(data, ...
+ 'Method', 'pearson', 'Channels', 1:3);
+
+ T = exploreFNIRS.export.connectivityToTable(result);
+
+ testCase.verifyTrue(istable(T));
+ testCase.verifyGreaterThan(height(T), 0);
+ testCase.verifyTrue(ismember('Coupling', T.Properties.VariableNames));
+ testCase.verifyTrue(ismember('ChannelA', T.Properties.VariableNames));
+ testCase.verifyTrue(ismember('ChannelB', T.Properties.VariableNames));
+ end
+
+ function testHyperscanningExport(testCase)
+ rng(42);
+ nDyads = 3;
+ data = createDyadDataWithSharedSignal(testCase, nDyads, 0.5);
+ pairs = exploreFNIRS.hyperscanning.pairSubjects(data);
+ result = exploreFNIRS.hyperscanning.computeGroup(data, pairs, ...
+ 'Method', 'pearson', 'Biomarker', 'HbO');
+ result.pairs = pairs;
+
+ T = exploreFNIRS.export.connectivityToTable(result, ...
+ 'IncludeDyads', true, 'IncludeGroup', true);
+
+ testCase.verifyTrue(istable(T));
+ testCase.verifyGreaterThan(height(T), 0);
+ testCase.verifyTrue(ismember('DyadID', T.Properties.VariableNames));
+ testCase.verifyTrue(ismember('Level', T.Properties.VariableNames));
+ end
+
+
+ %% Wavelet Coherence Tests
+
+ function testWcoherenceCorrelatedSignals(testCase)
+ % Two correlated signals should have high wavelet coherence
+ rng(42);
+ T = testCase.T;
+ fs = testCase.fs;
+ t = (0:T-1)' / fs;
+ shared = sin(2 * pi * 0.05 * t);
+ x = shared + randn(T, 1) * 0.2;
+ y = shared + randn(T, 1) * 0.2;
+
+ result = exploreFNIRS.coupling.wcoherence(x, y, fs, ...
+ 'FreqRange', [0.01, 0.2]);
+
+ testCase.verifyEqual(result.method, 'wcoherence');
+ testCase.verifyFalse(result.windowed);
+ testCase.verifyGreaterThan(result.value, 0.3);
+ testCase.verifyTrue(isnan(result.pvalue));
+ testCase.verifyTrue(isfield(result, 'wcoh'));
+ testCase.verifyTrue(isfield(result, 'freqs'));
+ testCase.verifyTrue(isfield(result, 'times'));
+ testCase.verifyTrue(isfield(result, 'coi'));
+ testCase.verifyEqual(size(result.wcoh, 2), T);
+ end
+
+ function testWcoherenceUncorrelatedSignals(testCase)
+ % Two independent signals should have low wavelet coherence
+ rng(42);
+ x = randn(testCase.T, 1);
+ y = randn(testCase.T, 1);
+
+ result = exploreFNIRS.coupling.wcoherence(x, y, testCase.fs, ...
+ 'FreqRange', [0.01, 0.5]);
+
+ testCase.verifyLessThan(result.value, 0.6);
+ end
+
+ function testWcoherencePhaseOutput(testCase)
+ % Phase output should be returned when requested
+ rng(42);
+ T = testCase.T;
+ fs = testCase.fs;
+ t = (0:T-1)' / fs;
+ x = sin(2 * pi * 0.05 * t);
+ y = cos(2 * pi * 0.05 * t);
+
+ result = exploreFNIRS.coupling.wcoherence(x, y, fs, ...
+ 'PhaseOutput', true);
+
+ testCase.verifyTrue(isfield(result, 'phase'));
+ testCase.verifyEqual(size(result.phase), size(result.wcoh));
+ end
+
+ function testWcoherenceInConnectivityMatrix(testCase)
+ % WCT should work through computeMatrix dispatch
+ rng(42);
+ data = createSyntheticSubject(testCase, 'correlated');
+ result = exploreFNIRS.connectivity.computeMatrix(data, ...
+ 'Method', 'wcoherence', 'Biomarker', 'HbO');
+
+ testCase.verifyTrue(isfield(result, 'matrix'));
+ nCh = testCase.nChannels;
+ testCase.verifyEqual(size(result.matrix), [nCh, nCh]);
+ testCase.verifyEqual(result.method, 'wcoherence');
+ end
+
+ function testWcoherenceInHyperscanning(testCase)
+ % WCT should work through computeDyad dispatch
+ rng(42);
+ data = createDyadDataWithSharedSignal(testCase, 2, 0.5);
+ pairs = exploreFNIRS.hyperscanning.pairSubjects(data);
+ result = exploreFNIRS.hyperscanning.computeGroup(data, pairs, ...
+ 'Method', 'wcoherence', 'Biomarker', 'HbO');
+
+ testCase.verifyTrue(isfield(result, 'Mean'));
+ testCase.verifyEqual(result.method, 'wcoherence');
+ end
+
+
+ %% Plot Tests (headless)
+
+ function testPlotWcoherence(testCase)
+ % plotWcoherence should create a figure without error
+ rng(42);
+ T = testCase.T;
+ fs = testCase.fs;
+ t = (0:T-1)' / fs;
+ x = sin(2 * pi * 0.05 * t) + randn(T, 1) * 0.3;
+ y = sin(2 * pi * 0.05 * t) + randn(T, 1) * 0.3;
+
+ result = exploreFNIRS.coupling.wcoherence(x, y, fs);
+ fig = exploreFNIRS.coupling.plotWcoherence(result, ...
+ 'Visible', 'off');
+
+ testCase.verifyTrue(ishandle(fig));
+ close(fig);
+ end
+
+ function testPlotWindowed(testCase)
+ % plotWindowed should create a figure for windowed results
+ rng(42);
+ [x, y] = generateCorrelatedPair(testCase.T, 0.5);
+ result = exploreFNIRS.coupling.pearson(x, y, testCase.fs, ...
+ 'WindowSize', 30);
+
+ testCase.verifyTrue(result.windowed);
+ fig = exploreFNIRS.coupling.plotWindowed(result, ...
+ 'Visible', 'off');
+
+ testCase.verifyTrue(ishandle(fig));
+ close(fig);
+ end
+
+ function testPlotGroupHyperscanning(testCase)
+ % plotGroup should create a figure for group results
+ rng(42);
+ data = createDyadDataWithSharedSignal(testCase, 3, 0.5);
+ pairs = exploreFNIRS.hyperscanning.pairSubjects(data);
+ result = exploreFNIRS.hyperscanning.computeGroup(data, pairs, ...
+ 'Method', 'pearson', 'Biomarker', 'HbO');
+
+ fig = exploreFNIRS.hyperscanning.plotGroup(result, ...
+ 'Visible', 'off');
+
+ testCase.verifyTrue(ishandle(fig));
+ close(fig);
+ end
+
+ %% ROI Mode Tests
+
+ function testConnectivityMatrixWithROI(testCase)
+ % computeMatrix with UseROI should use ROI data and return ROI labels
+ rng(42);
+ data = createSyntheticSubject(testCase, 'correlated');
+ data = addROIData(data, testCase.T, {'Left','Center','Right'});
+
+ result = exploreFNIRS.connectivity.computeMatrix(data, ...
+ 'UseROI', true, 'Method', 'pearson');
+
+ testCase.verifyEqual(size(result.matrix), [3, 3]);
+ testCase.verifyTrue(result.useROI);
+ testCase.verifyEqual(result.labels, {'Left';'Center';'Right'});
+ testCase.verifyEqual(result.channels, 1:3);
+ % Diagonal should be 1
+ testCase.verifyEqual(diag(result.matrix), ones(3,1), 'AbsTol', 1e-10);
+ end
+
+ function testConnectivityMatrixROISubsetChannels(testCase)
+ % computeMatrix with UseROI and Channels should subset ROIs
+ rng(42);
+ data = createSyntheticSubject(testCase, 'random');
+ data = addROIData(data, testCase.T, {'Left','Center','Right','Back'});
+
+ result = exploreFNIRS.connectivity.computeMatrix(data, ...
+ 'UseROI', true, 'Channels', [1 3]);
+
+ testCase.verifyEqual(size(result.matrix), [2, 2]);
+ testCase.verifyEqual(result.channels, [1 3]);
+ testCase.verifyEqual(result.labels, {'Left';'Right'});
+ end
+
+ function testDyadWithROI(testCase)
+ % computeDyad with UseROI should pair ROIs between subjects
+ rng(42);
+ roiNames = {'Left','Center','Right'};
+ dataA = createSyntheticSubject(testCase, 'random');
+ dataA = addROIData(dataA, testCase.T, roiNames);
+ dataA.info.DyadID = 'D01';
+ dataA.info.Role = 'Speaker';
+
+ dataB = createSyntheticSubject(testCase, 'random');
+ dataB = addROIData(dataB, testCase.T, roiNames);
+ dataB.info.DyadID = 'D01';
+ dataB.info.Role = 'Listener';
+
+ result = exploreFNIRS.hyperscanning.computeDyad(dataA, dataB, ...
+ 'UseROI', true, 'Method', 'pearson');
+
+ testCase.verifyEqual(length(result.values), 3);
+ testCase.verifyTrue(result.useROI);
+ testCase.verifyEqual(result.labelsA, {'Left';'Center';'Right'});
+ testCase.verifyEqual(result.labelsB, {'Left';'Center';'Right'});
+ end
+
+ function testDyadROIMissingError(testCase)
+ % computeDyad with UseROI should error when ROI data is missing
+ dataA = createSyntheticSubject(testCase, 'random');
+ dataB = createSyntheticSubject(testCase, 'random');
+
+ testCase.verifyError( ...
+ @() exploreFNIRS.hyperscanning.computeDyad(dataA, dataB, ...
+ 'UseROI', true), ...
+ 'exploreFNIRS:hyperscanning:computeDyad');
+ end
+
+ function testPlotMatrixWithROILabels(testCase)
+ % plotMatrix should display ROI labels when result has labels
+ rng(42);
+ data = createSyntheticSubject(testCase, 'correlated');
+ data = addROIData(data, testCase.T, {'Left','Center','Right'});
+
+ result = exploreFNIRS.connectivity.computeMatrix(data, ...
+ 'UseROI', true, 'Method', 'pearson');
+
+ fig = exploreFNIRS.connectivity.plotMatrix(result, ...
+ 'Visible', 'off');
+
+ testCase.verifyTrue(ishandle(fig));
+ % Check axis labels contain ROI names
+ ax = findobj(fig, 'Type', 'Axes');
+ labels = get(ax, 'XTickLabel');
+ testCase.verifyEqual(labels, {'Left';'Center';'Right'});
+ close(fig);
+ end
+
+ function testConnectivityMatrixROIMissingError(testCase)
+ % computeMatrix with UseROI should error when ROI data is missing
+ data = createSyntheticSubject(testCase, 'random');
+
+ testCase.verifyError( ...
+ @() exploreFNIRS.connectivity.computeMatrix(data, 'UseROI', true), ...
+ 'exploreFNIRS:connectivity:computeMatrix');
+ end
+
+ end
+
+end
+
+
+%% Helper functions
+
+function [x, y] = generateCorrelatedPair(T, target_r)
+ % Generate two signals with approximate Pearson r = target_r
+ x = randn(T, 1);
+ noise = randn(T, 1);
+ y = target_r * x + sqrt(1 - target_r^2) * noise;
+end
+
+
+function data = createSyntheticSubject(testCase, mode)
+ % Create a single fNIRS-like struct
+ T = testCase.T;
+ nCh = testCase.nChannels;
+ fs = testCase.fs;
+
+ data.time = (0:T-1)' / fs;
+ data.fs = fs;
+ data.fchMask = ones(1, nCh);
+
+ switch mode
+ case 'random'
+ data.HbO = randn(T, nCh);
+ data.HbR = randn(T, nCh) * 0.5;
+ case 'correlated'
+ shared = randn(T, 1);
+ data.HbO = repmat(shared, 1, nCh) + randn(T, nCh) * 0.3;
+ data.HbR = -data.HbO * 0.3 + randn(T, nCh) * 0.1;
+ end
+
+ data.info.SubjectID = 'TestSubject';
+end
+
+
+function data = createDyadData(testCase, nDyads)
+ % Create cell array of subjects paired by DyadID (independent signals)
+ data = cell(nDyads * 2, 1);
+ for d = 1:nDyads
+ for role = 1:2
+ idx = (d-1)*2 + role;
+ s = createSyntheticSubject(testCase, 'random');
+ s.info.SubjectID = sprintf('S%d%d', d, role);
+ s.info.DyadID = sprintf('D%02d', d);
+ if role == 1
+ s.info.Role = 'Speaker';
+ else
+ s.info.Role = 'Listener';
+ end
+ data{idx} = s;
+ end
+ end
+end
+
+
+function data = createDyadDataWithSharedSignal(testCase, nDyads, strength)
+ % Create dyad data where partners share a signal component
+ T = testCase.T;
+ nCh = testCase.nChannels;
+ data = cell(nDyads * 2, 1);
+
+ for d = 1:nDyads
+ shared = randn(T, nCh) * strength;
+ for role = 1:2
+ idx = (d-1)*2 + role;
+ s = createSyntheticSubject(testCase, 'random');
+ s.HbO = shared + randn(T, nCh) * (1 - strength);
+ s.info.SubjectID = sprintf('S%d%d', d, role);
+ s.info.DyadID = sprintf('D%02d', d);
+ if role == 1
+ s.info.Role = 'Speaker';
+ else
+ s.info.Role = 'Listener';
+ end
+ data{idx} = s;
+ end
+ end
+end
+
+
+function data = addROIData(data, T, roiNames)
+ % Add synthetic ROI data to an fNIRS struct
+ nROIs = length(roiNames);
+ data.ROI.HbO = randn(T, nROIs);
+ data.ROI.HbR = randn(T, nROIs) * 0.3;
+ data.ROI.info = table(repmat({[1,2,3]}, nROIs, 1), ...
+ 'VariableNames', {'Channels'}, ...
+ 'RowNames', roiNames);
+end
diff --git a/+pf2_base/+tests/+unit/DataManipulationTest.m b/+pf2_base/+tests/+unit/DataManipulationTest.m
index e01ba14d..8d398278 100644
--- a/+pf2_base/+tests/+unit/DataManipulationTest.m
+++ b/+pf2_base/+tests/+unit/DataManipulationTest.m
@@ -41,15 +41,15 @@ function loadSampleData(testCase)
% Create evenly spaced markers
testCase.dataWithMarkers.markers = [
- minT + duration*0.05, 49, 0; % Baseline marker
- minT + duration*0.10, 50, 0; % Task start 1
- minT + duration*0.20, 51, 0; % Task end 1
- minT + duration*0.30, 50, 0; % Task start 2
- minT + duration*0.40, 51, 0; % Task end 2
- minT + duration*0.50, 49, 0; % Baseline marker
- minT + duration*0.60, 50, 0; % Task start 3
- minT + duration*0.70, 51, 0; % Task end 3
- minT + duration*0.80, 52, 0; % Different marker
+ minT + duration*0.05, 49, 0, 1; % Baseline marker
+ minT + duration*0.10, 50, 0, 1; % Task start 1
+ minT + duration*0.20, 51, 0, 1; % Task end 1
+ minT + duration*0.30, 50, 0, 1; % Task start 2
+ minT + duration*0.40, 51, 0, 1; % Task end 2
+ minT + duration*0.50, 49, 0, 1; % Baseline marker
+ minT + duration*0.60, 50, 0, 1; % Task start 3
+ minT + duration*0.70, 51, 0, 1; % Task end 3
+ minT + duration*0.80, 52, 0, 1; % Different marker
];
end
end
diff --git a/+pf2_base/+tests/+unit/DataStructureTest.m b/+pf2_base/+tests/+unit/DataStructureTest.m
index d0c3029f..c5bd7316 100644
--- a/+pf2_base/+tests/+unit/DataStructureTest.m
+++ b/+pf2_base/+tests/+unit/DataStructureTest.m
@@ -114,16 +114,17 @@ function testTimeIsMonotonic(testCase)
'Time vector must be monotonically increasing');
end
- function testMarkersHasThreeColumns(testCase)
- % Verify markers matrix has exactly 3 columns when not empty
+ function testMarkersHasExpectedColumns(testCase)
+ % Verify markers matrix has at least 3 columns when not empty
%
- % Marker format: [time, value, duration]
+ % Marker format: [time, value, duration, amplitude]
+ % Columns 1-3 are required; column 4 (amplitude) defaults to 1
data = testCase.rawData;
if ~isempty(data.markers)
- testCase.verifyEqual(size(data.markers, 2), 3, ...
- 'Markers must have 3 columns: [time, value, duration]');
+ testCase.verifyGreaterThanOrEqual(size(data.markers, 2), 3, ...
+ 'Markers must have at least 3 columns: [time, value, duration]');
else
testCase.verifyTrue(true, 'Empty markers array is valid');
end
diff --git a/+pf2_base/+tests/+unit/GLMTest.m b/+pf2_base/+tests/+unit/GLMTest.m
new file mode 100644
index 00000000..3dc03ad1
--- /dev/null
+++ b/+pf2_base/+tests/+unit/GLMTest.m
@@ -0,0 +1,405 @@
+classdef GLMTest < matlab.unittest.TestCase
+ % GLMTEST Unit tests for GLM design matrix and solver
+ %
+ % Tests cover buildDesignMatrix and fitGLM for correctness of
+ % design matrix construction, OLS beta recovery, contrast testing,
+ % and AR-IRLS convergence using synthetic data.
+ %
+ % Example:
+ % results = runtests('pf2_base.tests.unit.GLMTest');
+ % disp(results);
+ %
+ % See also: pf2_base.fnirs.buildDesignMatrix, pf2_base.fnirs.fitGLM
+
+ properties
+ fs % Sampling frequency
+ time % Time vector
+ T % Number of samples
+ end
+
+ methods (TestClassSetup)
+ function setupSyntheticData(testCase)
+ testCase.fs = 10;
+ testCase.T = 3000; % 300 seconds at 10 Hz
+ testCase.time = (0:testCase.T-1)' / testCase.fs;
+ end
+ end
+
+ %% buildDesignMatrix Tests
+ methods (Test)
+
+ function testDesignMatrixDimensions(testCase)
+ % Design matrix should have T rows and correct number of columns
+ events(1).name = 'TaskA';
+ events(1).onsets = [10 40 70 100];
+ events(1).duration = 20;
+
+ [X, names] = pf2_base.fnirs.buildDesignMatrix(testCase.time, testCase.fs, events);
+
+ testCase.verifyEqual(size(X, 1), testCase.T, ...
+ 'Design matrix should have T rows');
+ % 1 condition + 4 drift (constant + linear + quad + cubic) = 5
+ testCase.verifyEqual(size(X, 2), 5, ...
+ 'Should have 1 stim + 4 drift columns');
+ testCase.verifyEqual(length(names), size(X, 2), ...
+ 'Names should match column count');
+ end
+
+ function testMultipleConditions(testCase)
+ events(1).name = 'TaskA';
+ events(1).onsets = [10 50 90];
+ events(1).duration = 15;
+ events(2).name = 'TaskB';
+ events(2).onsets = [30 70 110];
+ events(2).duration = 15;
+
+ [X, names] = pf2_base.fnirs.buildDesignMatrix(testCase.time, testCase.fs, events);
+
+ % 2 conditions + 4 drift = 6
+ testCase.verifyEqual(size(X, 2), 6);
+ testCase.verifyEqual(names{1}, 'TaskA');
+ testCase.verifyEqual(names{2}, 'TaskB');
+ end
+
+ function testWithDerivatives(testCase)
+ events(1).name = 'TaskA';
+ events(1).onsets = [10 50];
+ events(1).duration = 20;
+
+ [X, names] = pf2_base.fnirs.buildDesignMatrix(testCase.time, testCase.fs, events, ...
+ 'IncludeDerivative', true, 'IncludeDispersion', true);
+
+ % 1 condition * 3 (primary + deriv + disp) + 4 drift = 7
+ testCase.verifyEqual(size(X, 2), 7);
+ testCase.verifyTrue(any(contains(names, 'deriv')));
+ testCase.verifyTrue(any(contains(names, 'disp')));
+ end
+
+ function testShortChannelColumns(testCase)
+ events(1).name = 'TaskA';
+ events(1).onsets = [10];
+ events(1).duration = 20;
+
+ shortCh = randn(testCase.T, 3); % 3 short channels
+
+ [X, names] = pf2_base.fnirs.buildDesignMatrix(testCase.time, testCase.fs, events, ...
+ 'ShortChannels', shortCh);
+
+ % 1 condition + 4 drift + 3 short = 8
+ testCase.verifyEqual(size(X, 2), 8);
+ testCase.verifyEqual(sum(contains(names, 'short_ch')), 3);
+ end
+
+ function testNoDriftRegressors(testCase)
+ events(1).name = 'TaskA';
+ events(1).onsets = [10];
+ events(1).duration = 20;
+
+ [X, ~] = pf2_base.fnirs.buildDesignMatrix(testCase.time, testCase.fs, events, ...
+ 'DriftOrder', -1);
+
+ testCase.verifyEqual(size(X, 2), 1, ...
+ 'With DriftOrder=-1, only stimulus columns');
+ end
+
+ function testHRFConvolutionShape(testCase)
+ % Convolved regressor should peak after stimulus onset
+ events(1).name = 'TaskA';
+ events(1).onsets = 50;
+ events(1).duration = 0; % Impulse
+
+ [X, ~] = pf2_base.fnirs.buildDesignMatrix(testCase.time, testCase.fs, events, ...
+ 'DriftOrder', -1);
+
+ stimCol = X(:, 1);
+ [~, peakIdx] = max(stimCol);
+ peakTime = testCase.time(peakIdx);
+
+ testCase.verifyGreaterThan(peakTime, 50, ...
+ 'HRF peak should occur after stimulus onset');
+ testCase.verifyLessThan(peakTime, 60, ...
+ 'HRF peak should occur within ~10s of onset');
+ end
+
+ function testBoxcarConvolution(testCase)
+ % Boxcar stimulus should produce sustained response
+ events(1).name = 'TaskA';
+ events(1).onsets = 50;
+ events(1).duration = 30;
+
+ [X, ~] = pf2_base.fnirs.buildDesignMatrix(testCase.time, testCase.fs, events, ...
+ 'DriftOrder', -1);
+
+ stimCol = X(:, 1);
+ % Response should be sustained during block
+ midBlockIdx = round(65 * testCase.fs);
+ testCase.verifyGreaterThan(stimCol(midBlockIdx), 0, ...
+ 'Signal should be positive during block');
+ end
+
+ function testImpulseDesign(testCase)
+ % Impulse (duration=0) should produce single HRF response
+ events(1).name = 'TaskA';
+ events(1).onsets = [50 100 150];
+ events(1).duration = 0;
+
+ [X, ~] = pf2_base.fnirs.buildDesignMatrix(testCase.time, testCase.fs, events, ...
+ 'DriftOrder', -1);
+
+ % Should see 3 peaks
+ stimCol = X(:, 1);
+ testCase.verifyGreaterThan(max(stimCol), 0);
+ end
+
+ function testCustomHRF(testCase)
+ events(1).name = 'TaskA';
+ events(1).onsets = 50;
+ events(1).duration = 0;
+
+ customHRF = [0; 0.5; 1; 0.8; 0.3; 0];
+
+ [X, ~] = pf2_base.fnirs.buildDesignMatrix(testCase.time, testCase.fs, events, ...
+ 'HRF', customHRF, 'DriftOrder', -1);
+
+ testCase.verifyEqual(size(X, 2), 1);
+ testCase.verifyGreaterThan(max(X(:, 1)), 0);
+ end
+
+ end
+
+ %% fitGLM Tests
+ methods (Test)
+
+ function testOLSRecoversBetas(testCase)
+ % OLS should recover known betas from synthetic data
+ rng(42);
+
+ events(1).name = 'TaskA';
+ events(1).onsets = [20 60 100 140 180 220];
+ events(1).duration = 20;
+
+ [X, names] = pf2_base.fnirs.buildDesignMatrix(testCase.time, testCase.fs, events);
+
+ % Known betas: TaskA = 2.0, constant ~= 0, drifts ~= 0
+ trueBeta = [2.0; zeros(size(X, 2) - 1, 1)];
+ trueBeta(end-3) = 0; % constant
+
+ % Generate clean data
+ Y = X * trueBeta;
+ % Add small noise
+ noise = 0.05 * randn(testCase.T, 1);
+ Y = Y + noise;
+
+ results = pf2_base.fnirs.fitGLM(Y, X, names);
+
+ % TaskA beta should be close to 2.0
+ testCase.verifyEqual(results.beta(1), 2.0, 'AbsTol', 0.2, ...
+ 'OLS should recover known beta for TaskA');
+ end
+
+ function testOLSOutputStructure(testCase)
+ events(1).name = 'TaskA';
+ events(1).onsets = [20 60 100];
+ events(1).duration = 20;
+
+ [X, names] = pf2_base.fnirs.buildDesignMatrix(testCase.time, testCase.fs, events);
+ Y = randn(testCase.T, 4); % 4 channels
+
+ results = pf2_base.fnirs.fitGLM(Y, X, names);
+
+ testCase.verifySize(results.beta, [size(X, 2), 4]);
+ testCase.verifySize(results.tstat, [size(X, 2), 4]);
+ testCase.verifySize(results.pval, [size(X, 2), 4]);
+ testCase.verifySize(results.se, [size(X, 2), 4]);
+ testCase.verifySize(results.residuals, [testCase.T, 4]);
+ testCase.verifySize(results.R2, [1, 4]);
+ testCase.verifyEqual(results.method, 'OLS');
+ end
+
+ function testOLSSignificanceWithSignal(testCase)
+ % With strong signal, GLM should detect significance
+ rng(42);
+
+ events(1).name = 'TaskA';
+ events(1).onsets = [20 60 100 140 180 220];
+ events(1).duration = 20;
+
+ [X, names] = pf2_base.fnirs.buildDesignMatrix(testCase.time, testCase.fs, events);
+
+ trueBeta = zeros(size(X, 2), 1);
+ trueBeta(1) = 5.0; % Strong task effect
+
+ Y = X * trueBeta + 0.1 * randn(testCase.T, 1);
+ results = pf2_base.fnirs.fitGLM(Y, X, names);
+
+ testCase.verifyLessThan(results.pval(1), 0.001, ...
+ 'Strong signal should yield significant p-value');
+ end
+
+ function testOLSNoSignificanceWithNoise(testCase)
+ % Pure noise should not show significance
+ rng(42);
+
+ events(1).name = 'TaskA';
+ events(1).onsets = [20 60 100];
+ events(1).duration = 20;
+
+ [X, names] = pf2_base.fnirs.buildDesignMatrix(testCase.time, testCase.fs, events);
+ Y = randn(testCase.T, 1);
+
+ results = pf2_base.fnirs.fitGLM(Y, X, names);
+
+ % Not guaranteed to fail, but unlikely to be very significant
+ testCase.verifyGreaterThan(results.pval(1), 0.001, ...
+ 'Pure noise should not yield highly significant p-value');
+ end
+
+ function testContrastTesting(testCase)
+ % Contrast should test difference between conditions
+ rng(42);
+
+ events(1).name = 'TaskA';
+ events(1).onsets = [20 80 140 200];
+ events(1).duration = 20;
+ events(2).name = 'TaskB';
+ events(2).onsets = [50 110 170 230];
+ events(2).duration = 20;
+
+ [X, names] = pf2_base.fnirs.buildDesignMatrix(testCase.time, testCase.fs, events);
+
+ % TaskA = 3.0, TaskB = 1.0
+ trueBeta = zeros(size(X, 2), 1);
+ trueBeta(1) = 3.0;
+ trueBeta(2) = 1.0;
+
+ Y = X * trueBeta + 0.1 * randn(testCase.T, 1);
+
+ % Contrast: TaskA - TaskB
+ C = zeros(1, size(X, 2));
+ C(1) = 1; C(2) = -1;
+
+ results = pf2_base.fnirs.fitGLM(Y, X, names, ...
+ 'Contrasts', C, 'ContrastNames', {'A_vs_B'});
+
+ testCase.verifyTrue(isfield(results, 'contrast'));
+ testCase.verifyEqual(results.contrast.beta, 2.0, 'AbsTol', 0.3, ...
+ 'Contrast A-B should be ~2.0');
+ testCase.verifyLessThan(results.contrast.pval, 0.001, ...
+ 'Contrast should be significant');
+ testCase.verifyEqual(results.contrast.names{1}, 'A_vs_B');
+ end
+
+ function testMultipleContrasts(testCase)
+ events(1).name = 'TaskA';
+ events(1).onsets = [20 80 140];
+ events(1).duration = 20;
+ events(2).name = 'TaskB';
+ events(2).onsets = [50 110 170];
+ events(2).duration = 20;
+
+ [X, names] = pf2_base.fnirs.buildDesignMatrix(testCase.time, testCase.fs, events);
+ Y = randn(testCase.T, 2);
+
+ C = zeros(2, size(X, 2));
+ C(1, 1) = 1; C(1, 2) = -1; % A vs B
+ C(2, 1) = 1; C(2, 2) = 1; % A + B (mean)
+
+ results = pf2_base.fnirs.fitGLM(Y, X, names, 'Contrasts', C);
+
+ testCase.verifySize(results.contrast.beta, [2, 2]);
+ testCase.verifySize(results.contrast.pval, [2, 2]);
+ end
+
+ function testARIRLSConverges(testCase)
+ % AR-IRLS should converge and produce valid results
+ rng(42);
+
+ events(1).name = 'TaskA';
+ events(1).onsets = [20 60 100 140 180 220];
+ events(1).duration = 20;
+
+ [X, names] = pf2_base.fnirs.buildDesignMatrix(testCase.time, testCase.fs, events);
+
+ trueBeta = zeros(size(X, 2), 1);
+ trueBeta(1) = 3.0;
+
+ % Add AR(1) noise to simulate autocorrelated fNIRS
+ noise = zeros(testCase.T, 1);
+ noise(1) = randn;
+ for t = 2:testCase.T
+ noise(t) = 0.8 * noise(t-1) + randn;
+ end
+ Y = X * trueBeta + 0.3 * noise;
+
+ results = pf2_base.fnirs.fitGLM(Y, X, names, 'Method', 'AR-IRLS');
+
+ testCase.verifyEqual(results.method, 'AR-IRLS');
+ testCase.verifyEqual(results.beta(1), 3.0, 'AbsTol', 0.5, ...
+ 'AR-IRLS should recover beta with autocorrelated noise');
+ end
+
+ function testARIRLSMultiChannel(testCase)
+ rng(42);
+
+ events(1).name = 'TaskA';
+ events(1).onsets = [20 60 100];
+ events(1).duration = 20;
+
+ [X, names] = pf2_base.fnirs.buildDesignMatrix(testCase.time, testCase.fs, events);
+ Y = randn(testCase.T, 3);
+
+ results = pf2_base.fnirs.fitGLM(Y, X, names, 'Method', 'AR-IRLS');
+
+ testCase.verifySize(results.beta, [size(X, 2), 3]);
+ testCase.verifySize(results.residuals, [testCase.T, 3]);
+ end
+
+ function testR2Range(testCase)
+ % R2 should be between 0 and 1 for well-conditioned problems
+ rng(42);
+
+ events(1).name = 'TaskA';
+ events(1).onsets = [20 60 100 140 180];
+ events(1).duration = 20;
+
+ [X, names] = pf2_base.fnirs.buildDesignMatrix(testCase.time, testCase.fs, events);
+
+ trueBeta = zeros(size(X, 2), 1);
+ trueBeta(1) = 3.0;
+ Y = X * trueBeta + 0.5 * randn(testCase.T, 1);
+
+ results = pf2_base.fnirs.fitGLM(Y, X, names);
+
+ testCase.verifyGreaterThanOrEqual(results.R2, 0);
+ testCase.verifyLessThanOrEqual(results.R2, 1);
+ end
+
+ function testDOFCalculation(testCase)
+ events(1).name = 'TaskA';
+ events(1).onsets = [20 60];
+ events(1).duration = 20;
+
+ [X, names] = pf2_base.fnirs.buildDesignMatrix(testCase.time, testCase.fs, events);
+ Y = randn(testCase.T, 1);
+
+ results = pf2_base.fnirs.fitGLM(Y, X, names);
+
+ testCase.verifyEqual(results.dof, testCase.T - size(X, 2), ...
+ 'DOF should be T - P for OLS');
+ end
+
+ function testSizeMismatchError(testCase)
+ events(1).name = 'TaskA';
+ events(1).onsets = [20];
+ events(1).duration = 20;
+
+ [X, names] = pf2_base.fnirs.buildDesignMatrix(testCase.time, testCase.fs, events);
+ Y = randn(testCase.T + 10, 1); % Wrong size
+
+ testCase.verifyError(@() pf2_base.fnirs.fitGLM(Y, X, names), ...
+ 'pf2:fitGLM:sizeMismatch');
+ end
+
+ end
+
+end
diff --git a/+pf2_base/+tests/+unit/HierarchicalAverageTest.m b/+pf2_base/+tests/+unit/HierarchicalAverageTest.m
new file mode 100644
index 00000000..bb2961e3
--- /dev/null
+++ b/+pf2_base/+tests/+unit/HierarchicalAverageTest.m
@@ -0,0 +1,464 @@
+classdef HierarchicalAverageTest < matlab.unittest.TestCase
+ % HIERARCHICALAVERAGETEST Unit tests for pf2_base.hierarchicalAverage
+ %
+ % Tests the hierarchical (nested) averaging function used in group
+ % analysis to prevent pseudoreplication. Covers:
+ % - Two-level hierarchy (subject > trial)
+ % - Three-level hierarchy (group > subject > trial)
+ % - Single-level hierarchy (flat grouping)
+ % - Multi-column data
+ % - NaN handling
+ % - Custom averaging functions
+ % - Cell, numeric, and table input types
+ % - Edge cases and error conditions
+ %
+ % Run all tests:
+ % results = runtests('pf2_base.tests.unit.HierarchicalAverageTest');
+ %
+ % See also: pf2_base.hierarchicalAverage
+
+ %% Core Functionality Tests
+ methods (Test)
+ function testDocumentationExample(testCase)
+ % Verify the example from the function documentation
+ %
+ % Subject1 has two conditions with two trials each:
+ % Condition 1: [10, 10] -> mean = 10
+ % Condition 2: [5, 5] -> mean = 5
+ % Subject mean: mean([10, 5]) = 7.5
+ %
+ % Subject2 has one condition:
+ % Condition 1: [2, 2] -> mean = 2
+
+ arr = [10; 10; 5; 5; 2; 2];
+ hierarchy = cell(6, 2);
+ hierarchy(:,1) = {'Subject1';'Subject1';'Subject1';'Subject1';'Subject2';'Subject2'};
+ hierarchy(:,2) = {1; 1; 2; 2; 1; 1};
+
+ [avg, subjects] = pf2_base.hierarchicalAverage(arr, hierarchy);
+
+ testCase.verifyEqual(avg(1), 7.5, 'AbsTol', 1e-10, ...
+ 'Subject1 should average to 7.5');
+ testCase.verifyEqual(avg(2), 2, 'AbsTol', 1e-10, ...
+ 'Subject2 should average to 2');
+ testCase.verifyEqual(subjects, {'Subject1'; 'Subject2'}, ...
+ 'Highest tier should return subject labels');
+ end
+
+ function testTwoLevelHierarchyPreventsPseudoreplication(testCase)
+ % Verify that hierarchical averaging differs from flat averaging
+ %
+ % Subject1 has 4 observations, Subject2 has 2. A flat mean would
+ % weight Subject1 more heavily. Hierarchical averaging gives equal
+ % weight to each subject.
+
+ arr = [10; 10; 10; 10; 0; 0];
+ hierarchy = cell(6, 2);
+ hierarchy(:,1) = {'S1';'S1';'S1';'S1';'S2';'S2'};
+ hierarchy(:,2) = {1; 1; 2; 2; 1; 1};
+
+ [avg, ~] = pf2_base.hierarchicalAverage(arr, hierarchy);
+
+ % Hierarchical: S1=10, S2=0, grand mean would be [10; 0]
+ % Flat mean would be (10+10+10+10+0+0)/6 = 6.67 (wrong)
+ testCase.verifyEqual(avg(1), 10, 'AbsTol', 1e-10, ...
+ 'Subject1 mean should be 10');
+ testCase.verifyEqual(avg(2), 0, 'AbsTol', 1e-10, ...
+ 'Subject2 mean should be 0');
+ end
+
+ function testThreeLevelHierarchy(testCase)
+ % Three-level hierarchy: Group > Subject > Trial
+ %
+ % Group A:
+ % S1: trials [10, 20] -> mean 15
+ % S2: trials [30, 40] -> mean 35
+ % Group B:
+ % S3: trials [50, 60] -> mean 55
+
+ arr = [10; 20; 30; 40; 50; 60];
+ hierarchy = cell(6, 3);
+ hierarchy(:,1) = {'A';'A';'A';'A';'B';'B'};
+ hierarchy(:,2) = {'S1';'S1';'S2';'S2';'S3';'S3'};
+ hierarchy(:,3) = {1; 2; 1; 2; 1; 2};
+
+ [avg, groups] = pf2_base.hierarchicalAverage(arr, hierarchy);
+
+ % Group A: mean([15, 35]) = 25
+ % Group B: mean([55]) = 55
+ testCase.verifyEqual(avg(1), 25, 'AbsTol', 1e-10, ...
+ 'Group A should average to 25');
+ testCase.verifyEqual(avg(2), 55, 'AbsTol', 1e-10, ...
+ 'Group B should average to 55');
+ testCase.verifyEqual(groups, {'A'; 'B'}, ...
+ 'Highest tier should return group labels');
+ end
+
+ function testSingleLevelHierarchy(testCase)
+ % Single level: just group by one factor
+ %
+ % S1: [10, 20] -> mean 15
+ % S2: [30, 40] -> mean 35
+
+ arr = [10; 20; 30; 40];
+ hierarchy = {'S1'; 'S1'; 'S2'; 'S2'};
+
+ [avg, subjects] = pf2_base.hierarchicalAverage(arr, hierarchy);
+
+ testCase.verifyEqual(avg(1), 15, 'AbsTol', 1e-10, ...
+ 'S1 should average to 15');
+ testCase.verifyEqual(avg(2), 35, 'AbsTol', 1e-10, ...
+ 'S2 should average to 35');
+ testCase.verifyEqual(subjects, {'S1'; 'S2'}, ...
+ 'Should return subject labels');
+ end
+
+ function testMultiColumnData(testCase)
+ % Multiple data columns averaged independently
+
+ arr = [10, 100; 20, 200; 30, 300; 40, 400];
+ hierarchy = {'S1'; 'S1'; 'S2'; 'S2'};
+
+ [avg, ~] = pf2_base.hierarchicalAverage(arr, hierarchy);
+
+ testCase.verifyEqual(size(avg), [2, 2], ...
+ 'Output should have 2 rows x 2 columns');
+ testCase.verifyEqual(avg(1,1), 15, 'AbsTol', 1e-10, ...
+ 'S1 col1 should be 15');
+ testCase.verifyEqual(avg(1,2), 150, 'AbsTol', 1e-10, ...
+ 'S1 col2 should be 150');
+ testCase.verifyEqual(avg(2,1), 35, 'AbsTol', 1e-10, ...
+ 'S2 col1 should be 35');
+ testCase.verifyEqual(avg(2,2), 350, 'AbsTol', 1e-10, ...
+ 'S2 col2 should be 350');
+ end
+
+ function testAllUniqueRows(testCase)
+ % When all hierarchy rows are unique, no averaging needed
+
+ arr = [10; 20; 30];
+ hierarchy = {'S1'; 'S2'; 'S3'};
+
+ [avg, subjects] = pf2_base.hierarchicalAverage(arr, hierarchy);
+
+ testCase.verifyEqual(avg, [10; 20; 30], 'AbsTol', 1e-10, ...
+ 'No averaging should occur when all rows unique');
+ testCase.verifyEqual(numel(subjects), 3, ...
+ 'Should return 3 labels');
+ end
+
+ function testSingleObservation(testCase)
+ % Single observation returns itself
+
+ arr = [42];
+ hierarchy = {'S1'};
+
+ [avg, ~] = pf2_base.hierarchicalAverage(arr, hierarchy);
+
+ testCase.verifyEqual(avg, 42, 'AbsTol', 1e-10, ...
+ 'Single observation should return itself');
+ end
+ end
+
+ %% NaN Handling Tests
+ methods (Test)
+ function testNaNValuesIgnored(testCase)
+ % NaN values should be ignored in averaging (nanmean default)
+
+ arr = [10; NaN; 20; 30];
+ hierarchy = {'S1'; 'S1'; 'S2'; 'S2'};
+
+ [avg, ~] = pf2_base.hierarchicalAverage(arr, hierarchy);
+
+ testCase.verifyEqual(avg(1), 10, 'AbsTol', 1e-10, ...
+ 'S1 with one NaN should use only the valid value');
+ testCase.verifyEqual(avg(2), 25, 'AbsTol', 1e-10, ...
+ 'S2 should average normally');
+ end
+
+ function testAllNaNForOneGroup(testCase)
+ % All NaN for one group should produce NaN output
+
+ arr = [NaN; NaN; 20; 30];
+ hierarchy = {'S1'; 'S1'; 'S2'; 'S2'};
+
+ [avg, ~] = pf2_base.hierarchicalAverage(arr, hierarchy);
+
+ testCase.verifyTrue(isnan(avg(1)), ...
+ 'All-NaN group should produce NaN');
+ testCase.verifyEqual(avg(2), 25, 'AbsTol', 1e-10, ...
+ 'S2 should still average normally');
+ end
+ end
+
+ %% Custom Averaging Function Tests
+ methods (Test)
+ function testCustomFunctionHandle(testCase)
+ % Use @nanmedian instead of default @nanmean
+
+ arr = [1; 2; 100; 10; 20; 30];
+ hierarchy = cell(6, 2);
+ hierarchy(:,1) = {'S1';'S1';'S1';'S2';'S2';'S2'};
+ hierarchy(:,2) = {1; 2; 3; 1; 2; 3};
+
+ [avg_mean, ~] = pf2_base.hierarchicalAverage(arr, hierarchy);
+ [avg_median, ~] = pf2_base.hierarchicalAverage(arr, hierarchy, @nanmedian);
+
+ % S1: mean([1,2,100]) = 34.33, median([1,2,100]) = 2
+ testCase.verifyEqual(avg_mean(1), mean([1,2,100]), 'AbsTol', 1e-10, ...
+ 'Mean for S1 should match');
+ testCase.verifyEqual(avg_median(1), 2, 'AbsTol', 1e-10, ...
+ 'Median for S1 should be 2');
+ end
+
+ function testCustomFunctionString(testCase)
+ % Pass averaging function as a string name
+
+ arr = [10; 20; 30; 40];
+ hierarchy = {'S1'; 'S1'; 'S2'; 'S2'};
+
+ [avg, ~] = pf2_base.hierarchicalAverage(arr, hierarchy, 'nanmean');
+
+ testCase.verifyEqual(avg(1), 15, 'AbsTol', 1e-10, ...
+ 'String function name should work like handle');
+ end
+
+ function testInvalidFuncErrors(testCase)
+ % Invalid function argument should error
+
+ arr = [10; 20];
+ hierarchy = {'S1'; 'S1'};
+
+ testCase.verifyError(...
+ @() pf2_base.hierarchicalAverage(arr, hierarchy, 12345), ...
+ '', ...
+ 'Non-function non-string third argument should error');
+ end
+ end
+
+ %% Input Type Tests
+ methods (Test)
+ function testNumericHierarchy(testCase)
+ % Numeric hierarchy array
+
+ arr = [10; 20; 30; 40];
+ hierarchy = [1 1; 1 2; 2 1; 2 2];
+
+ [avg, ~] = pf2_base.hierarchicalAverage(arr, hierarchy);
+
+ testCase.verifyEqual(avg(1), 15, 'AbsTol', 1e-10, ...
+ 'Group 1 should average to 15');
+ testCase.verifyEqual(avg(2), 35, 'AbsTol', 1e-10, ...
+ 'Group 2 should average to 35');
+ end
+
+ function testCellArrayData(testCase)
+ % Cell array data input (converted to matrix internally)
+
+ arr = {10; 20; 30; 40};
+ hierarchy = {'S1'; 'S1'; 'S2'; 'S2'};
+
+ [avg, ~] = pf2_base.hierarchicalAverage(arr, hierarchy);
+
+ testCase.verifyEqual(avg(1), 15, 'AbsTol', 1e-10, ...
+ 'Cell array data should work like numeric');
+ end
+
+ function testTransposedDataAutoCorrects(testCase)
+ % Row vector data should be auto-transposed to match hierarchy
+
+ arr = [10, 20, 30, 40]; % 1x4 row vector
+ hierarchy = {'S1'; 'S1'; 'S2'; 'S2'}; % 4x1
+
+ [avg, ~] = pf2_base.hierarchicalAverage(arr, hierarchy);
+
+ testCase.verifyEqual(avg(1), 15, 'AbsTol', 1e-10, ...
+ 'Auto-transposed S1 should average to 15');
+ testCase.verifyEqual(avg(2), 35, 'AbsTol', 1e-10, ...
+ 'Auto-transposed S2 should average to 35');
+ end
+
+ function testCellHierarchyWithNumericValues(testCase)
+ % Cell hierarchy containing numeric values
+
+ arr = [10; 20; 30; 40];
+ hierarchy = cell(4, 2);
+ hierarchy(:,1) = {'S1';'S1';'S2';'S2'};
+ hierarchy(:,2) = {1; 2; 1; 2};
+
+ [avg, subjects] = pf2_base.hierarchicalAverage(arr, hierarchy);
+
+ testCase.verifyEqual(numel(avg), 2, ...
+ 'Should produce 2 output rows');
+ testCase.verifyEqual(subjects, {'S1'; 'S2'}, ...
+ 'Should return string labels from column 1');
+ end
+ end
+
+ %% Output Structure Tests
+ methods (Test)
+ function testOutputRowCount(testCase)
+ % Output rows should equal number of unique highest-tier groups
+
+ arr = [1; 2; 3; 4; 5; 6];
+ hierarchy = {'A';'A';'B';'B';'C';'C'};
+
+ [avg, subjects] = pf2_base.hierarchicalAverage(arr, hierarchy);
+
+ testCase.verifyEqual(size(avg, 1), 3, ...
+ 'Should have 3 output rows for 3 subjects');
+ testCase.verifyEqual(numel(subjects), 3, ...
+ 'Should have 3 labels');
+ end
+
+ function testHighestTierOutput(testCase)
+ % highestTier should contain unique labels from column 1
+
+ arr = [1; 2; 3; 4; 5; 6];
+ hierarchy = cell(6, 2);
+ hierarchy(:,1) = {'X';'X';'Y';'Y';'Z';'Z'};
+ hierarchy(:,2) = {1; 2; 1; 2; 1; 2};
+
+ [~, labels] = pf2_base.hierarchicalAverage(arr, hierarchy);
+
+ testCase.verifyEqual(numel(labels), 3, ...
+ 'Should have 3 unique top-level labels');
+ testCase.verifyTrue(all(ismember({'X';'Y';'Z'}, labels)), ...
+ 'Labels should contain X, Y, Z');
+ end
+
+ function testThirdOutputExists(testCase)
+ % outHarr (third output) should be returned for debugging
+
+ arr = [10; 20; 30; 40];
+ hierarchy = cell(4, 2);
+ hierarchy(:,1) = {'S1';'S1';'S2';'S2'};
+ hierarchy(:,2) = {1; 2; 1; 2};
+
+ [~, ~, outHarr] = pf2_base.hierarchicalAverage(arr, hierarchy);
+
+ testCase.verifyFalse(isempty(outHarr), ...
+ 'Third output (outHarr) should not be empty');
+ end
+ end
+
+ %% Unbalanced Design Tests
+ methods (Test)
+ function testUnbalancedTrialsPerSubject(testCase)
+ % Different number of trials per subject
+
+ arr = [10; 20; 30; 40; 50];
+ hierarchy = cell(5, 2);
+ hierarchy(:,1) = {'S1';'S1';'S1';'S2';'S2'};
+ hierarchy(:,2) = {1; 2; 3; 1; 2};
+
+ [avg, ~] = pf2_base.hierarchicalAverage(arr, hierarchy);
+
+ testCase.verifyEqual(avg(1), 20, 'AbsTol', 1e-10, ...
+ 'S1 with 3 trials: mean([10,20,30]) = 20');
+ testCase.verifyEqual(avg(2), 45, 'AbsTol', 1e-10, ...
+ 'S2 with 2 trials: mean([40,50]) = 45');
+ end
+
+ function testUnbalancedConditionsAcrossSubjects(testCase)
+ % Some subjects have conditions that others do not
+ %
+ % S1: Condition A [10,20] -> 15, Condition B [30] -> 30
+ % Subject mean: mean([15, 30]) = 22.5
+ % S2: Condition A [40,50] -> 45
+ % Subject mean: 45
+
+ arr = [10; 20; 30; 40; 50];
+ hierarchy = cell(5, 2);
+ hierarchy(:,1) = {'S1';'S1';'S1';'S2';'S2'};
+ hierarchy(:,2) = {'A'; 'A'; 'B'; 'A'; 'A'};
+
+ [avg, ~] = pf2_base.hierarchicalAverage(arr, hierarchy);
+
+ testCase.verifyEqual(avg(1), 22.5, 'AbsTol', 1e-10, ...
+ 'S1: mean([mean([10,20]), 30]) = 22.5');
+ testCase.verifyEqual(avg(2), 45, 'AbsTol', 1e-10, ...
+ 'S2: mean([40,50]) = 45');
+ end
+ end
+
+ %% Error Condition Tests
+ methods (Test)
+ function testNoInputErrors(testCase)
+ % No arguments should produce an error
+
+ testCase.verifyError(...
+ @() pf2_base.hierarchicalAverage(), ...
+ '', ...
+ 'No input should error');
+ end
+
+ function testMissingHierarchyErrors(testCase)
+ % Missing hierarchy should produce an error
+
+ testCase.verifyError(...
+ @() pf2_base.hierarchicalAverage([1; 2; 3]), ...
+ '', ...
+ 'Missing hierarchy should error');
+ end
+
+ function testMismatchedDimensionsErrors(testCase)
+ % Data rows not matching hierarchy rows should error
+
+ arr = [1; 2; 3];
+ hierarchy = {'S1'; 'S2'}; % 2 rows vs 3 data rows
+
+ testCase.verifyError(...
+ @() pf2_base.hierarchicalAverage(arr, hierarchy), ...
+ '', ...
+ 'Mismatched dimensions should error');
+ end
+ end
+
+ %% Realistic fNIRS-style Tests
+ methods (Test)
+ function testRealisticGroupAnalysis(testCase)
+ % Simulate a realistic fNIRS group analysis scenario
+ %
+ % 3 subjects, 2 conditions each, 2 trials per condition
+ % Data has 4 channels (columns)
+
+ nChannels = 4;
+ % S1: Cond1 trials
+ data = [
+ 1, 2, 3, 4; % S1, C1, T1
+ 3, 4, 5, 6; % S1, C1, T2
+ 5, 6, 7, 8; % S1, C2, T1
+ 7, 8, 9, 10; % S1, C2, T2
+ 10, 20, 30, 40; % S2, C1, T1
+ 12, 22, 32, 42; % S2, C1, T2
+ 14, 24, 34, 44; % S2, C2, T1
+ 16, 26, 36, 46; % S2, C2, T2
+ 100, 200, 300, 400; % S3, C1, T1
+ 100, 200, 300, 400; % S3, C1, T2
+ 100, 200, 300, 400; % S3, C2, T1
+ 100, 200, 300, 400; % S3, C2, T2
+ ];
+
+ hierarchy = cell(12, 3);
+ hierarchy(:,1) = {'S1';'S1';'S1';'S1';'S2';'S2';'S2';'S2';'S3';'S3';'S3';'S3'};
+ hierarchy(:,2) = {'C1';'C1';'C2';'C2';'C1';'C1';'C2';'C2';'C1';'C1';'C2';'C2'};
+ hierarchy(:,3) = {1; 2; 1; 2; 1; 2; 1; 2; 1; 2; 1; 2};
+
+ [avg, subjects] = pf2_base.hierarchicalAverage(data, hierarchy);
+
+ testCase.verifyEqual(size(avg, 1), 3, ...
+ 'Should have 3 subject rows');
+ testCase.verifyEqual(size(avg, 2), nChannels, ...
+ 'Should preserve 4 channels');
+
+ % S1 channel 1: mean([mean([1,3]), mean([5,7])]) = mean([2, 6]) = 4
+ testCase.verifyEqual(avg(1, 1), 4, 'AbsTol', 1e-10, ...
+ 'S1 channel 1 hierarchical average');
+
+ % S3 all channels should be 100, 200, 300, 400 (all identical)
+ testCase.verifyEqual(avg(3, :), [100, 200, 300, 400], 'AbsTol', 1e-10, ...
+ 'S3 with identical data should return same values');
+ end
+ end
+end
diff --git a/+pf2_base/+tests/+unit/NormalizeMarkersTest.m b/+pf2_base/+tests/+unit/NormalizeMarkersTest.m
new file mode 100644
index 00000000..60aeb55c
--- /dev/null
+++ b/+pf2_base/+tests/+unit/NormalizeMarkersTest.m
@@ -0,0 +1,69 @@
+classdef NormalizeMarkersTest < matlab.unittest.TestCase
+ % NORMALIZEMARKERSTEST Unit tests for pf2_base.normalizeMarkers
+ %
+ % Verifies that marker arrays are correctly padded to 4 columns
+ % with appropriate defaults (duration=0, amplitude=1).
+
+ methods (Test)
+ function testTwoColumnInput(testCase)
+ % 2-col input should get duration=0 and amplitude=1
+ mrk = [1.0, 50; 2.5, 51];
+ result = pf2_base.normalizeMarkers(mrk);
+
+ testCase.verifySize(result, [2, 4]);
+ testCase.verifyEqual(result(:,1:2), mrk);
+ testCase.verifyEqual(result(:,3), [0; 0], 'Duration should default to 0');
+ testCase.verifyEqual(result(:,4), [1; 1], 'Amplitude should default to 1');
+ end
+
+ function testThreeColumnInput(testCase)
+ % 3-col input should get amplitude=1
+ mrk = [1.0, 50, 5; 2.5, 51, 10];
+ result = pf2_base.normalizeMarkers(mrk);
+
+ testCase.verifySize(result, [2, 4]);
+ testCase.verifyEqual(result(:,1:3), mrk);
+ testCase.verifyEqual(result(:,4), [1; 1], 'Amplitude should default to 1');
+ end
+
+ function testFourColumnPassthrough(testCase)
+ % 4-col input should be returned as-is
+ mrk = [1.0, 50, 5, 0.8; 2.5, 51, 10, 1.5];
+ result = pf2_base.normalizeMarkers(mrk);
+
+ testCase.verifyEqual(result, mrk);
+ end
+
+ function testEmptyInput(testCase)
+ % Empty input should return zeros(0,4)
+ result = pf2_base.normalizeMarkers([]);
+
+ testCase.verifySize(result, [0, 4]);
+ testCase.verifyTrue(isempty(result));
+ end
+
+ function testAmplitudePreserved(testCase)
+ % Custom amplitude values should be preserved
+ mrk = [1.0, 50, 0, 2.0; 3.0, 51, 0, 0.5];
+ result = pf2_base.normalizeMarkers(mrk);
+
+ testCase.verifyEqual(result(:,4), [2.0; 0.5]);
+ end
+
+ function testSingleRow(testCase)
+ % Single-row marker should work
+ result = pf2_base.normalizeMarkers([5.0, 1]);
+
+ testCase.verifySize(result, [1, 4]);
+ testCase.verifyEqual(result, [5.0, 1, 0, 1]);
+ end
+
+ function testFiveColumnPassthrough(testCase)
+ % 5+ column input should be returned as-is
+ mrk = [1.0, 50, 5, 1, 99];
+ result = pf2_base.normalizeMarkers(mrk);
+
+ testCase.verifyEqual(result, mrk);
+ end
+ end
+end
diff --git a/+pf2_base/+tests/+unit/PlotHelpersTest.m b/+pf2_base/+tests/+unit/PlotHelpersTest.m
index de6e38cb..f9bf0662 100644
--- a/+pf2_base/+tests/+unit/PlotHelpersTest.m
+++ b/+pf2_base/+tests/+unit/PlotHelpersTest.m
@@ -11,7 +11,7 @@
methods (Test)
%% processMarkers tests
function testProcessMarkers_showAll(testCase)
- fNIR.markers = [1, 10, 0; 2, 20, 0; 3, 10, 0; 4, 30, 0];
+ fNIR.markers = [1, 10, 0, 1; 2, 20, 0, 1; 3, 10, 0, 1; 4, 30, 0, 1];
[codes, idx, data, counts] = pf2_base.plot.processMarkers(fNIR, true);
testCase.verifyEqual(length(codes), 3); % 10, 20, 30
@@ -20,7 +20,7 @@ function testProcessMarkers_showAll(testCase)
end
function testProcessMarkers_showNone(testCase)
- fNIR.markers = [1, 10, 0; 2, 20, 0];
+ fNIR.markers = [1, 10, 0, 1; 2, 20, 0, 1];
[codes, idx, data, counts] = pf2_base.plot.processMarkers(fNIR, false);
testCase.verifyEmpty(codes);
@@ -28,7 +28,7 @@ function testProcessMarkers_showNone(testCase)
end
function testProcessMarkers_showSpecific(testCase)
- fNIR.markers = [1, 10, 0; 2, 20, 0; 3, 10, 0; 4, 30, 0];
+ fNIR.markers = [1, 10, 0, 1; 2, 20, 0, 1; 3, 10, 0, 1; 4, 30, 0, 1];
[codes, idx, data, counts] = pf2_base.plot.processMarkers(fNIR, [10, 30]);
testCase.verifyEqual(length(codes), 2);
@@ -53,7 +53,7 @@ function testProcessMarkers_emptyMarkers(testCase)
end
function testProcessMarkers_allString(testCase)
- fNIR.markers = [1, 10, 0; 2, 20, 0];
+ fNIR.markers = [1, 10, 0, 1; 2, 20, 0, 1];
[codes, ~, ~, ~] = pf2_base.plot.processMarkers(fNIR, 'all');
testCase.verifyEqual(length(codes), 2);
diff --git a/+pf2_base/+tests/+unit/ProcessingContextTest.m b/+pf2_base/+tests/+unit/ProcessingContextTest.m
index fe3dd239..12ca65ae 100644
--- a/+pf2_base/+tests/+unit/ProcessingContextTest.m
+++ b/+pf2_base/+tests/+unit/ProcessingContextTest.m
@@ -404,6 +404,135 @@ function testCopyMethodForValueSemantics(testCase)
% ctx2 should be unchanged
testCase.verifyEqual(ctx2.subjectAge, 25);
end
+
+ %% Serialization Version Tests
+
+ function testToStructContainsVersion(testCase)
+ % TESTTOSTRUCTCONTAINSVERSION Verify contextVersion field exists
+ ctx = pf2_base.ProcessingContext();
+ s = ctx.toStruct();
+
+ testCase.verifyTrue(isfield(s, 'contextVersion'));
+ testCase.verifyEqual(s.contextVersion, '1.0');
+ end
+
+ function testToStructContainsCreated(testCase)
+ % TESTTOSTRUCTCONTAINSCREATED Verify created timestamp field
+ ctx = pf2_base.ProcessingContext();
+ s = ctx.toStruct();
+
+ testCase.verifyTrue(isfield(s, 'created'));
+ testCase.verifyClass(s.created, 'char');
+ testCase.verifyNotEmpty(s.created);
+ end
+
+ function testFromStructWithoutVersion(testCase)
+ % TESTFROMSTRUCTWITHOUTVERSION Legacy structs without version load correctly
+ s = struct();
+ s.dpfMode = 'Fixed';
+ s.dpfFixedValue = 6.0;
+
+ ctx = pf2_base.ProcessingContext.fromStruct(s);
+
+ testCase.verifyEqual(ctx.dpfMode, 'Fixed');
+ testCase.verifyEqual(ctx.dpfFixedValue, 6.0);
+ end
+
+ function testStructRoundtripPreservesVersion(testCase)
+ % TESTSTRUCTROUNDTRIPPRESERVESVERSION Version survives roundtrip
+ ctx1 = pf2_base.ProcessingContext();
+ s1 = ctx1.toStruct();
+
+ ctx2 = pf2_base.ProcessingContext.fromStruct(s1);
+ s2 = ctx2.toStruct();
+
+ testCase.verifyEqual(s2.contextVersion, '1.0');
+ testCase.verifyTrue(isfield(s2, 'created'));
+ end
+
+ %% fromStruct Edge Cases
+
+ function testFromStructEmpty(testCase)
+ % TESTFROMSTRUCTEMPTY Empty struct returns all defaults
+ s = struct();
+ ctx = pf2_base.ProcessingContext.fromStruct(s);
+
+ testCase.verifyEqual(ctx.dpfMode, 'Calc');
+ testCase.verifyEqual(ctx.dpfFixedValue, 5.93);
+ testCase.verifyEqual(ctx.subjectAge, 25);
+ testCase.verifyEqual(ctx.baselineLength, 10);
+ testCase.verifyEqual(ctx.rejectLevel, 0);
+ end
+
+ function testFromStructWithExtraFields(testCase)
+ % TESTFROMSTRUCTWITHEXTRAFIELDS Unknown fields silently ignored
+ s = struct();
+ s.dpfMode = 'Fixed';
+ s.unknownField = 42;
+ s.anotherExtra = 'hello';
+
+ ctx = pf2_base.ProcessingContext.fromStruct(s);
+
+ testCase.verifyEqual(ctx.dpfMode, 'Fixed');
+ end
+
+ function testFromStructPartialDPFOnly(testCase)
+ % TESTFROMSTRUCTPARTIALDPFONLY Partial struct sets DPF, defaults rest
+ s = struct();
+ s.dpfMode = 'Fixed';
+ s.dpfFixedValue = 7.5;
+
+ ctx = pf2_base.ProcessingContext.fromStruct(s);
+
+ testCase.verifyEqual(ctx.dpfMode, 'Fixed');
+ testCase.verifyEqual(ctx.dpfFixedValue, 7.5);
+ testCase.verifyEqual(ctx.baselineLength, 10); % default
+ testCase.verifyEqual(ctx.rejectLevel, 0); % default
+ end
+
+ function testFromStructPartialBaselineOnly(testCase)
+ % TESTFROMSTRUCTPARTIALBASELINEONLY Partial struct sets baseline, defaults rest
+ s = struct();
+ s.baselineStartTime = 5;
+ s.baselineLength = 20;
+
+ ctx = pf2_base.ProcessingContext.fromStruct(s);
+
+ testCase.verifyEqual(ctx.baselineStartTime, 5);
+ testCase.verifyEqual(ctx.baselineLength, 20);
+ testCase.verifyEqual(ctx.dpfMode, 'Calc'); % default
+ testCase.verifyEqual(ctx.subjectAge, 25); % default
+ end
+
+ %% rejectLevel Boundary Tests
+
+ function testRejectLevelBoundaryZero(testCase)
+ % TESTREJECTLEVELBOUNDARYZERO rejectLevel = 0 is valid
+ ctx = pf2_base.ProcessingContext();
+ ctx.rejectLevel = 0;
+ testCase.verifyEqual(ctx.rejectLevel, 0);
+ end
+
+ function testRejectLevelBoundaryOne(testCase)
+ % TESTREJECTLEVELBOUNDARYONE rejectLevel = 1 is valid
+ ctx = pf2_base.ProcessingContext();
+ ctx.rejectLevel = 1;
+ testCase.verifyEqual(ctx.rejectLevel, 1);
+ end
+
+ function testRejectLevelAboveRange(testCase)
+ % TESTREJECTLEVELABOVERANGE rejectLevel = 1.1 errors
+ ctx = pf2_base.ProcessingContext();
+ testCase.verifyError(@() setfield(ctx, 'rejectLevel', 1.1), ...
+ 'MATLAB:validators:mustBeInRange');
+ end
+
+ function testRejectLevelBelowRange(testCase)
+ % TESTREJECTLEVELBELOWRANGE rejectLevel = -0.1 errors
+ ctx = pf2_base.ProcessingContext();
+ testCase.verifyError(@() setfield(ctx, 'rejectLevel', -0.1), ...
+ 'MATLAB:validators:mustBeInRange');
+ end
end
end
diff --git a/+pf2_base/+tests/+unit/QualityControlTest.m b/+pf2_base/+tests/+unit/QualityControlTest.m
new file mode 100644
index 00000000..b1d76f81
--- /dev/null
+++ b/+pf2_base/+tests/+unit/QualityControlTest.m
@@ -0,0 +1,309 @@
+classdef QualityControlTest < matlab.unittest.TestCase
+ % QUALITYCONTROLTEST Unit tests for pf2.qc signal quality functions
+ %
+ % Tests cover:
+ % - SCI (Scalp Coupling Index): cardiac correlation, dead channels,
+ % threshold classification, explicit wavelength params, output dims
+ % - Power Spectrum: known sinusoid peaks, cardiac/respiratory detection,
+ % noise-only, HbO signal, channel subset, output dimensions
+ % - plotQuality: SCI bar chart, PSD overlay, PSD tiled
+ %
+ % Example:
+ % results = runtests('pf2_base.tests.unit.QualityControlTest');
+ % disp(results);
+
+ properties
+ dataWithHeart % Synthetic data with heartbeat
+ dataNoHeart % Synthetic data without heartbeat
+ dataWithResp % Synthetic data with heartbeat + respiration
+ end
+
+ methods (TestClassSetup)
+ function createTestData(testCase)
+ % Ensure functions/ is on the path (for bpf, etc.)
+ % mfilename path: +pf2_base/+tests/+unit/QualityControlTest
+ % Need 4 fileparts to get to repo root
+ rootPath = fileparts(fileparts(fileparts(fileparts(mfilename('fullpath')))));
+ addpath(rootPath);
+ addpath(fullfile(rootPath, 'functions'));
+
+ % Generate synthetic data with known properties
+ testCase.dataWithHeart = pf2_base.tests.synthetic.generateFNIRS( ...
+ 'duration', 60, ...
+ 'fs', 10, ...
+ 'nChannels', 4, ...
+ 'addHeartbeat', true, ...
+ 'heartRate', 70, ...
+ 'heartAmplitude', 0.01, ...
+ 'noiseLevel', 0.001, ...
+ 'seed', 42);
+
+ testCase.dataNoHeart = pf2_base.tests.synthetic.generateFNIRS( ...
+ 'duration', 60, ...
+ 'fs', 10, ...
+ 'nChannels', 4, ...
+ 'addHeartbeat', false, ...
+ 'noiseLevel', 0.02, ...
+ 'seed', 43);
+
+ testCase.dataWithResp = pf2_base.tests.synthetic.generateFNIRS( ...
+ 'duration', 120, ...
+ 'fs', 10, ...
+ 'nChannels', 4, ...
+ 'addHeartbeat', true, ...
+ 'heartRate', 70, ...
+ 'heartAmplitude', 0.01, ...
+ 'addRespiration', true, ...
+ 'respRate', 15, ...
+ 'respAmplitude', 0.01, ...
+ 'noiseLevel', 0.001, ...
+ 'seed', 44);
+ end
+ end
+
+
+ %% SCI Tests
+ methods (Test)
+
+ function testSCIWithCardiacSignal(testCase)
+ % Synthetic data with heartbeat should yield high SCI
+ result = pf2.qc.sci(testCase.dataWithHeart);
+
+ testCase.verifyGreaterThan(mean(result.sci), 0.7, ...
+ 'Mean SCI should be high (>0.7) when heartbeat is present.');
+ end
+
+ function testSCIWithoutCardiacSignal(testCase)
+ % Noise-only data should yield low SCI
+ result = pf2.qc.sci(testCase.dataNoHeart);
+
+ testCase.verifyLessThan(mean(result.sci), 0.5, ...
+ 'Mean SCI should be low (<0.5) without heartbeat.');
+ end
+
+ function testSCIDeadChannel(testCase)
+ % Set one channel to constant — SCI should be 0
+ data = testCase.dataWithHeart;
+ % Channel 1 uses columns 1 and 2 (alternating wavelengths)
+ data.raw(:, 1) = 1000;
+ data.raw(:, 2) = 1000;
+
+ result = pf2.qc.sci(data);
+
+ testCase.verifyEqual(result.sci(1), 0, ...
+ 'Dead channel (constant signal) should have SCI = 0.');
+ % Other channels should still have valid SCI
+ testCase.verifyGreaterThan(result.sci(2), 0, ...
+ 'Non-dead channels should have SCI > 0.');
+ end
+
+ function testSCIThresholdClassification(testCase)
+ % Verify isGood matches sci >= threshold
+ threshold = 0.6;
+ result = pf2.qc.sci(testCase.dataWithHeart, 'Threshold', threshold);
+
+ expected = result.sci >= threshold;
+ testCase.verifyEqual(result.isGood, expected, ...
+ 'isGood should match sci >= threshold.');
+ testCase.verifyEqual(result.threshold, threshold);
+ end
+
+ function testSCIExplicitWavelengthParams(testCase)
+ % Pass Wavelengths and ChannelNumbers manually
+ data = testCase.dataWithHeart;
+ nCh = data.info.synthetic.nChannels;
+ wl = repmat([730, 850], 1, nCh);
+ chNums = repelem(1:nCh, 2);
+
+ result = pf2.qc.sci(data, 'Wavelengths', wl, 'ChannelNumbers', chNums);
+
+ testCase.verifyEqual(numel(result.sci), nCh);
+ testCase.verifyGreaterThan(mean(result.sci), 0.7);
+ end
+
+ function testSCIOutputDimensions(testCase)
+ % Verify output shapes
+ result = pf2.qc.sci(testCase.dataWithHeart);
+ nCh = testCase.dataWithHeart.info.synthetic.nChannels;
+
+ testCase.verifySize(result.sci, [1, nCh]);
+ testCase.verifySize(result.isGood, [1, nCh]);
+ testCase.verifySize(result.channels, [1, nCh]);
+ testCase.verifyEqual(result.fs, testCase.dataWithHeart.fs);
+ end
+
+ end
+
+
+ %% Power Spectrum Tests
+ methods (Test)
+
+ function testPSDKnownSinusoid(testCase)
+ % Single-frequency signal should have peak at correct frequency
+ fs = 100;
+ t = (0:1/fs:30)';
+ targetFreq = 2.5;
+ sig = sin(2 * pi * targetFreq * t);
+
+ % Build minimal data struct
+ data = struct();
+ data.fs = fs;
+ data.HbO = sig;
+ data.fchMask = 1;
+
+ result = pf2.qc.powerSpectrum(data, 'Signal', 'HbO', ...
+ 'FreqRange', [0, fs/2], 'DetectPeaks', false);
+
+ % Find peak frequency in PSD
+ [~, peakIdx] = max(result.psd);
+ peakFreq = result.freqs(peakIdx);
+
+ testCase.verifyEqual(peakFreq, targetFreq, 'AbsTol', 0.2, ...
+ 'Peak frequency should match the injected sinusoid.');
+ end
+
+ function testPSDCardiacPeakDetection(testCase)
+ % Data with heartbeat should show cardiac peak near 1 Hz
+ result = pf2.qc.powerSpectrum(testCase.dataWithHeart, ...
+ 'Signal', 'raw', 'DetectPeaks', true);
+
+ expectedFreq = testCase.dataWithHeart.info.synthetic.heartRate / 60;
+
+ % At least some channels should have cardiac detected
+ testCase.verifyTrue(any(result.cardiac.detected), ...
+ 'Cardiac peak should be detected in at least one channel.');
+
+ % Detected frequency should be near the known heart rate
+ detectedIdx = find(result.cardiac.detected, 1);
+ if ~isempty(detectedIdx)
+ testCase.verifyEqual(result.cardiac.freq(detectedIdx), ...
+ expectedFreq, 'AbsTol', 0.7, ...
+ 'Detected cardiac frequency should be near heart rate.');
+ end
+ end
+
+ function testPSDNoCardiacInNoise(testCase)
+ % White noise should not reliably show cardiac peak
+ % Use high noise to drown out any structure
+ data = pf2_base.tests.synthetic.generateFNIRS( ...
+ 'duration', 60, 'fs', 10, 'nChannels', 4, ...
+ 'addHeartbeat', false, 'noiseLevel', 0.05, 'seed', 99);
+
+ result = pf2.qc.powerSpectrum(data, 'Signal', 'raw', ...
+ 'DetectPeaks', true);
+
+ % Most channels should not have cardiac detected
+ fractionDetected = sum(result.cardiac.detected) / numel(result.channels);
+ testCase.verifyLessThanOrEqual(fractionDetected, 0.75, ...
+ 'Noise-only data should not reliably show cardiac peaks.');
+ end
+
+ function testPSDRespiratoryPeakDetection(testCase)
+ % Signal with known respiratory frequency should show peak
+ fs = 100;
+ t = (0:1/fs:120)';
+ respFreq = 0.25; % 15 breaths/min
+ nCh = 2;
+ sig = repmat(sin(2 * pi * respFreq * t) + 0.1 * randn(size(t)), 1, nCh);
+
+ data = struct('fs', fs, 'HbO', sig, 'fchMask', ones(1, nCh));
+ result = pf2.qc.powerSpectrum(data, 'Signal', 'HbO', 'DetectPeaks', true);
+
+ testCase.verifyTrue(any(result.respiratory.detected), ...
+ 'Respiratory peak should be detected when respiration is present.');
+ detIdx = find(result.respiratory.detected, 1);
+ if ~isempty(detIdx)
+ testCase.verifyEqual(result.respiratory.freq(detIdx), ...
+ respFreq, 'AbsTol', 0.05, ...
+ 'Detected respiratory frequency should be near 0.25 Hz.');
+ end
+ end
+
+ function testPSDOnHbOSignal(testCase)
+ % Verify PSD works on processed (HbO) data
+ data = struct();
+ data.fs = 10;
+ t = (0:1/data.fs:60)';
+ nCh = 4;
+ data.HbO = randn(numel(t), nCh) * 0.01;
+ data.fchMask = ones(1, nCh);
+
+ result = pf2.qc.powerSpectrum(data, 'Signal', 'HbO');
+
+ testCase.verifyEqual(result.signal, 'HbO');
+ testCase.verifyEqual(size(result.psd, 2), nCh);
+ end
+
+ function testPSDChannelSubset(testCase)
+ % 'Channels' parameter should select only those channels
+ data = struct();
+ data.fs = 10;
+ t = (0:1/data.fs:60)';
+ data.HbO = randn(numel(t), 6);
+ data.fchMask = ones(1, 6);
+
+ result = pf2.qc.powerSpectrum(data, 'Signal', 'HbO', ...
+ 'Channels', [1, 3]);
+
+ testCase.verifyEqual(numel(result.channels), 2);
+ testCase.verifyEqual(result.channels, [1, 3]);
+ testCase.verifyEqual(size(result.psd, 2), 2);
+ end
+
+ function testPSDOutputDimensions(testCase)
+ % Verify [F x C] shape and freqs within FreqRange
+ freqRange = [0, 3];
+ result = pf2.qc.powerSpectrum(testCase.dataWithHeart, ...
+ 'Signal', 'raw', 'FreqRange', freqRange, 'DetectPeaks', false);
+
+ nCh = numel(result.channels);
+ nFreqs = numel(result.freqs);
+
+ testCase.verifySize(result.psd, [nFreqs, nCh]);
+ testCase.verifyGreaterThanOrEqual(min(result.freqs), freqRange(1));
+ testCase.verifyLessThanOrEqual(max(result.freqs), freqRange(2));
+ end
+
+ end
+
+
+ %% Plot Tests
+ methods (Test)
+
+ function testPlotSCI(testCase)
+ % plotQuality with SCI result should create a figure
+ result = pf2.qc.sci(testCase.dataWithHeart);
+ fig = pf2.qc.plotQuality(result, 'Visible', 'off');
+
+ testCase.addTeardown(@() close(fig));
+ testCase.verifyTrue(ishandle(fig), ...
+ 'plotQuality should return a valid figure handle for SCI.');
+ end
+
+ function testPlotPSDOverlay(testCase)
+ % plotQuality with PSD result in overlay mode
+ result = pf2.qc.powerSpectrum(testCase.dataWithHeart, ...
+ 'Signal', 'raw');
+ fig = pf2.qc.plotQuality(result, 'Visible', 'off', ...
+ 'Layout', 'overlay');
+
+ testCase.addTeardown(@() close(fig));
+ testCase.verifyTrue(ishandle(fig), ...
+ 'plotQuality should return a valid figure handle for PSD overlay.');
+ end
+
+ function testPlotPSDTiled(testCase)
+ % plotQuality with PSD result in tiled mode
+ result = pf2.qc.powerSpectrum(testCase.dataWithHeart, ...
+ 'Signal', 'raw');
+ fig = pf2.qc.plotQuality(result, 'Visible', 'off', ...
+ 'Layout', 'tiled');
+
+ testCase.addTeardown(@() close(fig));
+ testCase.verifyTrue(ishandle(fig), ...
+ 'plotQuality should return a valid figure handle for PSD tiled.');
+ end
+
+ end
+
+end
diff --git a/+pf2_base/+tests/+unit/SSRTest.m b/+pf2_base/+tests/+unit/SSRTest.m
new file mode 100644
index 00000000..2567b18f
--- /dev/null
+++ b/+pf2_base/+tests/+unit/SSRTest.m
@@ -0,0 +1,243 @@
+classdef SSRTest < matlab.unittest.TestCase
+ % SSRTEST Unit tests for short-channel regression
+ %
+ % Tests verify that shortChannelRegression correctly removes
+ % superficial signals while preserving brain signals, and that
+ % all three methods (nearest, pca, all) function correctly.
+ %
+ % Example:
+ % results = runtests('pf2_base.tests.unit.SSRTest');
+ % disp(results);
+ %
+ % See also: pf2_base.fnirs.shortChannelRegression, pf2_SSR
+
+ properties
+ testData % Synthetic fNIRS struct with short channels
+ end
+
+ methods (TestClassSetup)
+ function addFunctionsPath(~)
+ % Ensure functions/ directory is on path for pf2_SSR wrapper
+ projRoot = fileparts(fileparts(fileparts(fileparts(mfilename('fullpath')))));
+ funcDir = fullfile(projRoot, 'functions');
+ if isfolder(funcDir)
+ addpath(funcDir);
+ end
+ end
+
+ function buildSyntheticData(testCase)
+ % Create synthetic data with known brain and superficial signals
+ rng(42);
+
+ T = 1000;
+ fs = 10;
+ nLong = 8;
+ nShort = 2;
+ nOpt = nLong + nShort;
+
+ time = (0:T-1)' / fs;
+
+ % Superficial physiology (cardiac + respiration)
+ superficial = 0.5 * sin(2*pi*1.0*time) + 0.3 * sin(2*pi*0.25*time);
+
+ % Brain signal (task-related HRF response in channels 1-4)
+ brain = zeros(T, nLong);
+ hrf = pf2_base.fnirs.buildHRF(fs);
+ hrfVec = hrf(:, 2);
+ stim = zeros(T, 1);
+ stim(round([10 30 50 70] * fs)) = 1;
+ hrfResponse = conv(stim, hrfVec);
+ hrfResponse = hrfResponse(1:T);
+ brain(:, 1:4) = repmat(2 * hrfResponse, 1, 4);
+
+ % Build HbO: brain + superficial + noise
+ HbO = zeros(T, nOpt);
+ for ch = 1:nLong
+ HbO(:, ch) = brain(:, ch) + superficial + 0.1*randn(T, 1);
+ end
+ % Short channels: only superficial + noise
+ for ch = 1:nShort
+ HbO(:, nLong + ch) = superficial + 0.05*randn(T, 1);
+ end
+
+ HbR = -0.3 * HbO + 0.05*randn(T, nOpt);
+
+ % Build probe info with short-channel flags
+ probeInfo = struct();
+ probeInfo.Probe = cell(1, 1);
+
+ isShort = false(1, nOpt);
+ isShort(nLong+1:end) = true;
+ probeInfo.Probe{1}.IsShortSeparation = isShort;
+ probeInfo.Probe{1}.NumOptodes = nOpt;
+
+ % 3D positions: long channels spread out, short channels near
+ optX = [(1:nLong)*30, 15 75]';
+ optY = [zeros(1, nLong), 5 5]';
+ optZ = zeros(nOpt, 1);
+ probeInfo.Probe{1}.OptPosX = optX;
+ probeInfo.Probe{1}.OptPosY = optY;
+ probeInfo.Probe{1}.OptPosZ = optZ;
+ probeInfo.Probe{1}.OptPos3D = [optX, optY, optZ];
+
+ data = struct();
+ data.HbO = HbO;
+ data.HbR = HbR;
+ data.time = time;
+ data.fs = fs;
+ data.probeinfo = probeInfo;
+
+ % Store brain signal for verification
+ data.testBrain = brain;
+ data.testSuperficial = superficial;
+
+ testCase.testData = data;
+ end
+ end
+
+ %% Nearest method tests
+ methods (Test)
+
+ function testNearestReducesSuperficial(testCase)
+ % SSR with nearest method should reduce superficial component
+ data = testCase.testData;
+ corrected = pf2_base.fnirs.shortChannelRegression(data, 'Method', 'nearest');
+
+ % Correlation between corrected signal and superficial should decrease
+ longIdx = find(~data.probeinfo.Probe{1}.IsShortSeparation);
+ origCorr = abs(corr(data.HbO(:, longIdx(1)), data.testSuperficial));
+ corrCorr = abs(corr(corrected.HbO(:, longIdx(1)), data.testSuperficial));
+
+ testCase.verifyLessThan(corrCorr, origCorr, ...
+ 'Correlation with superficial signal should decrease after SSR');
+ end
+
+ function testNearestPreservesBrain(testCase)
+ % Brain signal should be largely preserved after SSR
+ data = testCase.testData;
+ corrected = pf2_base.fnirs.shortChannelRegression(data, 'Method', 'nearest');
+
+ % Channel 1 has brain signal - correlation should remain high
+ longIdx = find(~data.probeinfo.Probe{1}.IsShortSeparation);
+ brainCorr = corr(corrected.HbO(:, longIdx(1)), data.testBrain(:, 1));
+
+ testCase.verifyGreaterThan(brainCorr, 0.5, ...
+ 'Brain signal correlation should remain substantial after SSR');
+ end
+
+ function testNearestShortChannelsUnchanged(testCase)
+ % Short channel data should not be modified
+ data = testCase.testData;
+ corrected = pf2_base.fnirs.shortChannelRegression(data, 'Method', 'nearest');
+
+ shortIdx = find(data.probeinfo.Probe{1}.IsShortSeparation);
+ testCase.verifyEqual(corrected.HbO(:, shortIdx), data.HbO(:, shortIdx), ...
+ 'Short channel data should be unchanged');
+ end
+
+ function testNearestSSRInfoField(testCase)
+ data = testCase.testData;
+ corrected = pf2_base.fnirs.shortChannelRegression(data, 'Method', 'nearest');
+
+ testCase.verifyTrue(isfield(corrected, 'ssrInfo'));
+ testCase.verifyEqual(corrected.ssrInfo.method, 'nearest');
+ end
+
+ end
+
+ %% PCA method tests
+ methods (Test)
+
+ function testPCAReducesSuperficial(testCase)
+ data = testCase.testData;
+ corrected = pf2_base.fnirs.shortChannelRegression(data, 'Method', 'pca');
+
+ longIdx = find(~data.probeinfo.Probe{1}.IsShortSeparation);
+ origCorr = abs(corr(data.HbO(:, longIdx(1)), data.testSuperficial));
+ corrCorr = abs(corr(corrected.HbO(:, longIdx(1)), data.testSuperficial));
+
+ testCase.verifyLessThan(corrCorr, origCorr, ...
+ 'PCA method should reduce superficial correlation');
+ end
+
+ function testPCANumPCs(testCase)
+ data = testCase.testData;
+ corrected = pf2_base.fnirs.shortChannelRegression(data, ...
+ 'Method', 'pca', 'NumPCs', 2);
+
+ testCase.verifyTrue(isfield(corrected, 'ssrInfo'));
+ testCase.verifyEqual(corrected.ssrInfo.numPCs, 2);
+ end
+
+ end
+
+ %% All method tests
+ methods (Test)
+
+ function testAllReducesSuperficial(testCase)
+ data = testCase.testData;
+ corrected = pf2_base.fnirs.shortChannelRegression(data, 'Method', 'all');
+
+ longIdx = find(~data.probeinfo.Probe{1}.IsShortSeparation);
+ origCorr = abs(corr(data.HbO(:, longIdx(1)), data.testSuperficial));
+ corrCorr = abs(corr(corrected.HbO(:, longIdx(1)), data.testSuperficial));
+
+ testCase.verifyLessThan(corrCorr, origCorr, ...
+ 'All method should reduce superficial correlation');
+ end
+
+ end
+
+ %% Edge case tests
+ methods (Test)
+
+ function testNoShortChannelsWarning(testCase)
+ % Should warn when no short channels present
+ data = testCase.testData;
+ data.probeinfo.Probe{1}.IsShortSeparation = false(1, ...
+ data.probeinfo.Probe{1}.NumOptodes);
+
+ testCase.verifyWarning(@() pf2_base.fnirs.shortChannelRegression(data), ...
+ 'pf2:ssr:noShortChannels');
+ end
+
+ function testNoProbeInfoWarning(testCase)
+ % Should warn when no probe info present
+ data = struct('HbO', randn(100, 5), 'HbR', randn(100, 5));
+
+ testCase.verifyWarning(@() pf2_base.fnirs.shortChannelRegression(data), ...
+ 'pf2:ssr:noProbe');
+ end
+
+ function testHbRAlsoCorrected(testCase)
+ data = testCase.testData;
+ corrected = pf2_base.fnirs.shortChannelRegression(data, 'Method', 'nearest');
+
+ % HbR should also be modified
+ longIdx = find(~data.probeinfo.Probe{1}.IsShortSeparation);
+ testCase.verifyFalse(isequal(corrected.HbR(:, longIdx), data.HbR(:, longIdx)), ...
+ 'HbR should also be corrected');
+ end
+
+ function testCustomBiomarkers(testCase)
+ data = testCase.testData;
+ corrected = pf2_base.fnirs.shortChannelRegression(data, ...
+ 'Method', 'nearest', 'Biomarkers', {'HbO'});
+
+ % HbR should be unchanged when not in Biomarkers list
+ testCase.verifyEqual(corrected.HbR, data.HbR, ...
+ 'HbR should be unchanged when not specified in Biomarkers');
+ end
+
+ function testPf2SSRWrapper(testCase)
+ % Test the method-chain wrapper
+ data = testCase.testData;
+ corrected = pf2_SSR(data, 'nearest');
+
+ testCase.verifyTrue(isfield(corrected, 'ssrInfo'));
+ testCase.verifyEqual(corrected.ssrInfo.method, 'nearest');
+ end
+
+ end
+
+end
diff --git a/+pf2_base/+tests/+unit/SplitTest.m b/+pf2_base/+tests/+unit/SplitTest.m
new file mode 100644
index 00000000..8d25c8cc
--- /dev/null
+++ b/+pf2_base/+tests/+unit/SplitTest.m
@@ -0,0 +1,370 @@
+classdef SplitTest < matlab.unittest.TestCase
+ % SPLITTEST Unit tests for pf2.data.split
+ %
+ % Tests the time segmentation function including:
+ % - Basic time extraction (startTime, endTime)
+ % - segmentLength parameter
+ % - Relative vs absolute time modes
+ % - Baseline correction (blLength, blStartTime, blfNIR)
+ % - Marker filtering to extracted window
+ % - Field preservation (HbO, HbR, raw, etc.)
+ % - Edge cases and error conditions
+ %
+ % Run all tests:
+ % results = runtests('pf2_base.tests.unit.SplitTest');
+ %
+ % See also: pf2.data.split, pf2.data.resample, pf2.data.setT0
+
+ properties
+ processedData % Processed fNIRS sample data
+ rawData % Raw (unprocessed) fNIRS sample data
+ end
+
+ methods (TestClassSetup)
+ function loadSampleData(testCase)
+ testCase.rawData = pf2.import.sampleData.fNIR2000();
+ testCase.processedData = processFNIRS2(testCase.rawData, 'ShowGUI', false);
+ end
+ end
+
+ %% Basic Extraction Tests
+ methods (Test)
+ function testSplitStartTimeOnly(testCase)
+ % Providing only startTime extracts from start to end
+
+ timeVec = testCase.processedData.time;
+ midTime = mean(timeVec);
+
+ seg = pf2.data.split(testCase.processedData, midTime);
+
+ testCase.verifyGreaterThanOrEqual(min(seg.time), midTime, ...
+ 'Segment should start at or after startTime');
+ testCase.verifyEqual(max(seg.time), max(timeVec), 'AbsTol', 1/testCase.processedData.fs, ...
+ 'Segment should extend to end of recording');
+ end
+
+ function testSplitStartAndEndTime(testCase)
+ % Providing startTime and endTime extracts that window
+
+ timeVec = testCase.processedData.time;
+ startT = min(timeVec) + 50;
+ endT = min(timeVec) + 150;
+
+ seg = pf2.data.split(testCase.processedData, startT, endT);
+
+ testCase.verifyGreaterThanOrEqual(min(seg.time), startT, ...
+ 'Segment should start at or after startTime');
+ testCase.verifyLessThanOrEqual(max(seg.time), endT, ...
+ 'Segment should end at or before endTime');
+ end
+
+ function testSplitReducesSamples(testCase)
+ % Extracted segment should have fewer samples than original
+
+ timeVec = testCase.processedData.time;
+ startT = min(timeVec) + 50;
+ endT = min(timeVec) + 100;
+
+ seg = pf2.data.split(testCase.processedData, startT, endT);
+
+ testCase.verifyLessThan(length(seg.time), length(timeVec), ...
+ 'Segment should have fewer samples than original');
+ end
+
+ function testSplitApproximateSampleCount(testCase)
+ % Number of samples should match expected duration * fs
+
+ fs = testCase.processedData.fs;
+ timeVec = testCase.processedData.time;
+ startT = min(timeVec) + 20;
+ duration = 60; % 60 seconds
+ endT = startT + duration;
+
+ seg = pf2.data.split(testCase.processedData, startT, endT);
+
+ expectedSamples = duration * fs;
+ actualSamples = length(seg.time);
+ testCase.verifyEqual(actualSamples, expectedSamples, 'RelTol', 0.02, ...
+ 'Sample count should approximately match duration * fs');
+ end
+ end
+
+ %% segmentLength Parameter Tests
+ methods (Test)
+ function testSegmentLength(testCase)
+ % Using segmentLength instead of endTime
+
+ timeVec = testCase.processedData.time;
+ startT = min(timeVec) + 30;
+ segLen = 40;
+
+ seg = pf2.data.split(testCase.processedData, startT, nan, segLen);
+
+ segDuration = max(seg.time) - min(seg.time);
+ testCase.verifyEqual(segDuration, segLen, 'AbsTol', 1/testCase.processedData.fs, ...
+ 'Segment duration should match segmentLength');
+ end
+ end
+
+ %% Relative Time Mode Tests
+ methods (Test)
+ function testRelativeTimeMode(testCase)
+ % Relative mode equivalent: compute absolute times from offsets
+ %
+ % Since split() uses positional optionals, we test relative
+ % behavior by computing absolute times manually.
+
+ timeVec = testCase.processedData.time;
+ relStart = 20; % 20s from beginning
+ relEnd = 80; % 80s from beginning
+
+ absStart = min(timeVec) + relStart;
+ absEnd = min(timeVec) + relEnd;
+
+ seg = pf2.data.split(testCase.processedData, absStart, absEnd);
+
+ testCase.verifyGreaterThanOrEqual(min(seg.time), absStart, ...
+ 'Start should match expected absolute time');
+ testCase.verifyLessThanOrEqual(max(seg.time), absEnd, ...
+ 'End should match expected absolute time');
+
+ segDuration = max(seg.time) - min(seg.time);
+ testCase.verifyEqual(segDuration, relEnd - relStart, ...
+ 'AbsTol', 1/testCase.processedData.fs, ...
+ 'Duration should match relative time span');
+ end
+ end
+
+ %% Field Preservation Tests
+ methods (Test)
+ function testSplitPreservesOxyFields(testCase)
+ % HbO, HbR, HbDiff, HbTotal, CBSI should all be extracted
+
+ timeVec = testCase.processedData.time;
+ startT = min(timeVec) + 30;
+ endT = min(timeVec) + 90;
+
+ seg = pf2.data.split(testCase.processedData, startT, endT);
+
+ oxyFields = {'HbO', 'HbR', 'HbDiff', 'HbTotal', 'CBSI'};
+ for i = 1:length(oxyFields)
+ f = oxyFields{i};
+ testCase.verifyTrue(isfield(seg, f), ...
+ sprintf('%s field should be preserved', f));
+ testCase.verifyEqual(size(seg.(f), 1), length(seg.time), ...
+ sprintf('%s rows should match time vector length', f));
+ end
+ end
+
+ function testSplitPreservesChannelCount(testCase)
+ % Number of channels should be unchanged
+
+ origChannels = size(testCase.processedData.HbO, 2);
+ timeVec = testCase.processedData.time;
+ startT = min(timeVec) + 30;
+ endT = min(timeVec) + 90;
+
+ seg = pf2.data.split(testCase.processedData, startT, endT);
+
+ testCase.verifyEqual(size(seg.HbO, 2), origChannels, ...
+ 'Channel count should be preserved after split');
+ end
+
+ function testSplitPreservesMetadata(testCase)
+ % Non-timeseries fields should be preserved
+
+ timeVec = testCase.processedData.time;
+ startT = min(timeVec) + 30;
+ endT = min(timeVec) + 90;
+
+ seg = pf2.data.split(testCase.processedData, startT, endT);
+
+ metaFields = {'info', 'fchMask', 'fs', 'channels'};
+ for i = 1:length(metaFields)
+ f = metaFields{i};
+ if isfield(testCase.processedData, f)
+ testCase.verifyTrue(isfield(seg, f), ...
+ sprintf('%s should be preserved after split', f));
+ end
+ end
+ end
+
+ function testSplitPreservesFsUnchanged(testCase)
+ % Sampling rate should not change
+
+ timeVec = testCase.processedData.time;
+ startT = min(timeVec) + 10;
+ endT = min(timeVec) + 50;
+
+ seg = pf2.data.split(testCase.processedData, startT, endT);
+
+ testCase.verifyEqual(seg.fs, testCase.processedData.fs, ...
+ 'Sampling rate should be preserved after split');
+ end
+ end
+
+ %% Baseline Correction Tests
+ methods (Test)
+ function testBaselineCorrectionSubtractsMean(testCase)
+ % Baseline correction should subtract baseline mean from segment
+
+ timeVec = testCase.processedData.time;
+ startT = min(timeVec) + 20;
+ endT = min(timeVec) + 100;
+ blLen = 10;
+
+ % Split without baseline
+ segNoBL = pf2.data.split(testCase.processedData, startT, endT);
+
+ % Split with baseline (positional args: startTime, endTime, segLen, relative, blLength)
+ segLen = endT - startT;
+ segBL = pf2.data.split(testCase.processedData, startT, endT, segLen, false, blLen);
+
+ % Baseline-corrected data should differ from uncorrected
+ testCase.verifyNotEqual(segBL.HbO, segNoBL.HbO, ...
+ 'Baseline correction should change HbO values');
+ end
+
+ function testBaselineCorrectionUsesExternalBaseline(testCase)
+ % blfNIR parameter should use a separate struct for baseline
+
+ timeVec = testCase.processedData.time;
+ startT = min(timeVec) + 50;
+ endT = min(timeVec) + 100;
+
+ % Create a baseline from first 30s
+ blStart = min(timeVec);
+ blEnd = min(timeVec) + 30;
+ blData = pf2.data.split(testCase.processedData, blStart, blEnd);
+
+ % Split using external baseline
+ seg = pf2.data.split(testCase.processedData, startT, endT, ...
+ 'blfNIR', blData);
+
+ testCase.verifyTrue(isfield(seg, 'HbO'), ...
+ 'Should have HbO after split with external baseline');
+ testCase.verifyEqual(size(seg.HbO, 2), size(testCase.processedData.HbO, 2), ...
+ 'Channel count should be preserved with external baseline');
+ end
+ end
+
+ %% Marker Tests
+ methods (Test)
+ function testSplitFiltersMarkers(testCase)
+ % Markers outside the split window should be removed
+
+ % Add synthetic markers to data
+ dataWithMarkers = testCase.processedData;
+ timeVec = dataWithMarkers.time;
+ minT = min(timeVec);
+
+ dataWithMarkers.markers = [
+ minT + 10, 1, 0, 1; % Before split window
+ minT + 60, 2, 0, 1; % Inside split window
+ minT + 70, 3, 0, 1; % Inside split window
+ minT + 200, 4, 0, 1; % After split window
+ ];
+
+ startT = minT + 50;
+ endT = minT + 100;
+ seg = pf2.data.split(dataWithMarkers, startT, endT);
+
+ if isfield(seg, 'markers') && ~isempty(seg.markers)
+ if isnumeric(seg.markers)
+ markerTimes = seg.markers(:,1);
+ elseif isstruct(seg.markers) && isfield(seg.markers, 'data')
+ markerTimes = seg.markers.data(:,1);
+ else
+ markerTimes = [];
+ end
+
+ if ~isempty(markerTimes)
+ testCase.verifyGreaterThanOrEqual(min(markerTimes), startT, ...
+ 'Markers before start should be removed');
+ testCase.verifyLessThanOrEqual(max(markerTimes), endT, ...
+ 'Markers after end should be removed');
+ end
+ end
+ end
+ end
+
+ %% Raw Data Tests
+ methods (Test)
+ function testSplitRawData(testCase)
+ % Split should also work on raw (unprocessed) data
+
+ timeVec = testCase.rawData.time;
+ startT = min(timeVec) + 20;
+ endT = min(timeVec) + 80;
+
+ seg = pf2.data.split(testCase.rawData, startT, endT);
+
+ testCase.verifyTrue(isfield(seg, 'raw'), ...
+ 'Raw field should be preserved');
+ testCase.verifyTrue(isfield(seg, 'time'), ...
+ 'Time field should be preserved');
+ testCase.verifyEqual(size(seg.raw, 1), length(seg.time), ...
+ 'Raw rows should match time vector length');
+ end
+ end
+
+ %% Error and Edge Case Tests
+ methods (Test)
+ function testEndBeforeStartErrors(testCase)
+ % endTime before startTime should error
+
+ timeVec = testCase.processedData.time;
+ startT = min(timeVec) + 100;
+ endT = min(timeVec) + 50;
+
+ testCase.verifyError(...
+ @() pf2.data.split(testCase.processedData, startT, endT), ...
+ '', ...
+ 'End time before start time should error');
+ end
+
+ function testFullExtraction(testCase)
+ % Extracting the full time range should preserve all data
+
+ timeVec = testCase.processedData.time;
+ startT = min(timeVec);
+ endT = max(timeVec);
+
+ seg = pf2.data.split(testCase.processedData, startT, endT);
+
+ testCase.verifyEqual(length(seg.time), length(timeVec), ...
+ 'Full extraction should preserve all samples');
+ end
+
+ function testSplitTimeVectorIsMonotonic(testCase)
+ % Output time vector should be monotonically increasing
+
+ timeVec = testCase.processedData.time;
+ startT = min(timeVec) + 30;
+ endT = min(timeVec) + 90;
+
+ seg = pf2.data.split(testCase.processedData, startT, endT);
+
+ testCase.verifyTrue(all(diff(seg.time) > 0), ...
+ 'Output time vector should be monotonically increasing');
+ end
+
+ function testSplitDataValuesAreFromOriginal(testCase)
+ % Extracted values should be exact copies from original
+
+ timeVec = testCase.processedData.time;
+ startT = min(timeVec) + 30;
+ endT = min(timeVec) + 90;
+
+ seg = pf2.data.split(testCase.processedData, startT, endT);
+
+ % Find corresponding indices in original
+ idxStart = find(testCase.processedData.time >= startT, 1);
+ idxEnd = find(testCase.processedData.time <= endT, 1, 'last');
+ expectedHbO = testCase.processedData.HbO(idxStart:idxEnd, :);
+
+ testCase.verifyEqual(seg.HbO, expectedHbO, 'AbsTol', 1e-15, ...
+ 'Extracted HbO should exactly match original data slice');
+ end
+ end
+end
diff --git a/+pf2_base/+tests/generateGoldenFiles.m b/+pf2_base/+tests/generateGoldenFiles.m
new file mode 100644
index 00000000..28c506a7
--- /dev/null
+++ b/+pf2_base/+tests/generateGoldenFiles.m
@@ -0,0 +1,149 @@
+function generateGoldenFiles(varargin)
+% GENERATEGOLDENFILES Generate golden reference files for regression testing
+%
+% Creates golden .mat files from deterministic sample data using specific
+% processing configurations. Files are saved to golden/candidates/ for
+% review before promotion to golden/.
+%
+% Syntax:
+% pf2_base.tests.generateGoldenFiles()
+% pf2_base.tests.generateGoldenFiles('Promote', true)
+%
+% Options (Name-Value):
+% 'Promote' - If true, save directly to golden/ instead of candidates/
+% (default: false)
+%
+% Example:
+% pf2_base.tests.generateGoldenFiles();
+% pf2_base.tests.generateGoldenFiles('Promote', true);
+%
+% See also: pf2_base.tests.integration.GoldenFileTest,
+% pf2_base.tests.golden.computeHash
+
+p = inputParser;
+addParameter(p, 'Promote', false, @islogical);
+parse(p, varargin{:});
+
+promote = p.Results.Promote;
+
+% Get project root
+thisFile = mfilename('fullpath');
+projectRoot = fileparts(fileparts(fileparts(thisFile)));
+
+if promote
+ pipelineDir = fullfile(projectRoot, 'golden', 'processFNIRS2');
+ functionDir = fullfile(projectRoot, 'golden', 'functions');
+else
+ pipelineDir = fullfile(projectRoot, 'golden', 'candidates');
+ functionDir = fullfile(projectRoot, 'golden', 'candidates');
+end
+
+% Ensure directories exist
+if ~isfolder(pipelineDir), mkdir(pipelineDir); end
+if ~isfolder(functionDir), mkdir(functionDir); end
+
+% Load sample data
+fprintf('Loading sample data...\n');
+data = pf2.import.sampleData.fNIR2000();
+inputHash = pf2_base.tests.golden.computeHash(data.raw);
+
+[~, pf2ver] = pf2_base.pf2version();
+verStr = sprintf('%.1f', pf2ver);
+
+% --- Golden 1: Default processing (no raw method, no oxy method) ---
+fprintf('Generating: fNIR2000_default...\n');
+pf2.methods.raw.setMethod('None');
+pf2.methods.oxy.setMethod('None');
+processed = processFNIRS2(data, 'ShowGUI', false);
+
+golden = struct();
+golden.output = extractOutput(processed);
+golden.inputHash = inputHash;
+golden.version = verStr;
+golden.timestamp = char(datetime('now', 'Format', 'yyyy-MM-dd''T''HH:mm:ss'));
+golden.params = struct('rawMethod', 'None', 'oxyMethod', 'None');
+golden.matlabVersion = version();
+
+outFile = fullfile(pipelineDir, 'fNIR2000_default.mat');
+save(outFile, '-struct', 'golden');
+fprintf(' Saved: %s\n', outFile);
+
+% --- Golden 2: TDDR raw, no oxy ---
+fprintf('Generating: fNIR2000_TDDR_None...\n');
+
+% Check if x5_TDDR method exists
+global PF2
+rawMethods = PF2.myRawMethods.cfg.Sections;
+if ismember('x5_TDDR', rawMethods)
+ pf2.methods.raw.setMethod('x5_TDDR');
+ rawMethodName = 'x5_TDDR';
+else
+ % Create a temporary TDDR method
+ pf2.methods.raw.create('golden_TDDR', ...
+ {struct('f', 'pf2_MotionCorrectTDDR', 'args', {{'x', 'fs'}}, ...
+ 'argvals', {{'x', 'fs'}}, 'output', 'x')}, ...
+ 'Replace', true);
+ pf2.methods.raw.setMethod('golden_TDDR');
+ rawMethodName = 'golden_TDDR';
+end
+pf2.methods.oxy.setMethod('None');
+processed = processFNIRS2(data, 'ShowGUI', false);
+
+golden = struct();
+golden.output = extractOutput(processed);
+golden.inputHash = inputHash;
+golden.version = verStr;
+golden.timestamp = char(datetime('now', 'Format', 'yyyy-MM-dd''T''HH:mm:ss'));
+golden.params = struct('rawMethod', rawMethodName, 'oxyMethod', 'None');
+golden.matlabVersion = version();
+
+outFile = fullfile(pipelineDir, 'fNIR2000_TDDR_None.mat');
+save(outFile, '-struct', 'golden');
+fprintf(' Saved: %s\n', outFile);
+
+% Clean up temp method if created
+if strcmp(rawMethodName, 'golden_TDDR')
+ pf2.methods.raw.delete('golden_TDDR');
+end
+
+% --- Golden 3: TDDR function in isolation ---
+fprintf('Generating: pf2_TDDR_fNIR2000...\n');
+od = pf2_Intensity2OD(data.raw);
+tddrOutput = pf2_MotionCorrectTDDR(od, data.fs);
+
+golden = struct();
+golden.output = struct('corrected', tddrOutput);
+golden.inputHash = pf2_base.tests.golden.computeHash(od);
+golden.version = verStr;
+golden.timestamp = char(datetime('now', 'Format', 'yyyy-MM-dd''T''HH:mm:ss'));
+golden.params = struct('fs', data.fs, 'function', 'pf2_MotionCorrectTDDR');
+golden.matlabVersion = version();
+
+outFile = fullfile(functionDir, 'pf2_TDDR_fNIR2000.mat');
+save(outFile, '-struct', 'golden');
+fprintf(' Saved: %s\n', outFile);
+
+% Reset methods
+pf2.methods.raw.setMethod('None');
+pf2.methods.oxy.setMethod('None');
+
+fprintf('Golden file generation complete.\n');
+if ~promote
+ fprintf('Files saved to golden/candidates/. Review and promote with:\n');
+ fprintf(' movefile(''golden/candidates/file.mat'', ''golden/processFNIRS2/file.mat'')\n');
+end
+
+end
+
+
+function out = extractOutput(processed)
+% Extract key output fields for golden comparison
+out = struct();
+fields = {'HbO', 'HbR', 'HbTotal', 'HbDiff', 'CBSI', 'units', 'DPF_factor'};
+for i = 1:length(fields)
+ f = fields{i};
+ if isfield(processed, f)
+ out.(f) = processed.(f);
+ end
+end
+end
diff --git a/+pf2_base/+tests/testExperiment.m b/+pf2_base/+tests/testExperiment.m
new file mode 100644
index 00000000..c6fffab4
--- /dev/null
+++ b/+pf2_base/+tests/testExperiment.m
@@ -0,0 +1,750 @@
+classdef testExperiment < matlab.unittest.TestCase
+% TESTEXPERIMENT Unit tests for exploreFNIRS.core.Experiment
+%
+% Tests the scriptable Experiment container class including:
+% - Construction from processed data
+% - Selection and filtering
+% - Groupby operations
+% - Aggregation with preprocessing
+% - Info variable analysis (plotInfoVar)
+% - Export to long/wide format
+% - Headless plotting
+%
+% Run with:
+% results = runtests('pf2_base.tests.testExperiment');
+
+ properties (TestParameter)
+ end
+
+ properties
+ allData % Cell array of processed fNIRS structs
+ nSegments % Total number of segments
+ end
+
+ methods (TestClassSetup)
+ function buildTestData(tc)
+ % Process sample data and create multi-subject dataset
+ raw = pf2.import.sampleData.fNIR2000();
+ processed = processFNIRS2(raw, 'ShowGUI', false);
+
+ rng(42);
+ subjects = {'S01','S01','S01','S01', 'S02','S02','S02','S02', 'S03','S03','S03','S03'};
+ groups = {'Ctrl','Ctrl','Ctrl','Ctrl','Ctrl','Ctrl','Ctrl','Ctrl','Tx','Tx','Tx','Tx'};
+ conds = {'Easy','Hard','Easy','Hard','Easy','Hard','Easy','Hard','Easy','Hard','Easy','Hard'};
+ ages = [25, 25, 25, 25, 30, 30, 30, 30, 28, 28, 28, 28];
+ trials = [1, 1, 2, 2, 1, 1, 2, 2, 1, 1, 2, 2];
+
+ tc.nSegments = length(subjects);
+ tc.allData = cell(tc.nSegments, 1);
+ for i = 1:tc.nSegments
+ d = processed;
+ d.info.SubjectID = subjects{i};
+ d.info.Group = groups{i};
+ d.info.Condition = conds{i};
+ d.info.Age = ages(i);
+ d.info.Trial = trials(i);
+ d.info.reactionTime = 200 + 100*strcmp(conds{i},'Hard') + randn*20;
+ d.info.accuracy = 0.9 - 0.15*strcmp(conds{i},'Hard') + randn*0.03;
+
+ % Add synthetic Aux data
+ nSamples = length(d.time);
+ d.Aux.accelerometer.data = 0.01*randn(nSamples, 3);
+ d.Aux.accelerometer.time = d.time;
+ d.Aux.accelerometer.unit = 'g';
+ d.Aux.heartRate.data = 70 + 2*randn(nSamples, 1);
+ d.Aux.heartRate.time = d.time;
+ d.Aux.heartRate.unit = 'bpm';
+
+ tc.allData{i} = d;
+ end
+ end
+ end
+
+ methods (Test)
+
+ %% --- Construction ---
+
+ function testConstructor(tc)
+ ex = exploreFNIRS.core.Experiment(tc.allData);
+ tc.verifyEqual(length(ex.data), tc.nSegments);
+ tc.verifyEqual(height(ex.dataTable), tc.nSegments);
+ tc.verifyFalse(ex.isGrouped);
+ tc.verifyFalse(ex.isAggregated);
+ end
+
+ function testConstructorRejectsEmpty(tc)
+ tc.verifyError(@() exploreFNIRS.core.Experiment({}), ...
+ 'exploreFNIRS:core:Experiment');
+ end
+
+ function testConstructorRejectsNonCell(tc)
+ tc.verifyError(@() exploreFNIRS.core.Experiment(42), ...
+ 'exploreFNIRS:core:Experiment');
+ end
+
+ function testDataTableHasMissingFNIRS(tc)
+ ex = exploreFNIRS.core.Experiment(tc.allData);
+ tc.verifyTrue(ismember('missingFNIRS', ex.dataTable.Properties.VariableNames));
+ tc.verifyEqual(ex.dataTable.missingFNIRS, zeros(tc.nSegments, 1));
+ end
+
+ function testDataTableHasInfoFields(tc)
+ ex = exploreFNIRS.core.Experiment(tc.allData);
+ vars = ex.dataTable.Properties.VariableNames;
+ tc.verifyTrue(ismember('SubjectID', vars));
+ tc.verifyTrue(ismember('Group', vars));
+ tc.verifyTrue(ismember('Condition', vars));
+ tc.verifyTrue(ismember('reactionTime', vars));
+ tc.verifyTrue(ismember('accuracy', vars));
+ end
+
+ %% --- Selection ---
+
+ function testSelectByString(tc)
+ ex = exploreFNIRS.core.Experiment(tc.allData);
+ ex.select('Group', 'Ctrl');
+ sel = ex.getSelectedData();
+ tc.verifyEqual(length(sel), 8);
+ end
+
+ function testSelectByMultipleValues(tc)
+ ex = exploreFNIRS.core.Experiment(tc.allData);
+ ex.select('Condition', {'Easy', 'Hard'});
+ sel = ex.getSelectedData();
+ tc.verifyEqual(length(sel), tc.nSegments);
+ end
+
+ function testSelectNarrows(tc)
+ ex = exploreFNIRS.core.Experiment(tc.allData);
+ ex.select('Group', 'Ctrl');
+ ex.select('Condition', 'Easy');
+ sel = ex.getSelectedData();
+ tc.verifyEqual(length(sel), 4);
+ end
+
+ function testSelectInvalidVar(tc)
+ ex = exploreFNIRS.core.Experiment(tc.allData);
+ tc.verifyError(@() ex.select('Nonexistent', 'foo'), ...
+ 'exploreFNIRS:core:Experiment:select');
+ end
+
+ function testReset(tc)
+ ex = exploreFNIRS.core.Experiment(tc.allData);
+ ex.select('Group', 'Tx');
+ ex.reset();
+ sel = ex.getSelectedData();
+ tc.verifyEqual(length(sel), tc.nSegments);
+ end
+
+ %% --- Groupby ---
+
+ function testGroupbySingleVar(tc)
+ ex = exploreFNIRS.core.Experiment(tc.allData);
+ ex.groupby('Condition');
+ tc.verifyTrue(ex.isGrouped);
+ tc.verifyLength(ex.groups, 2);
+ end
+
+ function testGroupbyMultipleVars(tc)
+ ex = exploreFNIRS.core.Experiment(tc.allData);
+ ex.groupby({'Group', 'Condition'});
+ tc.verifyTrue(ex.isGrouped);
+ % Ctrl x {Easy,Hard} + Tx x {Easy,Hard} = 4 groups
+ tc.verifyLength(ex.groups, 4);
+ end
+
+ function testGroupbyWithSelection(tc)
+ ex = exploreFNIRS.core.Experiment(tc.allData);
+ ex.select('Condition', 'Easy');
+ ex.groupby('Group');
+ tc.verifyLength(ex.groups, 2);
+ % Each group should have only Easy segments
+ for g = 1:length(ex.groups)
+ conds = ex.groups(g).gbyTables.Condition;
+ tc.verifyTrue(all(conds == "Easy" | conds == 'Easy'));
+ end
+ end
+
+ function testGroupbyInvalidVar(tc)
+ ex = exploreFNIRS.core.Experiment(tc.allData);
+ tc.verifyError(@() ex.groupby('Nonexistent'), ...
+ 'exploreFNIRS:core:Experiment:groupby');
+ end
+
+ function testGroupbyHasLabel(tc)
+ ex = exploreFNIRS.core.Experiment(tc.allData);
+ ex.groupby('Group');
+ for g = 1:length(ex.groups)
+ tc.verifyNotEmpty(ex.groups(g).label);
+ end
+ end
+
+ %% --- Aggregate ---
+
+ function testAggregateRequiresGroupby(tc)
+ ex = exploreFNIRS.core.Experiment(tc.allData);
+ tc.verifyError(@() ex.aggregate(), ...
+ 'exploreFNIRS:core:Experiment:aggregate');
+ end
+
+ function testAggregateProducesGrandAverage(tc)
+ ex = exploreFNIRS.core.Experiment(tc.allData);
+ ex.groupby('Group');
+ ex.settings.resampleRate = 2;
+ ex.settings.useBaseline = false;
+ ex.aggregate();
+
+ tc.verifyTrue(ex.isAggregated);
+ for g = 1:length(ex.groups)
+ ga = ex.groups(g).gbyGrand;
+ tc.verifyNotEmpty(ga);
+ tc.verifyTrue(isfield(ga, 'HbO'));
+ tc.verifyTrue(isfield(ga, 'HbR'));
+ tc.verifyTrue(isfield(ga, 'time'));
+ tc.verifyTrue(isfield(ga.HbO, 'Mean'));
+ tc.verifyTrue(isfield(ga.HbO, 'SEM'));
+ end
+ end
+
+ function testAggregateWithPreprocessing(tc)
+ ex = exploreFNIRS.core.Experiment(tc.allData);
+ ex.groupby('Condition');
+ ex.settings.baseline = [-5, 0];
+ ex.settings.resampleRate = 1;
+ ex.settings.useBaseline = true;
+ ex.aggregate();
+
+ tc.verifyTrue(ex.isAggregated);
+ ga = ex.groups(1).gbyGrand;
+ tc.verifyNotEmpty(ga.time);
+ tc.verifyNotEmpty(ga.HbO.Mean);
+ end
+
+ function testAggregateBarFlat(tc)
+ ex = exploreFNIRS.core.Experiment(tc.allData);
+ ex.groupby('Group');
+ ex.settings.resampleRate = 1;
+ ex.settings.barBinSize = 10;
+ ex.settings.useBaseline = false;
+ ex.aggregate();
+
+ for g = 1:length(ex.groups)
+ tc.verifyNotEmpty(ex.groups(g).gbyGrandBarFlat);
+ tc.verifyTrue(isfield(ex.groups(g).gbyGrandBarFlat, 'HbO'));
+ end
+ end
+
+ %% --- Info Variable Plotting ---
+
+ function testPlotInfoVarBasic(tc)
+ ex = exploreFNIRS.core.Experiment(tc.allData);
+ ex.groupby('Condition');
+ fig = ex.plotInfoVar('reactionTime', 'Visible', 'off');
+ tc.verifyClass(fig, 'matlab.ui.Figure');
+ close(fig);
+ end
+
+ function testPlotInfoVarSavesFile(tc)
+ outPath = fullfile(tempdir, 'test_plotInfoVar.png');
+ if exist(outPath, 'file'), delete(outPath); end
+
+ ex = exploreFNIRS.core.Experiment(tc.allData);
+ ex.groupby('Condition');
+ fig = ex.plotInfoVar('reactionTime', ...
+ 'SavePath', outPath, 'Visible', 'off');
+ close(fig);
+
+ tc.verifyTrue(exist(outPath, 'file') > 0);
+ delete(outPath);
+ end
+
+ function testPlotInfoVarRequiresGroupby(tc)
+ ex = exploreFNIRS.core.Experiment(tc.allData);
+ tc.verifyError(@() ex.plotInfoVar('reactionTime'), ...
+ 'exploreFNIRS:core:Experiment:plotInfoVar');
+ end
+
+ function testPlotInfoVarRejectsNonNumeric(tc)
+ ex = exploreFNIRS.core.Experiment(tc.allData);
+ ex.groupby('Group');
+ tc.verifyError(@() ex.plotInfoVar('SubjectID', 'Visible', 'off'), ...
+ 'exploreFNIRS:core:Experiment:plotInfoVar');
+ end
+
+ function testPlotInfoVarMultiGroup(tc)
+ ex = exploreFNIRS.core.Experiment(tc.allData);
+ ex.groupby({'Group', 'Condition'});
+ fig = ex.plotInfoVar('reactionTime', 'Visible', 'off');
+ tc.verifyClass(fig, 'matlab.ui.Figure');
+ close(fig);
+ end
+
+ function testPlotInfoVarNoAggregateNeeded(tc)
+ % plotInfoVar should work without calling aggregate()
+ ex = exploreFNIRS.core.Experiment(tc.allData);
+ ex.groupby('Condition');
+ % Do NOT call aggregate
+ tc.verifyFalse(ex.isAggregated);
+ fig = ex.plotInfoVar('accuracy', 'Visible', 'off');
+ tc.verifyClass(fig, 'matlab.ui.Figure');
+ close(fig);
+ end
+
+ %% --- Scatter Plot ---
+
+ function testPlotScatterBasic(tc)
+ ex = exploreFNIRS.core.Experiment(tc.allData);
+ ex.groupby('Condition');
+ fig = ex.plotScatter('reactionTime', 'accuracy', 'Visible', 'off');
+ tc.verifyClass(fig, 'matlab.ui.Figure');
+ close(fig);
+ end
+
+ function testPlotScatterWithFitLine(tc)
+ ex = exploreFNIRS.core.Experiment(tc.allData);
+ ex.groupby('Group');
+ fig = ex.plotScatter('Age', 'reactionTime', ...
+ 'FitLine', true, 'Visible', 'off');
+ tc.verifyClass(fig, 'matlab.ui.Figure');
+ close(fig);
+ end
+
+ function testPlotScatterNoGrouping(tc)
+ ex = exploreFNIRS.core.Experiment(tc.allData);
+ % No groupby - should still work with single color
+ fig = ex.plotScatter('reactionTime', 'accuracy', 'Visible', 'off');
+ tc.verifyClass(fig, 'matlab.ui.Figure');
+ close(fig);
+ end
+
+ function testPlotScatterSavesFile(tc)
+ outPath = fullfile(tempdir, 'test_scatter.png');
+ if exist(outPath, 'file'), delete(outPath); end
+
+ ex = exploreFNIRS.core.Experiment(tc.allData);
+ fig = ex.plotScatter('reactionTime', 'accuracy', ...
+ 'SavePath', outPath, 'Visible', 'off');
+ close(fig);
+
+ tc.verifyTrue(exist(outPath, 'file') > 0);
+ delete(outPath);
+ end
+
+ function testPlotScatterRejectsNonNumeric(tc)
+ ex = exploreFNIRS.core.Experiment(tc.allData);
+ tc.verifyError(@() ex.plotScatter('SubjectID', 'accuracy', 'Visible', 'off'), ...
+ 'exploreFNIRS:core:Experiment:plotScatter');
+ end
+
+ %% --- InfoTable ---
+
+ function testInfoTable(tc)
+ ex = exploreFNIRS.core.Experiment(tc.allData);
+ T = ex.infoTable();
+ tc.verifyEqual(height(T), tc.nSegments);
+ tc.verifyTrue(ismember('reactionTime', T.Properties.VariableNames));
+ end
+
+ function testInfoTableRespectsSelection(tc)
+ ex = exploreFNIRS.core.Experiment(tc.allData);
+ ex.select('Group', 'Tx');
+ T = ex.infoTable();
+ tc.verifyEqual(height(T), 4);
+ end
+
+ %% --- Temporal Plot ---
+
+ function testPlotTemporal(tc)
+ ex = exploreFNIRS.core.Experiment(tc.allData);
+ ex.groupby('Group');
+ ex.settings.resampleRate = 2;
+ ex.settings.useBaseline = false;
+ ex.aggregate();
+
+ fig = ex.plotTemporal('Biomarkers', {'HbO'}, 'Channels', 1, ...
+ 'Visible', 'off');
+ tc.verifyClass(fig, 'matlab.ui.Figure');
+ close(fig);
+ end
+
+ function testPlotTemporalRequiresAggregate(tc)
+ ex = exploreFNIRS.core.Experiment(tc.allData);
+ ex.groupby('Group');
+ tc.verifyError(@() ex.plotTemporal('Visible', 'off'), ...
+ 'exploreFNIRS:core:Experiment:plotTemporal');
+ end
+
+ %% --- Bar Plot ---
+
+ function testPlotBar(tc)
+ ex = exploreFNIRS.core.Experiment(tc.allData);
+ ex.groupby('Condition');
+ ex.settings.resampleRate = 2;
+ ex.settings.useBaseline = false;
+ ex.aggregate();
+
+ fig = ex.plotBar('Biomarker', 'HbO', 'Channels', 1:3, ...
+ 'Visible', 'off');
+ tc.verifyClass(fig, 'matlab.ui.Figure');
+ close(fig);
+ end
+
+ %% --- Aux Plot ---
+
+ function testPlotAuxSingleChannel(tc)
+ ex = exploreFNIRS.core.Experiment(tc.allData);
+ ex.groupby('Condition');
+ ex.settings.resampleRate = 2;
+ ex.settings.useBaseline = false;
+ ex.aggregate();
+
+ fig = ex.plotAux('heartRate', 'Visible', 'off');
+ tc.verifyClass(fig, 'matlab.ui.Figure');
+ close(fig);
+ end
+
+ function testPlotAuxMultiChannel(tc)
+ ex = exploreFNIRS.core.Experiment(tc.allData);
+ ex.groupby('Group');
+ ex.settings.resampleRate = 2;
+ ex.settings.useBaseline = false;
+ ex.aggregate();
+
+ fig = ex.plotAux('accelerometer', 'Layout', 'grid', 'Visible', 'off');
+ tc.verifyClass(fig, 'matlab.ui.Figure');
+ close(fig);
+ end
+
+ function testPlotAuxOverlay(tc)
+ ex = exploreFNIRS.core.Experiment(tc.allData);
+ ex.groupby('Condition');
+ ex.settings.resampleRate = 2;
+ ex.settings.useBaseline = false;
+ ex.aggregate();
+
+ fig = ex.plotAux('accelerometer', 'Layout', 'overlay', 'Visible', 'off');
+ tc.verifyClass(fig, 'matlab.ui.Figure');
+ close(fig);
+ end
+
+ function testPlotAuxSelectChannels(tc)
+ ex = exploreFNIRS.core.Experiment(tc.allData);
+ ex.groupby('Condition');
+ ex.settings.resampleRate = 2;
+ ex.settings.useBaseline = false;
+ ex.aggregate();
+
+ fig = ex.plotAux('accelerometer', 'AuxChannels', [1, 3], 'Visible', 'off');
+ tc.verifyClass(fig, 'matlab.ui.Figure');
+ close(fig);
+ end
+
+ function testPlotAuxRequiresAggregate(tc)
+ ex = exploreFNIRS.core.Experiment(tc.allData);
+ ex.groupby('Group');
+ tc.verifyError(@() ex.plotAux('heartRate'), ...
+ 'exploreFNIRS:core:Experiment:plotAux');
+ end
+
+ function testPlotAuxSavesFile(tc)
+ outPath = fullfile(tempdir, 'test_aux.png');
+ if exist(outPath, 'file'), delete(outPath); end
+
+ ex = exploreFNIRS.core.Experiment(tc.allData);
+ ex.groupby('Condition');
+ ex.settings.resampleRate = 2;
+ ex.settings.useBaseline = false;
+ ex.aggregate();
+
+ fig = ex.plotAux('heartRate', 'SavePath', outPath, 'Visible', 'off');
+ close(fig);
+ tc.verifyTrue(exist(outPath, 'file') > 0);
+ delete(outPath);
+ end
+
+ %% --- AuxFields ---
+
+ function testAuxFieldsReturnsFields(tc)
+ ex = exploreFNIRS.core.Experiment(tc.allData);
+ ex.groupby('Group');
+ ex.settings.resampleRate = 2;
+ ex.settings.useBaseline = false;
+ ex.aggregate();
+
+ flds = ex.auxFields();
+ tc.verifyTrue(iscell(flds));
+ tc.verifyTrue(ismember('accelerometer', flds) || ismember('heartRate', flds));
+ end
+
+ function testAuxFieldsBeforeAggregate(tc)
+ ex = exploreFNIRS.core.Experiment(tc.allData);
+ flds = ex.auxFields();
+ tc.verifyTrue(isempty(flds));
+ end
+
+ %% --- Export ---
+
+ function testToLongTable(tc)
+ ex = exploreFNIRS.core.Experiment(tc.allData);
+ ex.groupby('Group');
+ ex.settings.resampleRate = 1;
+ ex.settings.barBinSize = 10;
+ ex.settings.useBaseline = false;
+ ex.aggregate();
+
+ T = ex.toLongTable({'HbO'}, 1:3);
+ tc.verifyClass(T, 'table');
+ tc.verifyGreaterThan(height(T), 0);
+ end
+
+ function testToLongTableWithAux(tc)
+ ex = exploreFNIRS.core.Experiment(tc.allData);
+ ex.groupby('Group');
+ ex.settings.resampleRate = 1;
+ ex.settings.barBinSize = 10;
+ ex.settings.useBaseline = false;
+ ex.aggregate();
+
+ T = ex.toLongTable({'HbO'}, 1:3, [], 'IncludeAux', true);
+ tc.verifyClass(T, 'table');
+ tc.verifyGreaterThan(height(T), 0);
+ % Check for aux columns
+ vars = T.Properties.VariableNames;
+ hasAux = any(startsWith(vars, 'aux_'));
+ tc.verifyTrue(hasAux, 'Expected aux_ columns in long table');
+ end
+
+ function testToWideTable(tc)
+ ex = exploreFNIRS.core.Experiment(tc.allData);
+ ex.groupby('Group');
+ ex.settings.resampleRate = 1;
+ ex.settings.barBinSize = 10;
+ ex.settings.useBaseline = false;
+ ex.aggregate();
+
+ T = ex.toWideTable({'HbO'}, 1:3);
+ tc.verifyClass(T, 'table');
+ tc.verifyGreaterThan(height(T), 0);
+ end
+
+ function testToWideTableWithAux(tc)
+ ex = exploreFNIRS.core.Experiment(tc.allData);
+ ex.groupby('Group');
+ ex.settings.resampleRate = 1;
+ ex.settings.barBinSize = 10;
+ ex.settings.useBaseline = false;
+ ex.aggregate();
+
+ T = ex.toWideTable({'HbO'}, 1:3, [], 'IncludeAux', true);
+ tc.verifyClass(T, 'table');
+ tc.verifyGreaterThan(height(T), 0);
+ vars = T.Properties.VariableNames;
+ hasAux = any(startsWith(vars, 'aux_'));
+ tc.verifyTrue(hasAux, 'Expected aux_ columns in wide table');
+ end
+
+ function testExportRequiresAggregate(tc)
+ ex = exploreFNIRS.core.Experiment(tc.allData);
+ ex.groupby('Group');
+ tc.verifyError(@() ex.toLongTable(), ...
+ 'exploreFNIRS:core:Experiment:toLongTable');
+ tc.verifyError(@() ex.toWideTable(), ...
+ 'exploreFNIRS:core:Experiment:toWideTable');
+ end
+
+ %% --- Summary ---
+
+ function testSummaryRuns(tc)
+ ex = exploreFNIRS.core.Experiment(tc.allData);
+ ex.groupby('Condition');
+ ex.summary(); % should not error
+ end
+
+ %% --- getGroupColors ---
+
+ function testGetGroupColorsSmall(tc)
+ colors = exploreFNIRS.core.getGroupColors(3);
+ tc.verifySize(colors, [3, 3]);
+ tc.verifyTrue(all(colors(:) >= 0 & colors(:) <= 1));
+ end
+
+ function testGetGroupColorsLarge(tc)
+ colors = exploreFNIRS.core.getGroupColors(20);
+ tc.verifySize(colors, [20, 3]);
+ end
+
+ %% --- ScatterFNIRS ---
+
+ function testScatterFNIRSBasic(tc)
+ ex = exploreFNIRS.core.Experiment(tc.allData);
+ ex.groupby('Condition');
+ ex.settings.resampleRate = 2;
+ ex.settings.barBinSize = 10;
+ ex.settings.useBaseline = false;
+ ex.aggregate();
+
+ [fig, stats] = ex.plotScatterFNIRS('reactionTime', ...
+ 'Biomarkers', {'HbO'}, 'Channels', 1, 'Visible', 'off');
+ tc.verifyClass(fig, 'matlab.ui.Figure');
+ tc.verifyNotEmpty(stats);
+ close(fig);
+ end
+
+ function testScatterFNIRSSpearman(tc)
+ ex = exploreFNIRS.core.Experiment(tc.allData);
+ ex.groupby('Group');
+ ex.settings.resampleRate = 2;
+ ex.settings.barBinSize = 10;
+ ex.settings.useBaseline = false;
+ ex.aggregate();
+
+ [fig, stats] = ex.plotScatterFNIRS('Age', ...
+ 'CorrType', 'Spearman', 'Channels', 1, 'Visible', 'off');
+ tc.verifyClass(fig, 'matlab.ui.Figure');
+ close(fig);
+ end
+
+ function testScatterFNIRSMultiChannel(tc)
+ ex = exploreFNIRS.core.Experiment(tc.allData);
+ ex.groupby('Condition');
+ ex.settings.resampleRate = 2;
+ ex.settings.barBinSize = 10;
+ ex.settings.useBaseline = false;
+ ex.aggregate();
+
+ [fig, ~] = ex.plotScatterFNIRS('reactionTime', ...
+ 'Channels', [1, 2, 3], 'Visible', 'off');
+ tc.verifyClass(fig, 'matlab.ui.Figure');
+ close(fig);
+ end
+
+ function testScatterFNIRSSavesFile(tc)
+ outPath = fullfile(tempdir, 'test_scatter_fnirs.png');
+ if exist(outPath, 'file'), delete(outPath); end
+
+ ex = exploreFNIRS.core.Experiment(tc.allData);
+ ex.groupby('Condition');
+ ex.settings.resampleRate = 2;
+ ex.settings.barBinSize = 10;
+ ex.settings.useBaseline = false;
+ ex.aggregate();
+
+ [fig, ~] = ex.plotScatterFNIRS('reactionTime', ...
+ 'Channels', 1, 'SavePath', outPath, 'Visible', 'off');
+ close(fig);
+ tc.verifyTrue(exist(outPath, 'file') > 0);
+ delete(outPath);
+ end
+
+ function testScatterFNIRSRequiresAggregate(tc)
+ ex = exploreFNIRS.core.Experiment(tc.allData);
+ ex.groupby('Condition');
+ tc.verifyError(@() ex.plotScatterFNIRS('reactionTime'), ...
+ 'exploreFNIRS:core:Experiment:plotScatterFNIRS');
+ end
+
+ function testScatterFNIRSWithFitLine(tc)
+ ex = exploreFNIRS.core.Experiment(tc.allData);
+ ex.groupby('Group');
+ ex.settings.resampleRate = 2;
+ ex.settings.barBinSize = 10;
+ ex.settings.useBaseline = false;
+ ex.aggregate();
+
+ [fig, stats] = ex.plotScatterFNIRS('Age', ...
+ 'FitLine', true, 'Channels', 1, 'Visible', 'off');
+ tc.verifyClass(fig, 'matlab.ui.Figure');
+ close(fig);
+ end
+
+ %% --- LME ---
+
+ function testLMEBasic(tc)
+ ex = exploreFNIRS.core.Experiment(tc.allData);
+ ex.groupby('Condition');
+ ex.settings.resampleRate = 2;
+ ex.settings.barBinSize = 10;
+ ex.settings.useBaseline = false;
+ ex.aggregate();
+
+ [fig, results] = ex.plotLME('Biomarkers', {'HbO'}, ...
+ 'Channels', 1, 'Visible', 'off');
+ tc.verifyNotEmpty(results);
+ tc.verifyNotEmpty(results.formula);
+ if ~isempty(fig)
+ close(fig);
+ end
+ end
+
+ function testLMEReturnsModel(tc)
+ ex = exploreFNIRS.core.Experiment(tc.allData);
+ ex.groupby('Condition');
+ ex.settings.resampleRate = 2;
+ ex.settings.barBinSize = 10;
+ ex.settings.useBaseline = false;
+ ex.aggregate();
+
+ [fig, results] = ex.plotLME('Biomarkers', {'HbO'}, ...
+ 'Channels', 1, 'Visible', 'off', 'ShowBar', false);
+ tc.verifyNotEmpty(results.models);
+ mdl = results.models{1, 1};
+ if ~isempty(mdl)
+ tc.verifyClass(mdl, 'LinearMixedModel');
+ end
+ if ~isempty(fig), close(fig); end
+ end
+
+ function testLMEMultiChannel(tc)
+ ex = exploreFNIRS.core.Experiment(tc.allData);
+ ex.groupby('Condition');
+ ex.settings.resampleRate = 2;
+ ex.settings.barBinSize = 10;
+ ex.settings.useBaseline = false;
+ ex.aggregate();
+
+ [fig, results] = ex.plotLME('Biomarkers', {'HbO'}, ...
+ 'Channels', [1, 2], 'Visible', 'off');
+ tc.verifyTrue(size(results.models, 2) >= 2);
+ if ~isempty(fig), close(fig); end
+ end
+
+ function testLMEAnovaTable(tc)
+ ex = exploreFNIRS.core.Experiment(tc.allData);
+ ex.groupby('Condition');
+ ex.settings.resampleRate = 2;
+ ex.settings.barBinSize = 10;
+ ex.settings.useBaseline = false;
+ ex.aggregate();
+
+ [fig, results] = ex.plotLME('Biomarkers', {'HbO'}, ...
+ 'Channels', 1, 'Visible', 'off', 'ShowBar', false);
+ tc.verifyClass(results.anova_pval, 'table');
+ tc.verifyGreaterThan(height(results.anova_pval), 0);
+ if ~isempty(fig), close(fig); end
+ end
+
+ function testLMERequiresAggregate(tc)
+ ex = exploreFNIRS.core.Experiment(tc.allData);
+ ex.groupby('Condition');
+ tc.verifyError(@() ex.plotLME(), ...
+ 'exploreFNIRS:core:Experiment:plotLME');
+ end
+
+ function testLMECustomFormula(tc)
+ ex = exploreFNIRS.core.Experiment(tc.allData);
+ ex.groupby('Condition');
+ ex.settings.resampleRate = 2;
+ ex.settings.barBinSize = 10;
+ ex.settings.useBaseline = false;
+ ex.aggregate();
+
+ [fig, results] = ex.plotLME('Biomarkers', {'HbO'}, ...
+ 'Channels', 1, 'Visible', 'off', 'ShowBar', false, ...
+ 'CustomFormula', 'Opt1_HbO ~ Condition + (1|SubjectID)');
+ tc.verifyNotEmpty(results.formula);
+ tc.verifyTrue(contains(results.formula, 'Condition'));
+ if ~isempty(fig), close(fig); end
+ end
+
+ end
+end
diff --git a/+pf2_base/ProcessingContext.m b/+pf2_base/ProcessingContext.m
index 0602a3fa..0b069077 100644
--- a/+pf2_base/ProcessingContext.m
+++ b/+pf2_base/ProcessingContext.m
@@ -216,6 +216,8 @@ function applyToGlobals(obj)
% save('settings.mat', '-struct', 's');
s = struct();
+ s.contextVersion = '1.0';
+ s.created = char(datetime('now', 'Format', 'yyyy-MM-dd''T''HH:mm:ss'));
s.dpfMode = obj.dpfMode;
s.dpfFixedValue = obj.dpfFixedValue;
s.subjectAge = obj.subjectAge;
@@ -367,6 +369,7 @@ function applyToGlobals(obj)
function mustBeInRange(value, minVal, maxVal)
% Custom validator for rejectLevel
if value < minVal || value > maxVal
- error('Value must be between %g and %g', minVal, maxVal);
+ error('MATLAB:validators:mustBeInRange', ...
+ 'Value must be between %g and %g', minVal, maxVal);
end
end
diff --git a/+pf2_base/applyLightTheme.m b/+pf2_base/applyLightTheme.m
new file mode 100644
index 00000000..c051c4d8
--- /dev/null
+++ b/+pf2_base/applyLightTheme.m
@@ -0,0 +1,72 @@
+function applyLightTheme(fig)
+% applyLightTheme Force light-mode colors on MATLAB GUIDE figures.
+%
+% pf2_base.applyLightTheme(fig) walks all uicontrol and uipanel children
+% of figure FIG and sets explicit foreground/background colors so that the
+% GUI remains readable on macOS dark mode.
+%
+% Only runs on macOS; returns immediately on other platforms.
+%
+% Usage:
+% % In any GUIDE OpeningFcn, after guidata(hObject, handles):
+% pf2_base.applyLightTheme(hObject);
+
+ if ~ismac
+ return;
+ end
+
+ bgLight = [0.94 0.94 0.94]; % standard MATLAB light gray
+ bgWhite = [1 1 1];
+ fgBlack = [0 0 0];
+ fgPanel = [0.2 0.2 0.2]; % slightly softer than pure black
+
+ % --- Figure ---
+ set(fig, 'Color', bgLight);
+
+ % --- Panels ---
+ panels = findall(fig, 'Type', 'uipanel');
+ for i = 1:numel(panels)
+ set(panels(i), 'BackgroundColor', bgLight, ...
+ 'ForegroundColor', fgPanel);
+ end
+
+ % --- Button groups ---
+ btnGroups = findall(fig, 'Type', 'uibuttongroup');
+ for i = 1:numel(btnGroups)
+ set(btnGroups(i), 'BackgroundColor', bgLight, ...
+ 'ForegroundColor', fgPanel);
+ end
+
+ % --- Controls ---
+ controls = findall(fig, 'Type', 'uicontrol');
+ for i = 1:numel(controls)
+ style = get(controls(i), 'Style');
+ switch style
+ case {'listbox', 'edit', 'popupmenu'}
+ set(controls(i), 'BackgroundColor', bgWhite, ...
+ 'ForegroundColor', fgBlack);
+ case {'text'}
+ set(controls(i), 'BackgroundColor', bgLight, ...
+ 'ForegroundColor', fgBlack);
+ case {'checkbox', 'radiobutton'}
+ set(controls(i), 'BackgroundColor', bgLight, ...
+ 'ForegroundColor', fgBlack);
+ case {'pushbutton', 'togglebutton'}
+ % Skip color-swatch buttons whose ForegroundColor is
+ % intentionally set to display a color (e.g. gui_color_N).
+ tag = get(controls(i), 'Tag');
+ if contains(tag, 'gui_color')
+ continue;
+ end
+ set(controls(i), 'BackgroundColor', bgLight, ...
+ 'ForegroundColor', fgBlack);
+ end
+ end
+
+ % --- Tables ---
+ tables = findall(fig, 'Type', 'uitable');
+ for i = 1:numel(tables)
+ set(tables(i), 'ForegroundColor', fgBlack, ...
+ 'BackgroundColor', bgWhite);
+ end
+end
diff --git a/+pf2_base/normalizeMarkers.m b/+pf2_base/normalizeMarkers.m
new file mode 100644
index 00000000..43651c13
--- /dev/null
+++ b/+pf2_base/normalizeMarkers.m
@@ -0,0 +1,43 @@
+function mrk = normalizeMarkers(mrk)
+% NORMALIZEMARKERS Pad marker array to 4 columns [time, value, duration, amplitude]
+%
+% Ensures marker arrays always have 4 columns. Missing columns are filled
+% with defaults: duration = 0, amplitude = 1. This standardizes the marker
+% format across all import functions and data generators.
+%
+% Syntax:
+% mrk = pf2_base.normalizeMarkers(mrk)
+%
+% Inputs:
+% mrk - Marker array in any of these formats:
+% [M x 2] - [time, value] (duration and amplitude added)
+% [M x 3] - [time, value, duration] (amplitude added)
+% [M x 4] - [time, value, duration, amplitude] (returned as-is)
+% [] - Empty input (returns zeros(0,4))
+%
+% Outputs:
+% mrk - Normalized marker array [M x 4]
+% Column 1: time (seconds)
+% Column 2: marker value/code
+% Column 3: duration (seconds), default 0
+% Column 4: amplitude/weight, default 1
+
+if isempty(mrk)
+ mrk = zeros(0, 4);
+ return;
+end
+
+nCols = size(mrk, 2);
+nRows = size(mrk, 1);
+
+if nCols < 3
+ % Add duration column (default 0)
+ mrk(:, 3) = zeros(nRows, 1);
+end
+
+if nCols < 4
+ % Add amplitude column (default 1)
+ mrk(:, 4) = ones(nRows, 1);
+end
+
+end
diff --git a/+pf2_base/pf2_plotArranged.m b/+pf2_base/pf2_plotArranged.m
index e2a7fda5..07fc577c 100644
--- a/+pf2_base/pf2_plotArranged.m
+++ b/+pf2_base/pf2_plotArranged.m
@@ -1,13 +1,41 @@
function [ figHandle ] = pf2_plotArranged(varargin)
-% pf2_plotArranged(fNIR,channels,showMarkers,wavelengths,ylimit,plotArranged,lineProps,rejectedLineProps)
-% This function plots and automatically arranges fNIRS data based on
-% the input values. Specify if raw or Oxy (or other) biomarkers
-% It expects fNIR (the data struct), the channels to plot:
-% (can be a number or a logical index)
-% a set of wavelengths (ie: 730 or logical index)
-% A specific ylimit ( to force all to use the same axes
-% a boolean to plot the arranged channels
-% Accepted line properties and for normal and rejected channels
+% PF2_PLOTARRANGED Plot fNIRS channels in device-arranged subplot layout
+%
+% Plots raw fNIRS data with channels arranged spatially according to the
+% device probe layout. Supports multi-probe devices, event markers,
+% wavelength selection, and visual indication of rejected channels.
+%
+% Syntax:
+% pf2_plotArranged(fNIR)
+% pf2_plotArranged(fNIR, channels)
+% pf2_plotArranged(fNIR, channels, showMarkers, signalNames, ylimit, ...
+% plotArranged, lineProps, rejectedLineProps, baseline, RawData)
+% figHandle = pf2_plotArranged(...)
+%
+% Inputs:
+% fNIR - fNIRS data struct (required)
+% channels - Channel indices to plot (default: all)
+% Numeric array or logical mask.
+% showMarkers - Display event markers (default: true)
+% signalNames - Signal/wavelength names to plot (default: all)
+% ylimit - Y-axis limits [min max] (default: auto from device)
+% plotArranged - Use spatial arrangement (default: false, true if all)
+% lineProps - Line properties for good channels (default: {'LineWidth', 1})
+% rejectedLineProps - Line properties for rejected channels (default: {'--', 'LineWidth', 1})
+% baseline - Show baseline indicator (default: false)
+% RawData - Plotting raw data flag (default: true)
+%
+% Outputs:
+% figHandle - Handle(s) to created figure(s)
+%
+% Example:
+% data = pf2.import.sampleData.fNIR2000();
+% pf2_base.pf2_plotArranged(data);
+%
+% % Plot specific channels
+% pf2_base.pf2_plotArranged(data, [1 3 5]);
+%
+% See also: pf2.data.plot.raw, pf2.data.plot.oxy
validFnirs = @(x) (iscell(x) || isstruct(x));
validChannels = @(x) (isnumeric(x) || ischar(x));
diff --git a/+pf2_base/pf2_unpackMethod.m b/+pf2_base/pf2_unpackMethod.m
index 139451b3..17a9bd16 100644
--- a/+pf2_base/pf2_unpackMethod.m
+++ b/+pf2_base/pf2_unpackMethod.m
@@ -91,7 +91,11 @@
for j=1:length(Fidx)
F_noarray.args{j}=Fidx(j).args;
F_noarray.argvals{j}=Fidx(j).argvals;
- F_noarray.default_argvals{j}=Fidx(j).default_argvals;
+ if isfield(Fidx, 'default_argvals')
+ F_noarray.default_argvals{j}=Fidx(j).default_argvals;
+ else
+ F_noarray.default_argvals{j}=Fidx(j).argvals;
+ end
if(isfield(Fidx(j),'output'))
F_noarray.output{j}=Fidx(j).output;
else
diff --git a/.gitignore b/.gitignore
index 31ca291e..3943f4ab 100644
--- a/.gitignore
+++ b/.gitignore
@@ -9,3 +9,5 @@ sampledata/sampleNIR_CH.mat
.vscode/settings.json
*.md
internal/
+benchmarks/data/
+benchmarks/results/
diff --git a/CHANGELOG.md b/CHANGELOG.md
index fb95c0b5..901171a3 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,72 @@
# Changelog
+## v1.0.0 (2026-02-06)
+Scriptable Group Analysis, Connectivity, Hyperscanning & Method CRUD
+
+### New Features
+
+**exploreFNIRS Scriptable API:**
+- **Experiment container class** (`+exploreFNIRS/+core/Experiment.m`) for fully scriptable group analysis without the GUI
+ - `select()`, `groupby()`, `aggregate()` methods for data organization
+ - `connectivity()` and `hyperscanning()` methods with block-wise support
+ - `plotTemporal()`, `plotBar()` wrappers that forward to headless plotting
+ - `exportLong()`, `exportWide()` for data export
+- **Headless plotting** (`+exploreFNIRS/+core/plotTemporal.m`, `plotBar.m`)
+ - Publication-ready temporal and bar chart plots without GUI
+ - ROI support via `'ROIs'` parameter (indices, names, or `'all'`)
+ - Configurable error bands (SEM/SD), layout (overlay/grid), and save options
+
+**Connectivity Analysis:**
+- **Connectivity module** (`+exploreFNIRS/+connectivity/`) with `computeMatrix`, `plotMatrix`, `plotBlockComparison`
+- **Coupling functions** (`+exploreFNIRS/+coupling/`) — Pearson, Spearman, cross-correlation, coherence, wavelet coherence
+- **Visualization** — `plotWcoherence` (time-freq heatmap), `plotWindowed` (windowed coupling time series)
+
+**Hyperscanning Analysis:**
+- **Hyperscanning module** (`+exploreFNIRS/+hyperscanning/`) with `pairSubjects`, `computeDyad`, `computeGroup`, `permutationTest`, `plotGroup`
+- Block-wise hyperscanning with struct array results
+
+**Block Definition & Extraction:**
+- `pf2.data.defineBlocks()` — convert markers to block struct array
+- `pf2.data.extractBlocks()` — extract fNIRS segments by block definitions
+- Positional API: `defineBlocks(data, [49,50], 30)` or name-value pairs
+
+**GLM & Short-Channel Regression:**
+- `pf2_base.fnirs.buildDesignMatrix()` — construct GLM design matrices from markers
+- `pf2_base.fnirs.fitGLM()` — fit general linear model to fNIRS data
+- `pf2_base.fnirs.shortChannelRegression()` — regress out short-channel signals
+- `functions/pf2_SSR.m` — short separation regression processing function
+
+**Method CRUD Operations:**
+- `pf2.methods.raw.create()` / `pf2.methods.oxy.create()` — create new methods programmatically
+- `pf2.methods.raw.delete()` / `pf2.methods.oxy.delete()` — delete methods
+- `pf2.methods.raw.editFunction()` / `removeFunction()` — modify method function chains
+- `pf2.methods.raw.exportMethod()` / `importMethod()` — portable method sharing
+- `pf2.methods.validateFunction()` — validate function compatibility
+
+**Other New Features:**
+- `pf2_base.normalizeMarkers()` — standardize marker codes across devices
+- `pf2_base.applyLightTheme()` — consistent light theme for figures
+- `+exploreFNIRS/+export/connectivityToTable.m` — export connectivity results as tables
+- `+exploreFNIRS/+core/getGroupColors.m` — consistent group coloring across plots
+
+### New Tests
+- **ConnectivityTest.m** — 31 tests covering all coupling functions, matrix computation, and block-wise connectivity
+- **BlockDefinitionTest.m** — 29 tests covering marker-to-block conversion and extraction
+- **GLMTest.m** — GLM design matrix and fitting tests
+- **SSRTest.m** — Short-channel regression tests
+- **HierarchicalAverageTest.m** — Hierarchical averaging validation
+- **NormalizeMarkersTest.m** — Marker normalization tests
+- **SplitTest.m** — Data splitting tests
+- **GoldenFileTest.m** — Regression testing with golden reference files
+- **testExperiment.m** — Experiment class integration tests
+- Test count: 225 → 300+ tests across 20+ test classes
+
+### Bug Fixes
+- Fixed `grandAvgFNIRS.m` crash when data lacks segmentTimes field
+- Fixed `plotTemporal.m` shadowing MATLAB builtins (`upper`/`lower` renamed to `upperBound`/`lowerBound`)
+
+---
+
## v0.9 (2026-01-23)
API Standardization, Testing Infrastructure & Context-Based Processing
diff --git a/ExploreFNIRS_README.md b/ExploreFNIRS_README.md
index 0d9508d7..3d8d302b 100644
--- a/ExploreFNIRS_README.md
+++ b/ExploreFNIRS_README.md
@@ -188,7 +188,65 @@ exploreFNIRS(allData);
% - Export statistics
```
-### Scriptable Analysis (Headless)
+### Scriptable Analysis (Headless) — New in v1.0.0
+
+The `Experiment` class enables complete group analysis without the GUI:
+
+```matlab
+% Create experiment from processed data
+ex = exploreFNIRS.core.Experiment(allData);
+
+% Filter and organize
+ex.select('Group', {'Control', 'Treatment'}, 'Condition', 'Task');
+ex.groupby({'Group', 'Condition'});
+ex.aggregate();
+
+% Headless temporal plot
+fig = ex.plotTemporal('Biomarkers', {'HbO'}, 'Channels', [1 5 10], ...
+ 'SavePath', 'temporal.png', 'SaveDPI', 300);
+
+% Headless bar chart with time window
+fig = ex.plotBar('Biomarker', 'HbO', 'TimeWindow', [5 25], ...
+ 'SavePath', 'bar.png');
+
+% ROI-based plotting
+fig = ex.plotTemporal('Biomarkers', {'HbO'}, 'ROIs', 'all');
+fig = ex.plotBar('Biomarker', 'HbO', 'ROIs', {'DLPFC_L', 'DLPFC_R'});
+
+% Export for external analysis
+longTable = ex.exportLong();
+wideTable = ex.exportWide();
+writetable(longTable, 'export_for_R.csv');
+```
+
+### Connectivity Analysis — New in v1.0.0
+
+```matlab
+% Compute connectivity matrices
+connResults = ex.connectivity('Method', 'pearson');
+fig = exploreFNIRS.connectivity.plotMatrix(connResults);
+
+% Block-wise connectivity
+connBlocks = ex.connectivity('Method', 'coherence', 'Blocks', blocks);
+fig = exploreFNIRS.connectivity.plotBlockComparison(connBlocks);
+
+% Export connectivity as table
+T = exploreFNIRS.export.connectivityToTable(connResults);
+```
+
+### Hyperscanning Analysis — New in v1.0.0
+
+```matlab
+% Pair subjects and compute inter-brain coupling
+hsResults = ex.hyperscanning('PairBy', 'Dyad', 'Method', 'wcoherence');
+
+% Group-level statistics with permutation testing
+groupStats = exploreFNIRS.hyperscanning.computeGroup(hsResults);
+pValues = exploreFNIRS.hyperscanning.permutationTest(hsResults, 1000);
+fig = exploreFNIRS.hyperscanning.plotGroup(groupStats);
+```
+
+### Other Scriptable Functions
```matlab
% Build segment info table
@@ -196,10 +254,6 @@ segmentTable = exploreFNIRS.dataset.buildSegmentInfoTable(allData);
% Standardize ROIs across subjects
allData = exploreFNIRS.dataset.standardizeROIs(allData);
-
-% Export for external analysis
-longData = exploreFNIRS.export.mergeGbyTablesLong(gbyData);
-writetable(longData, 'export_for_R.csv');
```
## Scriptable Functions
@@ -208,10 +262,24 @@ The following functions can be used outside the GUI:
| Package | Function | Purpose |
|---------|----------|---------|
+| `+core` | `Experiment` | Main experiment container class |
+| `+core` | `plotTemporal` | Headless temporal plots with ROI support |
+| `+core` | `plotBar` | Headless bar charts with ROI support |
+| `+core` | `getGroupColors` | Consistent group coloring |
+| `+connectivity` | `computeMatrix` | Channel-pair connectivity matrices |
+| `+connectivity` | `plotMatrix` | Matrix visualization |
+| `+connectivity` | `plotBlockComparison` | Block-wise comparison |
+| `+coupling` | `pearson`, `spearman`, `xcorr`, `coherence`, `wcoherence` | Coupling functions |
+| `+hyperscanning` | `pairSubjects` | Pair subjects by criteria |
+| `+hyperscanning` | `computeDyad` | Dyad-level coupling |
+| `+hyperscanning` | `computeGroup` | Group-level statistics |
+| `+hyperscanning` | `permutationTest` | Permutation significance testing |
+| `+hyperscanning` | `plotGroup` | Group visualization |
| `+dataset` | `buildSegmentInfoTable` | Create metadata table from structs |
| `+dataset` | `standardizeROIs` | Align ROI definitions across subjects |
| `+export` | `mergeGbyTablesLong` | Export to long format |
| `+export` | `mergeGbyTablesWide` | Export to wide format |
+| `+export` | `connectivityToTable` | Export connectivity results |
| `+fx` | `performFDR` | Benjamini-Hochberg FDR correction |
| `+fx` | `performFDR_twostep` | Adaptive two-step FDR |
| `+fx` | `autoContrast` | Generate post-hoc contrasts |
diff --git a/GUI/exploreFNIRS_browse.m b/GUI/exploreFNIRS_browse.m
index 773056f7..60a3549a 100644
--- a/GUI/exploreFNIRS_browse.m
+++ b/GUI/exploreFNIRS_browse.m
@@ -63,6 +63,8 @@ function exploreFNIRS_browse_OpeningFcn(hObject, eventdata, handles, varargin)
% Update handles structure
guidata(hObject, handles);
+pf2_base.applyLightTheme(hObject);
+
% UIWAIT makes exploreFNIRS_browse wait for user response (see UIRESUME)
% uiwait(handles.figure1);
global BrowseFNIRS
diff --git a/GUI/probeCheckGUI.m b/GUI/probeCheckGUI.m
index 2e643332..18ddbb6b 100644
--- a/GUI/probeCheckGUI.m
+++ b/GUI/probeCheckGUI.m
@@ -66,6 +66,8 @@ function probeCheckGUI_OpeningFcn(hObject, eventdata, handles, varargin)
% Update handles structure
guidata(hObject, handles);
+pf2_base.applyLightTheme(hObject);
+
global pf2ChannelCheck
if(isfield(pf2ChannelCheck,'autoscale'))
diff --git a/GUI/processFNIRS2_GUI.m b/GUI/processFNIRS2_GUI.m
index f58eb83c..8e7d48ba 100644
--- a/GUI/processFNIRS2_GUI.m
+++ b/GUI/processFNIRS2_GUI.m
@@ -60,6 +60,7 @@ function processFNIRS2_GUI_OpeningFcn(hObject, eventdata, handles, varargin)
global setF
global outputData
+pf2_base.applyLightTheme(hObject);
[~,pf2ver,~]=pf2_base.pf2version();
set(handles.uipanel_device_info,'Title',sprintf('ProcessFNIRS2 %s',pf2ver));
diff --git a/GUI/processFNIRS2_configureMethods.m b/GUI/processFNIRS2_configureMethods.m
index 12943580..38e665ba 100644
--- a/GUI/processFNIRS2_configureMethods.m
+++ b/GUI/processFNIRS2_configureMethods.m
@@ -55,6 +55,8 @@ function processFNIRS2_configureMethods_OpeningFcn(hObject, eventdata, handles,
% Choose default command line output for processFNIRS2_configureMethods
handles.output = hObject;
+pf2_base.applyLightTheme(hObject);
+
global configureMode
if(~isempty(varargin))
if(strcmp(varargin{1},'raw')||strcmp(varargin{1},'oxy'))
@@ -484,15 +486,15 @@ function pushbutton_rename_Callback(hObject, eventdata, handles)
case Numbers
case LowerCases
case UpperCases
- case {'À','?','Â','Ã','Ä','Å'}, Character = 'A';
+ case {'À','�?','Â','Ã','Ä','Å'}, Character = 'A';
case 'Æ', Character = 'AE';
case 'Ç', Character = 'C';
case {'È','É','Ê','Ë'}, Character = 'E';
- case {'Ì','?','Î','?'}, Character = 'I';
+ case {'Ì','�?','Î','�?'}, Character = 'I';
case 'Ñ', Character = 'N';
case {'Ò','Ó','Ô','Õ','Ö'}, Character = 'O';
case {'Ù','Ú','Û','Ü'}, Character = 'U';
- case '?', Character = 'Y';
+ case '�?', Character = 'Y';
case '²', Character = '2';
case '³', Character = '3';
case '¼', Character = '1_4';
@@ -1378,7 +1380,11 @@ function saveCurrentMethod()
for j=1:size(Fidx,2)
F_noarray.args{j}=Fidx(j).args;
F_noarray.argvals{j}=Fidx(j).argvals;
- F_noarray.default_argvals{j}=Fidx(j).default_argvals;
+ if isfield(Fidx, 'default_argvals')
+ F_noarray.default_argvals{j}=Fidx(j).default_argvals;
+ else
+ F_noarray.default_argvals{j}=Fidx(j).argvals;
+ end
if(isfield(Fidx(j),'output'))
F_noarray.output{1}=Fidx(1).output;
else
diff --git a/GUI/processFNIRS2_configureMethods_functionAddEdit.m b/GUI/processFNIRS2_configureMethods_functionAddEdit.m
index 9ef4d5ce..15f01902 100644
--- a/GUI/processFNIRS2_configureMethods_functionAddEdit.m
+++ b/GUI/processFNIRS2_configureMethods_functionAddEdit.m
@@ -55,6 +55,8 @@ function processFNIRS2_configureMethods_functionAddEdit_OpeningFcn(hObject, even
handles.output = hObject;
+pf2_base.applyLightTheme(hObject);
+
global compareMode
global curFunction
global outFunc
@@ -682,15 +684,15 @@ function checkbox_valid_raw_Callback(hObject, eventdata, handles)
case Numbers
case LowerCases
case UpperCases
- case {'À','?','Â','Ã','Ä','Å'}, Character = 'A';
+ case {'À','�?','Â','Ã','Ä','Å'}, Character = 'A';
case 'Æ', Character = 'AE';
case 'Ç', Character = 'C';
case {'È','É','Ê','Ë'}, Character = 'E';
- case {'Ì','?','Î','?'}, Character = 'I';
+ case {'Ì','�?','Î','�?'}, Character = 'I';
case 'Ñ', Character = 'N';
case {'Ò','Ó','Ô','Õ','Ö'}, Character = 'O';
case {'Ù','Ú','Û','Ü'}, Character = 'U';
- case '?', Character = 'Y';
+ case '�?', Character = 'Y';
case '²', Character = '2';
case '³', Character = '3';
case '¼', Character = '1_4';
diff --git a/OtherToolboxes/Wavelab850/InstallMEX.m b/OtherToolboxes/Wavelab850/InstallMEX.m
index cc62d338..6d688d53 100644
--- a/OtherToolboxes/Wavelab850/InstallMEX.m
+++ b/OtherToolboxes/Wavelab850/InstallMEX.m
@@ -30,6 +30,11 @@
% If not, install...
if ~MEX_OK,
disp('WaveLab detects that some or all of your MEX files are not installed,')
+ % Skip interactive prompt in batch/headless mode
+ if batchStartupOptionUsed
+ warning('WaveLab:MEXMissing', 'Some MEX files missing. Run InstallMEX interactively to compile.');
+ return;
+ end
R=input('do you want to install them now? [[Yes]/No] \n','s');
if strcmp(R,'') + strcmp(R,'Yes') | strcmp(R,'yes') | strcmp(R,'y') | strcmp(R,'Y') | strcmp(R,'YES'),
disp('INSTALLING MEX FILES, MAY TAKE A WHILE ...')
diff --git a/README.md b/README.md
index 35011126..d70bcfc3 100644
--- a/README.md
+++ b/README.md
@@ -1,4 +1,4 @@
-# processFNIRS2 v0.9
+# processFNIRS2 v1.0.0
## Overview
processFNIRS2 is a modular MATLAB toolbox designed for processing functional Near-Infrared Spectroscopy (fNIRS) data. The toolbox provides a flexible framework for importing, processing, analyzing, and visualizing fNIRS data from multiple device manufacturers.
@@ -59,7 +59,7 @@ exploreFNIRS(myprocesseddata);
## Processing Pipeline
### Data Import
-Use functions in the `pf2.Import` module to load data from various fNIRS devices:
+Use functions in the `pf2.import` module to load data from various fNIRS devices:
- `pf2.import.importNIR`: Import fNIR Devices/Biopac files
- `pf2.import.importHitachiMES`: Import Hitachi ETG-4000 files
- `pf2.import.importNIRX`: Import NIRx system files
@@ -74,16 +74,16 @@ The toolbox provides various functions for manipulating fNIRS data:
- `pf2.data.setT0`: Shift time to align with experiment start
- `pf2.data.split`: Split fNIRS segments based on time points
- `pf2.data.plot`: Visualize fNIRS data (Oxy, Raw, ROI, AuxData)
-- `pf2.Export`: Export data to NIR or SNIRF formats
+- `pf2.export`: Export data to NIR or SNIRF formats
### Method Configuration
processFNIRS2 uses a two-stage processing pipeline:
1. **Raw processing** (Raw → Optical Density)
- - Configure and select methods using `pf2.methods.Raw`
+ - Configure and select methods using `pf2.methods.raw`
- Common preprocessing includes: motion artifact correction, filtering, CAR, etc.
2. **Oxy processing** (Optical Density → Hemoglobin)
- - Configure and select methods using `pf2.methods.Oxy`
+ - Configure and select methods using `pf2.methods.oxy`
- Processing includes: Beer-Lambert conversion, filtering, ROI analysis, etc.
Methods can be configured through the GUI or programmatically:
@@ -93,12 +93,19 @@ pf2.methods.raw.configureMethods();
pf2.methods.oxy.configureMethods();
% List available methods
-pf2.methods.raw.List();
-pf2.methods.oxy.List();
+pf2.methods.raw.list();
+pf2.methods.oxy.list();
% Set methods programmatically
pf2.methods.raw.setMethod('MyRawMethod');
pf2.methods.oxy.setMethod('MyOxyMethod');
+
+% Create, modify, and share methods (new in v1.0.0)
+pf2.methods.raw.create('MyCustomMethod');
+pf2.methods.raw.editFunction('MyCustomMethod', 'pf2_lpf', struct('freq_cut', 0.08));
+pf2.methods.raw.exportMethod('MyCustomMethod', 'my_method.mat');
+pf2.methods.raw.importMethod('shared_method.mat');
+pf2.methods.raw.delete('OldMethod');
```
### Data Processing
@@ -178,15 +185,38 @@ allData = {subject1, subject2, subject3, ...};
% Launch exploreFNIRS GUI
exploreFNIRS(allData);
-% With configuration options
-exploreFNIRS(allData, 'timeShiftTo0', true, 'blStart', 0, 'blEnd', 5, ...
- 'blockStart', 5, 'blockEnd', 65, 'barSegmentLength', 60);
+% Or use the scriptable Experiment class (no GUI needed)
+ex = exploreFNIRS.core.Experiment(allData);
+ex.select('Group', {'Control', 'Treatment'});
+ex.groupby({'Group', 'Condition'});
+ex.aggregate();
+
+% Headless plots
+fig = ex.plotTemporal('Biomarkers', {'HbO'}, 'Channels', 1:5);
+fig = ex.plotBar('Biomarker', 'HbO', 'TimeWindow', [5, 25]);
+
+% ROI-based plotting
+fig = ex.plotTemporal('Biomarkers', {'HbO'}, 'ROIs', 'all');
+
+% Connectivity analysis
+connResults = ex.connectivity('Method', 'pearson');
+
+% Hyperscanning analysis
+hsResults = ex.hyperscanning('PairBy', 'Dyad', 'Method', 'coherence');
+
+% Export
+longTable = ex.exportLong();
```
exploreFNIRS features:
+- **Scriptable Experiment class** for complete headless group analysis
- Group-level analysis with hierarchical averaging
+- **Connectivity analysis** with 5 coupling methods (Pearson, Spearman, xcorr, coherence, wavelet coherence)
+- **Hyperscanning** with subject pairing, dyad/group computation, and permutation testing
+- **Block-wise analysis** for connectivity and hyperscanning
- Linear mixed-effects modeling with Satterthwaite degrees of freedom
- Visualization: temporal plots, bar charts, scatter plots, topographic maps
+- **Headless plotting** with ROI support (plotTemporal, plotBar)
- FDR correction: `exploreFNIRS.fx.performFDR()`
- Data export: `exploreFNIRS.export.mergeGbyTablesWide()` / `mergeGbyTablesLong()`
@@ -211,9 +241,9 @@ pf2.settings.selectDevice('fNIR_Devices_fNIR1200_16ch.cfg');
- `processFNIRS2.m`: Main function for processing fNIRS data
- `pf2.m`: Convenience wrapper for processFNIRS2
- `exploreFNIRS.m`: Group-level analysis GUI
-- `+pf2/`: User-facing API (Import, Export, Data, Methods, Settings, Probe)
-- `+pf2_base/`: Internal infrastructure and utilities
-- `+exploreFNIRS/`: Group analysis functions (plot, export, fx, dataset)
+- `+pf2/`: User-facing API (import, export, data, methods, settings, probe)
+- `+pf2_base/`: Internal infrastructure, utilities, and tests (300+ tests)
+- `+exploreFNIRS/`: Group analysis (core, connectivity, coupling, hyperscanning, plot, export, fx, dataset)
- `base_functions/`: Utility functions (legacy)
- `GUI/`: User interface components (legacy, GUIDE-based)
- `functions/`: Signal processing algorithms (filters, motion correction, etc.)
@@ -222,42 +252,42 @@ pf2.settings.selectDevice('fNIR_Devices_fNIR1200_16ch.cfg');
## Overall Structure
processFNIRS2 is laid out in the following manner:
-- **Data**: Functions to manipulate individual fNIRS segments
- - ApplyChannelMask: Set bad channels to nan
- - GetMarkers: Find timepoints of markers in a regex style
- - Resample: Resample or average fNIRS data
- - SetT0: Shift fNIRS time to match start of experiment
- - Split: Split fNIRS segment based on different input times
- - **Plot**: Functions to visualize fNIRS data
- - AuxData: Plot auxiliary data channels
- - Oxy: Plot oxygenation data
- - ROI: Plot Region of Interest data
- - Raw: Plot raw intensity data
- - **Export**: Functions to export fNIRS data
+- **data**: Functions to manipulate individual fNIRS segments
+ - applyChannelMask: Set bad channels to nan
+ - getMarkers: Find timepoints of markers in a regex style
+ - resample: Resample or average fNIRS data
+ - setT0: Shift fNIRS time to match start of experiment
+ - split: Split fNIRS segment based on different input times
+ - **plot**: Functions to visualize fNIRS data
+ - auxData: Plot auxiliary data channels
+ - oxy: Plot oxygenation data
+ - roi: Plot Region of Interest data
+ - raw: Plot raw intensity data
+ - **export**: Functions to export fNIRS data
- asNIR: Export to NIR file format
- asSNIRF: Export to SNIRF file format
-- **GUI**: Shortcut for accessing the GUI
-- **Help**: Access to help documentation
-- **Import**: Functions to import fNIRS files
+- **gui**: Shortcut for accessing the GUI
+- **help**: Access to help documentation
+- **import**: Functions to import fNIRS files
- importHitachiMES: Import Hitachi Probes
- importNIRX: Import NIRx files
- importNIR: Import fNIR Devices/Biopac files
- importSNIRF: Import SNIRF format files
- sampleData: Load sample data included with the toolbox
-- **Methods**: Functions to change and modify processing methods
- - Oxy: Oxy conversion pipeline methods
- - Raw: Raw domain pipeline methods
-- **Process**: Process fNIR segment data
+- **methods**: Functions to change and modify processing methods
+ - oxy: Oxy conversion pipeline methods
+ - raw: Raw domain pipeline methods
+- **process**: Process fNIR segment data
- processOxy: Run the Oxy Pipeline only
- processRaw: Run the Raw Pipeline only
-- **Settings**: Change settings related to processing
- - Baseline: Change baseline time settings
- - DPF: Change mode of Differential Path Length
+- **settings**: Change settings related to processing
+ - baseline: Change baseline time settings
+ - dpf: Change mode of Differential Path Length
- selectDevice: Reload device settings for FNIRS probe
## Troubleshooting Tips
- When importing data for the first time, verify that the probe configuration is correct
-- If you get errors about DPF factors, check the settings using `pf2.settings.DPF`
+- If you get errors about DPF factors, check the settings using `pf2.settings.dpf`
- For visualization issues, try running with default methods first
- If having trouble loading the software, check the MATLAB preference directory (`prefdir`) and delete any related settings files
- Remember that GUI settings are for visualization only and don't affect your data
@@ -270,7 +300,7 @@ Access this location using the MATLAB command: `prefdir`
For detailed function documentation, use MATLAB's `help` command:
```matlab
help processFNIRS2
-help pf2.methods.Raw
+help pf2.methods.raw
help pf2.import.importNIR
```
diff --git a/base_functions/grandAvgFNIRS.m b/base_functions/grandAvgFNIRS.m
index efdbdc18..e7991cbe 100644
--- a/base_functions/grandAvgFNIRS.m
+++ b/base_functions/grandAvgFNIRS.m
@@ -243,9 +243,9 @@
outGA.segmentTimes=unique(segmentTimesArr,'rows');
outGA.segmentTimes=sort(outGA.segmentTimes,1);
outGA.segmentTimes=round( outGA.segmentTimes,5);
+ outGA.time=outGA.time(ismember(outGA.time,outGA.segmentTimes(:,1)));
end
-outGA.time=outGA.time(ismember(outGA.time,outGA.segmentTimes(:,1)));
numSegs=length(outGA.time);
if(showProgress)
diff --git a/benchmarks/fresh/FreshData.csv b/benchmarks/fresh/FreshData.csv
new file mode 100644
index 0000000000000000000000000000000000000000..4cdc9a364854d8b3e835cff4106124aa61e655ef
GIT binary patch
literal 993660
zcmeFa>u#RMwdcvNY+yd^FWtayfTNAfrV~+?NFCG3G-w?x?UpP>Wa8dq7=m@QV_8;2
zIktk0xt83-U@m3uW@gr}{_mfw>N&kHMas6E2L#FYP){9Jt@B#7>i_$H_tsvoeYf^_
z?akWy+O4%~@%LuDe!2E-eE#LyleL4jui|%O?OgfhT)aAnZ@!3MKIiYFc=tR$|7-dG
zT|7e{w0u!|(fgr3w@b_ImHVXSvOaf8pF1n}Ny}w@o-ci#U%5|OF6;9`>GQ(MebRDS
zpBGD?7gz3+mdpCQRQkNMa-X!!J`=k(*4Vu<+}c>X7ufS^flprr{(TX>ypAWpyW{r)
zrX9Z*aO?QJfKA8m1^hXFFJR2^djUs|-^=9_>E+6a^m1$6{URKB6sa{sc-l<+z$KI)4R>$6{URKB6sa}qPy~w}gWA9Wit7Gp}FRNqk
zR4=Pz?^G|VWA9Wit7Gp}FRNqkR4=Pz?^G|VWA9WiN5Nh+t>a_wR4=Pz?^G|VWA9Wi
zt7Gp}FRNqkR4=Pz?^G|VWA9Wit7Gp}FRNqkR4>Bb-OygKiryUBPqlu8PCMN~Sob7s
zD=e&UqqmK-zi_u@lIN~&4#A=xL
zlbFZB+J9eZhC@$oto`wDjXCGPbXNPax7v1AtE{7j5AgaVAjSL-R+>9IgT;@H|1SZ<
zoAKYXdQu?1e23aT_F|X+-0S5|8OPYZ8|w!5fj?Bx-Q}aS?(Y1uXTjRUlT}*3irLoG
z3Y6P@8pocu?SC0gP50ZygXz{AYxiS)e{R;e9&`t0jPF?g8Ah#fp>6>#ZLDnuj_n6D
zzbY7au=d~LvwsQe^1FgFg7LxH2l45n82x4OThYrwwEw1Xr6<87uElrX2Y;wPrSF?-
ze~V|=*KS1r&@?nGJ~H%K7jODc(c-)4=i9*3r)Bnx@F3=vpRxit*UK1<*Kp0pW#r+R
zFwgP!*J8fF=}q+WUHliHTn*bcVvdaNr+HPcN!+adM0>IqV?GOL!UNU|CVUh9$lu3$
zSh>7~8*3j&n{Q*U#4qY652DrMn6DuYMZ5|)z)QXjh_Q!w8{z8U)<*?njSzW47A8^^qBLDJ%QE
z89m&HKET&!W$mn(9jW^oZ#CZY-39#F(SGmn`JBbFw%t+CT!3d^zSD9q=KnOl{XW{T
zn+gBn_QQKM-1ctl79%9~!97Gt%f@M++6yd0jxfKg;n&8y_#)mR)Aq{$ntOLKYWWjT
z#oh3qcb}YwonQeyJc@bk$2<-SoK2-*=d0jqsr{_SC~y|EA$XT|JE?IwRFX>yvVV(y
zHev+!P(8AXEvJW);L9x>#(DFyS?tFSFv>yvN9PmYf4TO%80UW${(L7^34g{D{JhK%
zN`}9KWkepJWK(!o`T`~KBOs(FzlmP(C_@V$MW1k-ujAX@;A5=&i=fEVC_gQ)A4jY2
zqYd-m`Hj*J{sbN^?*;jY9`Q7O8H2HJF8V}zy(r(nP2}qq1(`c<#MtP(U&grN_N?Y<
z@V{GuZF}(=ZfU;E9(HYxgIFPQ?s*xxtDOrn^&CVW@Q;H9{DE+HrfQWCZlO(UN#$&j
zvN~T4Kb7C-t5qP4uJt+~3?!|&!JTpr?Rc!Q%Ym~klp&tS`uiMitYypFrdyfnPez_g
zHo*tLCb$}$LvjiW~*ZWVq?G@YYtk~}SitR3}*zV$r
z?JljVeV<@p};21OjPlRW~QN{P6Ky*5=$1&AE1~otv-^7z}4k(yD
z;Cl3bJw}7Wd>emngs$T#?5$`GXF?)C<7Irv8I2Q=csv{$OnVXk@uPr82f_U%kH08B5u`r4FnExjlkK%L@#(F=srD%_%UVSD86xwR
zG-wvcQY71xHGB_Wl~<3~PRD1Gy&7#JTFP4^N%(clQPPpwBEgZLV61Es%?FGlivY>R
zdY*(tIwAg~gwnL9jFnvf_~JW-=PjPH=TjD_rlnzUkl?=i7xo`i)0mq4~K4ANE0k18twh+`)wY8a`V)AN+kI{+=&?iNv2@
zyIlTWiN6;jrpaHvy%0bAB?I6>X>(!i%9^0fKBEl|kFgQ$pdb1=SH?S6dOTNpJXd<$
zF5hm)_eA}f({^dI9c@03{f8Pqo3Q>gkF9tPo&jnt-$ndTEs3GUK2d9T~X^b$yMX~@w$?S*RbfWRP+SoB{6Mp
z0Mne|^8!~e1neef%lit?3&R9siaP1&LXK4nHV4Cn*ts`S&$;$Y7@KQq@QQO!
zOpTK6n!Pk6JPhK|7s8)-u5jvnUbeP$6_3>RB?_*oikVoxYsW3?p)_;}?~;3j&R)h3
z_=YvG8EwwQt22u|n3C|GYf9+0w+~7m9!dM@*cJXcpsYRj>}bOdv5(BZr5$ObmMY0P
zShL`|8@q-d8lE+r+jo9oPLerbAp0z@JDcus296kXC@iLOKp(
z>X{f1-P1qv{<5Fq5cG7We7hSh@igz2zvyHeWygjsGWWxt&MoAjpCKlGKty`RkJf`R
z!()cW%ukUzGj~Aqn5p4?poabM%KkC@)kG@&3J?j$&jwcm-xrgJibH$h}bdGdWi|-JAKxE)3#s!-BJ6f+xK64
zbpP+Sc)L$#KlSXVp8eFrSxVC$rP)tC`>FA`UF9^V@Gl>#1XDvzMuM&uuQOQ9Ji`t;aDyoEc_i9HIh0vf*mAA?K5D?
zDoE|pxe@NC{-0im=~rHH_y3TG0*b?Xfj=R8t$h{s>*X!#)`or8f5BLc>z%z&-t8BB
zC*nR-U&S%MI=-chGvo%bH{lVvO6^IRI=)Lp%%hU)!!Sg&VFo?vErK&xZl2E`&;Sk8EogP*fUxc%3!J?~%OFz{-@l|InV;kY!^e_3_*Pko%cLj0V5=QX?g>kjDLwO(
z_luuwg^n%_Sh-cQOQnO){KWj^e4!w+ygn#=0}Em^cnI%L|0uks=C$av;xn0@HpXkZ
zJOz$~#Q!GxK-$rlIkq%i@|ib_%$24~3wR9wcjFnj3$|;@7x*N+>RNmw-qX$&%rYP3
z^exu>qu^bfE)caO|MHyb^8?PS%UnX_9$J8&tIu?=dP?DFyfsfZ0?%b!6qzq>_;GQ^$O(!o@tbqRn{g#Z*%|G3d9Td#O4a9?nA4T`(M+{+
z@?I~Nncj~b_j!_cfe*rYsA~Kq3|NlVy&HQai~4TtA~f_W-k%CteYKpc-i%k$fze^n
z{AwMRGf>uH;U0gAm0?ec=ki^h*({w#`o?~=gzBI)D45*Q7X>a`MW;p%fagOxIyh;q{&bATX?f|&
zUF?DHtgroBP!aQ>zqvmbClMw||u>Uv`C#3y5(7{;*o(cFUfV5X`c
z2og^Vbgso{wH73rvm}xy53E7D$?vls?T7Vj+h3ATNIu2Zy)Qt9ta^N=-^L0#z5hIX
zz?zS^JzmB!|E1dD*@E6gTi6Te-SClb7Cj}^)%IMw&Q+{Oc07N&RI6{&YnZRv{Hpfd
zSgY2z$QrSJbUInNR@79%vzR@s`l{|bpV8^yA7U@%6{yE!9DY~->1Qq37pZUj&j0K7
z{>=aD?}ndtOaC8k)qdK(8-8+?HtmyoqW{;eW-pUv((GljE>Lr&*~?@vp=vi+a?W0+
zI=V7@neuG+`^HaZFH_C$Tw6|lU1^(G8?%>b_A+60tQP5>y-cT8@Lj#XUZ&d%`=j+S
z9oJigZGLnw754PI@-dw)UY)L&K%Nc!e^|r*&rKKUS;|*w?Pd*-z89tS9MyGCH5%6U
z>|;ug=1KUNQKj8_yZOT0F5|$1l+vI0&&lS85^Rs91$sX9O~{82%Qc
zEAxV<@nmMLAEJ*`@qlvo;x%(ecXlMwz2ZsbJI5OF_0vY5Gf%|~AC*DZ
zju$Il2L93qF*dv`JoBi$A`bmSU;_DEswubn{K(k>yPy=ZgPz0;zb$#)ITvLZF)Kyn
zGIHOYbwY@E%c}E+ul#%b&UIiEF#^lr5%30O<+3y22iQYQX>KvYW5$n}8fNY=^TJs+
z{4!-_T88C`*-r#3?^hrXFD-r3^C2K;F(|JmXtoA2FPynDeT
z;gdJR?|UcCx!x}vlYAp`3(42K8djic@#8Jf>1CY@<7^FnFMJ48x+f$4UcBelO7e2|
z!j`@tpY4aA!dbMs0qAs$h3CpSyZDa4vy43q#Y^aIC?|I#GDyi))trY~$fjS9F*yaV
z8L9>`+4T6s@OqpJA02twoJ}CBl+&|#<>Yaw?@K@`U5Ko3;PxiGDtH{A5%P}n95hd3
z9k8<3(H3a=+uuaqqQ0dC-*%%z#AdEXZZqCIL(efCMSNbD8sO-j83E4-4ByJ9`1@ED
z@L}x_I&&my&NJkV)8#>we_QP
z%5|1JAB#ZMG12L}t=MmFTad3q^#H~_hfo+12&h}x?y46qFJW878O(;e9RyW;8>{8D
zYryLpQwtnORG9ZT8SC3y2HmG925pgkp--v`)aT)8iO${~v+onFd%(ZByWmrQjb7z7
zg<9Yfc*k5ZLG+G7z+RUx$yYf#V^eZ
z#^dRK+SrRvh(Et5y$J?zQ%+T_lDSWOw*Em^*^3`O+by!Ku(p=}CCa`7~me>@y>^C$|)H)n4n#8oSNG&{m*|JZsK#XaeRL^$ej<
zQ9P}NpVjxj9{3i$%p6lOW32=2eEccS8d`Q3Q_;%43eLsqSU>WHo@=?N9oA{S++j|r
zrXEHcG<#0D`h+DXCRrU4k1#YlM@jXYiC}$tyT3t@-O}9Bn-SV*=KMR@zFDn@As)mWu8eY^E3xvN!=<9H53Pc4hUb7~RQJ$_g1O+O
zss~ff#o%z}fc9tH49iiJ0>;+-uCBz?w@*_y5_XDvCEkE-nHd5lfiVd@h^5r+$E)emmwB1L
zH{a4|^*Lom5Ik+YkV_$*6RssWloYFqCE};z8q^J8HCwTA{5G6#)cH-`Q~4+PmgycI
zVh-sM&_|V9k@+=(k0c8GNO?xmKR4I@ujn6q`X+wxU~H`Y-|^~BdHq3=$F!iX46>V+
zmx3y&!1E&7eX#cTLrG4>FDQmK-DfzZJM;`E>$l^X{}%n^GxsOQp|;n>1vdsV?mT+=S!$uE62v?dg#2m@*f?j-{RM
zsbpUhhh(Fw(v5e295bO>G$hfBa?aEPk`H>J53t}){H3}R_;4ov@5hQ+2fONuQ}iX9
zl`KOv8P?Xn=}t0!fHNKW-30tZ!P!QvxbFplGwRAYfkxhrgMA>?5uZBv<`PcN5Xf(x&))~U4b&?
z?GV2pXNsOP6HDI;LO`tDA21dVAXW#GY&-t5pN|fAPFjTU9F8S_2Rzd|C~IH(2zD9s
z7B5TP#Iv+dh840_p4BU$C0iN|wf5i$^VlNBgfJQGp%TE^un=7{N#4+Q_cjBn
zXX07ieWrRHxuf8P>XykD#u`j}>^3v2rJHG$(QUg5Qk#MQXfb$l#mV-90ro5&%$O$-cv7C|N`F5m-wF7|Ue-IrIN~R@cSVrryXn0W_DTO1
zA52bz=bSwPgV;r?d$1$)nzLpfs0VR5*2y>R3jFQ%;kV3E>+p9w@tyYYetcv8V0t&T
zpv_ZF0~snx(T#TEo^Zq0htPo+uo|e%_%LQ2^fA#R3MXp%JQ|*LcL{##+;DqHlN8?m#%K^zEJai5`kfI?`;C~Xfp_r9Dxe(!Lz9~PtX
z8o2HUW+%_2(ophF^kT>yK
z&}2n@IV3p(+|h>R{b|=ky&q*W!(AjPz|~>Py-S<2jAt8Vwo%4yq2t*o*1E+vi9P97
z^Dcj!Efn)?Xc&%>dV}TXY@vwX&KAnMwNR!t+_$kt&V-kI)%5R!
zEB-&@8$7`Ban*M4Z{iy&nEpr9sX{vKL?DfF%A3%?(z+oTXgxK<6uMlua#17
z&kcVa{y;L6pb&fPhWTU7e{L2pspfT3T*fl`W_TmTD>8n`{N(u%&9~fZh;MB~mD}#h
zejMJi7u<{e1xBl!&|8n(U)5y38an|-kT(Husx#%{+=sYDCEiXBNy^;(E?waBJX4UY%lGH78{F<
z0D9xbEGpu5wfxmTTMvv%cvBl&UNU4jm0N*?BMrnWEG3aQ8*Bf2{1P`>HL|1({Z|XW
zb{j29`A4<6Hz6B6E*?@e9Hf|hiOvsrS9#TrC9cOzGeSsiKX+dN1@_@a{OlLpxEY^;
zTSN}=SE60Cxauss)V0vq6npz3o|8UB_MvsL;qRR_yPJT9MXRA)Ym?V
z88a%LWNDIhORDRy2k&5HkLzmKH7!lwfyWr&HZ(=>%;wN4|Exk&Ko{N?8q>5pi4~KwOV8;2d2Xk;MQe!08<$AjY0m
z?uqrvn$e9#j7aQ&_O$G0^kGX80eTj{K+n$`zSqjCH6}At9N4&SUn6q8zc*O&!+D*G
zH;STwwJ8ZYm6Zbg>D93~t%_Nxl8Phh{TudpUUsSWi*R4TV5M$-H0!8gxDcBqj>jOi*^)z
z@ocqTb^shSo
zst5k8L@`HY(A_TZfinONs2)6Hhw4*ijR3FmQ#ch)&)qT)v#^FGLSn=?`6o8;DM{0EK$
z<2!JK9IC-r(LLpchYbr^9i1|B*FO;#cjNB1gyOOw^yx
zX`~y$hZQ;T3ASz~78;`>!)SpSZ^RCRW5cDM$LrliRIELpHD{upRF0TqglOZ$(dJB4
zbgCH(@na~yJw&Wu_eX>(3Hu99)
z^gp>2{alEivd!iAWT%Ydr>OjRE`Diq(m$u7=QDA3gg6=b&_u9_zT(lqzr*P)swaK5
z$V$E!R>^LB*G5Q?1$lbGCm6{Oq*5Eb-;3|R3ooqVx8wEod}|7HBL5O61$0V7XV~bK
zyc~M{DCSH|0?Vsky%DR``mN4=^X@`>stt76&4F~2D?JV$lVqc&
z@SB($Cy4Zv;#;zsSoOul+^)w=$X_SNUu$5mE(R`MF6*KX+EU+CpYZghiIugp&QlF|
z-|L!yK6T>N4zGsMsCUW^{jsdS&8YzghtH+d-e7hpBT}-^8LzH4f$xrXhCkQk>-
zF*DcGC2$`fQBwa52Kk(RB7r<o
z`HuN|@A`8drgyjyzD!dvvIH&!$AEh}s_nV+pO~snFRCsj&w5_;v>7LSa;AKeRkd{k
z5jP~|Inmwd4IC8YS%vX~Z@~1-i|U@J6*hzU#FS|XZGz$5^FgLzR#4NNc}YK_m#eC2
zsoovW8Emmc*^U0;+0Ib%*`9mx2{K7@V^=MmK8YD4A@9U@@P}6c0Z!|2)5ya^xX{0M
zcQvW7xftPPJSkamKYFu{c@$}3IWyMdF#m-%kTqk{MElLUyx&9X>dFM~ljihck%9MP
z6n~!F)H`bm&HYpikX%aQ4J|$luO;Vd{$hy&HAv4f=CS+GM0L3@B&Ia7%I;p8r3
z>9n4!sVDzrKfzCD^q~Q>Qu2krh^IO0_@?ZZtM;^WH*fB1ebSEbO8&0$SWm-cryW_W&u8B)WZ-q?xhHFFQdDo{KhFa&}Z
zm-a{{SKR|RwC40+bBdig%9HIIdWLp-JKC+X>a)>`6@x#UK^wA~-~w**fPM?2jEISmh50by0bo$eE?S{25ZL&rtWeLv{_UfDG`erpG`{D3n1_yujp##_!$
znD?b>bz>X~ng(K1qb=t)%y|`BxESrAtf|qKb0WClDs%!5gzBb7Yq?b7$?xJh`5m~?
zt;Q1|Io@u1w3dTqFTR0_I@;r7Sqgx=mk+UW3`bfHmSer`%7EfHMow`z;E3E3{*UoC
z`E|J0Wg3R7a6Rrq>dNPYKtdZR{_=QUuJm^+;&2z!$8c(l&k3Tye4OVP>UgK7C~}Oq
z`T3n9O^2j5y$)fQpR{xxTiqC+^Ya&bIF&kx+MpxtOs;N>yXkp&c2|G#Zh97Yk;Nm%
z)-lY}R>6gSo3Sy5Ia&Yt+>h=$;Lg+Q80KKBjL~m1wvI8rC1`KQ9LFV%+Q4h1xsQ!8
z&au?P`M?=rO&?m_mHhI!RhNd)g1a=(Z3Db?+xHmEr#aQRpyY9VMY{r>Nn(wSG0myK
zovsY(uH!NCXWT0-9o03i~PwI|A%x@XL0CQ!LjlpV~Cm|!2;T0>Lma@xu5^`+W
zj$k=lDEgBuM)_Is7ZVXDx`Fo*uaF|1L{aNIs)+`F5*{a=@3SvfkwU6Sa;uSQ**%U{
z_z+b;2|f6;g=T~8LOkQ#+UEsF)gx6F(+ed}g7$tdYMJ7ruKVKdEc^kE%aVuKkGv>q
zI`z%f;$3S!k@Wd`;a8SdT--%e&o@M(z7DUL=f}NL*Vq2ktl;+xF0mW!?+t$QB4GKa
z=mji>Qov#S8b3tbf<1T@-+4@;!7pO2FXA&rwzcwiRS6Tw67BLWdEFYA)c8~7!B=J8
z
zI-A6pR38g$;m&NlQg}jN#k|I=Zf4y~_7y*ETQM`^V0kjHjz`OP>q-b@4k${}*0#i7
zhKG#*S>5twe1_-zRg6iD?Kssd!9-aj%F_5S_y9CF96OU=xf`Ewh2#3{5ckQ$HHd9#
zhxz|`e6~3n9bv7i$4A6%Bi1F(p_3wYTr4vJ~oqxF`*Ljx8nE*!4nRr{u
z&&(j8-Op{RvyZG*%=uX8%IJ4n&76+~=Y~hf`wx!D-sHXV
zu`Ctw{-cAwPgz-5zj?2+JXG>eD*wSTyPcEcYDMKLB<9?WwrJ0^jCwg6(6}>GOQ?ySvZyv6I!%
z-&61o>{zzwUuYJ8sbnLK9->uc&*M)As|TDEqwm*|Ncw!Mf2d>sF$>7Gak
znHPl8dHAxaj%UVq3XYdpvW)7<
z|6HwKqKE7jyZXzZd-4LkqA7@;4bBKj1=HgK(pKOGL?Qs9n-;}9rh%KXkTD$G|j?Yg#V(QeL)H)SXe)TVZpD>E{3tvNd{0bU-&3u^=RnB}`Y@i6Rzd_*u}Ib7
z$P@Le0coSNxZ-FioUzgnIT^T6W&G4p(oOBdr;IW73F(I>c@9J(f35;}(B(YVc<*iZPK
zvJ%OJRp_Vtwq3s{xQ}2o7w<{z?miyW$%q4>ETG9e)hF
z(5c=3+hG^^AWzG@j9FwO_%5QG#pl`IjkeBPi=pQ8WMvB9ay8J+Fjk#E_Db4Wcr_=F
zRo`&ct3>b+ki5mqxb+nB{x%n4~i6*X2!
zoiFLO$g9^y8i}@9+xPM9)1uLmS)Up~>akg61nGnZg7yzgtc{T0u3N_2Q~#Msk8XYt
za{-ppQu>_C)sN6i>-u4;Cw}oh1YwGNauBPo5sy0dy1Q_D!sj`7jdmM#JMo5cgYJPsp{@drS|xi$gaR
zp>?yTV4v4aU($?JB~1E<<%DxJ_7|A@hrr>j3Io}1`g|CFe<LWv+8Cc>ZE>XEm
zP!z5x6Cb??NzR_PyP#G1L;2*c;Abb_u~I}yK_Y!$h!uT=7x8cUR9&5-Yxd3Ue&X*FIdaZqdio>EAs&2
zAe`t+iufsUlIS-!=Gh7-hjx5;FQwnYy}&zmP^Y1wtfOP-SQQ8GX1na@PWgMe{Jm2C
zl8KWXQ5h^rImd@Ab_6{Ai`Yf!&+H$b=`IbU$qqpYOaUI?4V?68{HQhzyb<}fD(;IH
zA}hN)1YR1;;i*G(=e_E}q`N-)c~Yo1;n-5PwQKPQ)(20wEPp$y=6TiDWPVx=oR@5D
zs7HPw{%U^XyjK~)wTylg&qLYb#Hz(PE`M{?_1Qr2InA<_li)?YdL#|fg)ybQ=ycXW
zp$M!UcoLi%gYeE&0JE=nJ_FJhloKg+qdaAI0ZBx2V%xwBd<7
z>)7R4a4Ti^!TD5uRNBhfwSU?`ZYcuvw|GVTLb-n4C#@5_0h`J8W?s{HXq$~Xl3&M@
z?6u%4{S~NcMQ9wht~o7+#1^(8y;}_}_hKBc;Zn@hGcd2WOcq=>*8cMbd+SP;zz0ms
zY++fO=$|?{@Lx8xvuC3BGle_3r&`mKUPKi~>|SvI{7N|=;4c4L4Nc@?D|I-Jhvp~Q
z0##BAko~c@E9F^R;jr6Nz4v{M(o&Dz;H)gzmzk5oMCnh$0Q849p{cVvY4eW*L-`KP
z%hFMM?d;O+a$17VEa$DUNt1e6a9UI-&kj5DFlGa_s6xDO!u$c~vaX?fw4uM;CH!yT
zIMC)4L3&X8eW8_TO-O_m4sEX&qEDo}_hi37;6=d;_bV^_cC;FXZL0^i{$**eqIgkE
zt)-<*1lN~)r>q@v;!gC2g@v4w7PDLSQrHik(f)euGw^3$p{K+N=G{Vr^QO|eH`e}8
zWHvWHiL-62{Xg;gR(yUxY7I*tgnFgj{yjd&(}%yi-=6(4rf4R|(e0TvRY%GC`Bn7X
z`eGgz$b7M&F?b-{yw>eJ6X*-A0lvmQ)iR>C$;9!}<}sIR$CZSTUX6tKr;Y9>NJF6Y
zI(8;^-dMI*_{PiFO?cHmZG2jI235k5K=Mvf2Ybl=hBvB;DeDlQ*E7OR>%GHXJM@%4
zm({;My*GUd7s|eoOo`eTu{#tl0`4ZBmvGBm58aH!|$nBb2dcJ7THO1BP#w)Gu3~VLj
z0USn}2CD`WlxwAUl53{^JbfAL-W2b$WbW2dOx&rdY4&hHlgH+vlCF#v$Hei{xaiM@z8%!h4Z2zg+)@Ea2n#H;ztq)&kp|?Aj65?!-lZ
zKE4lYxczmH;WE)8u=!kt&115n*2~`oZclkeoEqb0m+PYOA2Eta#Vn=w1AX?P->i@Bmn@QaU9
zK2U2F@++&d3Z76bA6g*!!9kv9{n8Q43A_6=Zo}r@&c2!1H#7TYeyP5hx6_a9K*+&t4hH%w-YR#Mks^4%;PIF26LtOy8F~ymn~*I5`w(^43Gr
zTcI5~ibqD8l03m=l>XC(&W5B18f{8Clzcs_>VWkOrN?KCKstaTQvZ2l7(xeYckNW%R^%>65na-QEe<~p~3{7WPo_sbBn4F$2Yt$!=OK*!o_NWcR7GJ)-
zz%A8axy-O7sy2*IQ0FvkcV^fI1l|qao}XcR0(w^kpB2F
z6xY^gmY4V9(|%-F5xF|of~YLM4|4H&0v%Vj8FgxN-SzB&!%sJVqzzo#kfGiq=Q7qK
zBFQNU&c4)Vu}&AyGJ2+N-sJhKWjBHKH|5+;dzR~d_~vwi%O2a!ctW{SU&a4N(ZjuX
zR{9t9FwltVJkopdo-7G&-*omHd2ou(k~4QTo^NyJZUx*GOXn$O^defwzqS|sffI_G
zlR3%!=!2X^)i&2S-HPRT=1{fXwW4k-pG#C-{69o&aW^VP-wKcXzeincB)Tw_tQ6Yf
zN2D_NdRSy_?DT^I|6%L{Z{FX$u++8{?H$eyYdz@LH|#6e+Gb^fh4xSOeU~}bG7n=L
zc}9yD-79MfojCOy=g&-M^ui0e*;eT7uExzSZkFm0tv5M)b%qOmRH%h0wqm_wn#2bM
zn+y-}BqHR>R1hbA72IgZ&$~LTq?yu7KJ&KT*mj?x1@lBa$LyNo;85|~!H5mV2ksm+QWnhcZ2=E-gz^{&VF=C_^~=m+j~5jCwcv(P89F=3D4V+ql!2?#NT~
zOg*QJl{)VwX#`$i@e`=|4a%WflB_C05URBeCN|D0mp
z3;2B(*bK%X!GT?!$C9U^PamF`Rz!zARxV23(wl8oJ5`NDB^#kVPkkz=2s*C4E5dHo
zn&i8@4?(^R=|$~r2XH~%b-*39WSx-5js)n&i?TQH3?7p(P!^praVWPYbt27zw$3r}
z%?dTO!PU`x$&h=s(4uNiPk%a*4n43ttavjpoBqjxglZngS~zveSk7763dqs#*JUrS
z|Kjq_ta0ib_V_7mP74$3xwNbUvVHq`dy@c=G$%k#qs?EBUnE-gF#W3TaqP8u6qZmu
z9_ObSYxN(z*FFU~{F)pu_6>=N@3Otmh-a3!s7IJ?J$?&8{mpbA?{-f0Jl4a$04iGF
zk2|a3-C7;eCk>_f{GxJ6r?tAYaOU%iSX0S{C8Oudy5Y-SYxdc^ulD-!Yt=&oBB6n7p@8FuEP5l#qmdn_xb_IGXpZgT=8jeraKJ$Krs-HEMSlU*eSS@ZfdpP=eD`VE=e10);?Opme=JShb
zp|9rqz^+GQ>T{bhT%OM_F871987GQ`&gU12B+lm-hor__TWYQ?g`cdOXCcm^{Iu3P
zWDU?!bytM?N2x9Kb{aBIp?}Qh7vaP|$BI+4S4O$=_&t|Jw^l#*gDyxOz%Qw;)O>!?
zH7shoUwS-|7kC`hgm$l-zo`>a&I-=cdE@|*mDjBwWk1(^e$g_xjVRCO7l-cDn9nbg
zU;nn$IiFv2#XNc6``(Io1>5tnYxB9=`yjI_oy?=5!==FQ{~2KMD~~;#0bN7XJyp6P;(odqz$L
zm8Hn`!4AQz=97|C<(7AuUh0TneMf_O$IG8q1-_a*-N>`}{?r5BeW28y?bP==rIjs*
zPFl)ttWU}@^ZBeKKGtqV`;T!}62842H*Ej3XC?FeclZ3@^ohm%cI?EoZglOQRQ&PI
zMRE@FQvCX~-nnwqJgbC$nca4ro7BlkeR8h6zZ~DW4SlTkx(mJIzCRU}QtFD&u9`sh
zFQxI?fK+hz}D^WMt|I{#{0W*UY2j~F8bH|yKDc);@j=`
zwtYui_099U0T1V3KuP9~yn<_Ioe9NIxCK0+XDyCsua%?sc&5N$^K1{I<-L
zn?vq~#Gr5AXTtcLX<)zJ9JbmI>3l6<#_W0OPTZyRQS6`Je-zNTQTA0QwSf<>1b?3R
zI9lI{|25ar+Bah~?)6#rfS8sA(mrw0FYu>FZdXjFa(UiF_==NrHVEI{jFo$R9Vn{(`}>R!Nv4Ej9(ugVg@
z;bvgQ%knK)@UI0+T(L?ptnYg4Qe%6QPkX=05al@rWq2FrCqs4P*eq{yfB!ZER$r7;
z80;^46fikoG(V_NH)cWkz>n|18EOLSwsF9H(v;l9-KE%GhsO{71C^Fn(Xhs}JKn@{P7o-min!
zxP|p{&`b)uzBnRhv;8bh&ThN4Sc0sHYx8EAUm4sqT*~=@$b~q;?Sd&K%U9(>OZ;
z4{%DiG-J0~3u)yB9%}X>sU*){N*6!^Bc0$dkIVSaOFcbKyFi-{4!>iD?2YH4)tN7D
z2A)5NRWhr;#O%;v;A^Mjuk7RcdwufdTaLtYLnFb9&MZa)nsCR@illPue~l_PAC^|S
z{jKilgMugUP2j3~*0__rf5)I{bUkpYu2jXE&&3bCHD7-ocqAPz|Kqi}9d%u($wrMr
zIKrvJ%4GdhZDYELbL&|xypXxO;uG9na!8-nlqYG>jwAm7_@HyUIy(y@VQzE!JR<8^!e;63ajEAd?p)aA)!_C>_
zJAvsNYyTF1dFNVUyuQEq?BnQ*no!i(;#p(?pU~!$MVtREUf+qG_#g58oAS-o(yw~A
zf6X^)??*pix8%H`|FG;kb3<=uCuCi~=g|7$e0;~JKJ!T{@kww`^2=I;D5<;GvRN?OtbI4us;=`zcYtdcn5veW4ELQ(P)3wKJ`nCOXYUOv`|I^=g
zTQ7TmyxykQGyU!C&jYva29?kLyxE^O`}3@CSgW1=dGH6kRkJ@&R@&^(dq4eow?qFy
zZ#q@Qrv4VuE^cU8c8f!O2Lzh0HCX@mCnxI*TwhpF_rhk}jXBDvh4q2{FJD+zVt){y
z+PiJrlRIqt?XYC}{=TF54{JZcaeRZa#^jNwGClt1PhxG@Ijq*6`5(iJhz|@{b*rd%
zp+)v6dN+0%S!TZV82ii`_K>&9o+i9Z%j%vl``-8_{oh#F{onN~P(cJ6dfYd^%)cN#
zULLveh>~nzaMU&PYulPR$j+`U{q@*+`C4sV4`smAF_I4owB?fr4$Cdc&q9vv#15Ye
zKA>H|%4e4B9Uco=`OGV$(e#bwsT$gXxc;sJ&;3{{5eFdfRq=nYhhFG1c3
ztq{M*$whDQ04<-Zd`j}l$-{u>YpeJ%6n$ieu*u>Tm61Dn;T9(
zWpme-^if}%y1G_3zF)FXWKq_BU_3}*wPFUsBk>I|3;z%hRg4_$QVXgM%hyG`El&ah
zW6_Yc&|4t^J_-*U9&2K0_|kL(HTwl0g7a&fwyL9hFtu)l9|6BAoEOg;_oT?zL;MOn
z{5t+pqpj@+PdT0|#=v`<*$a$|M@wfWzAo$B3>f{b$kkhm|Mz3m-SP?dSLoH%_)a+s
zyCXgKbN||?N$WRmA>7XfXjf#&?1q%pngl<5gv0Z>9V>+@YyP;#?%8hFjn
zmrdcV$CnEJ_I>K66;Zx+oc^Ii`Pu1>NWOi4Empu?b+Z?;=pgp)%NSqr{`$TxMgO{9
zbB%!i>Y(hMq6lLnOY4P8ft$vV&Bgri*HQNn9(-l(SBqVRPFSg8e0apC=jk(CL-*(a
z(f#6?gl6QaeNylVU1THHEj}w)LP1p7(VnnE$(DX>!am=jJpw~m*>^z?MCZU0^G8#I
zYWZLpFhl8vc-s8KoMkt7buXeP`e9Gm1I5~(1tg(luuoAJ+QBEicE5n8!!H2ZzQ
zP%8hR^XmkVK08-FG4>&1M!7|W@9^w{UFfK)BSe2hrGO!KnShZQ5B`18JJ2;0
z&3YWO_t|XSkj0biKhY!iK6I;*xYmp4pZpADuwur*M16oYNPA?qs{D#SOXfI&XI#r?
zB*T)gBg=9&U^j6o$3VzFqF+^vXI9{^a#5Izv7FDra#JEu0l#T`G4|nC@qZalS4I+C
zp_ZWV_eTL`-6g|!8CjgF${=b)Lj6AWj$J^`!9u%L`El?@
z=^KeR!a#BVJ{IU(^Z7bXn$&=H0BsJ@AQ032j;^%*QM6_jTOlQht|H;UA}B<-zgcA5
z(}0j-SQi(xL3=QTozVPs`OM2(ltmTQN7D
zFc(Jh!vB?lht>tnvm3GZM3`X2l#1JXcsDhu8BPY3+eNwK<`YV4`CJ3|o71IOu79czx
z&Qw#?>_M142*|71gYZl8AXpzjnyuogfh!%&<8U>efE!vzLBl#7)+W43XC5xbtD09!
zkCKlC34u;H&hy0S(0uwMsgm=RG8WmijEQeu|22nmF^lEn_u(YIIXo`s0~upG>^J?_
zOt!=K#QUAaTy~1a>Aq_})Ku@}!D$8Ro%4wrMP37LT-AH6bbA5q?XuDul7elk)%K-;
zEEi`yt|G081-5#wB%!oR^g78I@B$fxCM>zF*Va~$Lz4L)7Hze!J*1X0>}?<`)I959
zjN*F)&<+_L+ydV?L!P;{NO47gm~nkWL*BZJBu=dD$K2TgddBW0BX-!jsx@vb@;621
z?iKmV+GS%nBVYQGe6oyh%(1S3-*zY;1btANu67)2fHeXcBt1HYzS#^9;i;0z$9LrH
z-&p9|UxpW1Iiugk9Dj(p9>iZ-+e3ku;+PoB{w}H9J{9P>9B*pP^Ory*#Gu{nO$p
zf~LR&yheQ-ms~(z==@phoiOC4Jej}|=y;t(fPes}ro`jLeOL{%XMV^6d&6r_r2Bx=
za;F(EFg1A8!@N!
zsMNDhe+$r;_o?oI_4DiyU$K{@@7icC9W(jd7q3ev~Cl&eppcVc
zIGCk4ZH>2<;z%%`%6?nVyWk?Y3E!ZM_9_hF5&Oz6d7t10qTnga+UmrO*r}UkE_DP`_`ojYP7_}urYiZJyIseu
zyhC0$Q4QiMi~`@h8#7icbXY%RO{UW2vO1;k6kY+}mY)gY+I8VC`{N!9~c8xH02qQ(kj40ZNqG*jwmzKR7bRxY*lsV3;lVVyfFS&~@3Xfb>I%RL;GNzT*
z-eYm~@z{-QLC0qV^Xa%AoR#VVuR^{)D0=`zQnqdf{#;wcW9}4Q=a>N+nq(asZdZ~H
z!72F`>&TP)qNX98rH!Sv@tL*{&kEkpo=a1#*H+K|YBT<)v3D7i1IJM7S|j*+`4
ztOk1>yF+i2t^+U598ImXr3El2qa?~=%eWq~5+ZxKPY)s!Nci9}ZY)N?Cg^^Lg^}Yw
zh~2Dfp
zo-F>vywDEEBa2;}Tn6RDS;1Sj3iOV~thGVScmmMTQ@1=HuluVpU*doAb7WX8QHZZH
z2Ky>a@^;x@pnz5cWwPt&n461OtZ=a#fgR!z&Sx`D)VdGQ&5SmrjA42?UL!T614!0m
zjc9&?j!tmq-Vyykhj|_6B=@7`LD7>P0r)V!RlE&-TlyTFm=)4*t>_(8lV+b5uUcqLTw@IEG$uD_$2ecLEwFlS8;L6wL-UXsXW26=(1a-exBg;@au4)8!t6*`XEJRh!V@
zv|nx9P&iWOA>0UV0LMXE(g#?6GqB}G*?lxe?UuQStW_fYRHlFf_MrkU@5V1wrQ9)R
z7pv#wU8z^f;si6nP_NCo2)3N?)QJOwYcOp*Z!_kMR-oEx
zf+F)w-}!QyLUt!*9`(tij|m08-I#td$3*72qVG{xdh}=CW?auXA-wvye)erR5xLD^
zapJ6~R=zg*kfg^`Rl-yx{*)Zmm@I3+l6@I<)Y2q9RGAH6+&1qO^+>AA4=9T7@>j6N
z5^G#5f)cYYd$7@2_6nVG_C&t*+BYei6q(b_(XhM*e$v13@-RbWo%!n4Xl93`AJQ`Q
z1+OA?ns%+J_4ix8>QO(MHuKgF%lE`SZ$(S09!UnMPiq|QjwUA>^NuUi#(G>Y){_
zrmKi6mV
zsyx5MEF2BMqYHfFp*epLCpq;Jy_&Py!7x`294$v}_73c*?<*^cY3fo$mC6|b+>DkzvQ
z(+`#eePBC!SJD$QrcDSUw+ha{>Cxt8o8rZ8{dM-OYj?%T6NYt;Luw85l9-oL9~@Ec
z%PLaG@e{b0;fa(31ISr^UgoNbak?onttw!HHiKVMFlKhy5536VM?+$bpA_8tDqbn(
zP5#;~DqOC`0&!KG4TO@=25@bFSo$Y?p{Vc_g`?Ui;CYhs++%
znK{46S~q7lC-
zbAkT=M>O}0tBbbab?jKH-M5-1bT*EghGQN!XIHTj
zzlr&Q2af9Evz1j{a|Bj~Uhb$}?rt6X>|+=7<<+Z)q1+R}oSG@C!8pV>Dva?L{aK=d)?sWI$_<2`7n-glcgZ+e;Mn&8Y`o!HEx0Z
z#go!($i~i`25^OM!LF~P4@UP{U{;`cP~9S<@mW{nQsJ|Ty^pm~z1lyl{W1FG9)SJ0
z^I)v4`ujfS{e857dRQf669ZO7Uh7H9n($n}{u}X3#*uy>d1Azo(3GgNdZAQ1P+du%
zTKOfN-(z2E|HnZ&*{3QURLN`Odq68AV%=}Y>N2m#{{11e(Fd_Vb+k~spp$*Zd$PK2
z7p%J%Z4ct-Br6@`OL!ibfph|!RL5NN=+=2A7l;bH^%J%3m$*f(f%cEAjb@1@m0@(f_+Ob%7i74d)n_Op11K^kFu5u^G4dhHSd?T)*7|2MiSckWy+(j
z?u+CiSD5w6R@r7NYPn~<#2h!PL7?LxK^2WT>Drzt3A9>jspy`aq?z{B
zP{ok8erNkG+qa+BwLf&n&0avRt=vgXk9@vt$KPSC_)PxrYiwZTBQ`J?g6<;?mzws&&*V)!W8d@nbkDBu+3%tKrhS*KSL=Q3
z>G72R~0|O`W4^l+TOC9-X@o6I%^608`=6x)Mnrl
zH-mm6F<*Xsog<9A!wKbo${u8^r{&MS3J@
z(LYZ}J`wV58b_8pLbZ4_Uq1M|?eVH3Br}xf;QpfWq<uWrUG=NxlV?w3{1(B6Bb*<6hPIVwrn%>n4cX{uhpT#Z0+G_(d)-IBZ>iz`7_
zb2S3zauU0M|KO{W$4)V?He&=W$XIDjajr&SPlzFt9;)ss1iT(+g6rWo0@kfY-`1gQ
zbz30sD6u^0^0VhT$nID}5n-Kp*eiNEmWlHeu%_|V%9Ey<`RowZ0hah&jli6`&h)Ik
z6Q+pnS(G-W(WfBA=CCH_Y6P+&=V}D?DXoki;t9i|CX)b;FjphUGwSI_%=78uclLGE
zXDIsIH9a+RH3DvQY%LP$QsbD22S#$Krr;ro<&3Dg8iD<$mK(LCMH-+b+zNUh_pCZI
z3SS%A#9WP_&bSpOn$oeX;hAtye0lZxbNT)?AIi@f^-=%jT61N}L3r!CZ~NlDumHc5U+^i?7vHmH&nhhO;q4(QL4P
zu0~MTGmyr#8}LFKovRTbDX8#)?A{MdQk{_Wx|5HKJRzR~&riQL2v3^2>l1#{Px_4g
z(Z}?YIc9xwWqQ>9N{{~R+l=cu=W16U*Uuh3>jIl#u0|mFjBir${JjN6gX5fHue0&x
z>jcZ^Y6RrzxrUT$waQvT-$9l;r&8X^Znc!T8o_^BoD`_fiOtmrkRxb`I{7+Yw_&bE
zKsML-iB@uxt_B>*3BY0$Mb6a-rmKWFdWFv+wVb&c0kD{>5e!%7m)>wIC=H+BT#aC^
zMu5#TS0hk{>s*asW9@omBI93t9(XZVBS62aPmSUMBv*T`MlgN@%$ahg-{&dw4C*Iw
z7LwDKr{bK|-Kai$J$xXeB5~&
z_DsOgxo%UHmohVh_)n|3If6+pbJy!j#Br{Xjc?~DB)?6O1#p(CSI#O_
zjWL}!#nLh+qh;oq(7GxOCryb5wvh^J&UKZDJfmvd{xNH0qb*D7tw+WVu@+Tf;U2;9
zy86DAi)uZ}UQr)gv1YIFr>+!aXm%?DF<+g8(D`Pcj(~goI8`xEx@yozG4mH|e=hsQ
z>b{EDBDF-&=gCxm6?1mY$0@-C!+@Nc|F7(5+a7w&{A8%8F6e9OqfU
zlWtB7pL+LFJ2;tFrxrMzGd4b9_etZHova#`$WGxF)$FNqOSR{$F`!z^tAmaD;Z2$I|~ZfKqgA2G6AngI7Q9BroLe^KcPy#UV0
zOsV872z8;tsJpQe>Uu203jFAALNCw_=X;^?|30{nU=5Bi=8%Vw$g@{x454m!(d
z*Swj2@B+9t(0De9c4c{Qs`x`5SZZkILQqV~qM`Xn@0T}1J;?G>T?_m>sS7Rd4ey3(
zC-dZZPw>X|h@8k%z8gEqF1`uQ43@W%3Z14xmz36N9F`TjU!aM$h@A#pSy9^qV@Vq)
zW#D3L-nL+xG|29&M2ExHE_9m;9VmHQ;XF)Q7kh(D+l4y3(A#j1>{9plUAo%V`{_
z2Ylpn?ga;9`IaZZBeF$QIrv8T>@SP9WFvhS->IfCE0KnF+$zHt%ji#|@~>_UXI`!x
zy&k)toOfBJ;kIp^(LtZsW?h)2bwyQEpQ@qtXGPmr)oRxa_J3-0*IcK@3sLx1z#7>A
z1;Jn7R^Y>KU;?-WJ}3tap8Y&{Z9Ut~Xm=|t4tNI@@1Q*_CF88JHR}=5Zc3h<8p(d_
z8za71*d-66r?!Gv_l!ZxIJ8MuuQ%7Vya8iUDzJisSmV`LVar$73)NBKJn2k30^WTW
ze3QF!sqfl;e=hV@T0pZreXY!!o8@_0mCfs?hrZ7|?pu$nPh>!TP^
zJngTGG4tk_}NfK~1@yI?q>RFH(B|BPJ%BW;t8z1YL+=*|%cK8Tc#^BGZzzS*7
z$iVaqux>@mXv-QHN3|eiBaeO*$P?!A!JyRo*XF9Vpz15#@laK>xI%!
z`faA~->5Y>G-zrgBkAewKSJVu7{9OM)raw(`*zUwuz-FY8uOXZl`a=6mQi!g55faR
zP1jGNCp_}q)mp+RDll6Jt6K=jZ&|OW05E6xeOWyeLFOtt?vU=A{@>-AEvx>ibz3ak{?F>VZEBDC%6fPAKfV^Iy4GXi@zN~w
z`RP}+?pnXjr=SDu-toIjm<9Np{&B?|?#453!UB3&ZoPdSJPkV%TN*A!RA47~+NFZS
z*c*I$CO$bE|Ec`KT}IfVcy8<~>Ej{v^QPRUqqFxsdeNCSY`cH==Me^!fW8JPSJO+>q4v^oQF!(qvr-cUpBXMVJ}=PbCZV4
z-YW{>zPdBpiWZU@jGA-PlbVU=>-!+$2HFOJ}29yVKxv7v}`2?&^g5?bBBRhtC%{nHp*+5ot;48iRRv9550H
ziRKfE_?*#kHe&_3C-1298CD&id%nP_?)c2~v*1|;4^
zq8@O3UeQ8+Lse@vY@`?L#kbP+eS2xza5vXJE>dSNu=rZ+I=HRZqSw8k4L+yG`!Ofo
z3QF#XpM$Gk3pnXnz2dL<`d+*`xA>kua^zbfFWtvZX`!*znznait!=cKn}PR&(_M>M
zsWq=1zq(X?&TC;WGdNpX87emdO9Z_k$%$7{dKO;6-4P3nZ4r9o*)_?F`Wfz%0g3Yci>8)jQ2~YL!_b^MDi((pcZeC)~aQ
zu47#)&Y<-H6;602RS+%X9nXM?^a2f_XV;ori(`CS}OXF}8>1ymh5$~7fukCH{2#kQokkbs6!*~1?`#RJwe67un-i(&S
z#(+)>HIM&s;K{A{*4)L|1x(3!6|S^v;s$kmL1+eV$F(=kiY45_n^xZzZVn`BK_7e^PeTjpiNEdxOz+mq;d{JSL>IqW(C#g
zZ&H%vXaDAM#dC&d!q}`BY1@x*GH-oLuj1|EP@*`bj%jb$ZX@S}wE#z|D$BzGG@!Zh
zT~Upalid|?YX}N_?){QKLLA-v;wbyXTK}wB#}dyPmJZliy79-{4jsZ4kU2E1Y%fZWB(K3z$hrl5&b`sfSWwCvNxE1o(px)>`+
zS#~B*Dx?8IJ-
zV{P|0hpk!3>-e_ar(02*;C^JNb#qvOr1Frsx!&_
zx$7oj>s=25OWs7UEu6G_pT%>c3wus-CVYd|^L31spH627o-KH?8&AG$uv|TmNB?=*
z7x199Uo}2hCM?G~5)Q1o&eTNR2Uv(t0GyMZu6NpLaFX2?G>2w$CLpdK@@TMd%3Pjp?0g9-!(59Q+vdkU?D{kXYyNVwbT5AL
z6JBK$E?>17Z|%FWR&`_#$PpK6f1sqj#UA0CLi#$Ym%Y~g9Ot!n8-27MmHhasM-}F%g7p7d$LQ;fb5ue4^hQMW=BNT%
zBeusJRj6Z`{
z=9UNXmlIXM!%*wH_PI?*4yNJ`t3(yp6~`;&aXd;);bKUPj3_)gQY>LOlE7T>TI1Dx
zT>M&BV}|$x@hvJwmU(|ju8dmnE~~lWCB!f1hyl@~j9SS5`%}gXPDLKjRK$QN%~Zrd
zKGz%JVRY0$=VS3QDwEtz9gMx^2R(TL=qZ3N;}@dUq33<3-#`k-1d8U6YuxZ79p9hge~
z5%KE!c2VL8#7~5;e3x-Bav7-^Neevra6uws;(K7?yNi1$Kb#8Z#HIMj*vF3%k2rN0
z!>BWOi2!W|92H}2GkqOd0z!^MaqcJgqCTPH)@Zu^J@={+J
z8Zr6CaD*k}9r?;}sP6yOqA7|~DYB7$9w&lQKRXsp(YY-CDxzUIO#K+DQgDXwP01Ai
za%5iY#;4-rq7Bt}B>&fX<-f%qDmRsAI8W##A`xApX3igHe-w#>Qd(Mix`-9hE0p*=
zD2sik2aUsi^AtG7+@Lk$%+Q!Cw?f_UByI)KsVXQp^CXRX^h%Tz4}|CzoTcu^_3-Es
z=|YF$uSdyq5RA`^$Yu~0(hKy@H#f_F-alXaFB^Q4wUwZavunSK&&Wk!UeXr$PmDV6
znFT|cqcm-u_Mlz|x*0HF&c;!*=T8%27fFtCLIzG$#=$qKOGJaT*=Fiy+o`tT-PsuDLBF8`i
zu1$23H4=%}PO5Sg^OLVhQE+zPLi~*HLB?j8AtzuQOGi7Fk3v!G+*jtuIWT&4WEZXK
zg$rk8PWD~%%+tsvV>TKe<
z1+E--RA}k%hZ*Ac$}WOKU>e*5IEqUoEsL(4m8@3dlqIVMED^RsBi?^6tJd+9s8l;>*ur~wUfKsObhZRZ^+~KiwT#Js8|K}N@h~uj*}aLmnWHSD
zW%5~6VXconNzbOp)u_%?0j>&%lW$COqfZM|MKynqT_n=rdy$7|NL-`EGjXlFuTgMT
zOLmzVunS%45+^ZlyB&Er*P|6xoRBgvcLbe+_gu^ZFSz@Vrqbca85f%?bhzFkcNng!!jS5MLq
zT?0*$Gs(Y-|HE39d9ok^o$nRumHwjhGy6f4t#to5c&4g@vnFLM?}SFkILXb#^Y6q;
zQ&X7g3;innMp_Ru@=7>8gcgJrz__(!nc_R1(GW*mu2mpq1f5}>!Ru!fPW=LNTej9<
zX~S$C(?er0Hr9Nx<;><`;y~E6+b1OdYfd>ernL^A{yv_rBMttOz@{H!_3STr3}zv>
ztPM!=m>LoWsbwKcYG`|Z{Y-K|@=@;g7;3zQA2ArkB+my!SLFnqxa&BwJ$*u2x-{r2VkK%CVS0l
zh5cGT?Z>$zPe}7d2SIBduUDNIZOW(Leng}8a4HIAxI0W0#fZUkj;z6X7d
z`>((iW#m}5?MiZ;f3m%kkp=NR;gbIH-J7B>fi1~rjb$~TYwI=hynLtdbe)u-59FFS
zFn_zUg3s7J``2XmYaRHiPQqmdYhrpzfg$PvCRd{$`Ro1?
z`yp?3D}U98d`sxbDT~b2d|D`Be{>9J_gA$&@EZ!UeVmY|tKIQ$1J3NY`Bqob<{k;A
zs<(mOFBkw1Xea$q)V;z-@JY1b+ubrBje2w8>0su4-Y;#%5CSKZZ*X+(SMxaL?RqQT
zTeOWS3x?f9DlYRgwN}*-Uu1M^{=*(8=8C7E84*=MvJ2MftD=n{Xnt`!D3fYJwN7s?
zIP7h*1|q-h`Pe0TZGBFzSyua;;E~*s!Ctb<)=U&LI*w+tAFz@4n%D9zJh~Qx_Hj4-9b1vJaCd<*dSBnI^ta}F-Mh3E?eq2}w~(LnZrJzg
zQSU7U^f}LoZ?rzfA#R-da9@I~qRFgg56e41tr=6yQ1qkgRbAcpfG*-nI4wP%3%!pq
zYkE8z`*6OZNBAsfR};gQDTo&YChNGiPS}QGAbI2JG-&NjI(M`*oQ))gGm#
z*RYr0f_Q>vVIEL>X~8kek>nZdkz#V(}61{RQ|Z&89O
zbg}Qtd64kqPU!R33#A&$;1UsZaKFpHn2qGMtP;fse^|^xJ_Y7yxDQ*n=1WgP{;6Ui
zD-^{i^?_5lGWYn;8@2rZuN&|lxQ!sIh|QVU&A-O8d`FZ~*qL5MMX;zVO=b}^1pH*#
zf`hW4wYFzL1FTE@58V{Z;f@{p2SP0t3y)X@`#8jv1$VD!nEwM0tgqdeAHE1+imnA7
zqfJV-2qPE?43q?hKZ@IH@7mG+dBba`y6U&VH^4@%O0XmkNi{6%SRArQdOT7vEl@0B
z(O${`Nef~5ur6PF%j=vj>>H>>K4^LLoYTcBtz+ZuM$ZwJ%Qq}Oop=O%Sv!$7y`~4&
z(e5{Nrd*BOd1QjHPZSJx@m^d+UItsrj?}XQ!W+#}(p$8;2yViUJacZEyEfc*zr7ifSZV~Aw$8ylqPszk$3b79Vos)
zZe_-AgemZW-SVtPQXkS4fWA1KbULgGs$+R3(q|amyC7Z2b-7GeeeVI#_$qpXxAMAg
zi=}>{d4U7)9q7Z7T$(I9*PU;2d07QKXU#=(MW$)y@+{B;S`Iu*Rl|gBFUtIk9jn8&
zmNeo_;)do*;sKH?irPJo6|(QSSM1Fv@rh&T%rEt9bOYYG5*HB4b+_vowt)A`n(JK7
zYw-y(2hD-LCF!rlEAdq_jPt%;Xo=i&<*|dUH%d!9*H{r?m~sc1xSZt*g4OpvE1
zX(+&twXt_KjYE&LpmrbfPO@M6h~}&(sXeU-2lymj%?gwiT#vyV;MerVcxX4k7%D*b
zPu^j=xEY!U^n)y{;m_KPw~yl+&62(>E3`_qE5)Z@7hGo#z=4|gN%p=94F^4!8hxx0
zO@sGPhN&|lmvpMRwfLA^!x}4`!RyeO18PEguHfiatWo!uh-#jd_34HmOMh%dS-8R_
zIFwhj7td#YdfhX-lzHk&y)AZ2XMk%66SOWvk*A<)TD>ZI&5%?@GXfIU5oNWpckV|P
zw6eEEouUJJmnH0(*O=(1={tKDr{)e9tTApj5UwH5xX()2(8M8}!=Ky^JOA$m3#=d2
z?W7&hKc6N{vVUuXi3hbVS>s_Ba?;AczfWU^nsa7Vc%@H+OK}_2UbF!3vHjU?Rz>gH
z9VispgR*MnZ{a*YwI*Iyb>6*d+ZbW=Y+&
z<{dYTiFZ&R->d8Uclvj-Ns_QPE(3h){ew^q9m+{gWm)W^C?n>XnF3-skx*;A%WpBlI*%NMg
z%{^w%8BaLWwL-1W!5xHWw1lU1$L-G%zBFVzYKFD-Ftrn>V=VG_fR{#vo;-eQ2c<0z
z)^y{`imf8Am^)lfcNTahL*9~M+-M=c8$IzX{K@dr_`n1z%82&R7(-ApHy)4rvp&h$
z7B|TOzZ2g;8{+Pw64NwXoKK9KUcI$zpwnvCV6@@g=9bsFP+mFvRa@QSY~U-j{85a)
z?A~2bWZrOLUDb5!8MkLT3^8FT7=^S2>t$i*sTJW*M%Ue+^8rlxj`E9d%0obJo<~~n
zeJo)%Hv^9BlC(?UlDc?b@6LYVHuAxNXUINFf3y=(F4B95qYL_A26R4@ql@eUGrBno
zn?+v{W-to<3samUEs4>$#M1)__ZP?O!W10B3b(W*R_khj_N2)w)~evT&~{1!_jo5p
z07H=+@Lsqnt9I_$^bIDQhDCQ^b?|BQchTyTqKyDy;W|7K9?M@*ka)QK@9c=Q0d~5z
z^PviQ1@kweSJw4?S@}QWsuh=46w5ji^zu#oNKQf*P$P5#G}V^#pSuA&DEUoj@4C~{
zmN3@0I2xz0uEx`2mih;!A1osxkyH!CY9FeXEn9&$ry@T5B4#3gD%ua4knlwKX**iH
zSEn%h5?!%#$*-&~kyzcz5uy-XIbvx`Bu%Xq=S6?U3!ym2nUecJJXSSG=R0We?*pTg
zxPumaxE@Xzg}Tice@`2ch#~T}rky2Q&uc6BI_hbL1v`-n!*_JcEqeGkpzB|#@+%yV%&@x?N3-$LMf@+W{WD|!_5
zkAW&$o1V^ln5mGV`9jTshQyN+4O4`@xr@io99&&m$@s;$r6z8h=R)^^y5Ibm@R*XKsMXFv6m!=KhI
z>ety?cO;00@QUokUL#M@B2}N8F<1X_CgIrsPXe1V->}y^Xi1zge{h
z+JQr0D!P`mTAxvP9IKgIpf?~nps7{hF=ti5IpL*vb+M~tt@G)eRrPhDH7s~>cw~z+
zG2~T9Z>+3Nq{f_8rN}6`Za@j1Wo>27s!FW_4vxJ{bV)WrdVl7usyVAFbylcInfr5A
zmA!j&Ru%SEH{LU6Rq3p_eZ;orc6ej@<>`~Zkf-?+kBxGt=d3Du
z`L+WtSSR@9_=If3Ijaf`v;8|~Rn1vd(z%u4g3V8i@(w--2p{tDM{%#f(-
ztyr;q^U`|zvCsXGQPS@Y7CJR?*S1oh@QvySs*IpKapFltnALOE*7WJNIv)#r#!+4J
z0ptyrM}}|Mo9=m4=M{=Wg7MD3Gp`(aYB(an&Pz5EuTZwpeykI1pY_)Hddx~EB^2dV
zoH@@hs~@dMwgP7@TrXU6bR-+UtaWF_keM@b5ZwU(Hl77MyUH_C-`XW*?h9^2Z!^a<
zA<|t}W*8AGb_QbZeK%Ia{D$l$=Y%pJGKZD7jnvv)$SuJcs|njokn%|!vc(M#&5#@c
ztTLT;^3EzE>3z|8b!V&o|LmPxj}=Fj?|H5BIxo}Gj8sDEwTNl<0&7>2q3hajAl#^TAcljcEy8ddJ)
z%|KDPJm3VdIw(yLSpqk}C%S{N*WX0eoN@!GBcu9;T{2GcZs8)jm`|S24w|xng7S9Q
z;jH~IE-QkpzG}*Lgsu|qzHdm4%oII_IcvR!0J$0S@Shs#A->Zo$~k%9=gk<|SP9jT
zcY*Jl{nx6YG$?OuZ-FIP;hx#1{Wj-^${J>q_3ENnWu-Dl)vbtXiO6WZW6)^#Jb5Gd
z(BiY9<2?e}dZVI1$`kr(zlsa)CJ%`yM_nrf0p`IPd53e0>O2(2OFSIHgiP5n=ndBg
zuJ!X2Wf?bTC-N_eV$9v3NUksnvw?<8v8$2A4~8O@RA)=2*X=^*%I(76
zP{(t=i0?T@!Y?sT@zn9Be~3AdaR=Yftd*^92wB<#LAXoy3v8)kD^45U4vmA1HDfNF
znuLjY4Z3kYyhjwjYGu)l0y$Niy}kV56F%nXK+^Ieo%$_{iXcn$y
zL0a?eUG`)rMWfNtl}Y-gsah0+CR1lEi?i3dCiW|LypPqv2_x_yIsER{24{Zh;w(Y3
zEXldtdMPX3j`|HnhooUjD?1<9XDNm3r2;s0M5O++;=@ht{A}dbq7}3`-qa$ozsUOf
zBvIsmDmOj}h%hp`Igw7X&M(ARr(S06PE{N55#VziRU2eLPjwqp-3ApHUj;54u5P2&
zK35{sPP#Gj8~yySH5-dH^&x+eu8f6b&z@vyYdK>3$j9qSe|VTIQLIFBDOpc&awHyvF`7?-;9i@EMxNRq;?)Kg_Fn;C*L{6m~TM{LoeyFr%IokEiith;#S}
z1Y_CsXnJ_$xQ2zUJPh}uu`O?X5S=com1Ucc`3o`&jjk_HYucb|p)s9{H_0P(CO+}|
z{Zgw%MOgN_r3>_!T+Rw&=tka0&)6%L*Hx6{YOA_KpKNLQLlP!!b*SM%Mjcjp+TZ9H
znm?yUY5rKhypJBAKb_Dy3Rea?#o5?FbO!cOf7KV;kUkrpcvb$4?J>SU^cQ;|@a5s}
zk*AKGV9wly2LgYYYUD+A)XlW}YiQvefUDYbXWJMb8(KVhAvzILwSrtvXtdo&MbTsP
z60f1Q%KcpMKlVGf++2DXnF`t+OzaUp$gg315O$2ORX#;g3Z9D0^5=ebOY3*%)^=^|
zLOVLOhxVKUOZZssuxCHfFnDFGd|PQ4Yr*t-q1(yi$8U~S)u~@L50S7pORbtsPp=85
z3FqGhoIiWwQju45a|E-D|8y(R-lBP;d2I#0J+7!Y`FYCOaoq!oc@o+|E2I3TeA@Eu
zb!sL58Y^;b^t-E-r1oqb`RIR7Mr;5*3a_6%h3uIk3S<6%OI*6?C@uy^T?=gmT+&If
z_}Tghb^6{MkFI;685?*nnoQzL#)^?G;F#lxX>U0`DcO$XNo{C;-^lbf;S)#40v1^(
z=+!E&g(gAV=~1klxRQI*NTaL0u50ld=Uz%jV%F%CJohvz;qOK~MsaQFzrTb&c`J58
z-RrDRsBa*vVL!P4LC_a|GV-e1K7u+y>sHaZO>?KtZQ}^k~O*!FxsTpOt;({914UAL5-t2UteX!N$sAAxLy9)pU)|KH{O#cx-JQV(&SUtUIhJBfd
z$&7iqQWoe`3U*oC&~`ldO_ygd2dN}U!Fz4V)qcbDP(fjgrJM|Kj`xEnvV)x$0A^)>
zayo*NwII!DlA!kUncGYfR66D)LFE_3moiCEPH)0DBG1z#K^No1Bth}s{uU(Y7a?(x
z!el$jXN{dBJH~n(Tm)STeI5ye&aG?cK$78T(bg!@i}*(T`f0K3qywToNI#GkLG0hU
zG#Uf;CSDK4572HswyaNKB^r^nT6AV<_2>=sc4n679@aklj7jNl_PwMx?ci@fL&5q&
zcfyt(n?DwC##*>f_OjPyNCPouu30NfqRyBF>d(S9MuVeHs^7EpE7iLoc#iL!_JPj)
zhI#I#ASg)LXLNA*Xr*J&(JHQnE+>r2s)M%5f*CL?v`}f2+IO%EeO0Tj<27K8tQMfI
z{Anh>_ux>xX2zms1tyKot-o^Z}K^f=W|I5**3
zE1`x|IPp1T7CxC;>K{L!Bd!6DmgfL3E?klRobsKMm2oX}tQ*mtbf;vc;3?qD2YJbM
z<8HWc`f*oI{5iNcS%uueD9Sv+e}t})&*H^vpFI_KV*@Z=+w=2QfdtZ5J%8~KyXP-{
zuloB+=s{Gv%=a0s#F}cK7c`mPF=swo^O0m#bl<)eXKEmuQwwq>^jh&Xab}EkGGa^Q
zV(O&TFbqEzK0t77C4Q^@jPgK}qtp8unZ?08
zg*-9alZR)Kp4E#wE!UnmddyOP%87C0i$NV|OV7&{`LDTL38ztZE*H9gbWh(6n7PZZ
z`ui&9El15y^@!tF>pSYt