Converting for loops to parfor format gives error "Error: Unable to classify the variable 'f_est_1' in the body of the parfor-loop."

Hello, I am trying to convert the for-loop into parfor loop for faster computation (because of huge amount of data).

The code and data (sample data) is attached for the reference. The code works fine when compiled with for-loop (attached code is with for-loop).

When I change the for-loop to parfor I get the error "Error: Unable to classify the variable 'f_est_1' in the body of the parfor-loop".

The attached code segments the signal based on the adaptive thresholding. Then plots the segmented data and performs the fft of it and stores the dominant frequency from each fft plots in a vector. These tasks are done using the "for-loop". I wanted to convert the for-loops to parfor, I tried changing the the way variables to be used as per the parfor variable classification but I am not bale to crack it.

how to perform this task using parfor? any help is highly appreciable.

[Rajendra Rajak]

The following is the code:-

clear all;
close all;
clearvars;
tic
%% Read the data
load data2mw
figure(); 
plot(data(:,1), data(:,2)-mean(data(:,2))); 
xlabel('Time (s)','fontsize', 10); ylabel('Signal magnitude (V)', 'fontsize', 10); 
title('Raw Data signal', 'fontsize', 10);set(gca,'FontSize',10);
%% Extract the data for analysis
time=data(:,1);
voltage=data(:,2);
t1=2.1755
t2=2.1770
i1=find(time(:,1)>t1,1)
i2=find(time(:,1)>t2,1)
Time=time(i1:i2);
Volburst=voltage(i1:i2);
figure()
plot (Time,Volburst-mean(Volburst));xlabel('Time (seconds)');ylabel('Voltage (V)');

%% Denoise the signal 
Denoised_Burst_signal = wdenoise(Volburst,12, ...
    'Wavelet', 'sym8', ...
    'DenoisingMethod', 'UniversalThreshold', ...
    'ThresholdRule', 'Soft', ...
    'NoiseEstimate', 'LevelIndependent');
figure()
plot(Time,Denoised_Burst_signal)
%% calculate the analytical signal and get the envelope  %%
analytic = hilbert(Denoised_Burst_signal);
env = abs(analytic); 
figure(10);plot(env)
k = 20;   
%% moving average of the analytical signal
Smooth_window = 40;
env = conv(env,ones(1,Smooth_window)/Smooth_window);
env = env(:) - mean(env);
env = env/max(env); 
figure(11)
plot(env)

%% Threshold and initialize the buffers 
THR_SIG = 4*mean(env);
nois = mean(env)*(1/3);
threshold = mean(env);
thres_buf  = zeros(1,length(env)-k);
nois_buf = zeros(1,length(env)-k);
THR_buf = zeros(1,length(env));
h=1;
LogicalSignal =zeros(1,length(env));

for i = 1:length(env)-k
         %------ update threshold 10% of the maximum peaks found ------------%
  if env(i:i+k) > THR_SIG 
      LogicalSignal(i) = max(env);
      threshold = 0.1*mean(env(i:i+k)); 
      h=h+1;
  else
      % ----------- update noise --------------%
      if mean(env(i:i+k)) < THR_SIG
          nois = mean(env(i:i+k)); 
      else
          if ~isempty(nois_buf)
              nois = mean(nois_buf);
          end
      end
  end 
  
  thres_buf(i) = threshold;
  nois_buf(i) = nois;
  % ------------------------- update threshold -------------------------- %
  if h > 1
     THR_SIG = nois + 0.50*(abs(threshold - nois));            
  end
  THR_buf(i) = THR_SIG;
end 

     figure,ax(1)=subplot(211);
     plot(Volburst),hold on,plot(LogicalSignal.*max(Volburst),'r','LineWidth',0.5),
     hold on,plot(THR_buf,'--g','LineWidth',0.5);
     title('Raw Signal and detected Onsets of activity');
     legend('Raw Signal','Detected Activity in Signal',...
         'Adaptive Treshold',...
    'orientation','horizontal');
     grid on;axis tight;
     ax(2)=subplot(212);plot(env);
     hold on,plot(THR_buf,'--g','LineWidth',0.5),
     hold on,plot(thres_buf,'--r','LineWidth',0.5),
     hold on,plot(nois_buf,'--k','LineWidth',0.5),
     title('Smoothed Envelope of the signal(Hilbert Transform)');
     legend('Smoothed Envelope of the signal(Hilbert Transform)',...
         'Adaptive Treshold',...
    'Activity level','Noise Level','orientation','horizontal');
     linkaxes(ax,'x');
     zoom on;
     axis tight;
     grid on;

%% Calcuation of burst parameters and velocity prediction
digitalCell = bwconncomp(LogicalSignal);
%digitalCell.num = cellfun(@length,digitalCell.PixelIdxList)
CellLength = cellfun(@length,digitalCell.PixelIdxList);
CellLength2keep = CellLength>150;
N=  sum(CellLength2keep);
digitalCellNew = digitalCell.PixelIdxList(CellLength2keep);
%%
x=1:1:N;
for i = 1:length(x)
    %f_est_1 = zeros(1,N);
    yy = Volburst(digitalCellNew{1,i})
    tt = Time(digitalCellNew{1,i})
    figure()
    plot(tt,yy)
    %format long
    % FFT of the each burst
            % FFT of the Burst Voltage
    L=length(tt);                     % Length of signal
    nfft=2^nextpow2(L); %calculates the nearest power of 2 w.r.t. l
    %N=nfft;
    Ts=tt(end)-tt(1);                       % Total Observation time
    Fs=L/Ts;                            % Sampling frequency (will be same as fps)
    % F_nyquist=Fs/2;                     % Nyquist frequency
    yy_in_PRIME = yy - mean(yy);
    y1 = fft(yy_in_PRIME,L);            % Discrete Fourier transform (DFT)
    y2 = 2*((abs(y1))/L);               % Amplitude of the DFT
    [v,k] = max(y2);                      % find maximum value (v) and it's position in the array (k)
    %  f_totalscale = (0:(N-1))* Fs/N;       % total frequency scale
    f_scale=(0:(L-1)/2)* (Fs/L);            % Half the total frequency range
    y3=y2(1:round((L/2)));                       % half the power spectrum
    figure()
    plot(f_scale(1:length(y3))',y3);
    f_est_1(i)= f_scale(k);                   % dominant frequency estimate for each burst
    hold on;                              % Prevent image from being blown away
    plot(f_est_1,v,'ro');
    axis('auto'),grid('on');
    xlim([0 2e+07]);
    title(''),xlabel('Frequency (Hz)'),ylabel('Voltage (V)');
%     
    
end
%%
ElaspedTime1 = toc

Answer 1

[Walter Roberson]

    f_est_1(i)= f_scale(k);                   % dominant frequency estimate for each burst

Inside parfor i, at any one time, a (possibly discontinuous) subset of f_est_1 has had valid values written to it. parfor does not execute in sequential order (deliberately), so parfor will not write to f_est_1(1) then f_est_1(2) and so on. It might happen to write to f_est_1(17) then f_est_1(11) then f_est_1(20) then f_est_1(3) and so on. All you know about the order is that by the end of the loop, all of the entries have been filled in.

    hold on;                              % Prevent image from being blown away
    plot(f_est_1,v,'ro');

But there you try to plot using all of f_est_1 including the parts that have not been filled in yet.

[Rajendra Rajak]

I will modify the code and update back any errors.

[Walter Roberson]

I suggest you reconsider using parfor. Consider using parfeval() instead.

Initialize a parpool. Use parfeval() to start some workers, passing in values for i and values for THR_SIG .

Now, loop waiting for a "future" to be ready (finished computation.) Extract whatever values are needed such as thres_buf and nois_buf, and THR_SIG . Compare the retrieved THR_SIG to the "current" value, and choose the "best" of them. Now provided there is a remaining iteration to do, parfeval() the next iteration, passing in an i value and the "best" THR_SIG so far.

This assumes that THR_SIG can be decided as being better or not. You should not just use the "newest" value, because it is possible that the most recently returned value is from an older worker. But I guess in some cases it might make sense to use the value associated with the worker with highest index that has returned so far.

[Rajendra Rajak]

Thank you Walter for the valuable response.

I will incorporate these suggestions as you mentioned and get back to you if any errors persists.

One possible solution I am currently thinking is that I should avoid updating the "THR_SIG" as this is not a good thing to do in parfor, basically avoiding the adaptive thresholding. I'll choose a global threshold for whole signal ( which is basic method), then perform the parfor to extract the bursts in the signal.

PS: This parallel computing is something quite new to me.

[Rajendra Rajak]

Hi Walter, I made some minor changes in the code (avoided plotting as of now) and it works fine. I have followed your advice of not updating the "THR_SIG" (completely avoided for other buffer variables as well which I was updating over each iteration), and also haven't tried parfeval as of now. I have to further update code to inlcuded more feature extraction from the signal. "for-loop" take 6.8875 seconds for the sample signal. "parfor" takes 2.4215 seconds for the same sample signal. This would considerably reduce the computing time as my actual file is having close to 400 M data points. Loading it in MATLAB is another issue. I guess I may have to ask seperately that " how to load big amount of data".

Thank you.

[Walter Roberson]

In some cases you can improve performance by having each worker use matfile https://www.mathworks.com/help/matlab/ref/matlab.io.matfile.html instead of reading all the data in at the beginning. However, if the data is on hard disk (instead of on SSD) then this approach can end up slowing down computation, as the workers can end up competing for control over the disk resources.

[Rajendra Rajak]

Thank you Walter for the response.

I have a .bin data file. I used bin to ASCII conversion .m script, it took close to 80 minutes to load the data. After this, the data is processed using parfor. I think I may have to live with this.

I'll certainly follow the approach of not loading whole data at the beginning.

I'll update you once accomplished this task.

[Rajendra Rajak]

Hi Walter,

%%
parfor i = 1:length(LL)
             %------ update threshold 10% of the maximum peaks found ------------%
  if env(i:i+k) > THR_SIG 
      LogicalSignal(i) = max(env);
  else
            
  end 
  
  thres_buf(i) = threshold;
  % ------------------------- update threshold -------------------------- %
  
  THR_buf(i) = THR_SIG;
end 

This part of the code consumes more time and also it gives warning " the entire array or struct 'env' is a broadcast variable. This might result in unnecessary communication overhead".

Is there a way to get rid of this? (for 12 Mega data points it took 85 minutes only for this part)

Answer 2

[Raymond Norris]

I'm gathering the for-loop you're trying to re-write is

for i = 1:length(x)

Assuming you're able to re-write this, have you considered what will happen to plotting inside the parfor? Keep in mind, the entire body of code in the parfor is running elsewhere (i.e. the workers). Therefore, you won't see plotting. You have two options for plotting "parallelization" (I put this in quotes since it will need to happen outside of the parfor, where you are running serial code).

• DataQueue: look at creating a data queue to send data from the workers to the client, which then calls a function to update your plot(s).
• parfeval: you'll still need a parallel pool, but rather then using the built-in parfor, create tasks with parfeval and then update the plot(s) when the tasks come back.

This doesn't answer your parfor question, but if you implement parfeval, (A) you ought to be able to plot your data and (B) it'll probably solve your original issue.

[Rajendra Rajak]

I will try to modify the code as you suggested and update you for any further errors.

[Rajendra Rajak]

Even after globalization of the "THR_SIG" outside parfor, it throws this error "Error: The temporary variable 'THR_SIG' must be set before it is used."

THR_SIG = 4*mean(env);
nois = mean(env)*(1/3);
threshold = mean(env);
thres_buf  = zeros(1,length(env)-k);
nois_buf = zeros(1,length(env)-k);
THR_buf = zeros(1,length(env));
h=1;
LogicalSignal =zeros(1,length(env));
parfor i = 1:length(env)-k
         %------ update threshold 10% of the maximum peaks found ------------%
  if env(i:i+k) > THR_SIG 
      LogicalSignal(i) = max(env);
      threshold = 0.1*mean(env(i:i+k)); 
      h=h+1;
  else
      % ----------- update noise --------------%
      if mean(env(i:i+k)) < THR_SIG
          nois = mean(env(i:i+k)); 
      else
          if ~isempty(nois_buf)
              nois = mean(nois_buf);
          end
      end
  end 
  
  thres_buf(i) = threshold;
  nois_buf(i) = nois;

[Walter Roberson]

Do you still have

  if h > 1
     THR_SIG = nois + 0.50*(abs(threshold - nois));            
  end

You cannot have that inside the parfor loop.

[Rajendra Rajak]

Thank you Walter Roberson for answering. Yes, As my code try to adopt adaptive filtering using parfor, i guess the error is kept on popping as you said.

Is there a solution for adaptive filtering using parfor? Although the Parallel toolbox matlab documentation suggests that not all for oop can be coverted to parfor.

[Walter Roberson]

I feel I need to emphasize even more strongly: parfor will never execute the iterations "in order" ... well except for the case you are only doing one iteration.

Imagine if instead of using

for i = 1:length(env)-k

you had instead used

for i = randperm(length(env)-k)

If the order of iterations matters to you, then you should not be using parfor.

It is not completely impossible to update a value on a worker, but most of the time it is the wrong thing to do.

To get an adaptive value:

• create the parallel pool explicitly using parpool()
• use parallel.pool.pollableDataQueue() on the controller before starting the parfor. https://www.mathworks.com/help/parallel-computing/parallel.pool.pollabledataqueue.html . Call this queue Q_master for the moment
• Then use parfevalOnAll() to run a function on each of the workers that creates a pool, parallel.pool.pollableDataQueue and sends the pool variable back by writing to Q_master from each worker. The controller should store all of the received queues
• Now parfor()

There are more steps but I am falling asleep now.

[Rajendra Rajak]

Thank you Walter for the descriptive reposne to solution of the problem and sharing the link.

I am going to restructure the code and get back to you in the same thread to list out any discrepancies MATALB throws during compilation.