summaryrefslogtreecommitdiffstats
path: root/Wrappers
diff options
context:
space:
mode:
authorVaggelis Papoutsellis <22398586+epapoutsellis@users.noreply.github.com>2019-10-30 16:06:04 +0000
committerEdoardo Pasca <edo.paskino@gmail.com>2019-10-30 16:06:04 +0000
commit9bd69b755f90084db836479005733cefc170fcb8 (patch)
tree6625dd10e0d6bd1cd79de400f96db99e837aa40f /Wrappers
parenta4d3aa5d160f695acedef2d4fa39013ddd8d9bef (diff)
downloadframework-9bd69b755f90084db836479005733cefc170fcb8.tar.gz
framework-9bd69b755f90084db836479005733cefc170fcb8.tar.bz2
framework-9bd69b755f90084db836479005733cefc170fcb8.tar.xz
framework-9bd69b755f90084db836479005733cefc170fcb8.zip
fix MixedL21Norm (#417)
* fix MixedL21Norm * merged master * re added future imports
Diffstat (limited to 'Wrappers')
-rwxr-xr-xWrappers/Python/ccpi/optimisation/functions/MixedL21Norm.py162
1 files changed, 87 insertions, 75 deletions
diff --git a/Wrappers/Python/ccpi/optimisation/functions/MixedL21Norm.py b/Wrappers/Python/ccpi/optimisation/functions/MixedL21Norm.py
index 378cbda..55e6e53 100755
--- a/Wrappers/Python/ccpi/optimisation/functions/MixedL21Norm.py
+++ b/Wrappers/Python/ccpi/optimisation/functions/MixedL21Norm.py
@@ -1,22 +1,20 @@
# -*- coding: utf-8 -*-
-# Copyright 2019 Science Technology Facilities Council
-# Copyright 2019 University of Manchester
-#
-# This work is part of the Core Imaging Library developed by Science Technology
-# Facilities Council and University of Manchester
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-# http://www.apache.org/licenses/LICENSE-2.0.txt
-#
-# Unless required by applicable law or agreed to in writing, software
-# distributed under the License is distributed on an "AS IS" BASIS,
-# 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.
+# CCP in Tomographic Imaging (CCPi) Core Imaging Library (CIL).
+
+# Copyright 2017 UKRI-STFC
+# Copyright 2017 University of Manchester
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+# http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# 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.
from __future__ import absolute_import
from __future__ import division
@@ -25,17 +23,15 @@ from __future__ import unicode_literals
from ccpi.optimisation.functions import Function, ScaledFunction
from ccpi.framework import BlockDataContainer
+import numpy as np
import functools
-import numpy
class MixedL21Norm(Function):
- r'''MixedL21Norm: .. math:: f(x) = ||x||_{2,1} = \int \|x\|_{2} dx
-
- where x is a vector/tensor vield
-
+ '''
+ f(x) = ||x||_{2,1} = \sum |x|_{2}
'''
def __init__(self, **kwargs):
@@ -45,13 +41,15 @@ class MixedL21Norm(Function):
def __call__(self, x):
- '''Evaluates MixedL21Norm at point x
+ ''' Evaluates L2,1Norm at point x
- :param: x: is a BlockDataContainer
+ :param: x is a BlockDataContainer
+
'''
if not isinstance(x, BlockDataContainer):
raise ValueError('__call__ expected BlockDataContainer, got {}'.format(type(x)))
- tmp = x.get_item(0) * 0
+
+ tmp = x.get_item(0) * 0.
for el in x.containers:
tmp += el.power(2.)
return tmp.sqrt().sum()
@@ -62,11 +60,8 @@ class MixedL21Norm(Function):
def convex_conjugate(self,x):
- r'''Convex conjugate of of MixedL21Norm:
-
- Indicator function of .. math:: ||\cdot||_{2, \infty}
- which is either 0 if .. math:: ||x||_{2, \infty}<1 or \infty
-
+ ''' This is the Indicator function of ||\cdot||_{2, \infty}
+ which is either 0 if ||x||_{2, \infty} or \infty
'''
return 0.0
@@ -74,68 +69,59 @@ class MixedL21Norm(Function):
def proximal(self, x, tau, out=None):
- r'''Proximal operator of MixedL21Norm at x:
-
- .. math:: prox_{\tau * f(x)
- '''
- pass
-
- def proximal_conjugate(self, x, tau, out=None):
-
- r'''Proximal operator of the convex conjugate of MixedL21Norm at x:
-
- .. math:: prox_{\tau * f^{*}}(x)
+ if out is None:
+
+ tmp = sum([ el*el for el in x.containers]).sqrt()
+ res = (tmp - tau).maximum(0.0) * x/tmp
+ return res
+
+ else:
+
+ tmp = functools.reduce(lambda a,b: a + b*b, x.containers, x.get_item(0) * 0 ).sqrt()
+ res = (tmp - tau).maximum(0.0) * x/tmp
- '''
+ for el in res.containers:
+ el.as_array()[np.isnan(el.as_array())]=0
+ out.fill(res)
+
+
+ def proximal_conjugate(self, x, tau, out=None):
+
if out is None:
- # tmp = [ el*el for el in x.containers]
- # res = sum(tmp).sqrt().maximum(1.0)
- # frac = [el/res for el in x.containers]
- # return BlockDataContainer(*frac)
- tmp = x.get_item(0) * 0
- for el in x.containers:
- tmp += el.power(2.)
- tmp.sqrt(out=tmp)
- tmp.maximum(1.0, out=tmp)
- frac = [ el.divide(tmp) for el in x.containers ]
+ tmp = x.get_item(0) * 0
+ for el in x.containers:
+ tmp += el.power(2.)
+ tmp.sqrt(out=tmp)
+ tmp.maximum(1.0, out=tmp)
+ frac = [ el.divide(tmp) for el in x.containers ]
return BlockDataContainer(*frac)
-
-
+
+
else:
res1 = functools.reduce(lambda a,b: a + b*b, x.containers, x.get_item(0) * 0 )
- if False:
- res = res1.sqrt().maximum(1.0)
- x.divide(res, out=out)
- else:
- res1.sqrt(out=res1)
- res1.maximum(1.0, out=res1)
- x.divide(res1, out=out)
-
+ res1.sqrt(out=res1)
+ res1.maximum(1.0, out=res1)
+ x.divide(res1, out=out)
+
def __rmul__(self, scalar):
- '''Multiplication of MixedL21Norm with a scalar
-
- Returns: ScaledFunction
+ ''' Multiplication of MixedL21Norm with a scalar
+
+ Returns: ScaledFunction
'''
return ScaledFunction(self, scalar)
-def sqrt_maximum(x, a):
- y = numpy.sqrt(x)
- if y >= a:
- return y
- else:
- return a
#
if __name__ == '__main__':
- M, N, K = 2,3,5
- from ccpi.framework import BlockGeometry
+ M, N, K = 2,3,50
+ from ccpi.framework import BlockGeometry, ImageGeometry
import numpy
ig = ImageGeometry(M, N)
@@ -145,8 +131,9 @@ if __name__ == '__main__':
U = BG.allocate('random_int')
# Define no scale and scaled
+ alpha = 0.5
f_no_scaled = MixedL21Norm()
- f_scaled = 0.5 * MixedL21Norm()
+ f_scaled = alpha * MixedL21Norm()
# call
@@ -174,11 +161,36 @@ if __name__ == '__main__':
numpy.testing.assert_array_almost_equal(res_no_out[1].as_array(), \
res_out[1].as_array(), decimal=4)
-#
+ tau = 0.4
+ d1 = f_scaled.proximal(U, tau)
+
+ tmp = (U.get_item(0)**2 + U.get_item(1)**2).sqrt()
+
+ d2 = (tmp - alpha*tau).maximum(0) * U/tmp
+ numpy.testing.assert_array_almost_equal(d1.get_item(0).as_array(), \
+ d2.get_item(0).as_array(), decimal=4)
+ numpy.testing.assert_array_almost_equal(d1.get_item(1).as_array(), \
+ d2.get_item(1).as_array(), decimal=4)
+ out1 = BG.allocate('random_int')
+
+
+ f_scaled.proximal(U, tau, out = out1)
+
+ numpy.testing.assert_array_almost_equal(out1.get_item(0).as_array(), \
+ d1.get_item(0).as_array(), decimal=4)
+ numpy.testing.assert_array_almost_equal(out1.get_item(1).as_array(), \
+ d1.get_item(1).as_array(), decimal=4)
+#
+
+
+
+
+
+ \ No newline at end of file