Getting the coil-channel .tra matrix for the HCP-MEG preprocessed files about fieldtrip HOT 5 CLOSED

% Trying to rebuild missing from the HCP pre-processed files. 
% We want to reconstruct separate matrices from the .tra matrix, with the help of the .balance field:
%   - the sparse coil-channel correspondance matrix
%   - the noise-compensation matrix applied by 4D acquisition system
%   - the PCA/ICA cleaning procedures
%
% Authors: Francois Tadel, 2017

% Load example structure: Subject HCP/891667, resting state recordings pre-processed with FieldTrip
ftMat = load('891667_MEG_3-Restin_rmegpreproc_5trials.mat');
ftMat = ftMat.data;
grad = ftMat.grad;


% === FIX STRUCTURE ===
% Suggestion from JM: Fix glitches in the grad structure with ft_datatype_sens
% Reference: http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=3298#c6
% Problem: does not fix retroactively the history in the .balance field
grad = ft_datatype_sens(grad);


% === REVERT TO THE ORIGINAL .TRA MATRIX ===
% Suggestion from JM: Use ft_apply_montage to revert to the original .tra matrix
% Reference: http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=3298#c3
%
% Problems: There are lots of inconsistencies in the .balance fields of the HCP files:
%   1) The output channel units and types in transformation "pca" are set to "unknown"
%   2) Random naming chanunitorg/chanunitold fields
%   3) Channel types change randomly between "meg" (eg. pca) and "megmag" (eg. invcomp)
% Reference: http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=3298#c5

% Get the list of montages to undo
montageList = {grad.balance.current};
if ~isempty(grad.balance.previous)
    montageList = cat(2, montageList, grad.balance.previous{:});
end

% Undo montages one by one
for iMontage = 1:length(montageList)
    % Make sure the montage is defined in the structure
    mont = montageList{iMontage};
    if ~isfield(grad.balance, mont)
        continue;
    end
    
    % ISSUE #1: Trying to fix units and types  (with randomness of issue #2)
    if isfield(grad.balance.(mont), 'chanunitnew') && strcmpi(grad.balance.(mont).chanunitnew{1}, 'unknown')
        if isfield(grad.balance.(mont), 'chanunitorg')
            grad.balance.(mont).chanunitnew = grad.balance.(mont).chanunitorg;
        elseif isfield(grad.balance.(mont), 'chanunitold')
            grad.balance.(mont).chanunitnew = grad.balance.(mont).chanunitold;
        end
    end
    if isfield(grad.balance.(mont), 'chantypenew') && strcmpi(grad.balance.(mont).chantypenew{1}, 'unknown')
        if isfield(grad.balance.(mont), 'chantypeorg')
            grad.balance.(mont).chantypenew = grad.balance.(mont).chantypeorg;
        elseif isfield(grad.balance.(mont), 'chantypeold')
            grad.balance.(mont).chantypenew = grad.balance.(mont).chantypeold;
        end
    end
    
%     % ISSUE #3: Trying to fix the channel types
%     % DOES NOT WORK BECAUSE THE ORDER OF THE CHANNELS CHANGE
%     if isfield(grad.balance.(mont), 'chantypenew') && (length(grad.balance.(mont).chantypenew) == length(grad.chantype))
%         grad.balance.(mont).chantypenew = grad.chantype;
%     end
%     if isfield(grad.balance.(mont), 'chantypeorg') && (length(grad.balance.(mont).chantypeorg) == length(grad.chantype))
%         grad.balance.(mont).chantypeorg = grad.chantype;
%     end
%     if isfield(grad.balance.(mont), 'chantypeold') && (length(grad.balance.(mont).chantypeold) == length(grad.chantype))
%         grad.balance.(mont).chantypeold = grad.chantype;
%     end
    
    % Reverse transformation
    grad = ft_apply_montage(grad, grad.balance.(mont), 'inverse', 'yes', 'keepunused', 'yes');
end


% === REBUILD ORIGINAL COIL-CHANNEL MATRIX ===
% Remove small values (probably due to the multiple matrix inversions) to keep only the ones
grad.tra(abs(grad.tra) < 0.1) = 0;

% For each channel
for i = 1:length(ftMat.label)
    % MEG sensors
    if isfield(ftMat, 'grad') && ~isempty(ftMat.grad) && isfield(ftMat.grad, 'tra') && isfield(ftMat.grad, 'label') && ~isempty(ftMat.grad.label) && ismember(ftMat.label{i}, ftMat.grad.label)
        % Find channel index
        ichan = find(strcmpi(ftMat.label{i}, grad.label));
        % Find corresponding coils
        icoils = find(grad.tra(ichan,:));
        
        % Error: Two many coils
        if (length(icoils) > 2)
            error(['Wrong number of coils for channel ', ftMat.label{i}, ': Cannot import this file correctly...']);
        end
        % Get sensor properties
        Loc    = ftMat.grad.coilpos(icoils,:)';
        Orient = ftMat.grad.coilori(icoils,:)';
        Weight = ftMat.grad.tra(ichan,icoils);
    end
end