Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Spdhg unit tests - speed up #1500

Draft
wants to merge 4 commits into
base: master
Choose a base branch
from
Draft

Spdhg unit tests - speed up #1500

Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
079935b
initial testing
MargaretDuff Aug 23, 2023
43e3dc4
Sped up PDHG and SPDHG testing
MargaretDuff Aug 24, 2023
004ab2f
Removed timing statements
MargaretDuff Aug 24, 2023
ea4f114
Merge branch 'master' into SPDHG_unit_tests
MargaretDuff Nov 22, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
70 changes: 41 additions & 29 deletions Wrappers/Python/test/test_algorithms.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -755,7 +755,8 @@ class TestSPDHG(unittest.TestCase):

@unittest.skipUnless(has_astra, "cil-astra not available")
def test_SPDHG_vs_PDHG_implicit(self):
data = dataexample.SIMPLE_PHANTOM_2D.get(size=(128,128))

data = dataexample.SIMPLE_PHANTOM_2D.get(size=(16,16))

ig = data.geometry
ig.voxel_size_x = 0.1
Expand Down Expand Up @@ -798,14 +799,14 @@ def test_SPDHG_vs_PDHG_implicit(self):
sigma_tmp = 1.
tau = sigma_tmp / operator.adjoint(tau_tmp * operator.range_geometry().allocate(1.))
sigma = tau_tmp / operator.direct(sigma_tmp * operator.domain_geometry().allocate(1.))

# Setup and run the PDHG algorithm
pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma,
max_iteration = 1000,
update_objective_interval = 500)
max_iteration = 70,
update_objective_interval = 1000)
pdhg.run(verbose=0)
subsets = 10

subsets = 5
size_of_subsets = int(len(angles)/subsets)
# take angles and create uniform subsets in uniform+sequential setting
list_angles = [angles[i:i+size_of_subsets] for i in range(0, len(angles), size_of_subsets)]
Expand All @@ -831,10 +832,12 @@ def test_SPDHG_vs_PDHG_implicit(self):
G = alpha * TotalVariation(50, 1e-4, lower=0)

prob = [1/len(A)]*len(A)

spdhg = SPDHG(f=F,g=G,operator=A,
max_iteration = 1000,
update_objective_interval=200, prob = prob)
max_iteration = 320,
update_objective_interval=1000, prob = prob)
spdhg.run(1000, verbose=0)

qm = (mae(spdhg.get_output(), pdhg.get_output()),
mse(spdhg.get_output(), pdhg.get_output()),
psnr(spdhg.get_output(), pdhg.get_output())
Expand All @@ -849,7 +852,7 @@ def test_SPDHG_vs_PDHG_implicit(self):

@unittest.skipUnless(has_astra, "ccpi-astra not available")
def test_SPDHG_vs_PDHG_explicit(self):
data = dataexample.SIMPLE_PHANTOM_2D.get(size=(128,128))
data = dataexample.SIMPLE_PHANTOM_2D.get(size=(16,16))

ig = data.geometry
ig.voxel_size_x = 0.1
Expand Down Expand Up @@ -880,7 +883,7 @@ def test_SPDHG_vs_PDHG_explicit(self):
raise ValueError('Unsupported Noise ', noise)

#%% 'explicit' SPDHG, scalar step-sizes
subsets = 10
subsets = 5
size_of_subsets = int(len(angles)/subsets)
# create Gradient operator
op1 = GradientOperator(ig)
Expand Down Expand Up @@ -909,9 +912,11 @@ def test_SPDHG_vs_PDHG_explicit(self):

prob = [1/(2*subsets)]*(len(A)-1) + [1/2]
spdhg = SPDHG(f=F,g=G,operator=A,
max_iteration = 1000,
update_objective_interval=200, prob = prob)
max_iteration = 220,
update_objective_interval=220, prob = prob)

spdhg.run(1000, verbose=0)


#%% 'explicit' PDHG, scalar step-sizes
op1 = GradientOperator(ig)
Expand All @@ -928,10 +933,11 @@ def test_SPDHG_vs_PDHG_explicit(self):
f = BlockFunction(f1, f2)
# Setup and run the PDHG algorithm
pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma)
pdhg.max_iteration = 1000
pdhg.update_objective_interval = 200
pdhg.max_iteration = 180
pdhg.update_objective_interval =180

pdhg.run(1000, verbose=0)

#%% show diff between PDHG and SPDHG
# plt.imshow(spdhg.get_output().as_array() -pdhg.get_output().as_array())
# plt.colorbar()
Expand All @@ -950,7 +956,7 @@ def test_SPDHG_vs_PDHG_explicit(self):

@unittest.skipUnless(has_astra, "ccpi-astra not available")
def test_SPDHG_vs_SPDHG_explicit_axpby(self):
data = dataexample.SIMPLE_PHANTOM_2D.get(size=(128,128), dtype=numpy.float32)
data = dataexample.SIMPLE_PHANTOM_2D.get(size=(16,16), dtype=numpy.float32)

ig = data.geometry
ig.voxel_size_x = 0.1
Expand Down Expand Up @@ -986,7 +992,7 @@ def test_SPDHG_vs_SPDHG_explicit_axpby(self):
raise ValueError('Unsupported Noise ', noise)

#%% 'explicit' SPDHG, scalar step-sizes
subsets = 10
subsets = 5
size_of_subsets = int(len(angles)/subsets)
# create GradientOperator operator
op1 = GradientOperator(ig)
Expand Down Expand Up @@ -1018,17 +1024,20 @@ def test_SPDHG_vs_SPDHG_explicit_axpby(self):
prob = [1/(2*subsets)]*(len(A)-1) + [1/2]
algos = []
algos.append( SPDHG(f=F,g=G,operator=A,
max_iteration = 1000,
update_objective_interval=200, prob = prob.copy(), use_axpby=True)
max_iteration = 330,
update_objective_interval=330, prob = prob.copy(), use_axpby=True)
)

algos[0].run(1000, verbose=0)

algos.append( SPDHG(f=F,g=G,operator=A,
max_iteration = 1000,
update_objective_interval=200, prob = prob.copy(), use_axpby=False)
max_iteration = 330,
update_objective_interval=330, prob = prob.copy(), use_axpby=False)
)

algos[1].run(1000, verbose=0)



# np.testing.assert_array_almost_equal(algos[0].get_output().as_array(), algos[1].get_output().as_array())
qm = (mae(algos[0].get_output(), algos[1].get_output()),
Expand All @@ -1037,12 +1046,12 @@ def test_SPDHG_vs_SPDHG_explicit_axpby(self):
)
logging.info ("Quality measures {}".format(qm))
assert qm[0] < 0.005
assert qm[1] < 3.e-05
assert qm[1] < 5.e-05


@unittest.skipUnless(has_astra, "ccpi-astra not available")
def test_PDHG_vs_PDHG_explicit_axpby(self):
data = dataexample.SIMPLE_PHANTOM_2D.get(size=(128,128))
data = dataexample.SIMPLE_PHANTOM_2D.get(size=(16,16))
ig = data.geometry
ig.voxel_size_x = 0.1
ig.voxel_size_y = 0.1
Expand Down Expand Up @@ -1091,18 +1100,21 @@ def test_PDHG_vs_PDHG_explicit_axpby(self):
# Setup and run the PDHG algorithm

algos = []

algos.append( PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma,
max_iteration = 1000,
update_objective_interval=200, use_axpby=True)
max_iteration = 300,
update_objective_interval=1000, use_axpby=True)
)

algos[0].run(1000, verbose=0)

algos.append( PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma,
max_iteration = 1000,
update_objective_interval=200, use_axpby=False)
max_iteration = 300,
update_objective_interval=1000, use_axpby=False)
)
algos[1].run(1000, verbose=0)

algos[1].run(1000, verbose=0)

qm = (mae(algos[0].get_output(), algos[1].get_output()),
mse(algos[0].get_output(), algos[1].get_output()),
psnr(algos[0].get_output(), algos[1].get_output())
Expand Down