100 numpy exercises 后50

https://github.com/rougier/numpy-100/blob/master/100%20Numpy%20exercises.md

51.新建一个向量表示（x，y)和(r,g,b)

z = np.zeros(10,
[('position', [('x', float, 1),('y', float, 1)]),
('color', [('r', float, 1), ('g', float, 1), ('b', float, 1)])])

52.在一个10*2的矩阵求出两点间距离

z = np.random.random((10,2))
# atleast_2d查看至少2维的数组
x, y = np.atleast_2d(z[:,0]), np.atleast_2d(z[:,1])
D = np.sqrt((x-x.T)**2 + (y-y.T)**2)

#更快的版本
import scipy
import scipy.spatial
z = np.random.random((10, 2))
D = scipy.spatial.distance.cdist(z, z)

53.怎样原地把float的数组转换成int

z = np.arange(10, dtype=np.float32)
z = z.astype(np.int.32, copy=False)

54.怎样从文件读取(有问题)

from io import StringIO
#假设文件读入IO
s = StringIO("""1,2,3,4,5\n
6, , ,7,8\n
 , ,9,10,11\n""")
z = np.genfromtxt(s, delimiter=",", dtype=np.int)

55.和numpy枚举相同的做法

z = np.arange(9).reshape(3,3)
for index, value = in np.ndenumerate(z) :
    print(index, value)
for index in np.ndindex(z.shape) :
    print(index, z[index])

56.生成一个2维高斯分布

x, y = np.meshgrid(np.linspace(-1,1,10), np.linspace(-1,1,10))
D = np.sqrt(x**x + y**y)
sigma, mu = 1.0, 0.0
G = np.exp(-((D-mu)**2 / (2.0 *sigma**2)))

57.怎样在2维数组随机替换p个元素

n = 10
p = 3
z = np.zeros((n, n))
np.put(z, np.random.choice(range(n*n), p, replace=False), 1)

58.每行减平均数

y = x - x.mean(axis=1, keepdimas=True)

59.以第n列排序为重新排序矩阵

z = np.random.randint(0,10,(3,3))
z[z[:, 1].argsort()]

60.怎样发现给定的二维矩阵有空列

z = np.random.randint(0,3,(3,10))
print(~z.any(axis=0).any())

61.给定数组里的值找到最近的值

z = np.random.uniform(0,1,10)
x = 0.5
m = z.flat[np.abs(z-x).argmin()]

62.如何用生成器计算（1，3）和（3，1）数组的和

a = np.arange(3).reshape(3,1)
b = np.arange(3).reshape(1,3)
it = np.nditer([a, b, None])
for x, y, z in it:
    z[...] = x + y
print(it.operands[2])

63.新建一个数组类有数组属性

class NameArray(np.ndarray) :
    def __new__(cls, array, name="no name") :
        obj = np.asarray(array).view(cls)
        obj.name = name
        return obj
    def __array__finalize__(self, obj) :
        if obj is None : return 
        self.info = getattr(obj, 'name', "no name")

z = NameArray(np.arange(10), "range_10")
print(z.name)

64.新建一个向量，如何用第二个向量做下标向第一个元素加一

z = np.ones(10)
I = np.random.randint(0, len(z), 20)
z += np.bincount(I, minlength=len(z)) #bincount 返回每个数字出现次数

65.根据I，累计X到F

X = [1,2,3,4,5,6]
I = [1,3,9,3,4,1]
F = np.bincount(I, X) #相当于把X当权重付给I，然后计数

66.给定一个(w,h,3)的图片，计算不一样的颜色个数

w, h = 16, 16
I = np.random,randint(0,2,(h,w,3)).astype(np.ubyte)
F = I[...,0]*256*256 + I[...,1]*256 + I[...,2]
n = len(np.unique(F))

67.给定一个四维的数组，如果一次性求最后两维的总和

a = np.random.randint(0,10,(3,4,3,4))
#把最后两维前的保留，-1为根据需要压缩
sum = a.reshape(a.shape[:-2] + (-1,)).sum(axis=-1)

68.给定一个一维向量D，怎样用一个一样大小的向量S描述子集的索引计算D子集的平均值

例如 D=[0.1, 0.2, 0.3] S = [0, 1, 0] 即索引值为0的平均值为 (0.1+0.3)/2=0.2

D = np.random.uniform(0,1,100)
S = np.random.randint(0,10,100)
D_sums = np.bincount(S, weight=D)
D_count = np.bincount(S)
D_means = D_sums / D_counts

69.怎样获得矩阵乘法结果的对角线

a = np.random.uniform(0,1,(5,5))
b = np.random.uniform(0,1,(5,5))
#慢版本
np.diag(np.dot(a, b))
#快版本
np.sum(a * b.T, axis=1)
#更快版本，einsum貌似可以根据字符串来描述加的方法
np.einsum("ij, ji->i", a, b)

70.给定一个[1,2,3,4]，怎样新建一个新的向量两元素间有连续3个零

z = np.array([1,2,3,4,5])
nz = 3
z0 = np.zeros(len(z) + (len(z)-1)*(nz))
z0[::nz+1] = z

71.给定一个(5,5,3)的数组，怎样和一个(5,5)相乘

a = np.ones((5,5,3))
b = 2*np.ones((5,5))
a*b[:,:,None]

72.怎样交换两行

a = np.arange(25).reshape(5,5)
# 1赋值0，0赋值1。[0,1]代表两行
a[[0,1]] = a[[1,0]]

73.给定10个3元组来表述10个三角形，线段两两组合

faces = np.random.randint(0,100,(10,3))
# roll 滚动 repeat [1,2] => [1,1,2,2]
F = np.roll(faces.repeat(2, axis=1), -1, axis=-1)
F = F.reshape(len(F)*3, 2)
F = np.sort(F, axis=1)
G = F.view(dtype=[('p0', F.dtype), ('p1', F.dtype)])
G = np.unique(G)

# 先重复[1,2,3] => [1,1,2,2,3,3]
# 再滚动 [1,1,2,2,3,3] => [1,2,2,3,3,1]
# 这样就生成[1,2],[2,3],[3,1]

74.给定一个bincount数组C，怎样生成一个A，np.bincount == C

c = np.bincount([1,1,2,3,4,4,6])
[0,1,2,3,4,5,6] 0重复c[0]次 1重复c[1]次
a = np.repeat(np.arange(len(c)), c)

75.用计算滑窗内的平均值

def moving_average(a, n=3) :
# cumsum 累计和
    ret = np.cumsum(a, dtype=float)
    ret[n:] = ret[n:] - ret[:-n]
    return ret[n-1:] / n
z = np.arange(20)
moving_average(z, n=3)

76.给定一个一维向量，构建一个新二维向量，第一行z[0] z[1] z[2]，后面每行往后移一位

from numpy.lib import stride_tricks

def rolling(a, window) :
    shape = (a.size - window + 1, window)
    # srtides 代表跳过一行需要多少字节，跳过一列需要多少字节
    strides = (a.itemsize, a.itemsize)
    return stride_tricks.as_strided(a, shape=shapem strides=strides)
z = rolling(np.arange(10), 3)

77.布尔去反，或原地改浮点数符号

z = np.random.randint(0,2,100)
np.logical_not(z, out=z)

z = np.random.uniform(-1.0,1.0,100)
np.negative(z, out=z)

78.给定两个点集合P0, P1描述线和一个点p，怎样比较p点到每条线的距离

def distance(P0, P1, p) :
    T = P1 - P0
    L = (T**2).sum(axis=1)
    U = -((P0[:,0]-p[...,0]) * T[:,0] + (P0[:,1]-p[...,1]) * T[:,1]) / L
    U = U.reshape(len(U), 1)
    D = P0 + U*T - p
    return np.sqrt((D**2).sum(axis=1))

P0 = np.random.uniform(-10,10,(10,2))
P1 = np.random.uniform(-10,10,(10,2))
p = np.random.uniform(-10,10,(1,2))
ditance(P0,P1,p)

79.每个点到每条线的距离

P0 = np.random.uniform(-10,10,(10,2))
P1 = np.random.uniform(-10,10(10,2))
p = np.random.uniform(-10,10,(10,2))
np.array([distance(P0,P1,p_i) for p_i in p])

80.给定一个随机矩阵，写一个函数提取一个子集合，开始点位给定位置，其余用fill元素填充

Z = np.random.randint(0,10,(10,10))
shape = (5,5)
fill  = 0
position = (1,1)

R = np.ones(shape, dtype=Z.dtype)*fill
P  = np.array(list(position)).astype(int)
Rs = np.array(list(R.shape)).astype(int)
Zs = np.array(list(Z.shape)).astype(int)

R_start = np.zeros((len(shape),)).astype(int)
R_stop  = np.array(list(shape)).astype(int)
Z_start = (P-Rs//2)
Z_stop  = (P+Rs//2)+Rs%2

R_start = (R_start - np.minimum(Z_start,0)).tolist()
Z_start = (np.maximum(Z_start,0)).tolist()
R_stop = np.maximum(R_start, (R_stop - np.maximum(Z_stop-Zs,0))).tolist()
Z_stop = (np.minimum(Z_stop,Zs)).tolist()

r = [slice(start,stop) for start,stop in zip(R_start,R_stop)]
z = [slice(start,stop) for start,stop in zip(Z_start,Z_stop)]
R[r] = Z[z]
print(Z)
print(R)

81.给定一个z=[1:15]，怎样生成z[0] z[1] z[2] z[3],每行左移1的新矩阵

z = np.arange(1,15,dtype=np.unit32)
r = stride_tricks.as_strided(z, (11,4), (4,4))

82.计算矩阵的秩

z = np.random.uniform(0,1,(10,10))
# 分解矩阵z = u * np.diag(s) * v
u, s, v = np.linalg.svd(z)
# 大于0
rank = np.sum(s > 1e-10)

83.找众数

z = np.random.randint(0,10,50)
np.bincount(z).argmax()

84.在10*10的矩阵内用3*3的滑窗新建一个向量

z = np.random.randint(0,5,(10,10))
n = 3
i = 1+(z.shape[0]-3)
j = 1+(z.shape[1]-3)
c = stride_tricks.as_strided(z, shape(i,j,n,n), strides=z.strides+z.strides)

85.新建一个2维数组，z[i,j]=z[j,i]，修改一个值能自动更改

class Symetric(np.ndarray):
    def __setitem__(self, index, value):
        i,j = index
        super(Symetric, self).__setitem__((i,j), value)
        super(Symetric, self).__setitem__((j,i), value)

def symetric(Z):
    return np.asarray(Z + Z.T - np.diag(Z.diagonal())).view(Symetric)

S = symetric(np.random.randint(0,10,(5,5)))
S[2,3] = 42

86.假设p个(n,n)的矩阵，p个(n,1)的向量，怎样一次性计算和

p, n = 10, 20
M = np.ones((p,n,n))
V = np.ones((p,n,1))
# 0跟0，2跟1相乘
S = np.tensordot(M, V, axes=[[0,2], [0,1]])

87.给定一个16*16的数组，怎样获得模块和（模块大小4*4）

z = np.ones((16,16))
k = 4
s = np.reduceat(np.add.reduceat(z, np.arange(0,z.shape[0], k),axis=0),np.arange(0,z.shape[1],k), axis=1)

88. Game of Life

def iterate(Z):
    # Count neighbours
    N = (Z[0:-2,0:-2] + Z[0:-2,1:-1] + Z[0:-2,2:] +
         Z[1:-1,0:-2]                + Z[1:-1,2:] +
         Z[2:  ,0:-2] + Z[2:  ,1:-1] + Z[2:  ,2:])

    # Apply rules
    birth = (N==3) & (Z[1:-1,1:-1]==0)
    survive = ((N==2) | (N==3)) & (Z[1:-1,1:-1]==1)
    Z[...] = 0
    Z[1:-1,1:-1][birth | survive] = 1
    return Z
    Z = np.random.randint(0,2,(50,50))
for i in range(100): Z = iterate(Z)
print(Z)

89. N大数

z = np.arange(10000)
np.random.shuffle(z)
n = 5

#slow
z[np.argsort(z)[-n:]]

#fast argpartition 小于n的放在前面，大于的放在后面
z[np.argpartion(-z,n)[:n]]

90.给定一个随机向量，新建一个笛卡尔组合

def cartesian(arrays):
    arrays = [np.asarray(a) for a in arrays]
    shape = (len(x) for x in arrays)

    ix = np.indices(shape, dtype=int)
    ix = ix.reshape(len(arrays), -1).T

    for n, arr in enumerate(arrays):
        ix[:, n] = arrays[n][ix[:, n]]

    return ix

print (cartesian(([1, 2, 3], [4, 5], [6, 7])))

91.怎样从一个数组新建一个记录数组

z = np.array([("Hello" ,2.5, 3),("World", 3.6, 2)])
r = np.core.records.fromarray(z.T, names='col1, col2, col3', formats='S8,f8,i8')

92.给定一个很大的z，怎样求立方

x = np.random.rand(5e7)

%timeit np.power(x,3)
%timeit x*x*x
%time np.einsum('i,i,i->i', x,x,x)

93.给定A(8,3) B(2,2), 怎样找到A的行包括B每一行的元素不管顺序

A = np.random.randint(0,5,(8,3))
B = np.random.randint(0,5,(2,2))

C = (A[..., np.newaxis, np.newaxis] == B)
rows = (C.sum(axis=(1,2,3)) >= B.shape[1]).nonzero()[0]
print(rows)

94.给定一个10*3的矩阵，提取有不同值得行

z = np.random.randint(0,5,(10,3))
e = np.logical_and.reduce(z[:,1:]==z[:,:,-1], axis=-1)
u = z[~e]

95.转换一个整数的向量，变成二进制表达的矩阵

I = np.array([0, 1, 2, 3, 15, 16, 32, 64, 128])
B = ((I.reshape(-1,1) & (2**np.arange(8))) != 0).astype(int)
print(B[:,::-1])

I = np.array([0, 1, 2, 3, 15, 16, 32, 64, 128], dtype=np.uint8)
print(np.unpackbits(I[:, np.newaxis], axis=1))

96.给定一个2维数组，怎样提出唯一的行

Z = np.random.randint(0,2,(6,3))
T = np.ascontiguousarray(Z).view(np.dtype((np.void, Z.dtype.itemsize * Z.shape[1])))
_, idx = np.unique(T, return_index=True)
uZ = Z[idx]
print(uZ)

97.给定两个A，B, 用einsum写出inner outer sum mul功能

a = np.random.uniform(0,1,10)
b = np.random.uniform(0,1,10)

np.einsum('i->',A) #np.sum(a)
np.enisum('i,i->i',a,b) #a * b
np.enisum('i,i',a,b) #np.inner(a,b)
np.einsum('i,j',a,b) #np.outer(a,b)

98.给定一个用(x,y)描述的路，怎样等距抽样

phi = np.arange(0, 10*np.pi, 0.1)
a = 1
x = a*phi*np.cos(phi)
y = a*phi*np.sin(phi)

dr = (np.diff(x)**2 + np.diff(y)**2)**.5 # segment lengths
r = np.zeros_like(x)
r[1:] = np.cumsum(dr)                # integrate path
r_int = np.linspace(0, r.max(), 200) # regular spaced path
x_int = np.interp(r_int, r, x)       # integrate path
y_int = np.interp(r_int, r, y)

99.给定一个整数n和二维数组X，从x中选出符合n分布的行

X = np.asarray([[1.0, 0.0, 3.0, 8.0],
                [2.0, 0.0, 1.0, 1.0],
                [1.5, 2.5, 1.0, 0.0]])
n = 4
M = np.logical_and.reduce(np.mod(X, 1) == 0, axis=-1)
M &= (X.sum(axis=-1) == n)
print(X[M])

100.重新抽样，计算每个样本的平均值，计算平均值的百分位数的值

X = np.random.randn(100) # random 1D array
N = 1000 # number of bootstrap samples
idx = np.random.randint(0, X.size, (N, X.size))
means = X[idx].mean(axis=1)
confint = np.percentile(means, [2.5, 97.5])
print(confint)