2010-02-10 93 views
4

我不太確定這是否可行,但我對MATLAB的理解肯定會更好。在MATLAB中向量化循環

我有一些代碼我希望向量化,因爲它在我的程序中造成了相當多的瓶頸。它是優化程序的一部分,它具有許多可能的短期平均(STA),長期平均(LTA)和靈敏度(OnSense)配置。

時間是以矢量形式,FL2onSS是主數據(NX1雙),FL2onSSSTA是其STA(NxSTA雙),FL2onSSThresh是其閾值(NxLTAxOnSense雙)

的想法是,以計算紅色報警矩陣將爲4D - 在整個程序的其餘部分使用的alarmStatexSTAxLTAxOnSense。

Red = zeros(length(FL2onSS), length(STA), length(LTA), length(OnSense), 'double'); 
for i=1:length(STA) 
    for j=1:length(LTA) 
     for k=1:length(OnSense) 
      Red(:,i,j,k) = calcRedAlarm(Time, FL2onSS, FL2onSSSTA(:,i), FL2onSSThresh(:,j,k)); 
     end 
    end 
end 

我目前得到這個重複,企圖獲得更多一點的速度了它的功能,但顯然這將是更好的,如果整個事情可以矢量化。換句話說,如果有更好的解決方案,我不需要保留這個功能。

function [Red] = calcRedAlarm(Time, FL2onSS, FL2onSSSTA, FL2onSSThresh) 

% Calculate Alarms 
% Alarm triggers when STA > Threshold 

zeroSize = length(FL2onSS); 

%Precompose 
Red = zeros(zeroSize, 1, 'double'); 

for i=2:zeroSize 
    %Because of time chunks being butted up against each other, alarms can 
    %go off when they shouldn't. To fix this, timeDiff has been 
    %calculated to check if the last date is different to the current by 5 
    %seconds. If it isn't, don't generate an alarm as there is either a 
    %validity or time gap. 
    timeDiff = etime(Time(i,:), Time(i-1,:)); 
    if FL2onSSSTA(i) > FL2onSSThresh(i) && FL2onSSThresh(i) ~= 0 && timeDiff == 5 
     %If Short Term Avg is > Threshold, Trigger 
     Red(i) = 1; 
    elseif FL2onSSSTA(i) < FL2onSSThresh(i) && FL2onSSThresh(i) ~= 0 && timeDiff == 5 
     %If Short Term Avg is < Threshold, Turn off 
     Red(i) = 0; 
    else 
     %Otherwise keep current state 
     Red(i) = Red(i-1); 
    end 
end 
end 

代碼很簡單,所以我不會再解釋它。如果您需要澄清特定行的功能,請告訴我。

回答

5

訣竅是帶上所有的數據,以同樣的形式,主要使用repmat置換。那麼邏輯就是簡單的部分。

我需要一個討厭的技巧來實現最後一部分(如果沒有任何條件成立,請使用最後的結果)。通常這種邏輯是使用cumsum完成的。我不得不使用另一個2.^n矩陣來確保定義的值被使用(所以+ 1,+ 1,-1會真正給出1,1,0) - 只要看代碼:)

%// define size variables for better readability 
N = length(Time); 
M = length(STA); 
O = length(LTA); 
P = length(OnSense); 

%// transform the main data to same dimentions (3d matrices) 
%// note that I flatten FL2onSSThresh to be 2D first, to make things simpler. 
%// anyway you don't use the fact that its 3D except traversing it. 
FL2onSSThresh2 = reshape(FL2onSSThresh, [N, O*P]); 
FL2onSSThresh3 = repmat(FL2onSSThresh2, [1, 1, M]); 
FL2onSSSTA3 = permute(repmat(FL2onSSSTA, [1, 1, O*P]), [1, 3, 2]); 
timeDiff = diff(datenum(Time))*24*60*60; 
timeDiff3 = repmat(timeDiff, [1, O*P, M]); 
%// we also remove the 1st plain from each of the matrices (the vector equiv of running i=2:zeroSize 
FL2onSSThresh3 = FL2onSSThresh3(2:end, :, :); 
FL2onSSSTA3 = FL2onSSSTA3(2:end, :, :); 

Red3 = zeros(N-1, O*P, M, 'double'); 

%// now the logic in vector form 
%// note the chage of && (logical operator) to & (binary operator) 
Red3((FL2onSSSTA3 > FL2onSSThresh3) & (FL2onSSThresh3 ~= 0) & (timeDiff3 == 5)) = 1; 
Red3((FL2onSSSTA3 < FL2onSSThresh3) & (FL2onSSThresh3 ~= 0) & (timeDiff3 == 5)) = -1; 
%// now you have a matrix with +1 where alarm should start, and -1 where it should end. 

%// add the 0s at the begining 
Red3 = [zeros(1, O*P, M); Red3]; 

%// reshape back to the same shape 
Red2 = reshape(Red3, [N, O, P, M]); 
Red2 = permute(Red2, [1, 4, 2, 3]); 

%// and now some nasty trick to convert the start/end data to 1 where alarm is on, and 0 where it is off. 
Weights = 2.^repmat((1:N)', [1, M, O, P]); %// ' damn SO syntax highlighting. learn MATLAB already! 
Red = (sign(cumsum(Weights.*Red2))+1)==2; 

%// and we are done. 
%// print sum(Red(:)!=OldRed(:)), where OldRed is Red calculated in non vector form to test this. 
+0

太棒了!謝謝Ofri。 這肯定會讓我有點消化,我需要澄清幾點: 我假設行 FL2onSSThresh3 =重塑(FL2onSSThresh,[N,O * P]); 實際上應該與FL2onSSThresh2關聯? timeDiff = etime(Time(i,:),Time(i-1,:)); 在這種情況下是不是會工作?因爲我們不再使用我來迭代。時間是一個datevec,所以我假設我可以使用timeDiff = diff(Time)代替上面的代碼行,但是這會給出一個尺寸大不相同的數組,並且Red3操作失敗。 – Geodesic

+0

你說得對。我會編輯修復它。 –