2016-01-08 2 views
6

Мне нужно сгенерировать 1D-массив, где повторяющиеся последовательности целых чисел разделяются случайным числом нулей.Создать 1d numpy с кусками случайной длины

До сих пор я использую следующий код для этого:

from random import normalvariate 

regular_sequence = np.array([1,2,3,4,5], dtype=np.int) 
n_iter = 10 
lag_mean = 10 # mean length of zeros sequence 
lag_sd = 1 # standard deviation of zeros sequence length 

# Sequence of lags lengths 
lag_seq = [int(round(normalvariate(lag_mean, lag_sd))) for x in range(n_iter)] 

# Generate list of concatenated zeros and regular sequences 
seq = [np.concatenate((np.zeros(x, dtype=np.int), regular_sequence)) for x in lag_seq] 
seq = np.concatenate(seq) 

Это работает, но выглядит очень медленно, когда мне нужно много длинных последовательностей. Итак, как я могу его оптимизировать?

+0

Возможно, стоит сначала генерировать нулевой инициализированный массив, а затем помещать 'regular_sequence' в тщательно выбранные случайные места в этом массиве? –

+0

@Istrel, что вы будете делать, когда вы произнесете отрицательное число из своего 'N (10, 1)'? Шансы того, что это происходит, астрономически малы, но если эти параметры не исправлены, вы можете столкнуться с этой проблемой. Можно ли просто выбирать из однородного распределения по некоторому диапазону? –

+0

@N. Wouda благодарит вас за предложение. В моем случае любое значение, меньшее или равное 0, будет означать отсутствие отставания между последовательными «регулярными» последовательностями. Итак, будет еще одна строка кода, которая заменит любое отрицательное значение 'lag_seq' на 0. – Istrel

ответ

4

Вы можете предварительно вычислять индексы, где повторяющиеся regular_sequence элементов должны быть введены, а затем установить те с regular_sequence в векторизованных образом , Для предварительного вычисления этих индексов можно использовать np.cumsum, чтобы получить начало каждого такого chunk из regular_sequence, а затем добавить непрерывный набор целых чисел, расширяющийся до размера regular_sequence, чтобы получить все индексы, которые необходимо обновить.Таким образом, реализация будет выглядеть примерно так -

# Size of regular_sequence 
N = regular_sequence.size 

# Use cumsum to pre-compute start of every occurance of regular_sequence 
offset_arr = np.cumsum(lag_seq) 
idx = np.arange(offset_arr.size)*N + offset_arr 

# Setup output array 
out = np.zeros(idx.max() + N,dtype=regular_sequence.dtype) 

# Broadcast the start indices to include entire length of regular_sequence 
# to get all positions where regular_sequence elements are to be set 
np.put(out,idx[:,None] + np.arange(N),regular_sequence) 

времени выполнения тестов -

def original_app(lag_seq, regular_sequence): 
    seq = [np.concatenate((np.zeros(x, dtype=np.int), regular_sequence)) for x in lag_seq] 
    return np.concatenate(seq) 

def vectorized_app(lag_seq, regular_sequence): 
    N = regular_sequence.size  
    offset_arr = np.cumsum(lag_seq) 
    idx = np.arange(offset_arr.size)*N + offset_arr 
    out = np.zeros(idx.max() + N,dtype=regular_sequence.dtype) 
    np.put(out,idx[:,None] + np.arange(N),regular_sequence) 
    return out 

In [64]: # Setup inputs 
    ...: regular_sequence = np.array([1,2,3,4,5], dtype=np.int) 
    ...: n_iter = 1000 
    ...: lag_mean = 10 # mean length of zeros sequence 
    ...: lag_sd = 1 # standard deviation of zeros sequence length 
    ...: 
    ...: # Sequence of lags lengths 
    ...: lag_seq = [int(round(normalvariate(lag_mean, lag_sd))) for x in range(n_iter)] 
    ...: 

In [65]: out1 = original_app(lag_seq, regular_sequence) 

In [66]: out2 = vectorized_app(lag_seq, regular_sequence) 

In [67]: %timeit original_app(lag_seq, regular_sequence) 
100 loops, best of 3: 4.28 ms per loop 

In [68]: %timeit vectorized_app(lag_seq, regular_sequence) 
1000 loops, best of 3: 294 µs per loop 
0

Я думал, что ответ, отправленный по этому вопросу, имел хороший подход, используя двоичную маску и np.convolve, но ответ был удален, и я не знаю почему. Здесь речь идет о 2 проблемах.

def insert_sequence(lag_seq, regular_sequence): 
    offsets = np.cumsum(lag_seq) 
    start_locs = np.zeros(offsets[-1] + 1, dtype=regular_sequence.dtype) 
    start_locs[offsets] = 1 
    return np.convolve(start_locs, regular_sequence) 

lag_seq = np.random.normal(15,1,10) 
lag_seq = lag_seq.astype(np.uint8) 
regular_sequence = np.arange(1, 6) 
seq = insert_sequence(lag_seq, regular_sequence) 

print(repr(seq)) 
1

Не является тривиальной проблемой, поскольку данные являются несогласованными. Производительность зависит от длины . Возьмите пример квадратного проблемы: много, долго, регулярной и нули последовательность (n_iter==n_reg==lag_mean):

import numpy as np 
n_iter = 1000 
n_reg = 1000 
regular_sequence = np.arange(n_reg, dtype=np.int) 
lag_mean = n_reg # mean length of zeros sequence 
lag_sd = lag_mean/10 # standard deviation of zeros sequence length 
lag_seq=np.int64(np.random.normal(lag_mean,lag_sd,n_iter)) # Sequence of lags lengths 

Сначала ваше решение:

def seq_hybrid(): 
    seqs = [np.concatenate((np.zeros(x, dtype=np.int), regular_sequence)) for x in lag_seq] 
    seq = np.concatenate(seqs) 
    return seq 

Тогда чистый NumPy один:

def seq_numpy(): 
    seq=np.zeros(lag_seq.sum()+n_iter*n_reg,dtype=int) 
    cs=np.cumsum(lag_seq+n_reg)-n_reg 
    indexes=np.add.outer(cs,np.arange(n_reg)) 
    seq[indexes]=regular_sequence 
    return seq 

A для зольного золя социологическое загрязнение:

def seq_python(): 
    seq=np.empty(lag_seq.sum()+n_iter*n_reg,dtype=int) 
    i=0 
    for lag in lag_seq: 
     for k in range(lag): 
      seq[i]=0 
      i+=1 
     for k in range(n_reg): 
      seq[i]=regular_sequence[k] 
      i+=1  
    return seq 

И как раз во время компиляции с Numba:

from numba import jit 
seq_numba=jit(seq_python) 

Тесты теперь:

In [96]: %timeit seq_hybrid() 
10 loops, best of 3: 38.5 ms per loop 

In [97]: %timeit seq_numpy() 
10 loops, best of 3: 34.4 ms per loop 

In [98]: %timeit seq_python() 
1 loops, best of 3: 1.56 s per loop 

In [99]: %timeit seq_numba() 
100 loops, best of 3: 12.9 ms per loop 

Ваше гибридное решение является столь же скоростью, как и чистый NumPy один в этом случае потому что производительность зависит, главным образом, от внутреннего цикла. И ваш (нули и конкатенаты) - это числовое. Как и ожидалось, решение python медленнее с традиционным значением около 40x. Но numpy здесь не является оптимальным, потому что он использует фантастическую индексацию, необходимую с несогласованными данными. В этом случае numba может помочь вам: минимальные операции выполняются на уровне C, для 120x коэффициент усиления на этот раз по сравнению с решением python.

Для других значений n_iter,n_reg усилений фактора по сравнению с решением питона является:

n_iter= 1000, n_reg= 1000 : seq_numba 124, seq_hybrid 49, seq_numpy 44. 
n_iter= 10, n_reg= 100000 : seq_numba 123, seq_hybrid 104, seq_numpy 49. 
n_iter= 100000, n_reg= 10 : seq_numba 127, seq_hybrid 1, seq_numpy 42. 
+0

Будет ли назначение slice быстрее, чем for-loop здесь? –

2

Лучший подход, я думаю, было бы использовать свертку. Вы можете определить длину задержки, объединить это с длиной последовательности и использовать это, чтобы определить начальную точку каждой регулярной последовательности. Установите эти начальные точки на ноль, затем сверните с вашей обычной последовательностью, чтобы заполнить значения.

import numpy as np 

regular_sequence = np.array([1,2,3,4,5], dtype=np.int) 
n_iter = 10000000 
lag_mean = 10 # mean length of zeros sequence 
lag_sd = 1 # standard deviation of zeros sequence length 

# Sequence of lags lengths 
lag_lens = np.round(np.random.normal(lag_mean, lag_sd, n_iter)).astype(np.int) 
lag_lens[1:] += len(regular_sequence) 
starts_inds = lag_lens.cumsum()-1 

# Generate list of convolved ones and regular sequences 
seq = np.zeros(lag_lens.sum(), dtype=np.int) 
seq[starts_inds] = 1 
seq = np.convolve(seq, regular_sequence) 

Этот подход берет что-то вроде 1/20th время на больших последовательностей, даже после изменения версии, чтобы использовать Numpy генератор случайных чисел.

+0

let 'm, p = regular_sequence.size, seq.size'. 'np.convolve' дает 0 (mp) решение. Я думаю, что другое решение здесь 0 (p), так что асимптотически лучше. –

+0

Возможно, но свертка оптимизирована с использованием циклов C, а «фантазийное индексирование» - нет. Таким образом, в реальном коде он должен быть быстрее. – TheBlackCat

+0

Не понимаю. В моем тесте время исполнения исполняемого решения растет, как regular_sequence.size^1, когда свертка времени выполнения решения растет, как regular_sequence.size^2. для 'n_iter == n_reg == lag_mean = 1000', convolve в 100 раз медленнее. –