Merge branch 'master' into notimplementederror

CHANGELOG.md

@@ -28,8 +28,10 @@
     - Make Binner accept accelerated=False (#1887)
     - Added checks on memory allocations within `FiniteDifferenceLibrary.cpp` and verified the status of the return in `GradientOperator` (#1929)
     - Build release version of `cilacc.dll` for Windows. Previously was defaulting to the debug build (#1928)
+    - Armijo step size rule now by default initialises the search for a step size from the previously calculated step size (#1934)
   - Changes that break backwards compatibility:
     - CGLS will no longer automatically stop iterations once a default tolerance is reached. The option to pass `tolerance` will be deprecated to be replaced by `optimisation.utilities.callbacks` (#1892)
+
 
 * 24.1.0
   - New Features:
@@ -54,6 +56,7 @@
     - BlockOperator that would return a BlockDataContainer of shape (1,1) now returns the appropriate DataContainer. BlockDataContainer direct and adjoint methods accept DataContainer as parameter (#1802).
     - BlurringOperator: remove check for geometry class (old SIRF integration bug) (#1807)
     - The `ZeroFunction` and `ConstantFunction` now have a Lipschitz constant of 1. (#1768)
+    - Update dataexample remote data download to work with windows and use zenodo_get for data download (#1774)
   - Changes that break backwards compatibility:
     - Merged the files `BlockGeometry.py` and `BlockDataContainer.py` in `framework` to one file `block.py`. Please use `from cil.framework import BlockGeometry, BlockDataContainer` as before (#1799)
     - Bug fix in `FGP_TV` function to set the default behaviour not to enforce non-negativity (#1826).

Wrappers/Python/cil/optimisation/algorithms/FISTA.py

@@ -213,8 +213,19 @@ def update_objective(self):
         .. math:: f(x) + g(x)
 
         """
-        self.loss.append(self.f(self.x_old) + self.g(self.x_old))
+        self.loss.append(self.calculate_objective_function_at_point(self.x_old))
 
+    def calculate_objective_function_at_point(self, x):
+        """ Calculates the objective at a given point x
+
+        .. math:: f(x) + g(x)
+        
+        Parameters
+        ----------
+        x : DataContainer
+        
+        """
+        return self.f(x) + self.g(x)
 
 class FISTA(ISTA):
 

Wrappers/Python/cil/optimisation/algorithms/GD.py

@@ -84,7 +84,7 @@ def set_up(self, initial, objective_function, step_size, preconditioner):
         log.info("%s setting up", self.__class__.__name__)
 
         self.x = initial.copy()
-        self.objective_function = objective_function
+        self._objective_function = objective_function
 
         if step_size is None:
             self.step_size_rule = ArmijoStepSizeRule(
@@ -106,7 +106,7 @@ def set_up(self, initial, objective_function, step_size, preconditioner):
 
     def update(self):
         '''Performs a single iteration of the gradient descent algorithm'''
-        self.objective_function.gradient(self.x, out=self.gradient_update)
+        self._objective_function.gradient(self.x, out=self.gradient_update)
 
         if self.preconditioner is not None:
             self.preconditioner.apply(
@@ -117,7 +117,7 @@ def update(self):
         self.x.sapyb(1.0, self.gradient_update, -step_size, out=self.x)
 
     def update_objective(self):
-        self.loss.append(self.objective_function(self.solution))
+        self.loss.append(self._objective_function(self.solution))
 
     def should_stop(self):
         '''Stopping criterion for the gradient descent algorithm '''
@@ -132,3 +132,20 @@ def step_size(self):
         else:
             raise TypeError(
                 "There is not a constant step size, it is set by a step-size rule")
+
+    def calculate_objective_function_at_point(self, x):
+        """ Calculates the objective at a given point x
+
+        .. math:: f(x) + g(x)
+        
+        Parameters
+        ----------
+        x : DataContainer
+        
+        """
+        return self._objective_function(x)
+
+    @property
+    def objective_function(self):
+        warn('The attribute `objective_function` will be deprecated in the future. Please use `calculate_objective_function_at_point` instead.', DeprecationWarning, stacklevel=2)  
+        return self._objective_function

Wrappers/Python/cil/optimisation/utilities/StepSizeMethods.py

@@ -19,6 +19,9 @@
 from abc import ABC, abstractmethod
 import numpy
 from numbers import Number
+import logging
+
+log = logging.getLogger(__name__)
 
 class StepSizeRule(ABC):
     """
@@ -82,6 +85,9 @@ class ArmijoStepSizeRule(StepSizeRule):
         The amount the step_size is reduced if the criterion is not met
     max_iterations: integer, optional, default is numpy.ceil (2 * numpy.log10(alpha) / numpy.log10(2))
         The maximum number of iterations to find a suitable step size 
+    warmstart: Boolean, default is True
+        If `warmstart = True` the initial step size at each Armijo iteration is the calculated step size from the last iteration. If `warmstart = False` at each  Armijo iteration, the initial step size is reset to the original, large `alpha`. 
+        In the case of *well-behaved* convex functions, `warmstart = True` is likely to be computationally less expensive. In the case of non-convex functions, or particularly tricky functions, setting `warmstart = False` may be beneficial. 
 
     Reference
     ------------
@@ -91,21 +97,23 @@ class ArmijoStepSizeRule(StepSizeRule):
 
     """
 
-    def __init__(self, alpha=1e6, beta=0.5, max_iterations=None):
+    def __init__(self, alpha=1e6, beta=0.5, max_iterations=None, warmstart=True):
         '''Initialises the step size rule 
         '''
 
         self.alpha_orig = alpha
         if self.alpha_orig is None: # Can be removed when alpha and beta are deprecated in GD
             self.alpha_orig = 1e6 
-
+        self.alpha = self.alpha_orig
         self.beta = beta 
         if self.beta is None:  # Can be removed when alpha and beta are deprecated in GD
             self.beta = 0.5
 
         self.max_iterations = max_iterations
         if self.max_iterations is None:
             self.max_iterations = numpy.ceil(2 * numpy.log10(self.alpha_orig) / numpy.log10(2))
+
+        self.warmstart=warmstart
 
     def get_step_size(self, algorithm):
         """
@@ -117,26 +125,33 @@ def get_step_size(self, algorithm):
 
         """
         k = 0
-        self.alpha = self.alpha_orig
-        f_x = algorithm.objective_function(algorithm.solution)
+        if not self.warmstart:  
+            self.alpha = self.alpha_orig
+
+        f_x = algorithm.calculate_objective_function_at_point(algorithm.solution)
 
         self.x_armijo = algorithm.solution.copy()
-
+
+        log.debug("Starting Armijo backtracking with initial step size: %f", self.alpha)
+
         while k < self.max_iterations:
 
             algorithm.gradient_update.multiply(self.alpha, out=self.x_armijo)
             algorithm.solution.subtract(self.x_armijo, out=self.x_armijo)
 
-            f_x_a = algorithm.objective_function(self.x_armijo)
+            f_x_a = algorithm.calculate_objective_function_at_point(self.x_armijo)
             sqnorm = algorithm.gradient_update.squared_norm()
             if f_x_a - f_x <= - (self.alpha/2.) * sqnorm:
                 break
             k += 1.
             self.alpha *= self.beta
+
+        log.info("Armijo rule took %d iterations to find step size", k)
 
         if k == self.max_iterations:
             raise ValueError(
                 'Could not find a proper step_size in {} loops. Consider increasing alpha or max_iterations.'.format(self.max_iterations))
+
         return self.alpha
 
 

Wrappers/Python/cil/utilities/dataexample.py

@@ -25,10 +25,9 @@
 import os.path
 import sys
 from zipfile import ZipFile
-from urllib.request import urlopen
-from io import BytesIO
 from scipy.io import loadmat
 from cil.io import NEXUSDataReader, NikonDataReader, ZEISSDataReader
+from zenodo_get import zenodo_get
 
 class DATA(object):
     @classmethod
@@ -46,21 +45,15 @@ def get(cls, size=None, scale=(0,1), **kwargs):
 class REMOTEDATA(DATA):
 
     FOLDER = ''
-    URL = ''
-    FILE_SIZE = ''
+    ZENODO_RECORD = ''
+    ZIP_FILE = ''
 
     @classmethod
     def get(cls, data_dir):
         return None
 
     @classmethod
-    def _download_and_extract_from_url(cls, data_dir):
-        with urlopen(cls.URL) as response:
-            with BytesIO(response.read()) as bytes, ZipFile(bytes) as zipfile:
-                zipfile.extractall(path = data_dir)
-
-    @classmethod
-    def download_data(cls, data_dir):
+    def download_data(cls, data_dir, prompt=True):
         '''
         Download a dataset from a remote repository
 
@@ -71,14 +64,18 @@ def download_data(cls, data_dir):
 
         '''
         if os.path.isdir(os.path.join(data_dir, cls.FOLDER)):
-            print("Dataset already exists in " + data_dir)
+            print("Dataset folder already exists in " + data_dir)
         else:
-            if input("Are you sure you want to download " + cls.FILE_SIZE + " dataset from " + cls.URL + " ? (y/n)") == "y":
-                print('Downloading dataset from ' + cls.URL)
-                cls._download_and_extract_from_url(os.path.join(data_dir,cls.FOLDER))
-                print('Download complete')
-            else:
+            user_input = input("Are you sure you want to download {cls.ZIP_FILE} dataset from Zenodo record {cls.ZENODO_RECORD}? [Y/n]: ") if prompt else 'y'
+            if user_input.lower() not in ('y', 'yes'):
                 print('Download cancelled')
+                return False
+
+            zenodo_get([cls.ZENODO_RECORD, '-g', cls.ZIP_FILE, '-o', data_dir])
+            with ZipFile(os.path.join(data_dir, cls.ZIP_FILE), 'r') as zip_ref:
+                zip_ref.extractall(os.path.join(data_dir, cls.FOLDER))
+            os.remove(os.path.join(data_dir, cls.ZIP_FILE))
+            return True
 
 class BOAT(CILDATA):
     @classmethod
@@ -195,15 +192,21 @@ def get(cls, **kwargs):
 class WALNUT(REMOTEDATA):
     '''
     A microcomputed tomography dataset of a walnut from https://zenodo.org/records/4822516
+
+    Example
+    --------
+    >>> data_dir = 'my_PC/data_folder'
+    >>> dataexample.WALNUT.download_data(data_dir) # download the data
+    >>> dataexample.WALNUT.get(data_dir) # load the data
     '''
     FOLDER = 'walnut'
-    URL = 'https://zenodo.org/record/4822516/files/walnut.zip'
-    FILE_SIZE = '6.4 GB'
+    ZENODO_RECORD = '4822516'
+    ZIP_FILE = 'walnut.zip'
 
     @classmethod
     def get(cls, data_dir):
         '''
-        A microcomputed tomography dataset of a walnut from https://zenodo.org/records/4822516
+        Get the microcomputed tomography dataset of a walnut from https://zenodo.org/records/4822516
         This function returns the raw projection data from the .txrm file
 
         Parameters
@@ -227,15 +230,21 @@ def get(cls, data_dir):
 class USB(REMOTEDATA):
     '''
     A microcomputed tomography dataset of a usb memory stick from https://zenodo.org/records/4822516
+
+    Example
+    --------
+    >>> data_dir = 'my_PC/data_folder'
+    >>> dataexample.USB.download_data(data_dir) # download the data
+    >>> dataexample.USB.get(data_dir) # load the data
     '''
     FOLDER = 'USB'
-    URL = 'https://zenodo.org/record/4822516/files/usb.zip'
-    FILE_SIZE = '3.2 GB'
+    ZENODO_RECORD = '4822516'
+    ZIP_FILE = 'usb.zip'
 
     @classmethod
     def get(cls, data_dir):
         '''
-        A microcomputed tomography dataset of a usb memory stick from https://zenodo.org/records/4822516
+        Get the microcomputed tomography dataset of a usb memory stick from https://zenodo.org/records/4822516
         This function returns the raw projection data from the .txrm file
 
         Parameters
@@ -259,15 +268,21 @@ def get(cls, data_dir):
 class KORN(REMOTEDATA):
     '''
     A microcomputed tomography dataset of a sunflower seeds in a box from https://zenodo.org/records/6874123
+
+    Example
+    --------
+    >>> data_dir = 'my_PC/data_folder'
+    >>> dataexample.KORN.download_data(data_dir) # download the data
+    >>> dataexample.KORN.get(data_dir) # load the data
     '''
     FOLDER = 'korn'
-    URL = 'https://zenodo.org/record/6874123/files/korn.zip'
-    FILE_SIZE = '2.9 GB'
+    ZENODO_RECORD = '6874123'
+    ZIP_FILE = 'korn.zip'
 
     @classmethod
     def get(cls, data_dir):
         '''
-        A microcomputed tomography dataset of a sunflower seeds in a box from https://zenodo.org/records/6874123
+        Get the microcomputed tomography dataset of a sunflower seeds in a box from https://zenodo.org/records/6874123
         This function returns the raw projection data from the .xtekct file
 
         Parameters
@@ -279,6 +294,7 @@ def get(cls, data_dir):
         -------
         ImageData
             The korn dataset
+
         '''
         filepath = os.path.join(data_dir, cls.FOLDER, 'Korn i kasse','47209 testscan korn01_recon.xtekct')
         try:
@@ -293,10 +309,40 @@ class SANDSTONE(REMOTEDATA):
     '''
     A synchrotron x-ray tomography dataset of sandstone from https://zenodo.org/records/4912435
     A small subset of the data containing selected projections and 4 slices of the reconstruction
+
+    Example
+    --------
+    >>> data_dir = 'my_PC/data_folder'
+    >>> dataexample.SANDSTONE.download_data(data_dir) # download the data
+    >>> dataexample.SANDSTONE.get(data_dir) # load the data
     '''
     FOLDER = 'sandstone'
-    URL = 'https://zenodo.org/records/4912435/files/small.zip'
-    FILE_SIZE = '227 MB'
+    ZENODO_RECORD = '4912435'
+    ZIP_FILE = 'small.zip'
+
+    @classmethod
+    def get(cls, data_dir, filename):
+        '''
+        Get the synchrotron x-ray tomography dataset of sandstone from https://zenodo.org/records/4912435
+        A small subset of the data containing selected projections and 4 slices of the reconstruction
+        Parameters
+        ----------
+        data_dir: str
+           The path to the directory where the dataset is stored. Data can be downloaded with dataexample.SANDSTONE.download_data(data_dir)
+
+        file: str
+            The slices or projections to return, specify the path to the file within the data_dir
+
+        Returns
+        -------
+        ImageData
+            The selected sandstone dataset
+        '''
+        extension = os.path.splitext(filename)[1]
+        if extension == '.mat':
+            return loadmat(os.path.join(data_dir,filename))
+        raise KeyError(f"Unknown extension: {extension}")
+
 
 class TestData(object):
     '''Class to return test data