[PYTHON] My reverse numpy / scipy memo
Reverse numpy memo
numpy? It's ugly. I don't think you'll like it until you get used to it. But I won't lose. Well, num, I want to show you my python with numpy moving. ** * I think this article will be revised little by little. ** **
My minor idiom
The following is a writing style that I rarely see in commentary articles, I'll use it in this article because it's almost an idiom for me.
#Do something based on the gram matrix
# *Is a qualifier to expand when you hit a function
fuga = hoge(*np.meshgrid(x, y))
#Take out the hoge on the x-axis of the matrix A
#Something like a map. For some reason everyone uses a for statement, but isn't it easier to understand if numpy is written like a function?
fuga = np.apply_along_axis(hoge, x, A)
#It's an idiom, it's usually a dot product,
#Why not everyone use this np.Is it just dots? It's nice to write crisply, but ...
hoge @ fuga
Make i
I thought I didn't need this area, but I'm surprised that I'm googled.
np.eye(5)
Just make an arithmetic progression
Arithmetic progression of magnitudes from 0 to 100, with a difference of 1 from each other.
np.arange(0, 100, 1)
However, the last one will be set to 1 if you keep silent.
np.arange(100)
And. Well, it's basic.
Reverse order
np.flip
Linking
The one that connects in parallel or vertically is np.vstack The one that connects in series or horizontally np.hstack But this function is not np.vstack (hoge, fuga) Use it like np.vstack ((hoge, fuga)). There is also np.concatenate.
Make a gram matrix
Make a Gram matrix normally. In such a case, meshgrid is good. … But meshgrid itself makes me want to call it numpy art Am I the only one who has a lot of skill in writing?
xg, yg = np.meshgrid(x, y)
gram = xg * yg
But you can also write like this. In this case, in the sense that it is very similar to the kernel method below. I like this one because it's easy to understand. (The kernel method is a generalized dot product, isn't it?)
from operator import mul
gram = mul(*np.meshgrid(x, y))
Kernel method
For the time being, I wrote it with the Gaussian kernel. All you have to do is replace the mul in the above Gram matrix with the kernel.
def kernel(a: np.ndarray, b: np.ndarray, v: float) -> np.ndarray:
return np.exp(- v * np.square(a - b))
gram = kernel(*np.meshgrid(x, x))
With quadratic regularization
lmd = 0.1 #suitable
gram = kernel(*np.meshgrid(x, x)) + lmd * np.eye(x.shape[0])
Applying that, I used to playfully write kernel regression golf.
def solve_kernel(x: np.ndarray, y: np.ndarray,
kernel: Callable, lmd: float) -> np.ndarray:
K = kernel(*np.meshgrid(x, x)) + lmd * np.eye(x.shape[0])
alpha = np.linalg.solve(K, y)
return np.vectorize(lambda xs: np.sum(alpha * kernel(x, xs)))
Kernel regression can be written in only 4 lines even if PEP8 is protected. numpy is great ...
Try to make a mountain-like thing in three dimensions
Let's make a Gaussian mountain with 100 rows and 100 columns. You can also use meshgrid for this. It's sad that it's a little black magic, but it's not a big deal.
def mount(x: int, y: int) -> np.ndarray:
return np.exp(-(np.square(x) + np.square(y)) / 1000)
x_y = 100, 100
g = mount(*np.meshgrid(*(np.arange(-n/2, n/2, 1) for n in x_y)))
Why did you make this? This is the one I wrote for the clustering experiment. The combination of meshgrid and * may be golden, or it may be beaten for violating coding conventions.
When you don't like for statements
There are multiple ways. Well, it's okay because the map is new.
Use universal functions
You can create a function called universal function that behaves like numpy, so there is one way to use it.
A = [[2, 4, 8], [1, 1, 1]]
def test(x):
return x * 2
np.vectorize(test)(A)
>array([[4, 8, 16], [2, 2, 2]])
By the way, it can also be used as a decorator, but I don't know if it is allowed in the sense of coding conventions. It seems that there are various things such as the number of arguments must be just a number to become a function properly.
Use np.apply_along_axis
This method is quite applicable because the axis can be specified. Moreover, the official says it's fast. I don't know which is faster, vectorize.
#In the same way as above
np.apply_along_axis(test, 0, A)
It seems that there are some dimensional restrictions in order to move properly. There are several other functions. https://numpy.org/devdocs/reference/routines.functional.html
Solve simultaneous equations (if they can be solved)
Ordinary simultaneous equations. For example:
X = np.linalg.solve(A, Y)
In this case, the unknown number was the variable X. However, in practice, most of the time you want to know the constant A from the measured value X. Therefore, in the least squares method, A and X here are in opposite positions. It seems better to use this function without explicitly issuing the inverse matrix.
Solve the least squares method (if there is no solution)
For patterns that do not have a solution, the least squares method is used. "There is no solution" is, for example, the following.
There is no solution because these simultaneous equations are contradictory to each other. Since the solution of the least squares method is the following equation
A = np.linalg.solve(X.T @ X, X.T @ Y)
You can also use the reciprocal.
A = pinv(x) @ y
However, since there is already a function properly, the following is the correct answer instead of the above.
A = np.linalg.lstsq(X, Y)[0]
Least squares with regularization
I'm not sure what to do with regularization.
I googled, but I couldn't find the correct answer, so I wonder if I should write it myself.
lmd = 1
seisoku = lmd * np.eye(min(X.shape))
A = np.linalg.solve(X.T @ X + seisoku , X.T @ Y)
Solve the smallest norm solution (if the solution is infinite)
In such a case, the solution is infinite.
A = X.T @ np.linalg.inv(X @ X.T) @ Y
Isn't it? However, the inverse matrix of Moore Penrose is amazing, so I can go normally.
A = pinv(X) @ Y
However, the correct answer is as follows.
A = np.linalg.lstsq(X, Y)[0]
This function is amazing.
Least squares with regularization (when the solution is infinite)
I don't know how to attack this too. However, the formula is easy to make.
lmd = 1
seisoku = lmd * np.eye(min(X.shape))
A = X.T @ np.linalg.inv(X @ X.T + seisoku) @ Y
Hmm ... (･ ั ω ･ ั) Is it my lack of study that I couldn't find with regularization? It's okay because it's not too difficult ...
Other various optimization problems
https://docs.scipy.org/doc/scipy/reference/optimize.html#module-scipy.optimize
Make it faster
Operators such as `` `**are relatively slow in numpy. You shouldn't doe ** hoge```, it's faster to use np.exp.
That was the case with the Wavelet transform package I wrote. (So it's not readable ...)
hoge = np.arange(1000)
e ** hoge #slow
np.exp(hoge) #fast
hoge ** 2 #slow
np.square(hoge) #fast
hoge ** (0.5) #slow
np.sqart(hoge) #fast
hoge ** (3) #slow
np.power(hoge, 3) #fast
hoge ** (3.4) #slow
np.float_power(hoge, 3.4) #fast
Besides this, there are various things such as np.power, but it is faster to use basic and optimized functions.
Matrix multiplication
I'm not sure why everyone doesn't use the @ operator. I think it's easier to read than np.dot ... Are you late?
np.arange(10)[email protected](10)
Complex number related
Complex conjugate
np.conj
Real axis imaginary axis related
hoge.real, hoge.imag You can write like that, so when substituting this one
fuga.real, fuga.imag = hoge.real, hoge.imag
The reason why I want this is that it is surprising that some processing is done only on the imaginary axis and only on the real axis.
Fourier transform
I think that the Fourier transform often uses scipy. It's also in numpy.
scipy.fftpack.fft(wave)
Inverse conversion
comp #Let it be a complex matrix after Fourier transform
scipy.fftpack.ifft(comp)
cupy related
Convert to cupy
cp.asarray(numpy_array)
Convert to numpy
cp.asnumpy(cupy_array)
There are some differences between cupy itself and numpy, so be careful.
Divided case
Use np.where. If x matches the conditional expression, the second argument is issued, and if not, the third argument is issued.
x = np.arange(10)
np.where(x > 5, x * 10, x)
> array([0, 1, 2, 3, 4, 5, 60, 70, 80, 90])
Conditional statement related
When the conditional expression is inserted, it behaves like R.
x = np.arange(10)
x > 5
> array([False, False, False, False, False, False, True, True, True, True])
Math-ish guy
np.det(hoge) #Determinant
np.norm(hoge) #Norm
Zero padding
hoge is an array Next is to fill in the 3rd and earlier and 5th transitions. constant means to fill with constants. The last is a constant to fill with left and right. There are various other modes.
np.pad(hoge, [3, 5], 'constant' ,[0, 0])
Image hollow
Try to hollow out only the part with a large value. Create an array of True and False as before and hollow out with a conditional expression. In the example, after making something like a mountain with numpy, the top of the mountain is cut off. At this time, np.nan may be needed. This is because matplotlib's imshow ignores nan.
def mount(x: int, y: int) -> np.ndarray:
return np.exp(-(np.square(x) + np.square(y)) / 1000)
#A hollow function.
#A matrix where values is an ordinary matrix and cluster is bool
def make_single_cluster(values, cluster):
img = np.nan * np.zeros_like(values)
img[cluster] = tvalues[cluster]
return img
x_y = 100, 100
#This time I tried adding a little noise
X = np.random.rand(*x_y) / 100
#Make a mountain
X += mount(*np.meshgrid(*(np.arange(-n/2, n/2, 1) for n in x_y)))
plt.imshow(make_single_cluster(X, X > 0.5))
plt.show()
matlab read / write
scipy has matlab read / write function.
import scipy.io as sio
sio.loadmat('hoge.mat')
data = {'fuga', np.arange(100)}
sio.savemat(fuga.mat, data)
Unfortunately I've never used matlab so I'm not sure.
Integral
This is also scipy This is Majikayo.
from scipy.integrate import quad
def test(x):
return x ** 2 - x
print(quad, test, -5, 5)
Another thing is to integrate twice. It's dblquad
Special functions
Well, you like math, right? Look at this link. https://docs.scipy.org/doc/scipy/reference/special.html#module-scipy.special
Recommended Posts