Minimum of three iterations in robust re-reference?

[a-hurst]

Hi there,

I think I've discovered a small bug in the re-referencing logic of the PREP pipeline, but I wanted to make sure it isn't intentional behaviour before modifying anything on my end.

EEG-Clean-Tools/PrepPipeline/utilities/robustReference.m

Lines 58 to 87 in 3ed337e

	%% Remove reference from signal iteratively interpolating bad channels
	iterations = 0;
	noisyChannelsOld = [];
	while true % Do at least 1 iteration
	noisyStatistics = findNoisyChannels(signalTmp, referenceOut);
	referenceOut.badChannels = ...
	updateBadChannels(referenceOut.badChannels, noisyStatistics.noisyChannels);
	noisyChannels = referenceOut.badChannels.all(:)';
	if (iterations > 1 && (isempty(noisyChannels) ||...
	(isempty(setdiff(noisyChannels, noisyChannelsOld)) ...
	&& isempty(setdiff(noisyChannelsOld, noisyChannels))))) || ...
	iterations > referenceOut.maxReferenceIterations
	break;
	end
	noisyChannelsOld = noisyChannels;
	sourceChannels = setdiff(referenceOut.referenceChannels, noisyChannels);
	if length(sourceChannels) < 2
	error('robustReference:TooManyBad', ...
	'Could not perform a robust reference -- not enough good channels');
	end
	if ~isempty(noisyChannels)
	signalTmp = interpolateChannels(signal, noisyChannels, sourceChannels);
	else
	signalTmp = signal;
	end
	referenceSignal = nanmean(signalTmp.data(referenceChannels, :), 1);
	signalTmp = removeReference(signal, referenceSignal, referenceChannels);
	iterations = iterations + 1;
	fprintf('Iteration: %d\n', iterations);
	end

In the above loop, iterations starts at zero and only get incremented at the end, after the "check if we should break the loop" logic. As a consequence, it means that the code actually does a maximum of 5 iterations by default, even though the default maximum is documented to be 4. Additionally, it means that PREP needs to do at least three iterations of this loop, even if the detected bad channels are unchanged between the first and second iterations. Looking at the original PREP publication, the pseudocode for this loop doesn't say anything about a minimum of iterations, so I'm wondering whether this is a bug or an intentional behaviour I don't understand.

Thanks in advance, and thank you for creating and maintaining this software!

[VisLab]

I guess this is a case of lacking truth in advertising..... maxReferenceIterations really means how many "extra" ones.

The first call to findNoisyChannels (line 28) doesn't count as it only uses the globally bad channels to get things started. It computes the noisy channels based on the "faulty" average reference (line 62). Even if you don't find any bad channels, we don't exit because we want to compute at least one reference based on all of the methods (including Ransac) and then compute the bad channels from that. Thus there are a minimum of three calls to findNoisyChannels even when maxReferenceIterations is 0. This is because of the "chicken and egg" problem in getting the process started. We could have settled for fewer iterations, but some of our tests showed that for some datasets, the bad channel set really doesn't start to stabilize until after those three iterations. For most datasets this isn't necessary, but we thought it better to be on the safe side.

[a-hurst]

@VisLab Thanks for taking a look at this! I'm still somewhat confused though:

Even if you don't find any bad channels, we don't exit because we want to compute at least one reference based on all of the methods (including Ransac) and then compute the bad channels from that.

Unless I'm misreading something, the code still requires an additional 3 iterations of findNoisyChannels after the call on line 28 before it can break out of the loop. See the below simplified code to see what I mean:

 iterations = 0;                          
 while true  % Do at least 1 iteration 
     noisyStatistics = findNoisyChannels(signalTmp, referenceOut); 
     if (iterations > 1)
         break; 
     end  
     iterations = iterations + 1; 
 end 

Since iterations starts at 0 and increments at the end of the loop, and the iterations check happens after findNoisyChannels in the loop, iterations > 1 won't be true until the 3rd pass of the loop and thus the third iteration of findNoisyChannels. Is this still the expected behaviour?