2015-07-21 79 views
0

对于200000浮点数据集,此代码需要半个多小时。可以优化这个动态编程代码吗?

import numpy as np 
try: 
    import progressbar 
    pbar = progressbar.ProgressBar(widgets=[progressbar.Percentage(), 
     progressbar.Counter('%5d'), progressbar.Bar(), progressbar.ETA()]) 
except: 
    pbar = list 

block_length = np.loadtxt('bb.txt.gz') # get data file from http://filebin.ca/29LbYfKnsKqJ/bb.txt.gz (2MB, 200000 float numbers) 
N = len(block_length) - 1 

# arrays to store the best configuration 
best = np.zeros(N, dtype=float) 
last = np.zeros(N, dtype=int) 
log = np.log 

# Start with first data cell; add one cell at each iteration 
for R in pbar(range(N)): 
    # Compute fit_vec : fitness of putative last block (end at R) 
    #fit_vec = fitfunc.fitness(
    T_k = block_length[:R + 1] - block_length[R + 1] 
    #N_k = np.cumsum(x[:R + 1][::-1])[::-1] 
    N_k = np.arange(R + 1, 0, -1) 
    fit_vec = N_k * (log(N_k) - log(T_k)) 

    prior = 4 - log(73.53 * 0.05 * ((R+1) ** -0.478)) 
    A_R = fit_vec - prior #fitfunc.prior(R + 1, N) 

    A_R[1:] += best[:R] 

    i_max = np.argmax(A_R) 
    last[R] = i_max 
    best[R] = A_R[i_max] 

# Now find changepoints by iteratively peeling off the last block 
change_points = np.zeros(N, dtype=int) 
i_cp = N 
ind = N 
while True: 
    i_cp -= 1 
    change_points[i_cp] = ind 
    if ind == 0: 
     break 
    ind = last[ind - 1] 
    change_points = change_points[i_cp:] 

print edges[change_points] # show result 

第一个循环非常慢,因为在每次迭代中数组长度为R,即增加,导致N^2复杂度。

有什么办法可以进一步优化这个代码,例如:通过预先计算?我也很满意使用其他编程语言的解决方案。

+1

http://codereview.stackexchange.com/ – Mihai

+0

[跨代发表于代码评论](http://codereview.stackexchange.com/q/97539/9357) –

+0

这并没有得到太多的关注,因为它缺乏一个关键标签,'numpy'。 CR上没有那么多'numpy'知识丰富的海报。加CR对于问题格式更挑剔。但我同意CR的意见,这个问题需要更多的解释。它也应该有一个小的测试数据集。 'progressbar'是一个不必要的复杂因素。 – hpaulj

回答

0

我可以复制A_R(直至fit-prior步骤),为上三角NxN矩阵与:

def trilog(n): 
    nn = n[:-1,None]-n[None,1:] 
    nn[np.tril_indices_from(nn,-1)]=1 
    return nn 

T_k = trilog(block_length) 
N_k = trilog(-np.arange(N+1)) 
fit_vec = N_k * (np.log(N_k) - np.log(T_k)) 
R = np.arange(N)+1 
prior = 4 - log(73.53 * 0.05 * (R ** -0.478)) 
A_R = fit_vec - prior 
A_R = np.triu(A_R,0) 
print(A_R) 

我还没有经计算的逻辑和应用best工作。

我只对小数组做过这件事。对于你的完整问题,相应的矩阵对于我的记忆来说太大了。

B=np.ones((200000,200000),float) 

所以,仅仅从内存的考虑,你可能会坚持迭代for R in range(N)