標籤:

CodyNote012:再談動態長度子序列(Part.1)

CodyNote012:再談動態長度子序列(Part.1)

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

馬良,祁彬彬


序言

接連3個與索引多次定址有關的排序問題,Note009介紹索引構造幾種基本方法;Note010-011藉由排序系列問題,探討如何利用MATLAB的內置函數實現索引序列的構造,介紹包含字元、數值類型數組元素位置序列的構造、引用、變換、擴展的控制技巧,並在最後的問題求解代碼中,介紹一個少見的多重索引定址經典案例。自本篇開始,繼續深入探討另一個與索引密切關聯的主題:動態長度子序列。

這類問題一般要求某個序列中,按需提取具有某種特徵的子序列,所謂「某種特徵」通常不是靜態的,所提取子序列要麼元素集合的數量不定;要麼數組元素並非特定的值,且子序列元素間或與外部相鄰元素間有動態關聯,要根據實際狀況判斷。求解此類問題的代碼通常有一定難度,需要熟練掌握命令組合、元素的索引控制等,擁有比較好的代碼嗅覺。本篇以及後續幾篇,將通過多個動態子序列提取或判斷的實例,綜合理解和掌握MATLAB一些相對比較複雜的命令組合技巧,實際上這部分內容中,一些東西已經有了濃濃的演算法味道,很多情況下,並不是單純能夠用什麼生僻的命令或調用方法就能徹底解決的。

遞增子序列問題

Problem 42837. Increasing sub-sequence (Level 1)

【鏈接】:ww2.mathworks.cn/matlab

【原題】:Given a vector, v, of real numbers, return a positive integer, n, representing the longest contiguous increasing sub-sequence contained in v.

Example:

v = [2 18 9 6 11 20 25 3]n = 4

【解釋】:問題42837是本次探討的動態序列數組問題中最簡單的一個。求原序列中,最長連續遞增子序列的長度。印象中知乎之前某個帖子看到過,應該是微軟某年的程序員招聘題目還是什麼的...

Test

Code Input

1

%%v = [2 18 9 6 11 20 25 3];n_correct = 4;assert(isequal(subseq(v),n_correct))

2

%%v = [-6 -18 -9 -4 -11 -20 -25 -3];n_correct = 3;assert(isequal(subseq(v),n_correct))

3

%%v = zeros(1,30);n_correct = 1;assert(isequal(subseq(v),n_correct))

4

%%v = [exp(-1) sqrt(2) sqrt(3) exp(1) pi exp(2)];n_correct = 6;assert(isequal(subseq(v),n_correct))

5

%%v = 100:-10:-100;n_correct = 1;assert(isequal(subseq(v),n_correct))

6

%%v = [0:5 1:7 3:9 2:8];n_correct = 7;assert(isequal(subseq(v),n_correct))

【分析】:所得序列應滿足2個條件:從左到右,每個元素都比前一個要大,test suite中沒有相等的情況,所以不用管示例中某2個元素相等的狀況;其次,子序列必須在原序列中連續排布在一起,不能中斷跳過某個取得其他元素。實際上,Note005:Finding longest run中,對於如何在0-1序列中找到最長的0或者最長的1,展開了多種解法的比較與探討,其中有多種思路和命令組合與本問題有密切關聯。本問題和當時討論的情況略有區別,當時探討的是最長的0或者1,因此邏輯數組用起來得心應手,實在不行就取反,這個問題就不行了,Note005的代碼不能原樣照抄。

數值解法:max+diff+find組合

動態子序列長度問題的求解代碼中,max+diff+find屬於保留型經典解法,用的時候注意怎麼構造最裡層的邏輯數組判斷。

function ans = subseq(v) max(diff(find([1 diff(v)<=0 1])));end

假如和Note005中代碼做個比較,比如對0-1序列最長的0和1的長度序列計算可用如下方式:

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

比較兩代碼關鍵一句,發現區別在於對相鄰元素做差,計算非零元素的遞增,在最裡層的diff命令給一個不大於0的指示。這麼處理,其目的是什麼?最簡單的方法是直接用0-1序列測試運行結果,因為老代碼在Note005中已經介紹過,不再浪費篇幅重複:diff+find的連續組合本就是用來尋找數據發生變化的位置、而首尾補1目的是對應湊夠和原序列的位數對等,這兩點並不奇怪。

要說奇怪,或者用出新鮮味道的地方,只能是裡層diff命令上加的那個小於等於0,從下方代碼看出,通過這個小於等於0,恰好提取出那些diff之前,相鄰前後相等、或者發生前1後0變化的位置,也就是tx21變數中所有的邏輯真值1,而僅剩的兩個0,就是前0後1,相減之後大於0的位置。換句話說,遞增就發生這個位置!

>> x=[0 0 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 1 1 1];>> tx21=diff(x)<=0tx21 = 1×21 logical 數組 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1>> tx22=find([1 tx21 1])tx22 = 1 2 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 21 22 23>> tx23=diff(tx22)tx23 = 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1

既然遞增發生在此處,幹嘛不直接把小於等於0改成大於0,不就能直接算出起始位置嗎?因為要求的不是起始位置,而是遞增序列長度,因此這是為下一步所留的後手和布局:tx22變數是補1對齊後的索引查找,對應剛才兩個0的位置,發生索引號斷開的只有2和3兩個元素的索引2-4,以及18和19兩個元素的索引19-21,再做一次diff,兩處遞增序列的長度就用diff直接給減出來了。

再深想一步:如果去掉關鍵的大於等於0,最裡層得到的是只是發生0到1或者1到0變化的索引位置,顯然,這結果我題目要求就天差地遠了。

最後,不用1在首尾位置補,換成inf之類也是可以的,因為這反正是個補位對齊索引的操作:diff一次就少一個元素,補2個1或者inf對齊,兩次diff,正好恢復:

function ans = subseq(v) max(diff(find(diff([Inf v -Inf])<=0)));end

數值解法:cumsum+mode

動態序列問題中,用於累積求和的cumsum使用很頻繁,有興趣可以看看我和彬彬的書,裡面有不少例子用到;另一個mode命令獲得序列中出現頻數最高的元素以及出現了多少次,比如:

>> t=randi([1 3],1,12)t = 3 3 1 3 2 1 1 2 3 3 1 3>> [m,idx]=mode(t)m = 3idx = 6

數值為3的元素在t序列中出現6次,為最多。

現在的問題在於:需要獲得遞增序列長度。有點能肯定:既然要求「最長」遞增序列的長度,這是和mode命令建立關聯的關鍵點。

function n = subseq(v) [~,n] = mode(cumsum([0 diff(v)<=0]));end

想知道遞增序列究竟能有多長,其實關鍵在於遞增序列發生變化的起始和上一次變化或序列終止的位置,cumsum正好執行這個動作——使用diff(v)<=0獲得的是0-1序列,任何數加0不發生變化,0正好是在diff過程中,後數減前數大於0的,為什麼?不是要求遞增序列,為何還要變成0,搞成1不好嗎?看下面的累積求和結果:

>> cumsum([1 zeros(1,4) randi([0 3],1,4)])ans = 1 1 1 1 1 4 5 8 8

最開始那段連續1,就是所需要的,與原序列遞增索引相對應的位置,再用mode一加工...相信意思就很清楚了。方便理解,以test1為例,把完整變化過程放在下面:

>> v = [2 18 9 6 11 20 25 3];>> vx11=[0 diff(v)<=0]vx11 = 0 0 1 1 0 0 0 1>> vx12=cumsum(vx11)vx12 = 0 0 1 2 2 2 2 3>> [m,idx]=mode(vx12)m = 2idx = 4

最前面0顯然又是用來補位,更重要的是,diff是後數減前數,為前面負責,因此每個發生數字1的節點索引位置,都是遞增序列的開始,這樣,就變成,尋找索引累積增加序列中,出現頻數最大的值了,累加過程中還不用擔心重複,每段連續數值結果都獨一無二:畢竟遇到0,相加結果不變,遇到不是0的,結果一變,就成下一個序列了。

字元串類型解法:strsplit

之前說過這類計算辦法,把序列換成字元串,然後變成字元查找,首先是Tim的解法:

function ans=subseq(v)1+max(cellfun(@length,strsplit(char( +(diff(v)>0)))));end

cellfun里,用@length獲取指定子序列長度,其執行對象是文本序列滿足要求的所有連續長度;數值化字元型用了char命令,之前說過,char只轉換序列中的數字,其他分割數字的空格天然無需處理;再用strsplit一字排開,max遴選最大長度,加1又是給diff造成維度減小做補位。

關於strsplit,該命令和split有部分重疊,比如都用於分隔字元串,分隔字元都輸出cell類型,都能指定特定分隔標識符,例如空格、字元等等,但也有一定區別,最明顯是split函數為string類型配套函數,不過支持正則表達式等方面split相對strsplit有一定弱化,輸出cell是N×1,strsplit則是1×N。

Tim此處使用strsplit的主要構思是考慮空格字元與strsplit的天然契合,diff(v)>0得到0-1邏輯序列,真值1都是相鄰處兩個數據不同且遞增的,其他為0的都是原來相同或者遞減的序列,此時用diff得到的序列和空格ascii碼值相加,是0的通過char還是空格,但是1處與空格碼值32相加變成驚嘆號「!」的碼值33,通過在外圍由strsplit一處理,空格自動被過濾,得到了一系列驚嘆號,而驚嘆號個數肯定跟遞增序列長度有關,這是個巧妙的構思!由於strsplit輸出cell,天然外接一個cellfun對每段長度進行分別計算處理。

正則類解法:regexp

彬彬和我則不約而同,用正則表達式構造連續1序列的做法,在數值型換成字元型的方式上,他用了sprintf,我用的是num2str,二者並無區別,彬彬的做法:

function ans = subseq(v)sum(max(cellfun(length,regexp(sprintf(%d,(diff(v)>0)),1+,match)))) + 1;end

我的做法:

function ans = subseq(v) sum(max(cellfun(@numel,regexp(num2str(diff(v)>0,-6),1+,match))))+1;end

兩個代碼從數值轉為字元串,一個用sprintf指定』%d』參數消空格,一個用num2str,帶-6參數消空格,這倒無所謂,但前面不引人注意的長度命令卻發生一個小意外:彬彬用』length』,注意是字元串形式:inline命令時期的老用法,目的是降size,我用「@numel」,單就一維數組,注意:是一維數組中的元素多少,兩種做法等效。關鍵在於:如果仿效前者,用』numel』會提示出錯,這說明了:MATLAB對於新函數已經斷開了inline形式匿名函數的支持,這點需要注意。此外,我以前喝酒聊天好像聽吳鵬還是誰來過一段八卦:Mathworks公司內部對於字元型匿名函數的格式也不怎麼感冒,某次研討會上工程師也提出:內聯函數形式的語法結構會在以後的版本中慢慢取消,意味著匿名函數字元串書寫方式已慢慢不為新函數所支持,這個題目中,兩個做法一個提示出錯,另一個能通過,恐怕是明證之一,畢竟numel是V2006a之後的新...好吧,相對新一點的函數。

Problem 42838. Increasing sub-sequence (Level 2)

【鏈接】:ww2.mathworks.cn/matlab

【題目】:Given a vector, v, of real numbers, return a positive integer, n, representing the longest non-contiguous increasing sub-sequence contained in v.

Example:

v = [ 2 18 9 6 11 20 25 3]n = 5

【解釋】:還是子序列,還是遞增,不過呢,多了個小小的改動,不需要連續了,只要滿足子序列元素屬於原序列,跳著向前走也行。例如前面的例子,2在最前,但中間隔著18這元素,雖然18比2更大,但接受18,只能構造出相對比較短的子序列(2-18-20-25),而接受後面的9或者6,獲得的長度都是5。

Test

Code Input

1

%%v = [2 18 9 6 11 20 25 3];n_correct = 5;assert(isequal(subseq(v),n_correct))

2

%%v = [-2 -18 -9 -6 -11 -20 -25 -3];n_correct = 4;assert(isequal(subseq(v),n_correct))

3

%%v = zeros(1,30);n_correct = 1;assert(isequal(subseq(v),n_correct))

4

%%v = [exp(-1) sqrt(2) sqrt(3) exp(1) pi exp(2)];n_correct = 6;assert(isequal(subseq(v),n_correct))

5

%%v = 100:-10:-100;n_correct = 1;assert(isequal(subseq(v),n_correct))

6

%%v = [9 0:5 1:7 3:9 2:8];n_correct = 10;assert(isequal(subseq(v),n_correct))

7

%%v = [0 8 4 12 2 10 6 14 1 9 5 13 3 11 7 15];n_correct = 6;assert(isequal(subseq(v),n_correct))

【分析】:題意非常清楚,連test都和前面類似,但額外增加的test6通過難度比較大,還有test3的全0序列搞不好也會出問題,原因在於長度達到25以上,如果用nchoosek這類組合函數獲得某個元素後續的元素全排列方式,運行時間將難以忍受,也無法通過cody的伺服器。

非常低效且漏洞百出的nchoosek解法

之所以要專門開闢一節,而且是首節就提出反面例子。主因有二,一是強調不同的演算法對於實際問題的運行求解,效果可能判若雲泥,此類問題在cody比比皆是;二是因為寫出來不怕得罪人:畢竟這傻啦吧唧的解法是我自己提交的...

function ans = subseq(v)% better algorithm is needed for efficiently solving this problem. 1;if numel(v)>25 numel(unique(v));returnelse for k=1:min([numel(v) 10]) if nnz(all(diff(nchoosek(v,k))>0)),ans=k;end endendend

言歸正傳,貌似我這段代碼方案還沒那麼糟糕——如果序列都像test1那樣,元素數量很少(numel(v)<10),畢竟組合數在10以下序列中,低下的效率還不怎麼能顯示出來,但運行test3和test6,就會感到異常抓狂。所以函數前4行就根本在cheat,因為數量超過25根本無法通過算例3或者6,因此1用來應對test3、外層if則是針對test6,其結果正好是unique(v),真正的解法應該是裡層那個for循環,而如果去掉外層if和ans=1的包裝,運行效率將醜陋的難以想像,因為在這個循環中要計算28抽取10的組合數,官方命令在perms還是nchoosek的命令中給出提示:在10以上的序列中做排列組合,可能無法在短時間內得到結果。不過這個漏洞百出的解法,其意思倒是比較清楚,不停尋找某個元素後續的所有元素組合,看看是否是遞增,做一個極度低效的遍歷。

既然組合的方法行不通,還有沒有更好的思路呢?

幾種循環解法

動態數組序列難點之一是找不到合理的長度參照:它可能是不同元素的組合,所以合適的循環是必要的,下面提供幾種循環求解的思路:

function dmax = subseq(v)numlen = length(v);d = zeros(size(v));d(1) = 1;dmax = 1;for i = 1 : numlen submax = 0; for j = 1 : i if v(j) < v(i) submax = max(d(j), submax); end end d(i) = submax + 1; dmax = max(dmax, d(i));end

彬彬給出了看似普普通通,卻不易理解的求解思路,說不易理解,關鍵在子序列元素兩個大小比較條件,分別放在二重循環裡層if條件表達式(v(j) < v(i))和子序列長度計數表達式(submax = max(d(j), submax))中。他的代碼非常規範,看變數名稱很容易理解其意:首先,輸出dmax,顯然是變數d最大值,循環逐步進行元素比較,序列後面到底幾個比前面元素大,其數量就存儲在d中;但這個數量的存儲就是代碼亮點:if條件執行表達式submax = max(d(j), submax),蘊含著隱藏判斷條件:即使滿足判斷條件,靠前元素v(j)小於靠後元素v(i),但遞增子序列長度submax也要受d(j)元素值限制,因為submax在每個i循環中要重新置零,所以其增速滯後於d,以test1為例,經8次i循環,d數值如下:

>> dd = 1 2 2 2 3 4 5 2

循環運行結果說明:裡層j循環中獲得的子序列長度submax上限受控於d(j)的數值。裡面每個數值說明從第1個元素到第k(k=1:8)個元素,遞增子序列的最大長度,例如從v(1)到v(4),最多只能找到一個長度為2的遞增序列。

再看另一種循環解法,相對前一種,好理解多了:

function n = subseq(v)n = 0;for i=1:length(v) j=1; while(j<=n && T(j) < v(i)) j = j + 1; end T(j) = v(i); n = length(T);endend

外層i循環代表循環原序列第1:n個元素,裡層j循環則重新構造新序列T,每次存儲前一次i循環的v(i),即:T(j)=v(i),注意:第1次肯定存不上,因為第1次n=0,子序列沒開始構造,長度為0,j每次重置為1,肯定比此時n大,所以T(j)永遠存上一次的v(i),這就實現了前後元素錯位比較,如遞增,用while循環內的j=j+1跳過序號,T內新存一個數字,如不是遞增,例如對v=[2 18 9 6 11 20 25 3],當i=2時,v(2)=18,此時T=[2],因T(j)=T(1)=2<v(i)=v(2)=18,滿足遞增序列要求,T用j=j+1新增一位:T(j)=T(2)=v(2),即;T=[2 18];當i=3,j=1時,T(1)=2<v(3)=9,用j=j+1跳過1,此時T為[2 18],2的位置上用v(3)代替原來的T(2),因此T序列長度沒變,只是第2個元素從前次循環的18變成此次循環的9:T=[2 9];同理第4次循環用v(4)=6代替元素9:T=[2 6];第i=5次循環:j再被重置為1,while循環第1次因T(1)<v(5)=11,進入while循環得到j=j+1=2,再循環反覆一次時,因T(2)=6,仍小於v(5)=11,所以j=j+1=3新增序號3,如再while一次,因n=numel(T)=2,while循環第1條件不滿足,故跳出,下一句中T(3)=v(5)=11,T變成[2 6 11],以下類推,在i做一次循環找到整個序列v的最長遞增子序列。這同樣是個很棒的思路。

再來一個J.S.Kowontan的循環:

function ans = subseq(v)M = 0*v; 1; for k = v>v M(ans) = prod(max(M(k))+1); ans+1; end max(M);end

代碼外觀簡單,但幾處地方體現出作者J.S.Kowontan對循環和矢量化代碼透徹的理解。

(1). 第1句的M=0*v用於構造與v同維的全0向量,比zeros(size(v))要簡單;後面那個1相當於ans=1,這是for循環對新序列M計算的起點;

(2). 注意循環條件:這一步是代碼的關鍵,想必很少有人這麼寫循環條件,所以仔細看看v>v』的含義,仍以test1中的v=[2 18 9 6 11 20 25 3]為例:

>> v>vans = 8×8 logical 數組 0 1 1 1 1 1 1 1 0 0 0 0 0 1 1 0 0 1 0 0 1 1 1 0 0 1 1 0 1 1 1 0 0 1 0 0 0 1 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 0

暫不談邏輯數組的意義,先說說怎麼運行的,對2016b之前的老版本,因為普通矩陣運算不支持單一維隱式擴展,執行v>v』提示出錯,老版本用下面的代碼獲得同樣結果:

>> bsxfun(@gt,v,v)ans = 8×8 logical 數組 0 1 1 1 1 1 1 1 0 0 0 0 0 1 1 0 0 1 0 0 1 1 1 0 0 1 1 0 1 1 1 0 0 1 0 0 0 1 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 0

單一維擴展和矩陣基本操作的合併,初時令人感到不習慣,比如我自己,有時覺得不理解:既有bsxfun,何必讓隱式擴展的版圖混淆、甚至算是「侵犯」到普通加減乘除或者邏輯操作運算中去呢?這個觀點見仁見智,不過經過一段時間的練習和適應,現在對於類似v>v』這樣的操作已習慣。

言歸正傳,2個問題:結果怎麼算出來的,以及放在循環條件中的作用?先回答第1個問題

(1). 邏輯數組的結果怎麼算出來的?

因為v』的結果是:

>> vans = 2 18 9 6 11 20 25 3

因此8×8數組第1行:[2>2(false) 18>2(true) 9>2(true) ... 3>2(true)];同理,8×8邏輯數組第1列判斷v(i)>transpose(v),其中i=1,即:[2>2(false);2>18(false);...;2>3(false)],余類推。

(2). 把v>v』放在循環條件中的作用

首先還是再強調一遍:MATLAB循環中,每次循環執行的是當前列,例如本例中8×8邏輯數組,每次循環k相當於v>v』邏輯數組的第k列。下面說說把v>v』邏輯數組放循環條件的原因,以第2列為例,v(2)=18,按前一條敘述,第2列代表[18>2(1);18>18(0);18>9(1);...;18>25(0);18>3(1)],循環內部,此時ans=2,M=[1 0 ... 0],M(k)此時的最大值是1;按「prod(max(M(k))+1)」得到的結果,M(2)=2,代表v的前2個元素能構成的遞增序列最大是2個;循環到第3列,發現k第2-3個元素是0,而此時M除了前兩個元素不為0(1和2),後面全0,因此M(k)的最大值必然還是1,因此M(3)又是2,第4次循環的結果也是2,在4次循環後,M=[1 2 2 2 0 0 0 0],第5次則有變化,此時k=[1 0 1 1 0 0 0 1]』;M(k)結果為[1 0 2 2 0 0 0 0];最大值是2,執行prod(max(M(k))+1)得到3,換句話說,在前5個元素中,能形成的遞增序列,最大長度是3,此時M=[1 2 2 2 3 0 0 0];其餘類推。這是個出色的構思,代碼簡潔凝鍊,一重循環內除編號遞增外只有一句,就完成了前面兩個二重循環的工作。

同一種思路,還可以換個手法表達:

function ans = subseq(v) ; for i = triu(v<v) [ans max([ans(i)+1 1])]; end max(ans);end

J.S.Kowontan代碼用prod避免第1次循環時,M(k)為空而出現運算錯誤:因第1次循環k肯定全0;LY Cao的代碼使用max避免這一情況,二者殊途同歸;最明顯差異在於循環條件的處理:LY Cao發現每次循環都是從前往後,因此每次最大遞增序列肯定以:1:1,1:2,1:3,...的順序向後,顯然,遞增序列出現在邏輯數組的上三角陣內,因此代碼用triu控制所取得的最大遞增序列,每次處理長度依次是1,2,3,...,n個。

接著再看Alfonso的循環做法:

function ans = subseq(V) []; for v=V [ans max([0 ans(V(ans>0)<v)])+1]; end max(ans);end

上述代碼邏輯條件放在循環體內,這段代碼也表現出這位cody頂級大神在索引方面高超的應用技巧,把兩次索引定址、兩次邏輯索引構造集中在「ans(V(ans>0)<v)」方寸間,卻讓各自功能彼此相安無事,巧妙銜接。

循環體內,外層max當然還是起兩個作用:取最大長度和避免首次循環空矩陣時出現計算錯誤;代碼關鍵部分是循環體內的ans(V(ans>0)<v),從裡到外分析一下結果:

(1). 第1次循環結果不再多說,ans=1;

(2). 第2次循環從裡面開始ans>0即1>0,此為true,即結果是1,V(ans>0)=V(1)=2;在往外推一層:V(ans>0)<v,即:2<v=18,邏輯true(1),即:ans(1)=1;結果代入到外層,是max([0 ans(1)]+1)=max([1 2])=2,換句話說,前2個元素組成的序列,最大遞增子序列長度為2;第3次循環,ans經前2次循環變為ans=[1 2],因此ans>0得到邏輯數組結果應是[1 1],因此V([1 1])=[2 18],V(ans>0)<v=[2 18]<9=[1 0],代入ans(V(ans>0)<v)])=ans([1 0]),即:1,再到外圍變:max([0 1+1])=max([0 2])=2,這次循環過後,ans=[1 2 2];同理第4次循環結果2,ans=[1 2 2 2];第5次循環,ans>0勢必是4個邏輯真值[1 1 1 1],V([1 1 1 1])=[2 18 9 6],顯然[2 18 9 6]<11,11是第5次循環的v值,也就是V(5),結果變[1 0 1 1],ans([1 0 1 1])=[1 2 2],回外圍,max([0 1 2 2])+1=3,經這次循環,ans=[1 2 2 2 3];

V裡層的ans>0用於取得前面循環的遞增序列元素索引;

V(ans>0)相當於對該索引從原序列中取得實際數值;

V(ans>0)<v是第2次邏輯判斷,用於當擴維新增一個元素v之後,判定滿足遞增條件的索引值

外層再一次ans(V(ans>0)<v),則擴充為新的遞增序列長度值,換句話說ans從每次循環的[1]-[1 2]-[1 2 2]-[1 2 2 2]-[1 2 2 2 3]-[1 2 2 2 3 4]-[1 2 2 2 3 4 5]-[1 2 2 2 3 4 5 2],其中任意一次的結果,其中的第i個元素就是前i項元素能夠構造出最長的遞增序列長度,例如第6項為4,則前6個元素,最長遞增序列的長度為4.

結果和前面J.S.Kowontan以及LY Cao的一模一樣。從這個例子就可以看出Alfonso對於邏輯索引構造命令的組合,嫻熟到何等程度!

再看最後一個循環解法:

function ans = subseq(v)T = [];for i=v T(norm(find(T<i,1,last))+1) = i;endlength(T);end

這位ID:A.Sawas的夥計名不見經傳,做了幾道題很快從cody消失,這個問題他提交了幾個解法,這個版本我個人比較心水,代碼中用另一種方法避免第1次循環空矩陣運算錯誤:範數命令norm,因為norm([])=0,這是創新處;第2點比較喜歡他的序列組成方式:

第1次循環後,T結果為2,此時T=2;

第2次循環,i=18,find(T<i,1,』last』)=1,所以T(norm(1+1))=T(2)=18,因原來T小於i,允許擴維;

第3次循環i=9,T=[2 18],T<i結果一真一假:[1 0],通過find命令查找到最後一個1就在第1位,又是norm(1+1)=2,這次把T(2)從前一次的18替換為更小的9,但T總維度還是1×2;第4次循環把9替換為更小的6,同樣維度沒變還是2,第5次循環中,i=11,T=[2 6],這時T<11終於兩個都為真,代表T可以再度擴維,在遞增序列添加滿足遞增條件的新數,此時find最後一個1,應該是T第2個位置,因此find的結果為2,norm(2+1)=3,所以T(3)=11,其餘類推。這也是非常棒的循環遞進演算法!它和前面演算法不同處在於:通過循環直接給出實際遞增序列(T),而不是序列長度。

以上一共給出6種不同循環解法,各有其需仔細揣摩的精彩,體會後才會找到頂級高手依據問題特徵,選擇合適MATLAB變成的演算法,並進行代碼設計的思路。最重要是:得弄清楚他們為什麼會這麼考慮問題?為什麼會在這裡使用這個函數而不是那個?為什麼會選擇這樣的順序進行命令組合?揣摩多了,代碼水平自然就會有所提高,在擁有自己的特點之前,模仿的過程是多數天資並不高的人必須的步驟——比如我自己。

利用圖論中的有向圖函數:digraph

思路能從遞增子序列問題,轉到圖論連接矩陣構造,想到這樣思路的人,數學思維水平和跨度都不簡單!Jan Orwat就提交了這樣一個以圖論工具箱有向圖構造的digraph和有向圖內所有節點對間最短距離的distances兩個命令為基礎的代碼方案:

function n = subseq(v) n = 1-min(min(distances(digraph(-tril(bsxfun(@lt,v,v))))));end

劉鵬按照R2016b新版本歸併了的隱式擴展操作的方式修改和簡化了上述代碼。

function ans = subseq(v) 1-min(min(distances(digraph(-tril(v<v)))));end

寫出這樣的代碼,首先應理解圖論中,各點間連接矩陣的構造思想,尤其要清楚,連接矩陣構造在MATLAB中能用v<v』恰當地表述——我甚至都沒好意思提及該表述中,隱含對新版本中關於隱式擴展操作手法的熟悉:

>> v = [2 18 9 6 11 20 25 3];>> v<vans = 8×8 logical array 0 0 0 0 0 0 0 0 1 0 1 1 1 0 0 1 1 0 0 1 0 0 0 1 1 0 0 0 0 0 0 1 1 0 1 1 0 0 0 1 1 1 1 1 1 0 0 1 1 1 1 1 1 1 0 1 1 0 0 0 0 0 0 0

這是Test 1輸入向量v按v<v』獲得的邏輯矩陣,其中包含圖所有節點間通過錯序比較得到的8×8、共計64個判斷結果,邏輯矩陣中的元素1,即邏輯真值,意味著圖第(i,j)位置處,滿足v(i)<v』(j),即i→j節點對間處於連通狀態;外層接入一個tril提取下三角位置的元素,保證序列保持遞增順序,元素加負號意味著構成有向圖,所求結果是最大距離的相反數;digraph負責根據節點對之間的連通狀況,構造有向圖。

由於n個節點之間只能出n-1個距離,由兩個min獲得最大距離結果應該再加1。

哦,對了...

別忘記tril外面帶的負號,這意味著結果是用1減去最大值的相反數。

這個思路真是精巧,代碼實施無可挑剔,首先把遞增關係的概念轉化成鄰接矩陣,這是個極大的思維轉換幅度;其次,序列自左至右順序用下三角矩陣表示自動獲得遞增關係,最後,最大值的計算,用到的方法是倒過來:計算相反數的最小值。

...

WOW!!!

小結

從這兩道題目,可以看出動態序列問題的美妙所在,由於條件不定,哪個數據應當被選,條件非常靈活和隱蔽,同時還存在多解,用一般for+if很難處理。不過還好,幾位頂級高手給出角度不同的求解方法,並使用多個不同的函數命令組合,大開眼界同時,自身對於動態序列問題計算,也算懵懵懂懂有了還算不錯的開端,今後遇到類似問題,想必應該心裡略有底。

此外,演算法選擇的合適與否,計算效率判若雲泥,最開始我構思利用nchoosek遍歷顯得非常笨拙而低效,而按照錯序和對位這兩個動態子序列問題中最關鍵的關鍵詞構思出的思路,則簡潔而明快。

以上為動態序列問題第1個專題,後面還有幾個,再繼續探討其他規則、或者其他數據類型下的動態序列問題求解思路與方法。


推薦閱讀:

大家都來說說,matlab里有什麼函數,在python里是找不到的?
番外.飛秒脈衝診斷技術
RTLAB非同步通信案例
Matlab中的模態分析
第2節.超短脈衝在色散介質中的傳播模擬

TAG:MATLAB |