<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Sparse Linear Inverse Demo
In this demo, we illustrate how to use the vampyre package for a simple sparse linear inverse problem.   The problem is to estimate a sparse vector $z$ from linear measurements of the form $y=Az+w$ where $w$ is Gaussian noise and $A$ is a known linear transform -- a basic problem in compressed sensing.  By sparse, we mean that the vector $z$ has few non-zero values.  Knowing that the vector is sparse can be used for improved reconstruction if an appropriate sparse reconstruction algorithm is used.
There are a large number of algorithms for sparse linear inverse problems.  This demo uses the Vector Approximate Message Passing (VAMP) method, one of several methods that will be included in the vampyre package.  In going through this demo, you will learn to
Step1:  We will also load the other packages we will use in this demo.  This could be done before the above import.
Step2:  Generating Synthetic Data
We begin by generating synthetic data $z$ and measurements $y$ that we will use to test the algorithms. First, we set the dimensions and the shapes of the vectors we will use.
Step3:  To generate the synthetic data for this demo, we use the following simple probabilistic model.  For the input $z$, we will use Bernouli-Gaussian (BG) distribution, a simple  model in sparse signal processing.  In the BG model, the components $z_i$ are i.i.d. where each component $z_i=0$ with probability $1-\rho$ and $z_i \sim {\mathcal N}(0,1)$ with probability $\rho$.  The parameter $\rho$ is called the sparsity ratio and represents the average number of non-zero components.  When $\rho$ is small, the vector $z$ is sparse.  The components on which $z_i$ are non-zero are called the active components.  We set the parameters below.  We also set the SNR for the measurements.
Step4:  Using these parameters, we can generate random sparse z following this distribution with the following simple code.
Step5:  To illustrate the sparsity, we plot the vector z.  We can see from this plot that the majority of the components of z are zero.
Step6:  Now, we create a random transform $A$ and output $y_0 = Az$.
Step7:  Finally, we add noise at the desired SNR
Step8:  Creating the Vampyre estimators
Now that we have created the sparse data, we will use the vampyre package to recover z from y.  In vampyre the methods to perform this estimation are called solvers.  For this demo, we will use a simple solver called VAMP described in the paper
Step9:  We next use the vampyre class, MixEst, to describe a mixture of the two distributions.  This is done by creating a list, est_list, of the estimators and an array pz with the probability of each component.  The resulting estimator, est_in, is the estimator for the prior $z$, which is also the input to the transform $A$.  We give this a name Input since it corresponds to the input.  But, any naming is fine.  Or, you can let vampyre give it a generic name.
Step10:  Next, we describe the likelihood function, $p(y|z)$.  Since $y=Az+w$, we can first use the MatrixLT class to define a linear transform operator Aop corresponding to the matrix A.  Then, we use the LinEstim class to describe the likelihood $y=Az+w$.
Step11:  Finally, the VAMP method needs a message handler to describe how to perform the Gaussian message passing.  This is a more advanced feature.  For most applications, you can just use the simple message handler as follows.
Step12:  Running the VAMP Solver
Having described the input and output estimators and the variance handler, we can now construct a VAMP solver.   The construtor takes the input and output estimators, the variance handler and other parameters. The paramter nit is the number of iterations.  This is fixed for now.  Later, we will add auto-termination.  The other parameter, hist_list is optional, and will be described momentarily.
Step13:  We can print a summary of the model which indicates the dimensions and the estimators.
Step14:  We now run the solver by calling the solve() method.  For a small problem like this, this should be close to instantaneous.
Step15:  The VAMP solver estimate is the field zhat.  We plot one column of this (icol=0) and compare it to the corresponding column of the true matrix z.  You should see a very good match.
Step16:  We can measure the normalized mean squared error as follows.  The VAMP solver also produces an estimate of the MSE in the variable zhatvar.  We can extract this variable to compute the predicted MSE.  We see that the normalized MSE is indeed low and closely matches the predicted value from VAMP.
Step17:  Finally, we can plot the actual and predicted MSE as a function of the iteration number.  When solver was contructed, we passed an argument hist_list=['zhat', 'zhatvar'].  This indicated to store the value of the estimate zhat and predicted error variance zhatvar with each iteration.  We can recover these values from solver.hist_dict, the history dictionary.  Using the values we can compute and plot the normalized MSE on each iteartion.  We see that VAMP gets a low MSE in very few iterations, about 10.
<ASSISTANT_TASK:>
Python Code:
import os
import sys
vp_path = os.path.abspath('../../')
if not vp_path in sys.path:
    sys.path.append(vp_path)
import vampyre as vp
Explanation: Sparse Linear Inverse Demo
In this demo, we illustrate how to use the vampyre package for a simple sparse linear inverse problem.   The problem is to estimate a sparse vector $z$ from linear measurements of the form $y=Az+w$ where $w$ is Gaussian noise and $A$ is a known linear transform -- a basic problem in compressed sensing.  By sparse, we mean that the vector $z$ has few non-zero values.  Knowing that the vector is sparse can be used for improved reconstruction if an appropriate sparse reconstruction algorithm is used.
There are a large number of algorithms for sparse linear inverse problems.  This demo uses the Vector Approximate Message Passing (VAMP) method, one of several methods that will be included in the vampyre package.  In going through this demo, you will learn to:
* Load the vampyre package
* Create synthetic data for a sparse linear inverse problem
* Set up the VAMP method in the vampyre package to perform the estimation for the linear inverse problem
* Measure the mean squared error (MSE) and compare the value to the predicted value from the VAMP method.
* Using the hist_list feature to track variables per iteration of the algorithm.
Importing the Package
First we need to import the vampyre package.  Since python does not have relative imports, you need to add the path location for the vampyre package to the system path.  In this case, we have specified the path use a relative path location, but you can change this depending on where vampyre is located.
End of explanation
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
%matplotlib inline
Explanation: We will also load the other packages we will use in this demo.  This could be done before the above import.
End of explanation
# Parameters
nz = 1000   # number of components of z
ny = 500    # number of measurements y 
# Compute the shapes
zshape = (nz,)   # Shape of z matrix
yshape = (ny,)   # Shape of y matrix
Ashape = (ny,nz)   # Shape of A matrix
Explanation: Generating Synthetic Data
We begin by generating synthetic data $z$ and measurements $y$ that we will use to test the algorithms. First, we set the dimensions and the shapes of the vectors we will use.
End of explanation
sparse_rat = 0.1   # sparsity ratio
zmean1 = 0         # mean for the active components
zvar1 = 1          # variance for the active components
snr = 30           # SNR in dB
Explanation: To generate the synthetic data for this demo, we use the following simple probabilistic model.  For the input $z$, we will use Bernouli-Gaussian (BG) distribution, a simple  model in sparse signal processing.  In the BG model, the components $z_i$ are i.i.d. where each component $z_i=0$ with probability $1-\rho$ and $z_i \sim {\mathcal N}(0,1)$ with probability $\rho$.  The parameter $\rho$ is called the sparsity ratio and represents the average number of non-zero components.  When $\rho$ is small, the vector $z$ is sparse.  The components on which $z_i$ are non-zero are called the active components.  We set the parameters below.  We also set the SNR for the measurements.
End of explanation
# Generate the random input
z1 = np.random.normal(zmean1, np.sqrt(zvar1), zshape)
u = np.random.uniform(0, 1, zshape) < sparse_rat
z = z1*u
Explanation: Using these parameters, we can generate random sparse z following this distribution with the following simple code.
End of explanation
ind = np.array(range(nz))
plt.plot(ind,z)
Explanation: To illustrate the sparsity, we plot the vector z.  We can see from this plot that the majority of the components of z are zero.
End of explanation
A = np.random.normal(0, 1/np.sqrt(nz), Ashape)
y0 = A.dot(z)
Explanation: Now, we create a random transform $A$ and output $y_0 = Az$.
End of explanation
yvar = np.mean(np.abs(y0)**2)
wvar = yvar*np.power(10, -0.1*snr)
y = y0 + np.random.normal(0,np.sqrt(wvar), yshape)
Explanation: Finally, we add noise at the desired SNR
End of explanation
est0 = vp.estim.DiscreteEst(0,1,zshape)
est1 = vp.estim.GaussEst(zmean1,zvar1,zshape)
Explanation: Creating the Vampyre estimators
Now that we have created the sparse data, we will use the vampyre package to recover z from y.  In vampyre the methods to perform this estimation are called solvers.  For this demo, we will use a simple solver called VAMP described in the paper:
Rangan, Sundeep, Philip Schniter, and Alyson Fletcher. "Vector approximate message passing." arXiv preprint arXiv:1610.03082 (2016).
Similar to most of the solvers in the vampyre package, the VAMP solver needs precise specifications of the probability distributions of z and y.  The simplest way to use VAMP is to specify two densities:
* The prior $p(z)$; and
* The likelihood $p(y|z)$.
Each of the densities are described by estimators.  
We first describe the estimator for the prior $p(z)$.  The vampyre package will eventually have a large number of estimators to describe various densities.  In this simple demo, $p(z)$ is what is called a mixture distribution since $z$ is one distribution with probability $1-\rho$ and a second distribution with probability $\rho$.  To describe this mixture distribution in the vampyre package, we need to first create estimator classes for each component distribution.  The following code creates an estimator, est0, for a discrete distribution with a probability of 1 at a 0 and a second estimator, est1, for the Gaussian distribution with the active components.
End of explanation
est_list = [est0, est1]
pz = np.array([1-sparse_rat, sparse_rat])
est_in = vp.estim.MixEst(est_list, w=pz, name='Input')
Explanation: We next use the vampyre class, MixEst, to describe a mixture of the two distributions.  This is done by creating a list, est_list, of the estimators and an array pz with the probability of each component.  The resulting estimator, est_in, is the estimator for the prior $z$, which is also the input to the transform $A$.  We give this a name Input since it corresponds to the input.  But, any naming is fine.  Or, you can let vampyre give it a generic name.
End of explanation
Aop = vp.trans.MatrixLT(A,zshape)
est_out = vp.estim.LinEst(Aop,y,wvar,map_est=False, name='Output')
Explanation: Next, we describe the likelihood function, $p(y|z)$.  Since $y=Az+w$, we can first use the MatrixLT class to define a linear transform operator Aop corresponding to the matrix A.  Then, we use the LinEstim class to describe the likelihood $y=Az+w$.
End of explanation
msg_hdl = vp.estim.MsgHdlSimp(map_est=False, shape=zshape)
Explanation: Finally, the VAMP method needs a message handler to describe how to perform the Gaussian message passing.  This is a more advanced feature.  For most applications, you can just use the simple message handler as follows.
End of explanation
nit = 20  # number of iterations
solver = vp.solver.Vamp(est_in,est_out,msg_hdl,\
                        hist_list=['zhat', 'zhatvar'],nit=nit)
Explanation: Running the VAMP Solver
Having described the input and output estimators and the variance handler, we can now construct a VAMP solver.   The construtor takes the input and output estimators, the variance handler and other parameters. The paramter nit is the number of iterations.  This is fixed for now.  Later, we will add auto-termination.  The other parameter, hist_list is optional, and will be described momentarily.
End of explanation
solver.summary()
Explanation: We can print a summary of the model which indicates the dimensions and the estimators.
End of explanation
solver.solve()
Explanation: We now run the solver by calling the solve() method.  For a small problem like this, this should be close to instantaneous.
End of explanation
zhat = solver.zhat
ind = np.array(range(nz))
plt.plot(ind,z)
plt.plot(ind,zhat)
plt.legend(['True', 'Estimate'])
Explanation: The VAMP solver estimate is the field zhat.  We plot one column of this (icol=0) and compare it to the corresponding column of the true matrix z.  You should see a very good match.
End of explanation
zerr = np.mean(np.abs(zhat-z)**2)
zhatvar = solver.zhatvar
zpow = np.mean(np.abs(z)**2)
mse_act = 10*np.log10(zerr/zpow)
mse_pred = 10*np.log10(zhatvar/zpow)
print("Normalized MSE (dB): actual {0:f} pred {1:f}".format(mse_act, mse_pred))
Explanation: We can measure the normalized mean squared error as follows.  The VAMP solver also produces an estimate of the MSE in the variable zhatvar.  We can extract this variable to compute the predicted MSE.  We see that the normalized MSE is indeed low and closely matches the predicted value from VAMP.
End of explanation
# Compute the MSE as a function of the iteration
zhat_hist = solver.hist_dict['zhat']
zhatvar_hist = solver.hist_dict['zhatvar']
nit = len(zhat_hist)
mse_act = np.zeros(nit)
mse_pred = np.zeros(nit)
for it in range(nit):
    zerr = np.mean(np.abs(zhat_hist[it]-z)**2)
    mse_act[it] = 10*np.log10(zerr/zpow)
    mse_pred[it] = 10*np.log10(zhatvar_hist[it]/zpow)
    
plt.plot(range(nit), mse_act, 'o-', linewidth=2)
plt.plot(range(nit), mse_pred, 's', linewidth=1)
plt.xlabel('Iteration')
plt.ylabel('Normalized MSE (dB)')
plt.legend(['Actual', 'Predicted'])
plt.grid()
Explanation: Finally, we can plot the actual and predicted MSE as a function of the iteration number.  When solver was contructed, we passed an argument hist_list=['zhat', 'zhatvar'].  This indicated to store the value of the estimate zhat and predicted error variance zhatvar with each iteration.  We can recover these values from solver.hist_dict, the history dictionary.  Using the values we can compute and plot the normalized MSE on each iteartion.  We see that VAMP gets a low MSE in very few iterations, about 10.
End of explanation
<END_TASK>