AstraMC vs Channelwise operator

[epapoutsellis]

Below a comparison between AstraMC and Channelwise.
Channelwise is 1-2 sec slower.

Computing norm using MC operator: 0.0814945170423016sec
Computing norm using ChanWise operator: 0.08279416302684695sec
Computing direct out using MC operator: 7.874376090010628sec
Computing direct out using ChanWise operator: 8.888477805012371sec
Computing adjoint out using MC operator: 7.901970839011483sec
Computing adjoint out using ChanWise operator: 8.783782225975301sec
Computing direct (no out) using MC operator: 7.892078642034903sec
Computing direct (no out) using ChanWise operator: 8.707080838037655sec
Computing adjoint (no out) using MC operator: 7.789297449984588sec
Computing adjoint (no out) using ChanWise operator: 9.125933162984438sec

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wed May 13 10:01:00 2020

@author: evangelos
"""

from ccpi.framework import ImageGeometry, AcquisitionGeometry
from ccpi.optimisation.operators import ChannelwiseOperator
from ccpi.astra.operators import AstraProjectorMC, AstraProjectorSimple
from timeit import default_timer as timer
import numpy as np

chan = 50

#ig = ImageGeometry(voxel_num_x = 150, 
#                   voxel_num_y = 150,
#                   voxel_size_x = 1., 
#                   voxel_size_y = 1., 
#                   channels = chan,
#                   dimension_labels = ['channel','horizontal_y','horizontal_x'])
#
#ag = AcquisitionGeometry('cone',
#                         '2D',
#                         np.linspace(0,np.pi,150),
#                         pixel_num_h = 150,
#                         pixel_size_h = 1,                           
#                         dist_source_center = 10, 
#                         dist_center_detector = 10,
#                         channels = chan,
#                         dimension_labels = ['channel', 'angle', 'horizontal'])

ig = ImageGeometry(voxel_num_x = 50, 
                   voxel_num_y = 50,
                   voxel_size_x = 1., 
                   voxel_size_y = 1., 
                   channels = chan,
                   dimension_labels = ['channel','horizontal_y','horizontal_x'])

ag = AcquisitionGeometry('parallel',
                         '2D',
                         np.linspace(0,np.pi,50),
                         pixel_num_h = 50,
                         pixel_size_h = 1,                           
                         channels = chan,
                         dimension_labels = ['channel', 'angle', 'horizontal'])

Aop1 = AstraProjectorMC(ig, ag, 'cpu')


igtmp = ig.clone()
igtmp.shape = ig.shape[1:]
igtmp.dimension_labels = ['horizontal_y', 'horizontal_x']
igtmp.channels = 1

agtmp = ag.clone()
agtmp.shape = ag.shape[1:]
agtmp.dimension_labels = ['angle', 'horizontal']
agtmp.channels = 1

tmp_op = AstraProjectorSimple(igtmp, agtmp, 'cpu')
Aop2 = ChannelwiseOperator(tmp_op, chan)

x1 = ig.allocate('random')
tmp_x1a = ig.allocate()
tmp_x1b = ig.allocate()

x2 = ag.allocate('random')
tmp_x2a = ag.allocate()
tmp_x2b = ag.allocate()

t1 = timer()
a1 = Aop1.norm()
t2 = timer()

t3 = timer()
a2 = Aop2.norm()
t4 = timer()

np.testing.assert_almost_equal(a1, a2, decimal = 1)
print(" Computing norm using MC operator: {}sec ".format(t2 - t1))
print(" Computing norm using ChanWise operator: {}sec ".format(t4 - t3))

#%%

t5 = timer()
for i in range(100):
    Aop1.direct(x1, out = tmp_x1a)
t6 = timer()
print(" Computing direct out using MC operator: {}sec ".format(t6 - t5))

t7 = timer()
for i in range(100):
    Aop2.direct(x1, out = tmp_x1b)
t8 = timer()

#np.testing.assert_array_almost_equal(a1.as_array(), a2.as_array())

print(" Computing direct out using ChanWise operator: {}sec ".format(t8 - t7))


#%%

t9 = timer()
for i in range(100):
#    a1 = Aop1.adjoint(x2)
    Aop1.adjoint(x2, out = tmp_x2a)
t10 = timer()
print(" Computing adjoint out using MC operator: {}sec ".format(t10 - t9))

t11 = timer()
for i in range(100):
#    a2 = Aop2.adjoint(x2)
    Aop2.adjoint(x2, out = tmp_x2b)
t12 = timer()

#np.testing.assert_array_almost_equal(a1.as_array(), a2.as_array())

print(" Computing adjoint out using ChanWise operator: {}sec ".format(t12 - t11))


#%%

t13 = timer()
for i in range(100):
    a1 = Aop1.direct(x1)
t14 = timer()
print(" Computing direct (no out) using MC operator: {}sec ".format(t14 - t13))

t15 = timer()
for i in range(100):
    a2 = Aop2.direct(x1)
t16 = timer()

#np.testing.assert_array_almost_equal(a1.as_array(), a2.as_array())

print(" Computing direct (no out) using ChanWise operator: {}sec ".format(t16 - t15))


t17 = timer()
for i in range(100):
    a1 = Aop1.adjoint(x2)
t18 = timer()
print(" Computing adjoint (no out) using MC operator: {}sec ".format(t18 - t17))

t19 = timer()
for i in range(100):
    a2 = Aop2.adjoint(x2)
t20 = timer()

#np.testing.assert_array_almost_equal(a1.as_array(), a2.as_array())

print(" Computing adjoint (no out) using ChanWise operator: {}sec ".format(t20 - t19))

[epapoutsellis]

I suggest to uncomment these lines here as a temporary solution, because atm master branch does not work for multichannel data.