diff options
-rw-r--r-- | NOTICE.txt | 15 | ||||
-rwxr-xr-x | Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py | 30 | ||||
-rw-r--r-- | Wrappers/Python/ccpi/optimisation/functions/FunctionOperatorComposition.py | 64 |
3 files changed, 82 insertions, 27 deletions
diff --git a/NOTICE.txt b/NOTICE.txt new file mode 100644 index 0000000..c107329 --- /dev/null +++ b/NOTICE.txt @@ -0,0 +1,15 @@ +CCPi Core Imaging Library (CIL). +Copyright 2017 Rutherford Appleton Laboratory STFC +Copyright 2017 University of Manchester + +This software product is developed for the Collaborative Computational +Project in Tomographic Imaging CCPi (http://www.ccpi.ac.uk/) at RAL STFC (http://www.stfc.ac.uk), University of Manchester +and other contributing institutions. + +Main contributors in alphabetical order: +Evelina Ametova (UoM) +Jakob Jorgensen (UoM) +Daniil Kazantsev (Diamond Light Source) +Srikanth Nagella (STFC) +Edoardo Pasca (STFC) +Evangelos Papoutsellis (UoM) diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py b/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py index a14378c..0ba2c55 100755 --- a/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py @@ -16,7 +16,7 @@ # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. # See the License for the specific language governing permissions and # limitations under the License. -import time +import time, functools from numbers import Integral class Algorithm(object): @@ -148,14 +148,36 @@ class Algorithm(object): print ("Stop cryterion has been reached.") i = 0 + if verbose: + print ("{:>9} {:>10} {:>11} {:>20}".format('Iter', + 'Max Iter', + 's/Iter', + 'Objective Value')) for _ in self: if (self.iteration -1) % self.update_objective_interval == 0: if verbose: - print ("Iteration {}/{}, = {}".format(self.iteration-1, - self.max_iteration, self.get_last_objective()) ) + print (self.verbose_output()) if callback is not None: callback(self.iteration -1, self.get_last_objective(), self.x) i += 1 if i == iterations: break - + + def verbose_output(self): + '''Creates a nice tabulated output''' + timing = self.timing[-self.update_objective_interval-1:-1] + if len (timing) == 0: + t = 0 + else: + t = sum(timing)/len(timing) + el = [ self.iteration-1, + self.max_iteration, + "{:.3f} s/it".format(t), + self.get_last_objective() ] + + if type(el[-1] ) == list: + string = functools.reduce(lambda x,y: x+' {:>15.5e}'.format(y), el[-1],'') + out = "{:>9} {:>10} {:>11} {}".format(*el[:-1] , string) + else: + out = "{:>9} {:>10} {:>11} {:>20.5e}".format(*el) + return out diff --git a/Wrappers/Python/ccpi/optimisation/functions/FunctionOperatorComposition.py b/Wrappers/Python/ccpi/optimisation/functions/FunctionOperatorComposition.py index 8895f3d..e788d5a 100644 --- a/Wrappers/Python/ccpi/optimisation/functions/FunctionOperatorComposition.py +++ b/Wrappers/Python/ccpi/optimisation/functions/FunctionOperatorComposition.py @@ -50,41 +50,59 @@ class FunctionOperatorComposition(Function): self.operator.adjoint(tmp, out=out) - - - #TODO do not know if we need it - #def call_adjoint(self, x): - # - # return self.function(self.operator.adjoint(x)) - - - #def convex_conjugate(self, x): - # - # ''' convex_conjugate does not take into account the Operator''' - # return self.function.convex_conjugate(x) - - - if __name__ == '__main__': - from ccpi.framework import ImageGeometry + from ccpi.framework import ImageGeometry, AcquisitionGeometry from ccpi.optimisation.operators import Gradient from ccpi.optimisation.functions import L2NormSquared + from ccpi.astra.ops import AstraProjectorSimple + import numpy as np - M, N, K = 2,3 + M, N= 50, 50 ig = ImageGeometry(voxel_num_x=M, voxel_num_y = N) - G = Gradient(ig) + detectors = N + angles_num = N + det_w = 1.0 + + angles = np.linspace(0, np.pi, angles_num, endpoint=False) + ag = AcquisitionGeometry('parallel', + '2D', + angles, + detectors,det_w) + + + Aop = AstraProjectorSimple(ig, ag, 'cpu') + + u = ig.allocate('random_int') + u1 = ig.allocate('random_int') + b = Aop.direct(u1) + + +# G = Gradient(ig) alpha = 0.5 - f = L2NormSquared() - f_comp = FunctionOperatorComposition(G, alpha * f) - x = ig.allocate('random_int') - print(f_comp.gradient(x).shape + f1 = alpha * L2NormSquared(b=b) + + f_comp = FunctionOperatorComposition(f1, Aop) + + print(f_comp(u)) + + + z1 = Aop.direct(u) + tmp = 0.5 * ((z1 - b)**2).sum() + + + print(tmp) + + + + + + - ) |