Performance: keep the sources as cell array for as long a possible (#841)

Reduce the required RAM when browsing results (sources) which are saved as matrix decomposition (see #733)

---------

Co-authored-by: rcassani <raymundo.cassani@gmail.com>

Author: Edouard2laire

toolbox/core/bst_memory.m

@@ -2189,17 +2189,35 @@ function LoadResultsMatrix(iDS, iResult)
     % ===== GET RESULTS VALUES =====
     % === FULL RESULTS ===
     if ~isempty(GlobalData.DataSet(iDS).Results(iResult).ImageGridAmp)
-        % Get ImageGridAmp interesting sub-part
-        if isempty(iRows)
-            ResultsValues = double(GlobalData.DataSet(iDS).Results(iResult).ImageGridAmp(:, iTime));
-            if ~isempty(GlobalData.DataSet(iDS).Results(iResult).Std)
-                Std = double(GlobalData.DataSet(iDS).Results(iResult).Std(:, iTime, :, :));
+        % ImageGridAmp = [nSources, nTimes]
+        if isnumeric(GlobalData.DataSet(iDS).Results(iResult).ImageGridAmp)
+            % Get ImageGridAmp interesting sub-part
+            if isempty(iRows)
+                ResultsValues = double(GlobalData.DataSet(iDS).Results(iResult).ImageGridAmp(:, iTime));
+                if ~isempty(GlobalData.DataSet(iDS).Results(iResult).Std)
+                    Std = double(GlobalData.DataSet(iDS).Results(iResult).Std(:, iTime, :, :));
+                end
+            else
+                ResultsValues = double(GlobalData.DataSet(iDS).Results(iResult).ImageGridAmp(iRows, iTime));
+                if ~isempty(GlobalData.DataSet(iDS).Results(iResult).Std)
+                    Std = double(GlobalData.DataSet(iDS).Results(iResult).Std(iRows, iTime, :, :));
+                end
             end
-        else
-            ResultsValues = double(GlobalData.DataSet(iDS).Results(iResult).ImageGridAmp(iRows, iTime));
-            if ~isempty(GlobalData.DataSet(iDS).Results(iResult).Std)
-                Std = double(GlobalData.DataSet(iDS).Results(iResult).Std(iRows, iTime, :, :));
+        % ImageGridAmp = {[nSources,a], [a,b], [b, nTimes]}
+        elseif iscell(GlobalData.DataSet(iDS).Results(iResult).ImageGridAmp)
+            % Get ImageGridAmp interesting sub-part
+            ResultsValues = GlobalData.DataSet(iDS).Results(iResult).ImageGridAmp;
+            assert(isempty(GlobalData.DataSet(iDS).Results(iResult).Std), 'Storing Std a cell is not supported yet.')
+            ResultsValues{end}   = ResultsValues{end}(:, iTime);
+
+            if ~isempty(iRows)
+                ResultsValues{1} = ResultsValues{1}(iRows, :);
             end
+            
+            % Compute full results
+            ResultsValues = double(bst_multiply_cellmat(ResultsValues));
+            % Std should be empty if ImageGridAmp is a cell
+            Std = [];
         end
     % === KERNEL ONLY ===
     elseif ~isempty(GlobalData.DataSet(iDS).Results(iResult).ImagingKernel)
@@ -2586,13 +2604,29 @@ function LoadResultsMatrix(iDS, iResult)
         [maxGFP, iMax] = max(GFP);
         % Get the results values at this particular time point
         sources = GetResultsValues(iDS, iResult, [], iMax);
+        % Store minimum and maximum of displayed data
+        DataMinMax = [min(sources(:)), max(sources(:))];
     % Full results
     else
-        % Get the maximum on the full results matrix
-        sources = GlobalData.DataSet(iDS).Results(iResult).ImageGridAmp;
+        % Get the maximum on the full results matrix (process by time blocks)
+        nSamples = GlobalData.DataSet(iDS).Results(iResult).NumberOfSamples;
+        if isnumeric(GlobalData.DataSet(iDS).Results(iResult).ImageGridAmp)
+            nSources = size(GlobalData.DataSet(iDS).Results(iResult).ImageGridAmp, 1);
+        elseif iscell(GlobalData.DataSet(iDS).Results(iResult).ImageGridAmp)
+            nSources = size(GlobalData.DataSet(iDS).Results(iResult).ImageGridAmp{1}, 1);
+        end
+        ProcessOptions = bst_get('ProcessOptions');
+        MaxSizeDouble = ProcessOptions.MaxBlockSize;
+        blockSize  = max(floor(MaxSizeDouble / nSources), 1);
+        nBlocks    = ceil(nSamples / blockSize);
+        DataMinMax = [Inf, -Inf];
+        for iBlock = 1 : nBlocks
+            iTime = [((iBlock-1) * blockSize + 1) : min(iBlock * blockSize, nSamples)];
+            sources = GetResultsValues(iDS, iResult, [], iTime);
+            DataMinMax(1) = min(DataMinMax(1), min(sources(:)));
+            DataMinMax(2) = max(DataMinMax(2), max(sources(:)));
+        end
     end
-    % Store minimum and maximum of displayed data
-    DataMinMax = [min(sources(:)), max(sources(:))];
 end
 
 

toolbox/io/in_bst.m

@@ -89,8 +89,17 @@
         % FULL RESULTS
         if isfield(sMatrix, 'ImageGridAmp') && ~isempty(sMatrix.ImageGridAmp)
             iTime = GetTimeIndices(TimeBounds, sMatrix.Time);
-            sMatrix.ImageGridAmp  = sMatrix.ImageGridAmp(:,iTime);
             sMatrix.Time          = sMatrix.Time(iTime);
+            % Get results for iTime
+            if isnumeric(sMatrix.ImageGridAmp)
+                sMatrix.ImageGridAmp  = sMatrix.ImageGridAmp(:,iTime);
+            elseif iscell(sMatrix.ImageGridAmp)
+                % ImageGridAmp = {[nSources,a], [a,b], [b, nTimes]}
+                sMatrix.ImageGridAmp{end} = sMatrix.ImageGridAmp{end}(:, iTime);
+                if isLoadFull
+                    sMatrix.ImageGridAmp = bst_multiply_cellmat(sMatrix.ImageGridAmp);
+                end
+            end
             sMatrix.ImagingKernel = [];
             matName = 'ImageGridAmp';
         % KERNEL ONLY

toolbox/io/in_bst_results.m

@@ -261,13 +261,8 @@
         end
     end
 % If full results are saved as factor decomposition
-elseif isfield(Results,'ImageGridAmp') && iscell(Results.ImageGridAmp)
-    % ImageGridAmp = ImageGridAmp{1} * ImageGridAmp{2} * ... * ImageGridAmp{N}
-    tmp = Results.ImageGridAmp{1};
-    for iDecomposition = 2 : length(Results.ImageGridAmp)
-        tmp = tmp * Results.ImageGridAmp{iDecomposition};
-    end
-    Results.ImageGridAmp = full(tmp);
+elseif LoadFull &&  isfield(Results,'ImageGridAmp') && iscell(Results.ImageGridAmp)
+    Results.ImageGridAmp = bst_multiply_cellmat(Results.ImageGridAmp);
 end
 
 

toolbox/math/bst_multiply_cellmat.m

@@ -0,0 +1,103 @@
+function M = bst_multiply_cellmat(X)
+% BST_MULTIPLY_CELLMAT: Multiply all the matrices in cell array
+%
+% USAGE: M = bst_multiply_cellmat(X)
+%
+% INPUT:
+%    - X        : 1D cell array of 2D matrices to be multiplied
+%
+% OUTPUT:
+%    - M : Result of chain matrix multiplication
+%          M = X{1} * X{2} * ... * X{N}
+
+% @=============================================================================
+% This function is part of the Brainstorm software:
+% https://neuroimage.usc.edu/brainstorm
+% 
+% Copyright (c) University of Southern California & McGill University
+% This software is distributed under the terms of the GNU General Public License
+% as published by the Free Software Foundation. Further details on the GPLv3
+% license can be found at http://www.gnu.org/copyleft/gpl.html.
+% 
+% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
+% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
+% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
+% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
+% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
+%
+% For more information type "brainstorm license" at command prompt.
+% =============================================================================@
+%
+% Authors: Edouard Delaire, 2025
+%          Raymundo Cassani, 2025
+
+% Do nothing
+if ~iscell(X) || ~isvector(X)
+    M = X;
+    return
+end
+
+nMat = length(X);
+% One matrix
+if nMat == 1
+    M = X{1};
+    return
+end
+
+% Check sizes
+matSizes = cellfun(@size, X, 'UniformOutput', 0);
+if ~all(cellfun(@(x) length(x)==2, matSizes))
+    error('All matrices in cell array must be 2D')
+end
+dimDiff = diff([matSizes{1 : +1 : nMat}]);
+if ~all(dimDiff(2:2:nMat) == 0)
+    error('Matrices cannot be multiplied')
+end
+
+% Choose direction of association for multiplication
+if nMat == 3
+    % Cost is the number of scalar multiplications
+    cost1 = matSizes{1}(1) * matSizes{1}(2) * matSizes{2}(2) + matSizes{1}(1) * matSizes{2}(2) * matSizes{3}(2);  % (A*B)*C
+    cost2 = matSizes{2}(2) * matSizes{2}(1) * matSizes{3}(2) + matSizes{1}(1) * matSizes{2}(1) * matSizes{3}(2);  % A*(B*C)
+    fromRight = cost2 < cost1;
+else
+    % Heuristic rule
+    % Note, the optimal association approach can be obtained computing the cost recursively
+    fromRight = size(X{end}, 2) < size(X{1}, 1);
+end
+
+% Keep useful indices for inner dimensions
+if fromRight
+    for im = nMat : -1 : 2
+        usefulInnerIx       = any(X{im}, 2);
+        X{im}   = X{im}(usefulInnerIx, :);
+        X{im-1} = X{im-1}(:, usefulInnerIx);
+    end
+else
+    for im = 1 : nMat - 1
+        usefulInnerIx       = any(X{im}, 1);
+        X{im}   = X{im}(:, usefulInnerIx);
+        X{im+1} = X{im+1}(usefulInnerIx, :);
+    end
+
+end
+
+% Multiply matrices 
+if fromRight
+    % M = X{1} * ( X{2} * ( X{3} * X{4} ) )
+    M = X{end};
+    for iDecomposition = (nMat - 1) : -1 : 1
+        M = X{iDecomposition} * M;
+    end
+else
+    % M = ( ( X{1} * X{2} ) * X{3} ) * X{4}
+    M = X{1};
+    for iDecomposition = 2 : nMat
+        M = M * X{iDecomposition};
+    end
+end
+
+% Ensure results are full matrices
+M = full(M);
+end
+

toolbox/process/functions/process_extract_scout.m

@@ -665,7 +665,11 @@
             % FULL RESULTS
             if isfield(sResults, 'ImageGridAmp') && ~isempty(sResults.ImageGridAmp)
                 if nargout > 1
-                    matSourceValues = sResults.ImageGridAmp;
+                    if isnumeric(sResults.ImageGridAmp)
+                        matSourceValues = sResults.ImageGridAmp;
+                    elseif iscell(sResults.ImageGridAmp)
+                        matSourceValues = bst_multiply_cellmat(sResults.ImageGridAmp);
+                    end
                 end
                 % Drop large data field.
                 sResults = rmfield(sResults, 'ImageGridAmp');