In [1]:
# Author : Evangelos Papoutsellis
# Contact : http://epapoutsellis.github.io/
In [2]:
# Import libraries
from cil.framework import AcquisitionGeometry, VectorData
from cil.optimisation.functions import L2NormSquared, ZeroFunction, LeastSquares
from cil.optimisation.operators import MatrixOperator
from cil.optimisation.algorithms import FISTA, ISTA
from cil.plugins.astra.operators import ProjectionOperator
from cil.utilities.display import show2D
from cil.utilities import dataexample
from cil.plugins.astra.utilities import convert_geometry_to_astra
from cil.optimisation.functions import TotalVariation
from cil.plugins import TomoPhantom
import matplotlib.pyplot as plt
import cvxpy
import astra
import numpy as np
In [3]:
# Detectors
N = 64
detectors = N
# Angles
angles = np.linspace(0,180,90, dtype='float32')
# Setup acquisition geometry
ag = AcquisitionGeometry.create_Parallel2D()\
.set_angles(angles)\
.set_panel(detectors)
# Get image geometry
ig = ag.get_ImageGeometry()
# Create projection operator using Astra-Toolbox. Available CPU/CPU
A = ProjectionOperator(ig, ag, device = 'cpu')
# Get phantom
# phantom = dataexample.SIMPLE_PHANTOM_2D.get(size=(N, N))
phantom = TomoPhantom.get_ImageData(12, ig)
# Create an acqusition data
sino = A.direct(phantom)
# Simulate Gaussian noise for the sinogram
gaussian_var = 2.
gaussian_mean = 0
np.random.seed(10)
n1 = np.random.normal(gaussian_mean, gaussian_var, size = ag.shape)
noisy_sino = ag.allocate()
noisy_sino.fill(n1 + sino.array)
noisy_sino.array[noisy_sino.array<0]=0
# noisy_sino.fill(sino)
# Show numerical and noisy sinograms
show2D([phantom, sino, noisy_sino],
title = ['Ground Truth','Sinogram','Noisy Sinogram'],
num_cols=3, cmap = 'inferno')
show2D(A.adjoint(noisy_sino), cmap="inferno")
/opt/anaconda3/envs/cil_devel_epaps/lib/python3.7/site-packages/cil/utilities/display.py:358: UserWarning: FixedFormatter should only be used together with FixedLocator axes[i].set_yticklabels(labels_new)
Out[3]:
<cil.utilities.display.show2D at 0x7f4a0c746e90>
In [4]:
# convert to astra geometries
ig_astra, ag_astra = convert_geometry_to_astra(ig, ag)
# projection id
proj_id = astra.create_projector('line', ag_astra, ig_astra)
# matrix id
matrix_id = astra.projector.matrix(proj_id)
# Get the projection matrix as a Scipy sparse matrix.
W = astra.matrix.get(matrix_id)
sino_astra = (W*phantom.array.flatten()).reshape((len(ag.angles), N))
#show
show2D([sino_astra, sino.array, np.abs(sino.array-sino_astra)], num_cols=3, cmap="inferno")
# print("Norm of A_cil = {}".format(A.norm()))
mat_cil = MatrixOperator(W.T*W)
# print("Norm of astra mat = {}".format(np.sqrt(mat_cil.norm())))
In [5]:
alpha = 50.
In [6]:
N, M = ig.shape
u_cvx = cvxpy.Variable(N*M)
fidelity = cvxpy.sum_squares(W@u_cvx.flatten() - noisy_sino.array.flatten())
constraints = [u_cvx>=0]
regulariser = (alpha/ig.voxel_size_x) * cvxpy.tv(cvxpy.reshape(u_cvx, (N,M)))
obj = cvxpy.Minimize( fidelity + regulariser)
prob = cvxpy.Problem(obj, constraints)
res = prob.solve(solver = cvxpy.MOSEK, verbose = True)
# res = prob.solve(solver = SCS, verbose = True, eps=1e-3)
===============================================================================
CVXPY
v1.1.15
===============================================================================
(CVXPY) May 15 04:22:11 PM: Your problem has 4096 variables, 1 constraints, and 0 parameters.
(CVXPY) May 15 04:22:11 PM: It is compliant with the following grammars: DCP, DQCP
(CVXPY) May 15 04:22:11 PM: (If you need to solve this problem multiple times, but with different data, consider using parameters.)
(CVXPY) May 15 04:22:11 PM: CVXPY will first compile your problem; then, it will invoke a numerical solver to obtain a solution.
-------------------------------------------------------------------------------
Compilation
-------------------------------------------------------------------------------
(CVXPY) May 15 04:22:11 PM: Compiling problem (target solver=MOSEK).
(CVXPY) May 15 04:22:11 PM: Reduction chain: Dcp2Cone -> CvxAttr2Constr -> ConeMatrixStuffing -> MOSEK
(CVXPY) May 15 04:22:11 PM: Applying reduction Dcp2Cone
(CVXPY) May 15 04:22:11 PM: Applying reduction CvxAttr2Constr
(CVXPY) May 15 04:22:11 PM: Applying reduction ConeMatrixStuffing
(CVXPY) May 15 04:22:12 PM: Applying reduction MOSEK
(CVXPY) May 15 04:22:12 PM: Finished problem compilation (took 6.408e-01 seconds).
-------------------------------------------------------------------------------
Numerical solver
-------------------------------------------------------------------------------
(CVXPY) May 15 04:22:12 PM: Invoking solver MOSEK to obtain a solution.
Problem
Name :
Objective sense : max
Type : CONIC (conic optimization problem)
Constraints : 8066
Cones : 3970
Scalar variables : 21765
Matrix variables : 0
Integer variables : 0
Optimizer started.
Problem
Name :
Objective sense : max
Type : CONIC (conic optimization problem)
Constraints : 8066
Cones : 3970
Scalar variables : 21765
Matrix variables : 0
Integer variables : 0
Optimizer - threads : 20
Optimizer - solved problem : the primal
Optimizer - Constraints : 4097
Optimizer - Cones : 3970
Optimizer - Scalar variables : 21765 conic : 17669
Optimizer - Semi-definite variables: 0 scalarized : 0
Factor - setup time : 1.81 dense det. time : 0.00
Factor - ML order time : 0.69 GP order time : 0.00
Factor - nonzeros before factor : 7.62e+06 after factor : 8.39e+06
Factor - dense dim. : 2 flops : 2.30e+10
ITE PFEAS DFEAS GFEAS PRSTATUS POBJ DOBJ MU TIME
0 4.9e+01 6.6e+01 2.0e+05 0.00e+00 -1.000000000e+00 1.984500000e+05 1.0e+00 1.97
1 2.2e+01 3.0e+01 1.3e+05 -9.59e-01 -2.766915462e+03 1.908257894e+05 4.5e-01 2.27
2 1.5e+01 2.0e+01 1.0e+05 -9.10e-01 -5.180417132e+03 1.838848710e+05 3.0e-01 2.51
3 1.3e+01 1.7e+01 9.5e+04 -8.65e-01 -3.512507478e+03 1.840076679e+05 2.6e-01 2.71
4 7.7e+00 1.0e+01 6.9e+04 -8.46e-01 -1.046589718e+04 1.672586300e+05 1.6e-01 2.91
5 6.0e+00 8.1e+00 5.9e+04 -7.45e-01 -4.465478577e+03 1.675816914e+05 1.2e-01 3.11
6 3.3e+00 4.4e+00 3.7e+04 -6.83e-01 -1.106457642e+04 1.392270367e+05 6.7e-02 3.31
7 1.9e+00 2.5e+00 2.3e+04 -4.65e-01 -3.390794239e+02 1.261191156e+05 3.9e-02 3.50
8 6.2e-01 8.3e-01 6.8e+03 -1.78e-01 5.564239659e+03 7.032659851e+04 1.3e-02 3.69
9 8.0e-02 1.1e-01 3.6e+02 4.82e-01 3.706352281e+04 4.742101727e+04 1.6e-03 3.99
10 7.2e-02 9.7e-02 3.1e+02 1.02e+00 3.760792191e+04 4.688409187e+04 1.5e-03 4.19
11 4.8e-02 6.4e-02 1.6e+02 1.08e+00 4.080044122e+04 4.673216701e+04 9.7e-04 4.40
12 3.0e-02 4.1e-02 8.4e+01 1.02e+00 4.177797985e+04 4.552375390e+04 6.2e-04 4.58
13 7.9e-03 1.1e-02 1.1e+01 1.05e+00 4.411641083e+04 4.506652203e+04 1.6e-04 4.78
14 3.6e-04 4.8e-04 9.1e-02 1.01e+00 4.481875539e+04 4.486143967e+04 7.3e-06 5.09
15 2.4e-04 9.6e-01 5.2e-02 1.00e+00 4.483175010e+04 4.486056928e+04 4.9e-06 5.32
16 4.1e-05 1.6e-01 4.4e-03 1.00e+00 4.485298708e+04 4.485793316e+04 8.4e-07 5.65
17 2.3e-05 3.1e+00 1.8e-03 1.00e+00 4.485483002e+04 4.485753225e+04 4.6e-07 5.90
18 3.6e-06 8.1e-01 1.1e-04 1.00e+00 4.485670704e+04 4.485713350e+04 7.3e-08 6.19
19 2.4e-06 8.7e+00 6.2e-05 1.00e+00 4.485683101e+04 4.485711520e+04 4.9e-08 6.42
20 3.7e-07 2.6e+00 3.8e-06 1.00e+00 4.485703960e+04 4.485708339e+04 7.5e-09 6.68
21 6.3e-08 1.5e-01 3.2e-07 1.00e+00 4.485707007e+04 4.485707886e+04 1.3e-09 6.93
Optimizer terminated. Time: 8.10
Interior-point solution summary
Problem status : PRIMAL_AND_DUAL_FEASIBLE
Solution status : OPTIMAL
Primal. obj: 4.4857070070e+04 nrm: 1e+04 Viol. con: 4e-06 var: 0e+00 cones: 0e+00
Dual. obj: 4.4857078862e+04 nrm: 2e+04 Viol. con: 0e+00 var: 9e-08 cones: 4e-11
-------------------------------------------------------------------------------
Summary
-------------------------------------------------------------------------------
(CVXPY) May 15 04:22:20 PM: Problem status: optimal
(CVXPY) May 15 04:22:20 PM: Optimal value: 4.486e+04
(CVXPY) May 15 04:22:20 PM: Compilation took 6.408e-01 seconds
(CVXPY) May 15 04:22:20 PM: Solver (including time spent in interface) took 8.286e+00 seconds
ISTA algorithm¶
In [7]:
initial = ig.allocate()
f = LeastSquares(A, noisy_sino)
g = alpha * TotalVariation(max_iteration=50, tolerance=0, lower=0)
ista = ISTA(initial = initial, f = f, g = g, max_iteration=500, update_objective_interval=1)
ista.run(verbose=0)
FISTA algorithm¶
In [8]:
fista = FISTA(initial = initial, f = f, g = g, max_iteration=500, update_objective_interval=1)
fista.run(verbose=0)
In [9]:
plt.figure()
plt.semilogy(np.abs(fista.objective - obj.value), label="FISTA")
plt.semilogy(np.abs(ista.objective - obj.value), label="ISTA")
plt.legend()
plt.show()
In [10]:
print("Objective value of CVXpy = {}".format(obj.value))
print("Objective value of ISTA = {}".format(ista.objective[-1]))
print("Objective value of FISTA = {}".format(fista.objective[-1]))
# assert ista/fista vs cvxpy
cvx_sol = u_cvx.value.reshape(N,N)
show2D([ista.solution, fista.solution, cvx_sol,
np.abs(ista.solution.array-cvx_sol), np.abs(fista.solution.array-cvx_sol)],
cmap="inferno", origin="upper",
num_cols=3, title=['ISTA','FISTA','CVX',
'Diff: ISTA-CVX', 'Diff: FISTA-CVX'])
Objective value of CVXpy = 44857.07859587796 Objective value of ISTA = 44866.927001953125 Objective value of FISTA = 44865.864501953125
Out[10]:
<cil.utilities.display.show2D at 0x7f4a1407b950>
Solve liner system Ax=b (LeastSquares) with ISTA, FISTA, CVX¶
In [11]:
np.random.seed(10)
n = 100
m = 1000
A = np.random.uniform(0,1, (m, n)).astype('float32')
b = (A.dot(np.random.randn(n)) + 0.1*np.random.randn(m)).astype('float32')
In [12]:
Aop = MatrixOperator(A)
bop = VectorData(b)
f = LeastSquares(Aop, b=bop,c=0.5)
g = ZeroFunction()
ig = Aop.domain
initial = ig.allocate()
ista = ISTA(initial = initial, f = f, g = g, max_iteration=1000, update_objective_interval=1)
ista.run(verbose=0)
fista = FISTA(initial = initial, f = f, g = g, max_iteration=1000, update_objective_interval=1)
fista.run(verbose=0)
In [13]:
z = cvxpy.Variable(n)
objective = cvxpy.Minimize(0.5 * cvxpy.sum_squares(A @ z - b))
constraint = []
p = cvxpy.Problem(objective, constraint)
p.solve(verbose=True, solver=cvxpy.MOSEK)
===============================================================================
CVXPY
v1.1.15
===============================================================================
(CVXPY) May 15 04:23:20 PM: Your problem has 100 variables, 0 constraints, and 0 parameters.
(CVXPY) May 15 04:23:20 PM: It is compliant with the following grammars: DCP, DQCP
(CVXPY) May 15 04:23:20 PM: (If you need to solve this problem multiple times, but with different data, consider using parameters.)
(CVXPY) May 15 04:23:20 PM: CVXPY will first compile your problem; then, it will invoke a numerical solver to obtain a solution.
-------------------------------------------------------------------------------
Compilation
-------------------------------------------------------------------------------
(CVXPY) May 15 04:23:20 PM: Compiling problem (target solver=MOSEK).
(CVXPY) May 15 04:23:20 PM: Reduction chain: Dcp2Cone -> CvxAttr2Constr -> ConeMatrixStuffing -> MOSEK
(CVXPY) May 15 04:23:20 PM: Applying reduction Dcp2Cone
(CVXPY) May 15 04:23:20 PM: Applying reduction CvxAttr2Constr
(CVXPY) May 15 04:23:20 PM: Applying reduction ConeMatrixStuffing
(CVXPY) May 15 04:23:20 PM: Applying reduction MOSEK
(CVXPY) May 15 04:23:20 PM: Finished problem compilation (took 9.323e-02 seconds).
-------------------------------------------------------------------------------
Numerical solver
-------------------------------------------------------------------------------
(CVXPY) May 15 04:23:20 PM: Invoking solver MOSEK to obtain a solution.
Problem
Name :
Objective sense : max
Type : CONIC (conic optimization problem)
Constraints : 101
Cones : 1
Scalar variables : 1002
Matrix variables : 0
Integer variables : 0
Optimizer started.
Problem
Name :
Objective sense : max
Type : CONIC (conic optimization problem)
Constraints : 101
Cones : 1
Scalar variables : 1002
Matrix variables : 0
Integer variables : 0
Optimizer - threads : 20
Optimizer - solved problem : the primal
Optimizer - Constraints : 101
Optimizer - Cones : 1
Optimizer - Scalar variables : 1002 conic : 1002
Optimizer - Semi-definite variables: 0 scalarized : 0
Factor - setup time : 0.01 dense det. time : 0.00
Factor - ML order time : 0.00 GP order time : 0.00
Factor - nonzeros before factor : 5151 after factor : 5151
Factor - dense dim. : 0 flops : 1.05e+07
ITE PFEAS DFEAS GFEAS PRSTATUS POBJ DOBJ MU TIME
0 5.0e-01 2.7e+01 2.0e+00 0.00e+00 -1.000000000e+00 0.000000000e+00 1.0e+00 0.06
1 4.1e-02 2.2e+00 2.6e-01 -8.82e-01 4.040740509e+00 2.474629677e+00 8.1e-02 0.07
2 8.0e-04 4.3e-02 2.1e-04 5.06e-01 4.402794894e+00 4.419884400e+00 1.6e-03 0.08
3 1.8e-07 9.7e-06 5.0e-10 1.02e+00 4.365867996e+00 4.365872314e+00 3.6e-07 0.09
4 2.2e-11 1.2e-09 6.3e-16 1.00e+00 4.365863010e+00 4.365863010e+00 4.3e-11 0.10
Optimizer terminated. Time: 0.10
Interior-point solution summary
Problem status : PRIMAL_AND_DUAL_FEASIBLE
Solution status : OPTIMAL
Primal. obj: 4.3658630098e+00 nrm: 2e+00 Viol. con: 1e-10 var: 0e+00 cones: 0e+00
Dual. obj: 4.3658630103e+00 nrm: 1e+01 Viol. con: 0e+00 var: 1e-09 cones: 0e+00
-------------------------------------------------------------------------------
Summary
-------------------------------------------------------------------------------
(CVXPY) May 15 04:23:20 PM: Problem status: optimal
(CVXPY) May 15 04:23:20 PM: Optimal value: 4.366e+00
(CVXPY) May 15 04:23:20 PM: Compilation took 9.323e-02 seconds
(CVXPY) May 15 04:23:20 PM: Solver (including time spent in interface) took 1.677e-01 seconds
Out[13]:
4.365863010698794
In [14]:
plt.figure()
plt.semilogy(np.abs(fista.objective - p.value), label="FISTA")
plt.semilogy(np.abs(ista.objective - p.value), label="ISTA")
plt.legend()
plt.show()
print("Objective value of CVXpy = {}".format(p.value))
print("Objective value of ISTA = {}".format(ista.objective[-1]))
print("Objective value of FISTA = {}".format(fista.objective[-1]))
Objective value of CVXpy = 4.365863010698794 Objective value of ISTA = 4.6376729011535645 Objective value of FISTA = 4.3882551193237305
In [ ]: