CodyNote005：最大長度子序列問題

來自專欄 Cody習題6 人贊了文章

馬良，祁彬彬

序言

從1-D數組內提取滿足某種要求的動態長度子序列，算是有挑戰性的趣味代碼問題，首先明確挑戰性究竟在哪兒？所謂「某種要求」，通常和序列長度、兩序列內部元素兩相比較，內部各個元素排布特徵等，有密切關係，因此子序列提取往往具有明顯動態特徵，一般if判定流程描述複雜而冗長；其次，MATLAB中用邏輯索引、矩陣構造、函數命令組合等方式寫出合格的矢量化代碼，對代碼編寫人的函數運用理解能力也是個不小的考驗。

今天打算試試水，由其中最簡單的問題開始，初步探討如何從0-1一維數組中提取最長連續序列。關於動態長度子序列問題在後續內容中還會繼續結合實例做細緻討論。

子序列最大連續長度問題

Pro.42894：Find Longest Run

鏈接：https://ww2.mathworks.cn/matlabcentral/cody/problems/42894

【原題】：Write a function longest_run that takes as input an array of 0s and 1s and outputs the length and index of the longest consecutive run of either 0s or 1s. If there are multiple runs of the same length, it should output the first occurence.

Examples:

1. longest_run([1 1 0 1]) should output [2 1], since the longest consecutive run is [1 1], which starts at index 1

2. longest_run([1 1 0 1 0 0 0 0 1 1 1]) should ouptut [4 5], since the longest consecutive run is [0 0 0 0], which starts at index 5

3. longest_run([1 0]) should ouptut [1 1], since the first longest consecutive run is 1, which starts at index 1

【解釋】：給出一個0-1數組，即數組中所有元素全部由0和1組成，找到其中最長連續相同數字的起始位置和連續的長度值。例如對序列[1 1 0 1 0 0 0 0 1 1 1]而言，最長的應該是從第5位開始的連續4個0，所以結果是[4 5]。

Test

Code Input

1

%%binary_array = [1 1 0 1];expected_length = 2;expected_index = 1;[actual_length actual_index] = longest_run(binary_array);assert(isequal(expected_length, actual_length));assert(isequal(expected_index, actual_index));

2

%%binary_array = [1 1 0 1 0 0 0 0 1 1 1];expected_length = 4;expected_index = 5;[actual_length actual_index] = longest_run(binary_array);assert(isequal(expected_length, actual_length));assert(isequal(expected_index, actual_index));

3

%%binary_array = [1 0];expected_length = 1;expected_index = 1;[actual_length actual_index] = longest_run(binary_array);assert(isequal(expected_length, actual_length));assert(isequal(expected_index, actual_index));

4

%%binary_array = [1];expected_length = 1;expected_index = 1;[actual_length actual_index] = longest_run(binary_array);assert(isequal(expected_length, actual_length));assert(isequal(expected_index, actual_index));

5

%%binary_array = [1 1 1 1 1 1];expected_length = 6;expected_index = 1;[actual_length actual_index] = longest_run(binary_array);assert(isequal(expected_length, actual_length));assert(isequal(expected_index, actual_index));

【分析】：比較容易想到的思路是從數組內某元素開始，向右逐個循環查找相鄰元素是否與本身相同，不過難點也在這裡，因為這段最長的連續子序列，在數組中的起始位置未知，長度未知。

如果按照前述分析，從最左端開始，循環逐步移位，和前一元素相等則長度加1，把所有連續部分都提取出來，同時記錄起始位置，循環結束後，比較所找到最長的那段序列。這樣求解，可以想像程序會比較繁瑣；另一思路是把數組換成字元串進行文本搜索，是的，沒錯，典型的數組問題，也可以通過字元串命令求解，在CodyNote004中，也提到二者可以藉由Ascii碼切換。當輸入數組序列變成文本，可以使用各種文本查找命令；除此之外，卷積或濾波函數，結合分組函數splitapply等，也能得到簡約卻絕不簡單的矢量化求解思路。

劃線平推的循環移位

先說說多數人第一意識容易想到的劃線平推大法：循環右移，照著相等條件往長度上分段堆積木。

function [m, c] = longest_run(a)f=0;F=[];l=length(a);for i=1:l-1 if a(i)==a(i+1) f=f+1; else f=0; end F=[F f];end[k, p]=max(F);m=k+1;c=p-k+1;if isequal(2,l)c=c-1;elseif isequal(1,l)m=1; c=1;endend

能用for+if表達出這個動態序列的概念，其實不易，因為當代碼編寫人矢量化代碼思路不明確，對函數理解不透徹，硬想套用命令組合避免循環，會寫出東施效顰的恢弘視效：

function [run_length, start_index] = longest_run(binary_array)B_1=[0 binary_array 0]C_1=find(diff(B_1)==1)run_length_1=max(find(diff(B_1)==-1)-find(diff(B_1)==1))start_index_1=C_1(find(run_length_1==find(diff(B_1)==-1)-find(diff(B_1)==1))) B_0=[0 binary_array==0 0]C_0=find(diff(B_0)==1)run_length_0=max(find(diff(B_0)==-1)-find(diff(B_0)==1))start_index_0=C_0(find(run_length_0==find(diff(B_0)==-1)-find(diff(B_0)==1)))if run_length_0>run_length_1run_length=run_length_0start_index=start_index_0elserun_length=run_length_1start_index=start_index_1endend

這段程序的問題主要是沒處理好元素查找命令的設計安排：find被重複太多，多到影響對整體思路的理解。也有移位程序寫得相對比較成功的：

function [run_length, start_index] = longest_run(binary_array)varlength = 1;run_length = 1;start_index = 1; for I = 2:length(binary_array) if (binary_array(1,I-1) == binary_array(1,I)) varlength = varlength + 1 if varlength > run_length run_length = varlength start_index = I - varlength +1 end else varlength = 1 end endend

從第2個元素起開始判斷，如果相同，長度varlength就加1，移動中遇到不同元素，則停止本次計數累加，varlength重新賦值為1，重新開始前述過程。之所以說寫得不錯，因為沒太多糾結判斷究竟是等於0還是等於1上，因為代碼最終只要最長序列，沒說是全0子序列還是全1子序列。

字元類型解法1：num2str+strfind

我本想提交正則版本，想了想估計用num2str+strfind做這題的人不會多，於是寫了下面的代碼：

function [r,s] = longest_run(x)t=num2str(x,-6);s=[];t2=[];r=7;while isempty(s)&&isempty(t2) r=r-1; s=strfind(t,repelem(0,r)); t2=strfind(t,repelem(1,r));endif ~isempty(t2) s=t2;endend

首先num2str後置參數「-6」，這是沒有在MATLAB幫助中正式提及的載入參數，從下面例子里看看區別：

>> num2str(ones(1,4))ans = 1 1 1 1>> +ansans = 49 32 32 49 32 32 49 32 32 49>> num2str(ones(1,4),-6)ans = 1111>> +ansans = 49 49 49 49

如不加後綴參數「-6」，轉化為字元串後，多了許多空格，運行結果顯示ASCII碼值為32；有參數-6，則自動去掉這些字元元素間的空格。本問題中，無疑帶參數更合適，關於字元串內去空格辦法較多，比如在前者基礎上，用ans(ans==』 』)=』』來去掉元素間的空格，更多轉換類型時去掉空格的方法在下面將繼續介紹。

言歸正傳，生成類似下面無空格字元串類型數組時，可對接各種字元串查找方法尋找最長序列：

>> x=randi([0 1],1,8)x = 0 0 1 1 1 0 1 1>> num2str(x,-6)ans = 00111011

我用的是strfind，字元串查找符合pattern子串位置的利器之一。原因是命令如：strcmp、contains等等，只能返回是否存在子字元串的真假判斷，不返回索引位置。

字元類型解法2：sprintf+regexp

正則表達式用得順了，看到字元串總習慣性朝這條思路上拐，但與其強大而靈活的搜索、匹配替代功能相匹配的是regexp和regexprep超級複雜的各類底層語法規則、輸入輸出參數。不過MATLAB幫助寫得很細緻，示例和語法之間一一對應，但除語言障礙外，缺乏有效的實例訓練(自己不找其他問題反覆練，幫助例子再典型也白搭)是多數人不能很好掌握正則文本處理的主要困難之一。因為層次豐富的正則語法規則及諸多後綴參數、輸入輸出的多變設定，沒下一定功夫，精通不過是種妄想。我和彬彬的書中，專門用一章通過大量實例描述這兩個命令的用法，今後其他筆記中，這兩個命令也無疑會是常客。

下面就是把字元串的正則搜索命令用在數組問題中的案例：

function [run_length, start_index] = longest_run(binary_array) str = sprintf(%d, binary_array); [s, idx] = regexp(str, (1+)|(0+), match, start); [run_length,id] = max(cellfun(@length,s)); start_index = idx(id);end

首先討論上述代碼第1行，sprintf函數以不帶空格的方式，將數組轉字元串的方法。這是前一節num2str的「-6」參數講解的延續。sprintf在構造動態字元串時使用頻繁，工程實驗或者數據讀寫，在構造輸出文本文件時，其中包含計算、模擬、試驗數據和文字的混合，用文本說明數據的用途、單位等等。sprintf可以在固定的文本字元串中，動態替換數據，另外，有時也用於動態構造文件名、變數名等，或在正則表達式中構造動態搜索或替換表達式等。

>> sprintf(the price of this book is $%6.2f,randi([4 8]))ans = the price of this book is $ 6.00>> sprintf(the price of this book is $%6.2f,randi([4 8]))ans = the price of this book is $ 8.00>> sprintf(the price of this book is $%6.2f,randi([4 8]))ans =the price of this book is $ 4.00

顯然，字元串中數字在隨機數的不同調用時發生改變，同時，用』%6.2f』約定了其默認保留小數點後兩位有效數字。本問題中生成無空格字元型數字序列其實只能算sprintf小試牛刀：

>> sprintf(%d,x)ans = 00111011

接下來就是把它放進正則搜索命令「regexp」：

注意正則搜索表達式中的邏輯「或」操作，這意味著優先搜索前一個表達式』(1+)』，括弧中的內容代表搜索連續由字元』1』組成的字串，如果搜索到就不執行後面第2個』0+』，也就是連續1個以上的字元』0』。輸出：』[s,idx]』不但需要得到所匹配的多個連續1或者0，主要為後面統計連續數量，idx用於判定起始索引位，用regexp搜索到一系列不定數量字元串，並放進cell數組，繼續用cellfun統計數量，套一層max獲最大值，這個解法是值得提倡的，因為把搜索任務交給正則表達式，而不是自己費半天勁自定義循環套路，cell數組搜索輸出結果也沒浪費，轉入cellfun調用「@length」句柄，自動遍歷所有前述連續長度cell數組，得到各自長度大小。

這是個常用的MATLAB矢量化編程套路。

大神Alfonso Nieto-Castanon也用了類似方法，但手法更加不落俗套：

function [run_length, start_index] = longest_run(x)[a,b]=regexp(char(a+x),a+|b+,start,end);[run_length,idx]=max(b-a+1);start_index=a(idx);end

第一句調用regexp命令，選擇兩個不很常用的參數重載：』start』,』end』，分別獲得搜索字元串起始和結束位置；正則搜索表達式中匹配一段連續的』a+』，或者連續』b+』(+符號代表至少搜索到1個以上的a(b)才算匹配成功)，估計有人奇怪為什麼問題中的0-1數組突然變成搜索連續字元串』a』，和』b』了？可能他感覺字元形式的』0』或者』1』看起來不夠「文本化」吧，用char把01變成了』ab』，注意代碼char(『a』+x)的輸入中包含著兩個截然不同的數據類型：char形式的』a』和數值類型的0-1數組，也就是輸入x，基於Ascii碼的兩者相加，再用char換成char字元類型，這種用char天然去空格的「數值-文本」轉換法，我記得LY Cao、劉鵬都常用，彬彬有次跟我說某個題目就偷學了這招，弄過一個leading：

>> d=randi([0 1],1,10)d = 1 1 0 1 1 0 0 1 1 1>> a+dans = 98 98 97 98 98 97 97 98 98 98>> char(ans)ans = bbabbaabbb

順勢接外層regexp正則搜索，返回的a和b變數依次為搜索到的多個a或者多個b的始末索引端點。當然，正則搜索搜ab也是搜，搜01也是搜，所以如下代碼是等效的：

[a,b]=regexp(num2str(x,-6),0+|1+,start,end)

或者：

[a,b]=regexp(sprintf(%d,x),0+|1+,start,end)

總體來講，Alfonso用基本Ascii碼值把數組轉成文本，接入regexp正則搜索，銜接起來行雲流水，細微處體現了Alfonso在使用MATLAB處理具體問題時的精湛功力。

但這還不是regexp篇章終結：如果不知道start和end這兩個重載參數，還能否通過regexp命令獲得比較簡潔的代碼？

J.S.Kowotan給出了答案：

function [run_length,ans] = longest_run(x) [regexp(num2str(x),(1+|0+)) numel(x)+1]; [run_length,ind] = max(diff(ans)); ans(ind);end

上述程序第1句中，先構造貌似突兀的「正則搜索索引+長度」數組，Test.2執行結果：

x = [1 1 0 1 0 0 0 0 1 1 1];[regexp(num2str(x),(1+|0+)) numel(x)+1]ans = 1 3 4 5 9 12

ans前5項代表連續1或者連續0的搜索起始索引位，最後一項是輸入數組長度+1.這是為獲得最後一段連續相同數字的長度[x(9)-x(end)]，為下一行命令diff相鄰做差準備，Kowontan讓regexp純粹成為搜索索引的工具，處理索引長度的工作則交給了max+diff的數值型命令組合。

分析到這裡，就該進入使用純粹數值命令組合完成工作的時候了，畢竟問題本身是數值序列，要說沒有純數值解決方案，那是說不過去的。

數值解法1：cumsum+mode

cumsum對數組元素正反向累加求和，例如：

>> cumsum(1:10)ans = 1 3 6 10 15 21 28 36 45 55>> cumsum(1:10,r)ans = 55 54 52 49 45 40 34 27 19 10

mode用於統計數組中出現最頻繁的數，以及出現的頻次。2個命令都是2006a版本出現的。如下是個堪稱漂亮的組合解決方案：

function [r, s] = longest_run(x)c = cumsum([2,~x,2,x]);[m,r]=mode(c);r = r-1;s = rem(find(c==m,1),length(x)+1);end

最難想到的，是利用cumsum構造備統計數組c的第1句，邏輯值~x和原數組x是為把最長連續的0和最長連續的1統計出來，原因很簡單，cumsum累加時，連續0的位置處，累加值不發生變化，而連續1則梯次加1，不利於mode統計，故用邏輯反來翻轉數組，乾脆統計兩個連續0的長度，此外，由於cumsum的累加，使得每段連續長度的數值都是獨一無二的，便於後續一句mode的元素數量統計；同時這兩個2和最後一句用rem求余命令相互呼應，動機實際上是為求余而試湊。

第2句mode返回兩個輸出，m是數組c中出現最頻繁的數，r則是m在數組c中的頻數，減1得輸入變數x中最長連續1或者0的長度。

代碼短，思路卻不易想得到，綜合求余試湊、構造0-1數組統一處理、頻數統計等多個技巧。

數值解法2：filter2+splitapply

如下卷積+分組構造算是比之前mode+cumsum組合更難想到的方案。

function [len, idx] = longest_run(x) counts = splitapply(@numel, x, cumsum(filter2(eye(1, 3), x)~=x)); [len, idx] = max(repelem(counts, counts));end

利用二維濾波filter2+splitapply+cumsum+repelem的複雜做法，來自大神Jan Orwat的手筆。

從內到外，逐層分析

第1步是二維濾波函數執行結果，窗函數選擇[1 0 0]，結果：

>> filter2(eye(1, 3), x)ans = 0 1 1 0 1 0 0 0 0 1 1

二維濾波命令和卷積有著千絲萬縷瓜葛，二者數學上同源，只是窗函數位置和構造有變動，一維數組用filter2，第1參數eye(1,3)，即：[1 0 0]，如果用conv2，第2，注意，是第2參數，就得用rot90(eye(1,3),2)，也就是旋轉180°代替，且放在第2參數位上，加』same』參數保證維數不變。

>> conv2(x,[0 0 1],same)ans = 0 1 1 0 1 0 0 0 0 1 1

卷積命令conv結果相同：

>> conv(x,[0 0 1],same)ans = 0 1 1 0 1 0 0 0 0 1 1

對卷積計算基本原理如有興趣，參考我和彬彬寫的書第1.6.3節，或任一本和數字信號處理有關聯的教材。下面繼續討論正題：後接和原數組x的不等比較，得到：

>> conv(x,[0 0 1],same)~=xans = 1×11 logical 數組 1 0 1 1 1 0 0 0 1 0 0

估計絕大多數人看不出什麼跡象——就像我剛開始一樣，那乾脆把x和這個邏輯數組放一塊兒對比一下：

>> [x;conv(x,[0 0 1],same)~=x]ans = 1 1 0 1 0 0 0 0 1 1 1 1 0 1 1 1 0 0 0 1 0 0

上邊是x，下邊是卷積+邏輯操作的返回結果：所有原x數組中連續n個相同數構成的子數組，下邊一行的對應操作結果都變成了諸如：[1 0 0 … 0]這種形式，使用cumsum時，這一段累加值就不會發生變化；此外，原數組在相鄰兩段連續數字拼接處，對應位置不是由0變1(前一段連續數值[1 0 0 … 0]的末尾0)，就是數字1不發生變化，這意味著當使用累積求和命令cumsum時，累加值將向上增加1，也就是統計下一段連續數，方便比較起見，再次把分組結果和原數組x放在一起：

>> [x;cumsum(conv(x,[0 0 1],same)~=x)]ans = 1 1 0 1 0 0 0 0 1 1 1 1 1 2 3 4 4 4 4 5 5 5

第2行的累積求和結果實際上為最外層的splitapply提供了統計分組的依據，剩下的就是讓splitapply去按照指定分組，調用「@numel」或者「@length」這類函數，統計各個分組的長度，最後一句也很有意思：

>> counts=splitapply(@numel, x, cumsum(filter2(eye(1, 3), x)~=x))counts = 2 1 1 4 3>> repelem(counts, counts)ans = 2 2 1 1 4 4 4 4 3 3 3>> [t,it]=max(ans)t = 4it = 5

再度把長度向量擴展成和原輸入數組同維向量，用max命令找到最大值及其索引位，存在多個相等max值時，max函數默認返回第1個索引位(按從左到右)，圓滿解決問題。

整個看下來，問題難點在濾波函數filter2的eye(1,3)構造，感覺幾乎神來之筆！

感嘆：人和人智商方面的差距，有時真TM比人和狗還大！

數值解法3：max+diff+find

diff命令人人知道，真懂怎麼用的實際鳳毛麟角：

function [len, idx] = longest_run(x) diff(find([1 diff(x) 1])); [len, idx] = max(repelem(x, x));end

這個看似簡單而不簡單的動態數組序列，其實只需2句代碼，真正核心的部分，是第1句的連續兩次diff。

第1句內圈diff做差，返回結果中不為0的地方，都是數字發生變化的位置，但不能光diff，因為做1次diff維度就減1,2次diff意味著必須做2次補位，所以頭尾各加1，但這裡存在一個問題：輸入x首尾元素數值是0還是1並不確定，首尾都加1會不會影響結果的變化？看兩個例子：

>> x=randi([0 1],1,10)x = 0 1 1 0 1 0 0 1 1 1>> find([1 diff(x) 1])ans = 1 2 4 5 6 8 11>> x=randi([0 1],1,10)x = 1 0 1 1 1 1 1 0 1 0>> find([1 diff(x) 1])ans = 1 2 3 8 9 10 11

顯然，首尾元素對結果是沒有影響的，因為結果要的是0-1的開關變化條件，因此find([1 diff(x) 1])找到的不為0的索引位中，永遠存在1和numel(x)+1。最外層再加diff，得到所有連續數字長度值，剩下的就是用max+repelem組合尋求起始位置和長度值了。

這個代碼方案，我感覺是解決整數動態序列數組最大連續長度問題的leading solution。

Pro.672.Longest run of consecutive numbers

前一問題中的序列元素只有0和1兩種，下面再擴展一下，在全部由整數元素構成的數組中，最長連續子序列的元素值。

【鏈接】：https://ww2.mathworks.cn/matlabcentral/cody/problems/672

【原題】：Given a vector a, find the number(s) that is/are repeated consecutively most often. For example, if you have

a = [1 2 2 2 1 3 2 1 4 5 1]

The answer would be 2, because it shows up three consecutive times.

If your vector is a row vector, your output should be a row vector. If your input is a column vector, your input should be a column vector. You can assume there are no Inf or NaN in the input. Super (albeit non-scored) bonus points if you get a solution that works with these, though.

【解釋】：這次不找索引位置，也不找最大長度值，專找最大連續子序列中的元素值本身。

Test

Code Input

1

%%a = [1 2 2 2 1 3 2 1 4 5 1];y_correct = 2;assert(isequal(longrun(a),y_correct))

2

%%a = [1 1 1 2 2 2 3 3 3];y_correct = [1 2 3];assert(isequal(longrun(a),y_correct))

3

%%a = [-2 -2 -2 -2 -1 0 3];y_correct = -2;assert(isequal(longrun(a),y_correct))

4

%%a=[0 1 1 1 0 2 2 0 1 1 1 0];y_correct = [1 1];assert(isequal(longrun(a),y_correct))

5

%%a=[3 3 3 2 2 1 6];y_correct=3;assert(isequal(longrun(a),y_correct))

6

%%a=[3 3 3 2 2 2 1 6];y_correct=[3 2];assert(isequal(longrun(a),y_correct))

7

%%a=[1 2 3 4 5];y_correct=a;assert(isequal(longrun(a),y_correct))

【分析】：看起來好像簡單了，其實是錯覺：因為長度相等的最大子序列若不止1個，各自的元素值即使是相等，都要依次輸出，例如下面的輸入變數：

>> repelem(1:3,3)ans = 1 1 1 2 2 2 3 3 3

最大序列長度都是3，因此返回[1 2 3]，再如：

>> [3 3 0 2 2]ans = 3 3 0 2 2

最大序列長度2，返回結果是[3 2]。此外返回結果的順序和自左至右的出現順序相同。

問題要求返回最大連續子序列的元素值本身，需滿足2個要求：

1. 返回元素值t必須同時包含於輸入序列a和其分屬的子序列；
2. 返回元素值t所在數組中最長子序列m滿足unique(m)=t。

事實上Pro.672突破口就在前一問題各類思路里，它同樣能利用find+max+diff組合尋找連續相同元素。根據上述分析，尋找最長連續子序列長度的問題可以直接套用max+find+diff命令：

function ans=longrun(a)find(diff([nan,a(:),nan]));a(ans(find(diff(ans)==max(diff(ans)))));end

Alfonso這2行代碼中，第1行用於標識原數組中元素改變的索引位:

>> a = [1 2 2 2 1 3 2 1 4 5 1]a = 1 2 2 2 1 3 2 1 4 5 1>> find(diff([nan,a(:),nan]))ans = 1 2 5 6 7 8 9 10 11 12

代碼返回結果中需要注意2點：由於暫時只用1次diff，但首尾均添加了nan元素，所以索引數組ans最大值(12)比原序列長度大1；其次，從ans看出，如果再來一次diff得到的就是各個子序列的長度值，顯然第2和第3個元素做差的結果就是最大子序列長度5-2=3。那麼最大值自然需要再用一次find命令，以diff構造輸入邏輯數組返回所有，注意是所有和最大長度相等的元素索引位，因為又用了一次diff，此時的邏輯數組diff(ans)長度和原輸入序列a相同，所找到的索引位置正好是a序列中，最長子序列的第1位，用a(...)檢索到該元素返回即可。

同類思路做法大同小異，而且首尾元素，Alfonso用了nan，Tim用元素1做端點：

function a=longrun(a)k=find([1 diff(a(:)) 1]);a=a(k(diff(k)==max(diff(k))));end

我用的是inf添加首尾元素：

function ans=longrun(a)find(diff([inf a(:) inf]));a(ans(find(diff(ans)==max(diff(ans)))));end

無論選擇何種首尾添加數值，返回結果都一樣。當然，在MATLAB靈活的語法規則下，第2個find可以省略：

function ans=longrun(a)find([1 diff(a(:)) 1]);a(ans(diff(ans)==max(diff(ans))));end

劉鵬和LY Cao發現如果把「[1 diff(a(:)) 1]」寫成「vertcat(1,diff(a(:)),1)」可以進一步減小size，看似沒什麼區別，但後者結構中的「1」、「diff(a(:))」、「1」是vercat的3個調用參數，而前者則相當於3個並列的函數，於是：

function val=longrun(a)find(vertcat(1,diff(a(:)),1));val = a(ans(diff(ans)==max(diff(ans))));end

的確是小細節，但體現了頂級高手對所寫代碼的完美程度要求，達到近乎於殘忍和苛刻的態度。

這個問題用regexp做正則搜索，就比較麻煩，例如端點索引搜索方案：

function val=longrun(a) [vals,vale]=regexp(sprintf(%d,logical(diff(a))),0*,start,end); val = a(vals(vale - vals == max(vale - vals))); if isempty(val) val = a; endend

regexp後置參數start和end實際就是索引查找，且增加字元→數值類型的中間步驟，還要做非空判定，因此肯定繁瑣一點，此外還有結合regexp的token存儲所需索引的正則搜索方法：

function ans=longrun(a)regexp(num2str(a(:)),(?<n>-?d)(?<len>s*k<n>)*,names);l = cellfun(length,{ans.len})cellfun(@str2double,{ans(l==max(l)).n})if iscolumn(a), ans; endend

正則搜索表達式用小括弧表示為2部分，前一部分token名稱為n，後一部分token名為len，前一部分是連續長度子序列的第1個元素，len則對應式子序列除第1個元素之外的其他部分。然後再通過接入「@length」，用cellfun判斷每個len的長度並找到最大值。創意很好，只是實在太麻煩了。

因此文本搜索方法效率感覺略吃虧。

Cody論壇的維護工程師Ned Gully給出用hankel矩陣構造基本結構的新奇做法：

function val = longrun(a) r = sum(cumprod(hankel([diff(a(:))==0;0]))); val = a(unique(find(r==max(r))))end

hankel矩陣每條副對角線上的元素都相等：

>> v=randi([0 3],1,5)v = 0 3 3 1 3>> hankel(v)ans = 0 3 3 1 3 3 3 1 3 0 3 1 3 0 0 1 3 0 0 0 3 0 0 0 0>> cumprod(ans)ans = 0 3 3 1 3 0 9 3 3 0 0 9 9 0 0 0 27 0 0 0 0 0 0 0 0

容易看出每行(列)視為原數組v循環從數組首位抽取i個元素，再在尾部補i個0元素獲得，重點是為什麼用hankel+cumprod+diff組合？這需要從裡層[diff(a(:))==0;0]說起：裡層代表輸入變數做差之後等於0的位置，意味著通過做差查找連續相等元素的索引，尾部加0之前說過，用一個diff要在端點補1位對齊，此時等於0的位置，因為邏輯值為真(==0)，全部變成1，因此hankel矩陣接受的輸入又是個0-1數組：

>> v=randi([0 3],1,8)v = 2 0 2 1 2 2 2 1>> v1=[diff(v)==0 0]v1 = 0 0 0 0 1 1 0 0

比較v和v1，後者起始元素1的位置和連續數組起始位置相同，由此hankel矩陣：

>> v2=hankel(v1)v2 = 0 0 0 0 1 1 0 0 0 0 0 1 1 0 0 0 0 0 1 1 0 0 0 0 0 1 1 0 0 0 0 0 1 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

因為是對稱矩陣，從行或者列方向看都一樣，不同的行(列)相當於讓連續1在矩陣中從原位置逐行(列)上(左)移，這麼做的目的是什麼呢？

>> v3=cumprod(v2)v3 = 0 0 0 0 1 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0>> v4=sum(v3)v4 = 0 0 0 0 2 1 0 0

通過cumprod命令執行的結果應該就很清楚了：如果在1的位置之前存在任何一個0元素，相乘的結果都是0，從原輸入v而言，意味著之前的這段序列並不是連續相同元素(別忘了邏輯操作「==0」，有0說明元素和v(v1==1))不相同，而一段連續相同元素(連續1)之後接任何0元素，同樣後續全部為0，連續的1用矩陣思想，做一次整體運算獲得，該演算法思想隱隱包含著動態長度子序列問題兩個本質性的關鍵詞:錯序和對位，所以我說這是個異常精巧的矢量化代碼構思！

延伸討論

如果轉回本篇內容第1題，發現題意要求得到的是最長的連續0或者連續1，如果把這個要求稍微修改一下：指定找到最長的連續1，怎麼改代碼？以下方一維0-1數組t為例：

>> t = [1 0 1 1 1 1 1 0 1 0]t = 1 0 1 1 1 1 1 0 1 0

如果是查找最長連續1的個數，照例前後補0，仍用diff，但是改成先查找非1元素(元素1的邊界)，再做差：

>> t1=diff(find([0 t 0]~=1))t1 = 2 6 2 1

返回結果t1意味著用前後端點0以及數組t內部的0，把各段元素1間隔成3段，連續非1的索引差為1，中間有n個1的，差值為n+1,因此減1並提取非0元素，即為所有連續1的長度值：

>> nonzeros(t1-1)ans = 1 5 1

上述過程再推廣，合寫成一行代碼，即序列x中所有連續y的個數:

f = @(x,y)nonzeros(diff(find([~y,x,~y]~=y))-1);>> f(t,1)ans = 1 5 1>> f(t,0)ans = 1 1 1>> t=randi([0 2],1,12)t = 0 1 0 2 0 0 0 0 1 0 2 1>> f(t,0)ans = 1 1 4 1>> t=randi([0 2],1,12)t = 0 2 2 1 0 0 1 1 0 1 2 0>> f(t,1)ans = 1 2 1

呃...

世界終於清凈了...

小結

(1). MATLAB中的數組和字元串命令間以Ascii碼為橋樑，二者間的轉換無縫銜接。通過恰當的類型轉換，正則搜索命令regexp、字元串命令char、strfind等，都可以在數組中發揮作用，順便提一句：strfind從名稱上看似是字元串搜索函數，其實輸入變數也接受數組：

>> strfind([1 1 0 0 1 1],[1 1])ans = 1 5

(2). 最大長度連續子序列問題中，包含著做差、索引端點查找、矩陣構造等基本內容。一些數組命令，見過未必代表會用，很多命令看似不起眼，只有深入分析問題背景，深挖後置參數的潛力，盡心安排多個命令的銜接搭配，才能真正透徹理解它們在思路上的不平凡，例如前述cumsum+mode，再如我和彬彬書中以及本篇中提到的max+diff+find；

(3). 盡量拓展和推演新老函數間的組合變化，例如splitapply、repelem、filter2、convn等，這裡面組合變化的維度可謂無窮無盡。問題解決在Cody從來不是目的，而是從中獲得什麼。例如本問題的解決方案中filter2中的索引序列構造思路，就是令人眼界大開的做法。

最後是CodyNote005主要函數一覽：