From 824663d20ea31d49ee5771c7c3f7925f8487ac0e Mon Sep 17 00:00:00 2001 From: vais-ral Date: Tue, 19 Mar 2019 16:50:33 +0000 Subject: Set theme jekyll-theme-slate --- _config.yml | 1 + 1 file changed, 1 insertion(+) create mode 100644 _config.yml diff --git a/_config.yml b/_config.yml new file mode 100644 index 0000000..c741881 --- /dev/null +++ b/_config.yml @@ -0,0 +1 @@ +theme: jekyll-theme-slate \ No newline at end of file -- cgit v1.2.3 From 91ac62b7069222c7cdb7e6777a3bdbda65383757 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Tue, 19 Mar 2019 16:55:02 +0000 Subject: initial sphinx doc --- docs/Makefile | 19 +++++ docs/make.bat | 35 ++++++++ docs/source/conf.py | 187 +++++++++++++++++++++++++++++++++++++++++++ docs/source/framework.rst | 20 +++++ docs/source/index.rst | 23 ++++++ docs/source/optimisation.md | 113 ++++++++++++++++++++++++++ docs/source/optimisation.rst | 92 +++++++++++++++++++++ 7 files changed, 489 insertions(+) create mode 100755 docs/Makefile create mode 100755 docs/make.bat create mode 100755 docs/source/conf.py create mode 100644 docs/source/framework.rst create mode 100755 docs/source/index.rst create mode 100755 docs/source/optimisation.md create mode 100644 docs/source/optimisation.rst diff --git a/docs/Makefile b/docs/Makefile new file mode 100755 index 0000000..4801716 --- /dev/null +++ b/docs/Makefile @@ -0,0 +1,19 @@ +# Minimal makefile for Sphinx documentation +# + +# You can set these variables from the command line. +SPHINXOPTS = +SPHINXBUILD = sphinx-build +SOURCEDIR = source +BUILDDIR = docs/build + +# Put it first so that "make" without argument is like "make help". +help: + @$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) + +.PHONY: help Makefile + +# Catch-all target: route all unknown targets to Sphinx using the new +# "make mode" option. $(O) is meant as a shortcut for $(SPHINXOPTS). +%: Makefile + @$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) \ No newline at end of file diff --git a/docs/make.bat b/docs/make.bat new file mode 100755 index 0000000..f88c76c --- /dev/null +++ b/docs/make.bat @@ -0,0 +1,35 @@ +@ECHO OFF + +pushd %~dp0 + +REM Command file for Sphinx documentation + +if "%SPHINXBUILD%" == "" ( + set SPHINXBUILD=sphinx-build +) +set SOURCEDIR=source +set BUILDDIR=docs/build + +if "%1" == "" goto help + +%SPHINXBUILD% >NUL 2>NUL +if errorlevel 9009 ( + echo. + echo.The 'sphinx-build' command was not found. Make sure you have Sphinx + echo.installed, then set the SPHINXBUILD environment variable to point + echo.to the full path of the 'sphinx-build' executable. Alternatively you + echo.may add the Sphinx directory to PATH. + echo. + echo.If you don't have Sphinx installed, grab it from + echo.http://sphinx-doc.org/ + exit /b 1 +) + +%SPHINXBUILD% -M %1 %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% +goto end + +:help +%SPHINXBUILD% -M help %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% + +:end +popd diff --git a/docs/source/conf.py b/docs/source/conf.py new file mode 100755 index 0000000..406877b --- /dev/null +++ b/docs/source/conf.py @@ -0,0 +1,187 @@ +# -*- coding: utf-8 -*- +# +# Configuration file for the Sphinx documentation builder. +# +# This file does only contain a selection of the most common options. For a +# full list see the documentation: +# http://www.sphinx-doc.org/en/master/config + +# -- Path setup -------------------------------------------------------------- + +# If extensions (or modules to document with autodoc) are in another directory, +# add these directories to sys.path here. If the directory is relative to the +# documentation root, use os.path.abspath to make it absolute, like shown here. +# +# import os +# import sys +# sys.path.insert(0, os.path.abspath('.')) + + +# -- Project information ----------------------------------------------------- + +project = 'CCPi-Framework' +copyright = '2019, Edoardo Pasca' +author = 'Edoardo Pasca' + +# The short X.Y version +version = '' +# The full version, including alpha/beta/rc tags +release = '19.02' + + +# -- General configuration --------------------------------------------------- + +# If your documentation needs a minimal Sphinx version, state it here. +# +# needs_sphinx = '1.0' + +# Add any Sphinx extension module names here, as strings. They can be +# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom +# ones. +extensions = [ + 'sphinx.ext.autodoc', + 'sphinx.ext.doctest', + 'sphinx.ext.todo', + 'sphinx.ext.coverage', + 'sphinx.ext.mathjax', + 'sphinx.ext.viewcode', +] + +# Add any paths that contain templates here, relative to this directory. +templates_path = ['docstemplates'] + +# The suffix(es) of source filenames. +# You can specify multiple suffix as a list of string: +# +source_suffix = ['.rst', '.md'] +#source_suffix = '.rst' + +# The master toctree document. +master_doc = 'index' + +# The language for content autogenerated by Sphinx. Refer to documentation +# for a list of supported languages. +# +# This is also used if you do content translation via gettext catalogs. +# Usually you set "language" from the command line for these cases. +language = 'en' + +# List of patterns, relative to source directory, that match files and +# directories to ignore when looking for source files. +# This pattern also affects html_static_path and html_extra_path. +exclude_patterns = [] + +# The name of the Pygments (syntax highlighting) style to use. +pygments_style = None + + +# -- Options for HTML output ------------------------------------------------- + +# The theme to use for HTML and HTML Help pages. See the documentation for +# a list of builtin themes. +# +html_theme = 'classic' + +# Theme options are theme-specific and customize the look and feel of a theme +# further. For a list of options available for each theme, see the +# documentation. +# +# html_theme_options = {} + +# Add any paths that contain custom static files (such as style sheets) here, +# relative to this directory. They are copied after the builtin static files, +# so a file named "default.css" will overwrite the builtin "default.css". +html_static_path = ['docsstatic'] + +# Custom sidebar templates, must be a dictionary that maps document names +# to template names. +# +# The default sidebars (for documents that don't match any pattern) are +# defined by theme itself. Builtin themes are using these templates by +# default: ``['localtoc.html', 'relations.html', 'sourcelink.html', +# 'searchbox.html']``. +# +# html_sidebars = {} + + +# -- Options for HTMLHelp output --------------------------------------------- + +# Output file base name for HTML help builder. +htmlhelp_basename = 'CCPi-Frameworkdoc' + + +# -- Options for LaTeX output ------------------------------------------------ + +latex_elements = { + # The paper size ('letterpaper' or 'a4paper'). + # + # 'papersize': 'letterpaper', + + # The font size ('10pt', '11pt' or '12pt'). + # + # 'pointsize': '10pt', + + # Additional stuff for the LaTeX preamble. + # + # 'preamble': '', + + # Latex figure (float) alignment + # + # 'figure_align': 'htbp', +} + +# Grouping the document tree into LaTeX files. List of tuples +# (source start file, target name, title, +# author, documentclass [howto, manual, or own class]). +latex_documents = [ + (master_doc, 'CCPi-Framework.tex', 'CCPi-Framework Documentation', + 'Edoardo Pasca', 'manual'), +] + + +# -- Options for manual page output ------------------------------------------ + +# One entry per manual page. List of tuples +# (source start file, name, description, authors, manual section). +man_pages = [ + (master_doc, 'ccpi-framework', 'CCPi-Framework Documentation', + [author], 1) +] + + +# -- Options for Texinfo output ---------------------------------------------- + +# Grouping the document tree into Texinfo files. List of tuples +# (source start file, target name, title, author, +# dir menu entry, description, category) +texinfo_documents = [ + (master_doc, 'CCPi-Framework', 'CCPi-Framework Documentation', + author, 'CCPi-Framework', 'One line description of project.', + 'Miscellaneous'), +] + + +# -- Options for Epub output ------------------------------------------------- + +# Bibliographic Dublin Core info. +epub_title = project + +# The unique identifier of the text. This can be a ISBN number +# or the project homepage. +# +# epub_identifier = '' + +# A unique identification for the text. +# +# epub_uid = '' + +# A list of files that should not be packed into the epub file. +epub_exclude_files = ['search.html'] + + +# -- Extension configuration ------------------------------------------------- + +# -- Options for todo extension ---------------------------------------------- + +# If true, `todo` and `todoList` produce output, else they produce nothing. +todo_include_todos = True diff --git a/docs/source/framework.rst b/docs/source/framework.rst new file mode 100644 index 0000000..4da0d6a --- /dev/null +++ b/docs/source/framework.rst @@ -0,0 +1,20 @@ +DataContainers and Geometry +======================================== +| + +.. autoclass:: ccpi.framework.DataContainer + :members: +.. autoclass:: ccpi.framework.ImageData + :members: +.. autoclass:: ccpi.framework.AcquisitionData + :members: +.. autoclass:: ccpi.framework.AcquisitionGeometry + :members: +.. autoclass:: ccpi.framework.ImageGeometry + :members: +.. autoclass:: ccpi.framework.DataProcessor + :members: + +| + +:ref:`Return Home ` diff --git a/docs/source/index.rst b/docs/source/index.rst new file mode 100755 index 0000000..4342bec --- /dev/null +++ b/docs/source/index.rst @@ -0,0 +1,23 @@ +.. CCPi-Framework documentation master file, created by + sphinx-quickstart on Tue Mar 19 15:12:44 2019. + You can adapt this file completely to your liking, but it should at least + contain the root `toctree` directive. + +Welcome to CCPi-Framework's documentation! +========================================== + +.. toctree:: + :maxdepth: 2 + :caption: Contents: + :name: mastertoc + + + framework + optimisation + +Indices and tables +================== + +* :ref:`genindex` +* :ref:`modindex` +* :ref:`search` diff --git a/docs/source/optimisation.md b/docs/source/optimisation.md new file mode 100755 index 0000000..a3b9039 --- /dev/null +++ b/docs/source/optimisation.md @@ -0,0 +1,113 @@ +Optimisation framework +====================== + +This package allows rapid prototyping of optimisation-based +reconstruction problems, i.e. defining and solving different +optimization problems to enforce different properties on the +reconstructed image. + +Firstly, it provides an object-oriented framework for defining +mathematical operators and functions as well a collection of useful +example operators and functions. Both smooth and non-smooth functions +can be used. + +Further, it provides a number of high-level generic implementations of +optimisation algorithms to solve genericlly formulated optimisation +problems constructed from operator and function objects. + +The fundamental components are: + +- Operator: A class specifying a (currently linear) operator +- Function: A class specifying mathematical functions such as a least + squares data fidelity. +- Algorithm: Implementation of an iterative optimisation algorithm to + solve a particular generic optimisation problem. Algorithms are + iterable Python object which can be run in a for loop. Can be + stopped and warm restarted. + +Algorithm +--------- + +A number of generic algorithm implementations are provided including +Gradient Descent CGLS and FISTA. An algorithm is designed for a +particular generic optimisation problem accepts and number of Functions +and/or Operators as input to define a specific instance of the generic +optimisation problem to be solved. They are iterable objects which can +be run in a for loop. The user can provide a stopping criterion +different than the default max\_iteration. + +New algorithms can be easily created by extending the Algorithm class. +The user is required to implement only 4 methods: set\_up, \_\_init\_\_, +update and update\_objective. + +- `set_up` and `__init__` are used to configure the algorithm +- `update` is the actual iteration updating the solution +- `update_objective` defines how the objective is calculated. + +For example, the implementation of the update of the Gradient Descent +algorithm to minimise a Function will only be: + +The `Algorithm` provides the infrastructure to continue iteration, to +access the values of the objective function in subsequent iterations, +the time for each iteration. + +::: {.autoclass members="" private-members="" special-members=""} +ccpi.optimisation.algorithms.Algorithm +::: + +::: {.autoclass members=""} +ccpi.optimisation.algorithms.GradientDescent +::: + +::: {.autoclass members=""} +ccpi.optimisation.algorithms.CGLS +::: + +::: {.autoclass members=""} +ccpi.optimisation.algorithms.FISTA +::: + +Operator +-------- + +The two most important methods are `direct` and `adjoint` methods that +describe the result of applying the operator, and its adjoint +respectively, onto a compatible `DataContainer` input. The output is +another `DataContainer` object or subclass hereof. An important special +case is to represent the tomographic forward and backprojection +operations. + +::: {.autoclass members=""} +ccpi.optimisation.operators.Operator +::: + +::: {.autoclass members=""} +ccpi.optimisation.operators.LinearOperator +::: + +::: {.autoclass members=""} +ccpi.optimisation.operators.ScaledOperator +::: + +Function +-------- + +A `Function` represents a mathematical function of one or more inputs +and is intended to accept `DataContainers` as input as well as any +additional parameters. + +Fixed parameters can be passed in during the creation of the function +object. The methods of the function reflect the properties of it, for +example, if the function represented is differentiable the function +should contain a method `gradient` which should return the gradient of +the function evaluated at an input point. If the function is not +differentiable but allows a simple proximal operator, the method +`proximal` should return the proximal operator evaluated at an input +point. The function value is evaluated by calling the function itself, +e.g. `f(x)` for a `Function f` and input point `x`. + +::: {.autoclass members=""} +ccpi.optimisation.functions.Function +::: + +`Return Home `{.interpreted-text role="ref"} diff --git a/docs/source/optimisation.rst b/docs/source/optimisation.rst new file mode 100644 index 0000000..7c1eda6 --- /dev/null +++ b/docs/source/optimisation.rst @@ -0,0 +1,92 @@ +Optimisation framework +********************** +This package allows rapid prototyping of optimisation-based reconstruction problems, i.e. defining and solving different optimization problems to enforce different properties on the reconstructed image. + +Firstly, it provides an object-oriented framework for defining mathematical operators and functions as well a collection of useful example operators and functions. Both smooth and non-smooth functions can be used. + +Further, it provides a number of high-level generic implementations of optimisation algorithms to solve genericlly formulated optimisation problems constructed from operator and function objects. + +The fundamental components are: + ++ Operator: A class specifying a (currently linear) operator ++ Function: A class specifying mathematical functions such as a least squares data fidelity. ++ Algorithm: Implementation of an iterative optimisation algorithm to solve a particular generic optimisation problem. Algorithms are iterable Python object which can be run in a for loop. Can be stopped and warm restarted. + +Algorithm +========= + +A number of generic algorithm implementations are provided including +Gradient Descent CGLS and FISTA. An algorithm is designed for a +particular generic optimisation problem accepts and number of +Functions and/or Operators as input to define a specific instance of +the generic optimisation problem to be solved. +They are iterable objects which can be run in a for loop. +The user can provide a stopping criterion different than the default max_iteration. + +New algorithms can be easily created by extending the Algorithm class. The user is required to implement only 4 methods: set_up, __init__, update and update_objective. + ++ :code:`set_up` and :code:`__init__` are used to configure the algorithm ++ :code:`update` is the actual iteration updating the solution ++ :code:`update_objective` defines how the objective is calculated. + +For example, the implementation of the update of the Gradient Descent +algorithm to minimise a Function will only be: + +.. code-block :: python + + def update(self): + self.x += -self.rate * self.objective_function.gradient(self.x) + def update_objective(self): + self.loss.append(self.objective_function(self.x)) + +The :code:`Algorithm` provides the infrastructure to continue iteration, to access the values of the objective function in subsequent iterations, the time for each iteration. + +.. autoclass:: ccpi.optimisation.algorithms.Algorithm + :members: + :private-members: + :special-members: +.. autoclass:: ccpi.optimisation.algorithms.GradientDescent + :members: +.. autoclass:: ccpi.optimisation.algorithms.CGLS + :members: +.. autoclass:: ccpi.optimisation.algorithms.FISTA + :members: + +Operator +======== +The two most important methods are :code:`direct` and :code:`adjoint` +methods that describe the result of applying the operator, and its +adjoint respectively, onto a compatible :code:`DataContainer` input. +The output is another :code:`DataContainer` object or subclass +hereof. An important special case is to represent the tomographic +forward and backprojection operations. + +.. autoclass:: ccpi.optimisation.operators.Operator + :members: +.. autoclass:: ccpi.optimisation.operators.LinearOperator + :members: +.. autoclass:: ccpi.optimisation.operators.ScaledOperator + :members: + +Function +======== + +A :code:`Function` represents a mathematical function of one or more inputs +and is intended to accept :code:`DataContainers` as input as well as any +additional parameters. + +Fixed parameters can be passed in during the creation of the function object. +The methods of the function reflect the properties of it, for example, if the function +represented is differentiable the function should contain a method :code:`gradient` +which should return the gradient of the function evaluated at an input point. +If the function is not differentiable but allows a simple proximal operator, +the method :code:`proximal` should return the proximal operator evaluated at an +input point. The function value is evaluated by calling the function itself, +e.g. :code:`f(x)` for a :code:`Function f` and input point :code:`x`. + + +.. autoclass:: ccpi.optimisation.functions.Function + :members: + + +:ref:`Return Home ` -- cgit v1.2.3 From cdacb8f7421acfa7df6178843478ba23fe5d840f Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Tue, 16 Apr 2019 16:29:28 +0100 Subject: added KullbackLeibler out implementation and test --- .../ccpi/optimisation/functions/KullbackLeibler.py | 52 +++++++++++++++----- Wrappers/Python/test/test_functions.py | 55 ++++++++++------------ 2 files changed, 65 insertions(+), 42 deletions(-) diff --git a/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py b/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py index 18af154..a925f6d 100644 --- a/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py +++ b/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py @@ -43,20 +43,35 @@ class KullbackLeibler(Function): return numpy.sum( x.as_array() - self.b.as_array() * numpy.log(self.sum_value.as_array())) - def gradient(self, x): - - #TODO Division check - return 1 - self.b/(x + self.bnoise) + def gradient(self, x, out=None): + if out is None: + #TODO Division check + return 1 - self.b/(x + self.bnoise) + else: + x.add(self.bnoise, out=out) + self.b.divide(out, out=out) + out *= -1 + out += 1 def convex_conjugate(self, x, out=None): pass def proximal(self, x, tau, out=None): - - z = x + tau * self.bnoise - return (z + 1) - ((z-1)**2 + 4 * tau * self.b).sqrt() - - + if out is None: + z = x + tau * self.bnoise + return (z + 1) - ((z-1)**2 + 4 * tau * self.b).sqrt() + else: + z_m = x + tau * self.bnoise - 1 + self.b.multiply(4*tau, out=out) + z_m.multiply(z_m, out=z_m) + out += z_m + out.sqrt(out=out) + # z = z_m + 2 + z_m.sqrt(out=z_m) + z_m += 2 + out *= -1 + out += z_m + def proximal_conjugate(self, x, tau, out=None): pass @@ -67,15 +82,28 @@ if __name__ == '__main__': N, M = 2,3 ig = ImageGeometry(N, M) - data = ImageData(numpy.random.randint(-10, 100, size=(M, N))) - x = ImageData(numpy.random.randint(-10, 100, size=(M, N))) + #data = ImageData(numpy.random.randint(-10, 100, size=(M, N))) + #x = ImageData(numpy.random.randint(-10, 100, size=(M, N))) + #bnoise = ImageData(numpy.random.randint(-100, 100, size=(M, N))) + data = ig.allocate(ImageGeometry.RANDOM_INT) + x = ig.allocate(ImageGeometry.RANDOM_INT) + bnoise = ig.allocate(ImageGeometry.RANDOM_INT) - bnoise = ImageData(numpy.random.randint(-100, 100, size=(M, N))) + out = ig.allocate() f = KullbackLeibler(data, bnoise=bnoise) print(f.sum_value) print(f(x)) + grad = f.gradient(x) + f.gradient(x, out=out) + numpy.testing.assert_array_equal(grad.as_array(), out.as_array()) + + prox = f.proximal(x,1.2) + #print(grad.as_array()) + f.proximal(x, 1.2, out=out) + numpy.testing.assert_array_equal(prox.as_array(), out.as_array()) + diff --git a/Wrappers/Python/test/test_functions.py b/Wrappers/Python/test/test_functions.py index b428c12..7d68019 100644 --- a/Wrappers/Python/test/test_functions.py +++ b/Wrappers/Python/test/test_functions.py @@ -9,7 +9,7 @@ Created on Sat Mar 2 19:24:37 2019 import numpy as np #from ccpi.optimisation.funcs import Function -from ccpi.optimisation.functions import Function +from ccpi.optimisation.functions import Function, KullbackLeibler from ccpi.framework import DataContainer, ImageData, ImageGeometry from ccpi.optimisation.operators import Identity from ccpi.optimisation.operators import BlockOperator @@ -341,32 +341,27 @@ class TestFunction(unittest.TestCase): f_no_scaled.proximal_conjugate(U, 1, out=z3) self.assertBlockDataContainerEqual(z3,z1) -# -# f1 = L2NormSq(alpha=1, b=noisy_data) -# print(f1(noisy_data)) -# -# f2 = L2NormSq(alpha=5, b=noisy_data).composition_with(op2) -# print(f2(noisy_data)) -# -# print(f1.gradient(noisy_data).as_array()) -# print(f2.gradient(noisy_data).as_array()) -## -# print(f1.proximal(noisy_data,1).as_array()) -# print(f2.proximal(noisy_data,1).as_array()) -# -# -# f3 = mixed_L12Norm(alpha = 1).composition_with(op1) -# print(f3(noisy_data)) -# -# print(ImageData(op1.direct(noisy_data).power(2).sum(axis=0)).sqrt().sum()) -# -# print( 5*(op2.direct(d) - noisy_data).power(2).sum(), f2(d)) -# -# from functions import mixed_L12Norm as mixed_L12Norm_old -# -# print(mixed_L12Norm_old(op1,None,alpha)(noisy_data)) - - - # - - + def test_KullbackLeibler(self): + N, M = 2,3 + ig = ImageGeometry(N, M) + #data = ImageData(numpy.random.randint(-10, 100, size=(M, N))) + #x = ImageData(numpy.random.randint(-10, 100, size=(M, N))) + #bnoise = ImageData(numpy.random.randint(-100, 100, size=(M, N))) + data = ig.allocate(ImageGeometry.RANDOM_INT) + x = ig.allocate(ImageGeometry.RANDOM_INT) + bnoise = ig.allocate(ImageGeometry.RANDOM_INT) + + out = ig.allocate() + + f = KullbackLeibler(data, bnoise=bnoise) + print(f.sum_value) + + print(f(x)) + grad = f.gradient(x) + f.gradient(x, out=out) + numpy.testing.assert_array_equal(grad.as_array(), out.as_array()) + + prox = f.proximal(x,1.2) + #print(grad.as_array()) + f.proximal(x, 1.2, out=out) + numpy.testing.assert_array_equal(prox.as_array(), out.as_array()) \ No newline at end of file -- cgit v1.2.3 From 38cbce611ebe35cf803ccfc269d13e874fcf7bec Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Tue, 16 Apr 2019 16:45:40 +0100 Subject: fix L2NormSquared --- Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py | 4 ++-- Wrappers/Python/test/test_functions.py | 3 ++- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py b/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py index 1946d67..2bb4ca7 100644 --- a/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py +++ b/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py @@ -93,8 +93,8 @@ class L2NormSquared(Function): if self.b is None: return x/(1+2*tau) else: - tmp = x - tmp -= self.b + tmp = x.subtract(self.b) + #tmp -= self.b tmp /= (1+2*tau) tmp += self.b return tmp diff --git a/Wrappers/Python/test/test_functions.py b/Wrappers/Python/test/test_functions.py index 7d68019..530cd21 100644 --- a/Wrappers/Python/test/test_functions.py +++ b/Wrappers/Python/test/test_functions.py @@ -99,6 +99,7 @@ class TestFunction(unittest.TestCase): def test_L2NormSquared(self): # TESTS for L2 and scalar * L2 + print ("Test L2NormSquared") M, N, K = 2,3,5 ig = ImageGeometry(voxel_num_x=M, voxel_num_y = N, voxel_num_z = K) @@ -324,7 +325,7 @@ class TestFunction(unittest.TestCase): a1 = f_no_scaled(U) a2 = f_scaled(U) - self.assertNumpyArrayAlmostEqual(a1.as_array(),a2.as_array()) + self.assertNumpyArrayAlmostEqual(a1,a2) tmp = [ el**2 for el in U.containers ] -- cgit v1.2.3 From a8bdbf487d144fa5c9e549a0acea44aeadb41739 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Tue, 16 Apr 2019 16:52:41 +0100 Subject: removed comments --- Wrappers/Python/test/test_functions.py | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/Wrappers/Python/test/test_functions.py b/Wrappers/Python/test/test_functions.py index 530cd21..33a69d7 100644 --- a/Wrappers/Python/test/test_functions.py +++ b/Wrappers/Python/test/test_functions.py @@ -343,11 +343,9 @@ class TestFunction(unittest.TestCase): self.assertBlockDataContainerEqual(z3,z1) def test_KullbackLeibler(self): + print ("test_KullbackLeibler") N, M = 2,3 ig = ImageGeometry(N, M) - #data = ImageData(numpy.random.randint(-10, 100, size=(M, N))) - #x = ImageData(numpy.random.randint(-10, 100, size=(M, N))) - #bnoise = ImageData(numpy.random.randint(-100, 100, size=(M, N))) data = ig.allocate(ImageGeometry.RANDOM_INT) x = ig.allocate(ImageGeometry.RANDOM_INT) bnoise = ig.allocate(ImageGeometry.RANDOM_INT) @@ -355,14 +353,11 @@ class TestFunction(unittest.TestCase): out = ig.allocate() f = KullbackLeibler(data, bnoise=bnoise) - print(f.sum_value) - print(f(x)) grad = f.gradient(x) f.gradient(x, out=out) numpy.testing.assert_array_equal(grad.as_array(), out.as_array()) prox = f.proximal(x,1.2) - #print(grad.as_array()) f.proximal(x, 1.2, out=out) numpy.testing.assert_array_equal(prox.as_array(), out.as_array()) \ No newline at end of file -- cgit v1.2.3 From 956fcf2544f7cbafba7ceb139fdcbd3ec2ef13fe Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Wed, 17 Apr 2019 10:34:40 +0100 Subject: add implementation of KullbackLeibler with tests and memopt --- .../ccpi/optimisation/functions/KullbackLeibler.py | 109 +++++++++++++++------ Wrappers/Python/test/test_functions.py | 6 +- 2 files changed, 82 insertions(+), 33 deletions(-) diff --git a/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py b/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py index a925f6d..e7e41f7 100644 --- a/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py +++ b/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py @@ -19,44 +19,94 @@ import numpy from ccpi.optimisation.functions import Function -from ccpi.optimisation.functions.ScaledFunction import ScaledFunction -from ccpi.framework import DataContainer, ImageData, ImageGeometry +from ccpi.optimisation.functions.ScaledFunction import ScaledFunction +from ccpi.framework import ImageData, ImageGeometry +import functools class KullbackLeibler(Function): - def __init__(self,data,**kwargs): + ''' Assume that data > 0 + + ''' + + def __init__(self,data, **kwargs): super(KullbackLeibler, self).__init__() self.b = data self.bnoise = kwargs.get('bnoise', 0) - self.sum_value = self.b + self.bnoise - if (self.sum_value.as_array()<0).any(): - self.sum_value = numpy.inf def __call__(self, x): + # TODO check + + self.sum_value = x + self.bnoise + if (self.sum_value.as_array()<0).any(): + self.sum_value = numpy.inf + if self.sum_value==numpy.inf: return numpy.inf else: - return numpy.sum( x.as_array() - self.b.as_array() * numpy.log(self.sum_value.as_array())) + tmp = self.sum_value + #tmp.fill( numpy.log(tmp.as_array()) ) + self.log(tmp) + return (x - self.b * tmp ).sum() + +# return numpy.sum( x.as_array() - self.b.as_array() * numpy.log(self.sum_value.as_array())) + def log(self, datacontainer): + '''calculates the in-place log of the datacontainer''' + if not functools.reduce(lambda x,y: x and y>0, + datacontainer.as_array().ravel(), True): + raise ValueError('KullbackLeibler. Cannot calculate log of negative number') + datacontainer.fill( numpy.log(datacontainer.as_array()) ) def gradient(self, x, out=None): + + #TODO Division check if out is None: - #TODO Division check return 1 - self.b/(x + self.bnoise) else: x.add(self.bnoise, out=out) self.b.divide(out, out=out) + out.subtract(1, out=out) out *= -1 - out += 1 - - def convex_conjugate(self, x, out=None): - pass + + def convex_conjugate(self, x): + + tmp = self.b/( 1 - x ) + self.log(tmp) + return (self.b * ( tmp - 1 ) - self.bnoise * (x - 1)).sum() +# return self.b * ( ImageData(numpy.log(self.b/(1-x)) - 1 )) - self.bnoise * (x - 1) def proximal(self, x, tau, out=None): + + if out is None: + return 0.5 *( (x - self.bnoise - tau) + ( (x + self.bnoise - tau)**2 + 4*tau*self.b ) .sqrt() ) + else: + tmp = 0.5 *( (x - self.bnoise - tau) + + ( (x + self.bnoise - tau)**2 + 4*tau*self.b ) .sqrt() + ) + x.add(self.bnoise, out=out) + out -= tau + out *= out + tmp = self.b * (4 * tau) + out.add(tmp, out=out) + out.sqrt(out=out) + + x.subtract(self.bnoise, out=tmp) + tmp -= tau + + out += tmp + + out *= 0.5 + + + + def proximal_conjugate(self, x, tau, out=None): + + if out is None: z = x + tau * self.bnoise return (z + 1) - ((z-1)**2 + 4 * tau * self.b).sqrt() @@ -71,9 +121,17 @@ class KullbackLeibler(Function): z_m += 2 out *= -1 out += z_m - - def proximal_conjugate(self, x, tau, out=None): - pass + + + def __rmul__(self, scalar): + + ''' Multiplication of L2NormSquared with a scalar + + Returns: ScaledFunction + + ''' + + return ScaledFunction(self, scalar) @@ -82,29 +140,16 @@ if __name__ == '__main__': N, M = 2,3 ig = ImageGeometry(N, M) - #data = ImageData(numpy.random.randint(-10, 100, size=(M, N))) - #x = ImageData(numpy.random.randint(-10, 100, size=(M, N))) - #bnoise = ImageData(numpy.random.randint(-100, 100, size=(M, N))) - data = ig.allocate(ImageGeometry.RANDOM_INT) - x = ig.allocate(ImageGeometry.RANDOM_INT) - bnoise = ig.allocate(ImageGeometry.RANDOM_INT) + data = ImageData(numpy.random.randint(-10, 100, size=(M, N))) + x = ImageData(numpy.random.randint(-10, 100, size=(M, N))) - out = ig.allocate() + bnoise = ImageData(numpy.random.randint(-100, 100, size=(M, N))) f = KullbackLeibler(data, bnoise=bnoise) print(f.sum_value) print(f(x)) - grad = f.gradient(x) - f.gradient(x, out=out) - numpy.testing.assert_array_equal(grad.as_array(), out.as_array()) - - prox = f.proximal(x,1.2) - #print(grad.as_array()) - f.proximal(x, 1.2, out=out) - numpy.testing.assert_array_equal(prox.as_array(), out.as_array()) - - \ No newline at end of file + diff --git a/Wrappers/Python/test/test_functions.py b/Wrappers/Python/test/test_functions.py index 33a69d7..af419c7 100644 --- a/Wrappers/Python/test/test_functions.py +++ b/Wrappers/Python/test/test_functions.py @@ -360,4 +360,8 @@ class TestFunction(unittest.TestCase): prox = f.proximal(x,1.2) f.proximal(x, 1.2, out=out) - numpy.testing.assert_array_equal(prox.as_array(), out.as_array()) \ No newline at end of file + numpy.testing.assert_array_equal(prox.as_array(), out.as_array()) + + proxc = f.proximal_conjugate(x,1.2) + f.proximal_conjugate(x, 1.2, out=out) + numpy.testing.assert_array_equal(proxc.as_array(), out.as_array()) \ No newline at end of file -- cgit v1.2.3