others - python - 给定一个向量快速查找类似的所有向量的方法

通过"相似的矢量",我定义了一个矢量,该矢量在一个位置上与给定的矢量相差-1或1。但是,如果给定的元素为零,则只能差1。例如:


similar_vectors(np.array([0,0,0]))



array([[ 1., 0., 0.],


 [ 0., 1., 0.],


 [ 0., 0., 1.]])



similar_vectors(np.array([1,0,2,3,0,0,1]))



array([[ 0., 0., 2., 3., 0., 0., 1.],


 [ 2., 0., 2., 3., 0., 0., 1.],


 [ 1., 1., 2., 3., 0., 0., 1.],


 [ 1., 0., 1., 3., 0., 0., 1.],


 [ 1., 0., 3., 3., 0., 0., 1.],


 [ 1., 0., 2., 2., 0., 0., 1.],


 [ 1., 0., 2., 4., 0., 0., 1.],


 [ 1., 0., 2., 3., 1., 0., 1.],


 [ 1., 0., 2., 3., 0., 1., 1.],


 [ 1., 0., 2., 3., 0., 0., 0.],


 [ 1., 0., 2., 3., 0., 0., 2.]])



我想获得上面提到的similar_vectors(vector)函数的最大快速实现, 在我的模拟中, 它至少运行了数百万次vector ,所以,速度是至关重要的, 我对纯numpy解决方案以及其他语言的包装器都感兴趣, 代码可以是并行的。

我当前的实现是:


def singleOne(length,positionOne): #generates a vector of len length filled with zeros apart from a single one at positionOne


 arr=np.zeros(length)


 arr[positionOne]=1


 return arr



def similar_vectors(state):


 connected=[]


 for i in range(len(state)):


 if(state[i]!=0):


 connected.append(state-singleOne(state.shape[0],i)) 


 connected.append(state+singleOne(state.shape[0],i))


 return np.array(connected)



由于循环,我不能轻松地解除这个问题,这是非常慢的,

为了参考,我附加了similar_vectors(vector)的1000次执行配置文件:


 37003 function calls in 0.070 seconds



 Ordered by: cumulative time



 ncalls tottime percall cumtime percall filename:lineno(function)


 1 0.000 0.000 0.070 0.070 {built-in method builtins.exec}


 1 0.003 0.003 0.070 0.070 <string>:2(<module>)


 1000 0.035 0.000 0.064 0.000 <ipython-input-19-69431411f902>:6(similar_vectors)


 11000 0.007 0.000 0.021 0.000 <ipython-input-19-69431411f902>:1(singleOne)


 11000 0.014 0.000 0.014 0.000 {built-in method numpy.core.multiarray.zeros}


 2000 0.009 0.000 0.009 0.000 {built-in method numpy.core.multiarray.array}


 11000 0.002 0.000 0.002 0.000 {method 'append' of 'list' objects}


 1000 0.000 0.000 0.000 0.000 {built-in method builtins.len}


 1 0.000 0.000 0.000 0.000 {method 'disable' of '_lsprof.Profiler' objects}



时间:

可以创建1s和另一个 -1s的对角矩阵,然后在原始阵列中将前者添加到原始阵列, 然后使用 np.concatenate或者ndarrays来连接两个 :


def similar_vectors(a):


 ones = np.ones(len(a))


 w = np.flatnonzero(a!=0)


 return np.concatenate([np.diag(-ones)[w]+a, np.diag(ones)+a])



运行示例


a = np.array([1,0,2,3,0,0,1])


similar_vectors(a)



array([[0., 0., 2., 3., 0., 0., 1.],


 [1., 0., 1., 3., 0., 0., 1.],


 [1., 0., 2., 2., 0., 0., 1.],


 [1., 0., 2., 3., 0., 0., 0.],


 [2., 0., 2., 3., 0., 0., 1.],


 [1., 1., 2., 3., 0., 0., 1.],


 [1., 0., 3., 3., 0., 0., 1.],


 [1., 0., 2., 4., 0., 0., 1.],


 [1., 0., 2., 3., 1., 0., 1.],


 [1., 0., 2., 3., 0., 1., 1.],


 [1., 0., 2., 3., 0., 0., 2.]])



a = np.array([0,0,0])


similar_vectors(a)



array([[1., 0., 0.],


 [0., 1., 0.],


 [0., 0., 1.]])



只是yatu代码的注释, 您可以就地添加,并且更改dtype,性能提升约33%,


def similar_vectors_v2(a):


 eye = np.eye(len(a), dtype=a.dtype)


 neg_eye = -(eye[a!=0])


 eye += a


 neg_eye +=a


 return np.concatenate([neg_eye, eye])



另外,对于非常大的输出,concatenate 可能会有帮助,

方法#1

这里有一个面具


def similar_vectors_masking(a):


 n = len(a)


 m = n*2+1


 ar = np.repeat(a[None],len(a)*2,0)


 ar.ravel()[::m] -= 1


 ar.ravel()[n::m] += 1


 mask = np.ones(len(ar),dtype=bool)


 mask[::2] = a!=0


 out = ar[mask]


 return out



运行示例-


In [142]: similar_vectors_masking(np.array([0,0,0]))


Out[142]: 


array([[1, 0, 0],


 [0, 1, 0],


 [0, 0, 1]])



In [143]: similar_vectors_masking(np.array([1,0,2,3,0,0,1]))


Out[143]: 


array([[0, 0, 2, 3, 0, 0, 1],


 [2, 0, 2, 3, 0, 0, 1],


 [1, 1, 2, 3, 0, 0, 1],


 [1, 0, 1, 3, 0, 0, 1],


 [1, 0, 3, 3, 0, 0, 1],


 [1, 0, 2, 2, 0, 0, 1],


 [1, 0, 2, 4, 0, 0, 1],


 [1, 0, 2, 3, 1, 0, 1],


 [1, 0, 2, 3, 0, 1, 1],


 [1, 0, 2, 3, 0, 0, 0],


 [1, 0, 2, 3, 0, 0, 2]])



方法#2

对于零填充数组,我们可以使用np.repeat更好地使用进行复制,并且使用一些屏蔽来增加和减少 1s,


def similar_vectors_repeat(a):


 mask = a!=0


 ra = np.arange(len(a))


 r = mask+1


 n = r.sum()


 ar = np.repeat(a[None],n,axis=0)


 add_idx = r.cumsum()-1


 ar[add_idx,ra] += 1


 ar[(add_idx-1)[mask],ra[mask]] -= 1


 return ar



基准测试

在大数组上采用所有发布方法的计时 -


In [414]: # Setup input array with ~80% zeros


 . ..: np.random.seed(0)


 . ..: a = np.random.randint(1,5,(5000))


 . ..: a[np.random.choice(range(len(a)),int(len(a)*0.8),replace=0)] = 0



In [415]: %timeit similar_vectors(a) # Original soln


 . ..: %timeit similar_vectors_flatnonzero_concat(a) # @yatu's soln


 . ..: %timeit similar_vectors_v2(a) # @Brenlla's soln


 . ..: %timeit similar_vectors_masking(a)


 . ..: %timeit similar_vectors_repeat(a)


1 loop, best of 3: 195 ms per loop


1 loop, best of 3: 234 ms per loop


1 loop, best of 3: 231 ms per loop


1 loop, best of 3: 238 ms per loop


10 loops, best of 3: 82.8 ms per loop



使用 Numba 两个解决方案


@nb.njit()


def similar_vectors_nb(data):


 n=data.shape[0]


 out=np.empty((n*2,n),dtype=data.dtype)


 ii=0


 for i in range(n):


 if data[i]==0:


 out[ii,:]=data[:]


 out[ii,i]+=1


 ii+=1


 else:


 out[ii,:]=data[:]


 out[ii,i]+=1


 ii+=1



 out[ii,:]=data[:]


 out[ii,i]-=1


 ii+=1


 return out[0:ii,:]



@nb.njit()


def similar_vectors_nb_2(data):


 n=data.shape[0]



 #Determine the final array size


 num=0


 for i in range(n):


 if data[i]==0:


 num+=1


 else:


 num+=2



 out=np.empty((num,n),dtype=data.dtype)


 ii=0


 for i in range(n):


 if data[i]==0:


 out[ii,:]=data[:]


 out[ii,i]+=1


 ii+=1


 else:


 out[ii,:]=data[:]


 out[ii,i]+=1


 ii+=1



 out[ii,:]=data[:]


 out[ii,i]-=1


 ii+=1


 return out



基准测试器


#https://github.com/MSeifert04/simple_benchmark


from simple_benchmark import benchmark



np.random.seed(0)


data=[]


for i in range(10,1000,10):


 a = np.random.randint(1,5,(i))


 a[np.random.choice(range(len(a)),int(len(a)*0.5),replace=0)] = 0


 data.append(a)



arguments = arguments = {10*i: data[i] for i in range(len(data))}



b = benchmark([similar_vectors_nb,similar_vectors_nb_2,similar_vectors_yatu,similar_vectors_Brenlla,similar_vectors_repeat_Divakar,similar_vectors_masking_Divakar], arguments, warmups=[similar_vectors_nb,similar_vectors_nb_2])



%matplotlib notebook


b.plot()



Timings

...