Added LSQR algorithm - #1975

# Load components needed
from cil.framework import ImageGeometry, AcquisitionGeometry, ImageData
from cil.optimisation.algorithms import LSQR
from cil.plugins.astra import ProjectionOperator
from cil.utilities.display import show2D

# Load third-party inputs (using phantominator because the TomoPhantom plug-in does not work on mac for me)
import numpy as np    
import matplotlib.pyplot as plt
from phantominator import shepp_logan

# Set up toy problem
# Geometry
n = 256
ig = ImageGeometry(voxel_num_x=n, 
                   voxel_num_y=n, 
                   voxel_size_x=2/n, 
                   voxel_size_y=2/n)

num_angles = 180
ag = AcquisitionGeometry.create_Parallel2D()  \
                   .set_angles(np.linspace(0, 180, num_angles, endpoint=False))  \
                   .set_panel(n, 2/n)

# Shepp-Logan phantom
phantom2DArray = shepp_logan(n)
phantom2D = ImageData(array = phantom2DArray, deep_copy=False, geometry=ig)

# Load operator
device = "cpu"
A = ProjectionOperator(ig, ag, device)

# Create noise-less sinogram
sinogram = A.direct(phantom2D)

# Add noise to sinogram
background_counts = 10000 
counts = background_counts * np.exp(-sinogram.as_array())
noisy_counts = np.random.poisson(counts)
sino_out = -np.log(noisy_counts/background_counts)
sinogram_noisy = ag.allocate()
sinogram_noisy.fill(sino_out)

# Solve problem with LSQR
maxit = 25
initial = ig.allocate(0)
lsqr_simple = LSQR(initial=initial, operator=A, data=sinogram_noisy)
lsqr_simple.run(maxit, verbose=True)

# Plot the solution
plots = [lsqr_simple.solution, lsqr_simple.solution - phantom2D]
titles = ["LSQR solution","Difference from ground truth" ]
show2D(plots, titles, fix_range=[(-0.2,1.2),(-0.2,0.2)])

# Include Tikhonov regularisation without building a block operator
alpha = 0.1
lsqr_reg = LSQR(initial=initial, operator=A, data=sinogram_noisy, alpha = alpha)
lsqr_reg.run(maxit, verbose=True)

# Plot the solution
plots = [lsqr_reg.solution, lsqr_reg.solution - phantom2D]
titles = ["Regularised LSQR solution","Difference from ground truth" ]
show2D(plots, titles, fix_range=[(-0.2,1.2),(-0.2,0.2)])