Potential Bug? Freezing while making plots iteratively

Hi community,

I am at a loss for why my matlab code is freezing up, I am not using that much of matlab's capabilities...I am attempting to read through a directory, and create a tiledlayout 2x2 plot for each folder in the directory. In this, I create an exponential trendline to fit some of my data, and once I try plotting all of the stuff in the tiled layout, it can cause the system to freeze up right before saving the picture, preventing the for loop from reaching the next file. It shouldn't be a code failure, as this code has worked in pieces before, but causing the entire app become entirely unresponsive for hours requiring a force quit on a 2022 computer makes me relatively surprised.

For the purposes of ease of viewing, I will attach the code in a .m, and paste the file's full code here, and I will also attach a link to a compressed folder of a few files that will allow one to run this code.

To be clear the freeze and crash happens sporadically, sometimes on the first file, sometimes on the 5 file, which makes me believe there is a bug of some sorts! I see it more frequently occur when I have a debug run-stop at the last end line after the saveas function, so that I can view the figure for any issues before moving to the next folder. I appreciate the help!

Link to the compressed file folder: https://drive.google.com/file/d/1_8fpljROfMsGG5oMc1oIvaM0xDyrKlj0/view?usp=sharing

clear all 
close all
clc

%% Change filepath and timestep frequency

filepath = '/Users/nickscott/Library/CloudStorage/Box-Box/Parekh Lab Public/Nick/CellProfiler/ProperFPS_NikonExp';
cd '/Users/nickscott/Library/CloudStorage/Box-Box/Parekh Lab Public/Nick/CellProfiler/ProperFPS_NikonExp/'
topLevelFolder = pwd;
files = dir(topLevelFolder);
dirFlags = [files.isdir];
% Extract only those that are directories.
subFolders = files(dirFlags); % A structure with extra info.
% Get only the folder names into a cell array.
folderList = {subFolders(3:end).name}; % Start at 3 to skip . and ..

for i= 1: length(folderList)
    thisFolder = folderList{i};
    fprintf('Processing folder %s\n', thisFolder);
    % Get txt files.
    filePattern = sprintf('%s/*.txt', thisFolder);
    baseFileNames = dir(filePattern);
    bFN = struct2cell(baseFileNames)
    dt{i}= bFN{1}
end

%% Constants

FigList = {'Median Trendline of Perinuclear Intensity', 'Median Trendline of Nuclei Intensity'...
    'Best Fit Using Exponential Function', 'Median Trendline of Nuclei/Perinuclear Intensity'};

PeriRing = '/ShNucDiffPipe_Perinuclear_Ring.csv'; %Perinuclear Ring file name
Nuc = '/ShNucDiffPipe_ShrunkenNuclei.csv'; %Stained Nuclei file name

%% Giant 'for loop' for every folder in the batch

for p= 1:length(dt)  
%% Naming Extraction

    newstr= extractBetween([filepath '/' folderList{p} '/'], "NikonExp/","/"); %Pull Metadata folder name
    FigTitle = regexprep(newstr, '_', ' '); %Change _ to spaces


%% Read Table, Extract Variables of Interest
%Automatically extract the name

    TPR = readtable([filepath '/' folderList{p} PeriRing], 'Delimiter',','); %perinuclear table
    TN = readtable([filepath '/' folderList{p} Nuc], 'Delimiter',','); %nuclei table
    TT = readtable([filepath '/' folderList{p} '/' char(dt{1,p})], 'Delimiter', 'tab');

% Modify Time frame to start out at 0 instead of being image set index from
% ... Cellprofiler, may need to change if pipeline can be adjusted
    TPR{:,1}=TPR{:,1}-TPR{1,1};
    TN{:,1}=TN{:,1}-TN{1,1};
    


%Extract Unique nuclei objects (must be in column 2)
    [~,~,j]=unique(TPR(:,2));
    ArrayPR = accumarray(j,(1:length(j)).', [],  @(k) {TPR(k,:)});
    [~,~,l]=unique(TN(:,2));
    ArrayN = accumarray(l,(1:length(l)).', [],  @(k) {TN(k,:)});

    
    table_heightsPR = cellfun(@height,ArrayPR);
    table_heightsN = cellfun(@height,ArrayN);
    height_thresholdPR = 0.75*max(table_heightsPR);
    height_thresholdN = 0.75*max(table_heightsN);
    LineNumb = 0;

    if contains(TT.Properties.VariableNames{2}, 'Time_s_ms_')
        for i= 1:length(table_heightsN)
            if height(ArrayN{i}) >= height_thresholdN
                ArrayN{i,1}.Timestamp = zeros(table_heightsN(i),1);
                [row,col] = size(ArrayN{i,1});
                LineNumb = LineNumb + 1;
                for j=1:row
                    FrameNum = ArrayN{i,1}.ImageNumber(j);
                    FrameNum = FrameNum+1;
                    ind1 = find(TT.Index == FrameNum);
                    ArrayN{i,1}.Timestamp(j) = sum(sscanf(num2str(TT.Time_s_ms_(ind1)),'%f.%f').*[1;0.0001]);         
                end
            else
            end
        end
    
        for i= 1:length(table_heightsPR)
            if height(ArrayPR{i}) >= height_thresholdPR
                ArrayPR{i,1}.Timestamp = zeros(table_heightsPR(i),1);
                [row,col] = size(ArrayPR{i,1});
                for j=1:row
                    FrameNum = ArrayPR{i,1}.ImageNumber(j);
                    FrameNum = FrameNum+1;
                    ind2 = find(TT.Index == FrameNum);
                    ArrayPR{i,1}.Timestamp(j) = sum(sscanf(num2str(TT.Time_s_ms_(ind2)),'%f.%f').*[1;0.0001]);   %sum(sscanf(char(TT.Time_m_s_ms_(ind1)),'%f:%f').*[60;1]);
                end
            else
            end
        end
    
    elseif contains(TT.Properties.VariableNames{2}, 'Time_m_s_ms_')
        for i= 1:length(table_heightsN)
            if height(ArrayN{i}) >= height_thresholdN
                ArrayN{i,1}.Timestamp = zeros(table_heightsN(i),1);
                [row,col] = size(ArrayN{i,1});
                LineNumb = LineNumb + 1;
                for j=1:row
                    FrameNum = ArrayN{i,1}.ImageNumber(j);
                    FrameNum = FrameNum+1;
                    ind1 = find(TT.Index == FrameNum);
                    ArrayN{i,1}.Timestamp(j) = sum(sscanf(char(TT.Time_m_s_ms_(ind1)),'%f:%f').*[60;1]);
                end
            else
            end
        end
    
        for i= 1:length(table_heightsPR)
            if height(ArrayPR{i}) >= height_thresholdPR
                ArrayPR{i,1}.Timestamp = zeros(table_heightsPR(i),1);
                [row,col] = size(ArrayPR{i,1});
                for j=1:row
                    FrameNum = ArrayPR{i,1}.ImageNumber(j);
                    FrameNum = FrameNum+1;
                    ind2 = find(TT.Index == FrameNum);
                    ArrayPR{i,1}.Timestamp(j) = sum(sscanf(char(TT.Time_m_s_ms_(ind2)),'%f:%f').*[60;1]);
                end
            else
            end
        end

    elseif contains(TT.Properties.VariableNames{2}, 'Time_s_')
        for i= 1:length(table_heightsN)
            if height(ArrayN{i}) >= height_thresholdN
                ArrayN{i,1}.Timestamp = zeros(table_heightsN(i),1);
                [row,col] = size(ArrayN{i,1});
                LineNumb = LineNumb + 1;
                for j=1:row
                    FrameNum = ArrayN{i,1}.ImageNumber(j);
                    FrameNum = FrameNum+1;
                    ind1 = find(TT.Index == FrameNum);
                    ArrayN{i,1}.Timestamp(j) = sum(sscanf(num2str(TT.Time_s_(ind1)),'%f').*[1]);
                end
            else
            end
        end
    
        for i= 1:length(table_heightsPR)
            if height(ArrayPR{i}) >= height_thresholdPR
                ArrayPR{i,1}.Timestamp = zeros(table_heightsPR(i),1);
                [row,col] = size(ArrayPR{i,1});
                for j=1:row
                    FrameNum = ArrayPR{i,1}.ImageNumber(j);
                    FrameNum = FrameNum+1;
                    ind2 = find(TT.Index == FrameNum);
                    ArrayPR{i,1}.Timestamp(j) = sum(sscanf(num2str(TT.Time_s_(ind2)),'%f').*[1]);
                end
            else
            end
        end
    
    end
 

    
%% Plot to the tiled subplots *if* the nuclei doesn't float in and out of frame significantly 
    Tile= tiledlayout(2,2);
    set(gcf, 'WindowState', 'maximized');

    nexttile
    hold on
        for i= 1:length(ArrayPR) %i=Nucleus number
            if height(ArrayPR{i}) >= height_thresholdPR
                plot(ArrayPR{i,1}.Timestamp,ArrayPR{i,1}.Intensity_MedianIntensity_Alexa488,'color',[0,0,0]+0.75)
            else
            end
        end
        
    % median intensity for each image number calculation:
    % first, concatenate all the tables into one big table:
    % t = vertcat(ArrayPR{:});

    % or concatenate just the tables whose height is >= height_threshold:
    t = vertcat(ArrayPR{table_heightsPR >= height_thresholdPR});

    % get all the image numbers:
    all_image_time = t.Timestamp;
   
    % get all the intensities:
    all_intensity = t.Intensity_MedianIntensity_Alexa488;
    mask = isnan(all_intensity);
    all_intensity(mask) = 0;
    % get the unique image numbers:
    u_image_time = unique(all_image_time);

    % calculate the median intensity for each image number:
    % number of unique image numbers:
    n_images = numel(u_image_time);

    % pre-allocating the median intensity array:
    median_intensity_pr = zeros(1,n_images);

    % for each unique image number:
    for f = 1:n_images
        % get a logical index saying whether each element of all_image_number
        % is equal to this image number (idx is true where 
        % all_image_number == u_image_number(i) and false elsewhere):
        idx_pr = all_image_time == u_image_time(f);
        % use that logical index to calculate the median intensity for this
        % image number:
        median_intensity_pr(f) = median(all_intensity(idx_pr));
    end
    % plot
    plot(u_image_time,median_intensity_pr,'k','LineWidth',2); %Black median line
    
    title(FigList{1})
    hold off


    nexttile
    hold on
        for i= 1:length(ArrayN) %i=Nucleus number
            if height(ArrayN{i}) >= height_thresholdN
            plot(ArrayN{i,1}.Timestamp,ArrayN{i,1}.Intensity_MedianIntensity_Alexa488,'color',[0,0,0]+0.75)
            else
            end
        end
        
    % median intensity for each image number calculation:
    % first, concatenate all the tables into one big table:
    % t = vertcat(ArrayPR{:});

    % or concatenate just the tables whose height is >= height_threshold:
    N = vertcat(ArrayN{table_heightsN >= height_thresholdN});

    % get all the image numbers:
    all_image_time_N = N.Timestamp;

    % get all the intensities:
    all_intensity_N = N.Intensity_MedianIntensity_Alexa488;

    % get the unique image numbers:
    u_image_time_N = unique(all_image_time_N);

    % calculate the median intensity for each image number:
    % number of unique image numbers:
    n_images_N = numel(u_image_time_N);

    % pre-allocating the median intensity array:
    median_intensity_N = zeros(1,n_images);

    % for each unique image number:
    for f = 1:n_images
        % get a logical index saying whether each element of all_image_number
        % is equal to this image number (idx is true where 
        % all_image_number == u_image_number(i) and false elsewhere):
        idx = all_image_time_N == u_image_time_N(f);
        % use that logical index to calculate the median intensity for this
        % image number:
        median_intensity_N(f) = median(all_intensity_N(idx));
    end
    % plot
    plot(u_image_time_N,median_intensity_N,'k','LineWidth',2); %Black median line
    title(FigList{2})
    hold off

    
    %plot(t(indexOfInterest),signal(indexOfInterest)) % plot on new axes
    %if height(ArrayN{i}) >= height_thresholdN
        indexOfInterest = (ArrayN{1,1}.Timestamp <= 12.5) & (ArrayN{1,1}.Timestamp >= 2.4);
        IoI1 = find(indexOfInterest,1,'first');
        nexttile
        hold on
                
        xdata = u_image_time(IoI1:127);
        ydata = transpose(median_intensity_N(IoI1:127)./median_intensity_pr(IoI1:127));
        fun = @(cs,xdata) cs(1)-cs(2)*exp(-cs(3)*xdata);
        x0 = [1.5, 2.5, 0.3];
        lb = [-Inf, 0,0];
        ub = [Inf, Inf, Inf];
        x = lsqcurvefit(fun,x0,xdata,ydata,lb,ub)
        times = linspace(xdata(1),xdata(end));
        plot(xdata,ydata,'ko',times,fun(x,times),'b-')
        txt = {'a-b*exp(-c*t)', ['a= ' num2str(x(1))], ['b= ' num2str(x(2))], ['c= ' num2str(x(3))]};
        text(max(u_image_time(IoI1:127))*2/3, (max(median_intensity_N(IoI1:127)./median_intensity_pr(IoI1:127))-min(median_intensity_N(IoI1:127)./median_intensity_pr(IoI1:127)))/5+min(median_intensity_N(IoI1:127)./median_intensity_pr(IoI1:127)), txt)
        title(FigList{3})
 
        hold off

    nexttile
    hold on
        for i= 1:length(ArrayN) %i=Nucleus number
            if height(ArrayN{i}) >= 0.75*max(table_heightsN)
            plot(ArrayN{i,1}.Timestamp,ArrayN{i,1}.Intensity_MedianIntensity_Alexa488./ArrayPR{i,1}.Intensity_MedianIntensity_Alexa488,'color',[0,0,0]+0.75)
            else
            end
        end

    plot(u_image_time,(median_intensity_N./median_intensity_pr),'k','LineWidth',2); %Black median line
    title(FigList{4})
    hold off

    a = string(LineNumb);
    tite = strcat(FigTitle, '   N=', a);

    title(Tile,tite, 'fontsize', 16)
    xlabel(Tile,'Time Scale (sec)', 'fontsize', 16)
    ylabel(Tile,'Median Fluorescence Intensity', 'fontsize', 16)
    Tile.TileSpacing = 'compact';
    Tile.Padding = 'normal';
    set(findall(Tile, '-property', 'fontsize'), 'fontsize', 16)

    picName = char(strcat(filepath, '/', folderList{p}, '/NormalizedTrendFit', newstr, '_', datestr(now, 'yyyymmddHHMM'), '.png'))
    saveas(gcf, picName);

end

[dpb]

That's a lot of code; not that it's all that complicated, just a lot for somebody to look through carefully enough to find issues...

My first suspicion is the very long file path and the reference therein to "cloud" -- have you tried moving a group of files to the local file system and seeing if it still hangs? Or, even first just making the write to local file.

[Nicholas Scott]

I can confirm the freeze is still happening locally. Since it is sporadic, this could be a red herring, but I tried without the saveas function, it worked without a hitch, then when using the saveas it froze on the first iteration. Maybe the saveas is slowing down the process, causing an issue? Kind of like the hands are wanting to move faster than the brain?

I even thought, 'well what if it is because the saveas function is still saving to the cloud, and the uploading is the hookup?' I changed the save and read path simultaneously to local, and it still hangs up and freezes. The freeze up screen shows all 4 graphs fully formed, but right before the tiledlayout goes to full screen before saving (aka still stuck in the default widow sizing), it freezes.

Thoughts?

[dpb]

Those details would tend to make one think along the line @Matt J suggests below...

Which renderer are you using; hardware I suppose? If so, just for grins, see if setting to 'software' makes any difference (at the expense, probably, of a visible slowdown, but if it were something in the graphics renderer, it just might make a difference).

[Nicholas Scott]

how does one change the graphics renderer? I can definitely check this out!

[dpb]

See what

rendererinfo

returns to see what it says.

Then

opengl software|hardware

swaps back/forth. See the doc for all the gory details, of course...

Answer 1

[Matt J]

Does the freeze always happen right before the saveas? If so, it might be worth issuing a drawnow() command right before the save step.

    picName = char(strcat(filepath, '/', folderList{p}, '/NormalizedTrendFit', newstr, '_', datestr(now, 'yyyymmddHHMM'), '.png'))
    
    drawnow
    saveas(gcf, picName);

My suspicion is that the figure cannot be saved because the drawing of the figure is not fully complete, resulting in some sort of tug of war between saveas and the renderer.

[Nicholas Scott]

it always occurs right before the saveas and the figure is not at full size. I'll add that and see if it alleviates the issue. Thanks Matt J!