#!/usr/bin/env python
import time
import numpy as np
import tomograpy
import lo
# object
obj = tomograpy.centered_cubic_map(10, 64)
obj[:] = tomograpy.phantom.shepp_logan(obj.shape)
# data 
radius = 200
a = tomograpy.fov(obj, radius)
data = tomograpy.centered_stack(a, 128, n_images=60, radius=radius, max_lon=np.pi)
# model
kwargs = {"pb":"pb", "obj_rmin":1.5, "data_rmin":1.5}
P, D, obj_mask, data_mask = tomograpy.models.thomson(data, obj, u=.5, **kwargs)
# projection
t = time.time()
data[:] = (P * obj.ravel()).reshape(data.shape)
print("projection time : " + str(time.time() - t))
# data
# backprojection
t = time.time()
x0 = P.T * data.ravel()
bpj = x0.reshape(obj.shape)
print("backprojection time : " + str(time.time() - t))
# inversion using scipy.sparse.linalg
t = time.time()
sol = lo.acg(P, data.ravel(), D, 1e-3 * np.ones(3),  maxiter=100, tol=1e-8)
sol = sol.reshape(obj.shape)
print("inversion time : " + str(time.time() - t))

[(7, 'tomograpy.centered_cubic_map', 'tomograpy.centered_cubic_map', False, 'import tomograpy\n'), (8, 'tomograpy.phantom.shepp_logan', 'tomograpy.phantom.shepp_logan', False, 'import tomograpy\n'), (11, 'tomograpy.fov', 'tomograpy.fov', False, 'import tomograpy\n'), (12, 'tomograpy.centered_stack', 'tomograpy.centered_stack', False, 'import tomograpy\n'), (15, 'tomograpy.models.thomson', 'tomograpy.models.thomson', False, 'import tomograpy\n'), (17, 'time.time', 'time.time', False, 'import time\n'), (22, 'time.time', 'time.time', False, 'import time\n'), (27, 'time.time', 'time.time', False, 'import time\n'), (28, 'numpy.ones', 'np.ones', True, 'import numpy as np\n'), (19, 'time.time', 'time.time', False, 'import time\n'), (25, 'time.time', 'time.time', False, 'import time\n'), (30, 'time.time', 'time.time', False, 'import time\n')]