Re: [討論] 簡單的計算 想不到暴力解之外的方法

[celestialgod (天)]

%% function rle
% to comput the lengths and values along a vector
function [len, values] = rle(x)
i = [find(x(2:end) ~= x(1:end-1)); length(x)];
len = diff([0; i]);
if nargout > 1
    values = x(i);
end
end

%% function inverse_rle
% produce a vector with lengths and coresponding values
function out = inverse_rle(len, values)
res = arrayfun(@(x, y) y*ones(x, 1), len, values, 'UniformOutput', false);
out = cat(1, res{:});
end

%% main
% data generation
N = 10000;
occurPoints = unique(randi(N, 100, 1));
if mod(length(occurPoints), 2) == 1
    occurPoints = [occurPoints; N];
end
dat = zeros(N, 1);
for i = 1:2:length(occurPoints)
    dat((occurPoints(i)+1):occurPoints(i+1)) = 1;
end

%% data generation - method 2
% len = diff([0; occurPoints; N]);
% values = repmat(0:1, 1, ceil(length(len)/2));
% values = values(1:length(len))';
% v = inverse_rle(len, values);
% isequal(v, dat) % 1

% part 1
[len, values] = rle(dat);
valuesCusum = cumsum(values);
values(values==1) = valuesCusum(values==1);
out_part1 = inverse_rle(len, values);

% part 2
[len, values] = rle(dat);
tmp = values(values == 0);
tmp(len(values==0) < 100) = 1;
values(values == 0) = tmp;
outTmp = inverse_rle(len, values);
[len2, values2] = rle(outTmp);
valuesCusum = cumsum(values2);
values2(values2==1) = valuesCusum(values2==1);
out_part2 = inverse_rle(len2, values2);

% part 2 - method 2 (faster, about 2 times)
[len, values] = rle(dat);
[valuesLen, valuesValues] = rle(values);
tmp = [0; cumsum(valuesLen)];
valuesLenNew = zeros(length(valuesValues), 1);
k = 1;
for i = 1:(length(tmp)-1)
    valuesLenNew(i) = sum(len((tmp(i)+1):tmp(i+1)));
end
valuesCusum = cumsum(valuesValues);
valuesValues(valuesValues==1) = valuesCusum(valuesValues==1);
out_part2_2 = inverse_rle(valuesLenNew, valuesValues);

isequal(out_part2, out_part2_2) % 1

第三部分就利用類似第二部分的概念去寫吧，當作練習吧XDD


※ 引述《JamesChen ( )》之銘言：
: 問題很簡單，分兩個部分
: 一串 0 1 數列
: 大致長得像是 0 0 0 0 0 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 1 1 1 1 0 0 0
: 簡單來說 1 代表 A 事件發生 0 代表沒有 然後這是每個時間點的紀錄
: 所以 A 一旦發生不會只發生 1 個點就結束 一定是一串
: 我想要做的是創一個新的數列 第 N 次 一串 1 變成一串 N
: 以上面的例子來說就是後面 4 個  1 都變成 2
: 我只想到 for loop + if 硬幹的方法
: 但實際資料很長  又有好多受試者  感覺很耗時間
: 第二個部份是
: 如果兩串 1 之前的 0 少於 200 個  需要把這兩串 1 合併 (中間的 0 都當作 1)
: 我一樣只想到硬解
: 我猜是小弟我不夠熟 Matlab 平常都在用一些自己習慣的 function
: 沒有做過類似的事情  但應該都有速解
: 希望高手可以幫個忙
: 甚至提點一下就好

--
※ 發信站: 批踢踢實業坊(ptt.cc), 來自: 123.205.27.107
※ 文章網址: https://www.ptt.cc/bbs/MATLAB/M.1441424461.A.57D.html
推 JamesChen: 回應超快  感謝... 好神
第二部分可以再加速，直接用values跟lens做新的lens跟values

不過這部分要再想一想怎麼做XD (已新增)
推 profyang: 其實c大你這樣寫以演算法而言也算是暴力法就是...只是
→ profyang: 你向量化得漂亮
→ profyang: 不過我也想不到如何將演算法變成非暴力法就是...感覺
→ profyang: 就ㄧ定要暴力法
對我來說這只是一種資料整理 我沒有想過要用一個演算法去解決這個問題

資料整理只要有效率就好了拉XDD

況且這個方法 1e6筆、100個change點的資料也只要0.2秒

畢竟我是統計出身...我並不在意是不是暴力，只在乎是否能夠整理成我要的資料

演算法就交給其他人做吧XDDDDDD


可能建議原PO標題改成如何加速這段程式，就不需要著眼在是不是暴力法這件事

另外就是原PO可能使用的方式是 如果資料有N筆，有M個change(0->1 or 1->0)的點

那麼for-loop + if 就是 O(N)，我這樣計算就只是O(M)

(如果原PO不是一個個算就不是O(N))

就計算來說，我降低了計算的複雜度，而達到加速的功效

根據我自己讀過的內容，廣義而言，演算法就是一種降低計算複雜度的方法

因此，廣義而言，我這也是在提供一種演算法XD

推 profyang: 沒啦 因為原原PO說他不會非暴力法 讓我直覺他想改進演
→ profyang: 算法 但看來他應該只是要比loop快的方法 那對matlab而
→ profyang: 言當然就是向量化了
我覺得我太認真了 哈
推 profyang: 不是喔~以part1來說 你的O(N)的部分等於是在cumsum中
cumsum那裏也是O(M)才對

因為出來的值是change point的值
→ profyang: 你指的O(M)是在後面的values(values==1)中對吧? 這個其
那裏應該也是O(M)
→ profyang: 實也還是matlab內建一個個元素去search吧...所以應該還
→ profyang: 是O(N)才對?
你說的是rle的find嗎？那部分確實是O(N)
→ profyang: 等下我是搞錯什麼 cumsum不就是把一個個元素疊加起來?
→ profyang: 這樣應該是O(N)阿~然後values(values==1)的部分也是一個
→ profyang: 個元素去看你values是否=1 是的話就回傳index給你 這樣
可能有誤會 cumsum確實是O(length(input))

但是我這裡的values 只會有M個值 所以我寫O(M)

N是指原本的資料長度
→ profyang: 還是O(N)阿~ 也就是matlab內建這功能還是有掃個N次的for
推 profyang: 確實有誤會 你的O(N)應該在rle裡面就做完了
對，我的O(N)在rle就做完了
→ profyang: sorry我沒細看XD 不過這樣複雜度就還是O(N)
是阿(茶，我想錯了

所以我還是沒有降低計算的複雜度QQ
※ 編輯: celestialgod (123.205.27.107), 09/05/2015 14:35:47
推 profyang: 沒差 matlab就是向量化就快