<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
Step2:  Linear inverse solutions in NumPy
Step3:  Mauricio's 1D problem in 2D
We first define the model.
Step4:  Form the discrete kernel, G.
Step5:  Compute the data; this is the forward problem.
Step6:  Minimum least squares solution
Step7:  We can also just use NumPy's solver directly.
Step8:  More complex models
Step9:  I changed the shape of m so we have to make a new G.
Step10:  Forward model the data.
Step11:  Noise-free
Step12:  Solve with LAPACK
Step13:  With noise
Step14:  With noise
Step15:  Now we solve
<ASSISTANT_TASK:>
Python Code:
import numpy as np
import numpy.linalg as la
import matplotlib.pyplot as plt
import matplotlib.image as mpimg
%matplotlib inline
def plot_all(m, d, m_est, d_pred):
    
    Helper function for plotting. You can ignore this.
    
    fig = plt.figure(figsize=(10,6))
    ax0 = fig.add_subplot(2,2,1)
    ax0.imshow(m, cmap='viridis', aspect='auto')
    ax0.set_title("$\mathrm{Model}\ m$")
    ax1 = fig.add_subplot(2,2,2)
    ax1.imshow(d, cmap='viridis', aspect='auto')
    ax1.set_title("$\mathrm{Data}\ d$")
    ax2 = fig.add_subplot(2,2,3)
    ax2.imshow(m_est, cmap='viridis', aspect='auto')
    ax2.set_title("$\mathrm{Estimated\ model}\ m_\mathrm{est}$")
    ax3 = fig.add_subplot(2,2,4)
    ax3.imshow(d_pred, cmap='viridis', aspect='auto')
    ax3.set_title("$\mathrm{Predicted\ data}\ d_\mathrm{pred}$")
    plt.show()
import scipy.linalg
def convmtx(h, n):
    
    Equivalent of MATLAB's convmtx function, http://www.mathworks.com/help/signal/ref/convmtx.html.
    
    Makes the convolution matrix, C. The product C.x is the convolution of h and x.
    
    Args
        h (ndarray): a 1D array, the kernel.
        n (int): the number of rows to make.
        
    Returns
        ndarray. Size m+n-1
    
    col_1 = np.r_[h[0], np.zeros(n-1)]
    row_1 = np.r_[h, np.zeros(n-1)]
    return scipy.linalg.toeplitz(col_1, row_1)
Explanation: Linear inverse solutions in NumPy: 2D
Reproducing Mauricio Sacchi's GEOPH 431 and 531 tutorial on linear inversion, which was written in MATLAB. Mauricio's version is also recorded at SubSurfWiki on the Linear inversion chrestomathy page.
This version of the notebook solves the same problem, but in 2D. It then goes on to look at a different model, also in 2D. The math is all exactly the same.
NB For the @ notation to work (matrix mulitplication) you need Python 3.5 and Numpy 1.11.
First, the usual preliminaries.
End of explanation
M = 50
m = np.zeros((M+1, 1))  # Have to do +1 because we're going to lose one in computing RC series.
m[10:15,:] = 1.0
m[15:27,:] = -0.3
m[27:35,:] = 2.1
m = (m[1:] - m[:-1]) / (m[1:] + m[:-1] + 1e-9)  # Small number avoid division by zero.
m = np.repeat(m, 50, axis=-1)
plt.imshow(m, cmap='viridis')
Explanation: Mauricio's 1D problem in 2D
We first define the model.
End of explanation
N = 20
L = 100
alpha = 0.08
x = np.arange(0, M, 1) * L/(M-1)
dx = L/(M-1)
r = np.arange(0, N, 1) * L/(N-1)
G = np.zeros((N, M))
for j in range(M):
    for k in range(N):
        G[k,j] = dx * np.exp(-alpha * np.abs(r[k] - x[j])**2)
Explanation: Form the discrete kernel, G.
End of explanation
d = G @ m
plt.imshow(d, cmap='viridis')
Explanation: Compute the data; this is the forward problem.
End of explanation
m_est = G.T @ la.inv(G @ G.T) @ d
d_pred = G @ m_est
plot_all(m, d, m_est, d_pred)
Explanation: Minimum least squares solution:
End of explanation
m_est = la.lstsq(G, d)[0]
d_pred = G @ m_est
plot_all(m, d, m_est, d_pred)
Explanation: We can also just use NumPy's solver directly.
End of explanation
from io import BytesIO
import requests
from urllib.parse import quote
text = '+ ='
url = "https://chart.googleapis.com/chart"
params = {'chst': 'd_text_outline',
          'chld': '000000|36|h|000000|_|{}'.format(text),
         }
r = requests.get(url, params)
b = BytesIO(r.content)
img = mpimg.imread(b)
m = np.pad(img[...,3], 20, 'constant')
plt.imshow(m, cmap='viridis')
m = (m[1:] - m[:-1]) / (m[1:] + m[:-1] + 1e-9)  # Small number avoid division by zero.
Explanation: More complex models
End of explanation
M = m.shape[0]
N = 20
L = 100
alpha = 0.08
x = np.arange(0, M, 1) * L/(M-1)
dx = L/(M-1)
r = np.arange(0, N, 1) * L/(N-1)
G = np.zeros((N, M))
for j in range(M):
    for k in range(N):
        G[k,j] = dx * np.exp(-alpha * np.abs(r[k] - x[j])**2)
Explanation: I changed the shape of m so we have to make a new G.
End of explanation
d = G @ m
plt.imshow(d, cmap='viridis')
Explanation: Forward model the data.
End of explanation
m_est = G.T @ la.inv(G @ G.T) @ d
d_pred = G @ m_est
plot_all(m, d, m_est, d_pred)
Explanation: Noise-free: minimum norm
End of explanation
m_est = la.lstsq(G, d)[0]
d_pred = G @ m_est
plot_all(m, d, m_est, d_pred)
Explanation: Solve with LAPACK
End of explanation
s = 1
d += s * np.random.random(d.shape)
I = np.eye(N)
µ = 2.5
m_est = G.T @ la.inv(G @ G.T + µ * I) @ d
d_pred = G @ m_est
plot_all(m, d, m_est, d_pred)
Explanation: With noise: damped least squares
End of explanation
W = convmtx([1,-1], M)[:,:-1]  # Skip last column
Explanation: With noise: damped least squares with first derivative regularization
End of explanation
m_est = la.inv(G.T @ G + µ * W.T @ W) @ G.T @ d
d_pred = G @ m_est
plot_all(m, d, m_est, d_pred)
Explanation: Now we solve:
$$ \hat{\mathbf{m}} = (\mathbf{G}^\mathrm{T} \mathbf{G} + \mu \mathbf{W}^\mathrm{T} \mathbf{W})^{-1} \mathbf{G}^\mathrm{T} \mathbf{d} \ \ $$
End of explanation
<END_TASK>