diff --git a/demos/laplace_testing.py b/demos/laplace_testing.py
index be878e2feddcdf991bfef9f325ee3c098bea124b..3577720c8dcb69f6ed494a639dd0befb6a843a8b 100644
--- a/demos/laplace_testing.py
+++ b/demos/laplace_testing.py
@@ -58,13 +58,13 @@ if __name__ == "__main__":
# Setting up power space
p_space = PowerSpace(h_space, logarithmic=False,
- distribution_strategy=distribution_strategy, nbin=16)
+ distribution_strategy=distribution_strategy, nbin=70)
# Choosing the prior correlation structure and defining correlation operator
pow_spec = (lambda k: (.05 / (k + 1) ** 2))
# t = Field(p_space, val=pow_spec)
t= Field.from_random("normal", domain=p_space)
- lap = LogLaplaceOperator(p_space)
+ lap = LaplaceOperator(p_space)
T = SmoothnessOperator(p_space,sigma=1.)
test_energy = TestEnergy(t,T)
diff --git a/nifty/library/operator_library/__init__.py b/nifty/library/operator_library/__init__.py
index 64e43c39c6c21c4925e45a10fa295bfb8f91c3de..386f038783339c3582ac46ae9bd007f1d57bbebb 100644
--- a/nifty/library/operator_library/__init__.py
+++ b/nifty/library/operator_library/__init__.py
@@ -1,5 +1,4 @@
from wiener_filter_curvature import WienerFilterCurvature
from critical_power_curvature import CriticalPowerCurvature
from laplace_operator import LaplaceOperator
-from laplace_operator import LogLaplaceOperator
from smoothness_operator import SmoothnessOperator
\ No newline at end of file
diff --git a/nifty/library/operator_library/laplace_operator.py b/nifty/library/operator_library/laplace_operator.py
index 98934775da199f8bc2a4a4a1028af20b687a0f8b..b914cd7105271b6020b624cfbf54ecd343d7b572 100644
--- a/nifty/library/operator_library/laplace_operator.py
+++ b/nifty/library/operator_library/laplace_operator.py
@@ -3,98 +3,40 @@ import numpy as np
from nifty import Field,\
EndomorphicOperator,\
PowerSpace
-class LaplaceOperator(EndomorphicOperator):
- def __init__(self, domain,
- default_spaces=None):
- super(LaplaceOperator, self).__init__(default_spaces)
- if (domain is not None):
- if (not isinstance(domain, PowerSpace)):
- raise TypeError("The domain has to live over a PowerSpace")
- self._domain = self._parse_domain(domain)
-
- """
- input parameters:
- domain- can only live over one domain
- to do:
- correct implementation of D20 object
- """
-
- @property
- def target(self):
- return self._domain
-
- @property
- def domain(self):
- return self._domain
-
- @property
- def unitary(self):
- return False
-
- @property
- def symmetric(self):
- return False
-
- @property
- def self_adjoint(self):
- return False
-
- def _times(self, t, spaces):
- if t.val.distribution_strategy != 'not':
- # self.logger.warn("distribution_strategy should be 'not' but is %s"
- # % t.val.distribution_strategy)
- array = t.val.get_full_data()
- else:
- array = t.val.get_local_data(copy=False)
- ret = 2 * array
- ret -= np.roll(array, 1)
- ret -= np.roll(array, -1)
- ret[0:2] = 0
- ret[-1] = 0
- return Field(self.domain, val=ret)
-
- def _adjoint_times(self, t, spaces):
- t = t.copy().weight(1)
- t.val[1:] /= self.domain[0].kindex[1:]
- if t.val.distribution_strategy != 'not':
- # self.logger.warn("distribution_strategy should be 'not' but is %s"
- # % t.val.distribution_strategy)
- array = t.val.get_full_data()
- else:
- array = t.val.get_local_data(copy=False)
- ret = np.copy(array)
- ret[0:2] = 0
- ret[-1] = 0
- ret = 2 * ret - np.roll(ret, 1) - np.roll(ret, -1)
- result = Field(self.domain, val=ret).weight(-1)
- return result
+# def _irregular_nabla(x,k):
+# #applies forward differences and does nonesense at the edge. Thus needs cutting
+# y = -x
+# y[:-1] += x[1:]
+# y[1:-1] /= - k[1:-1] + k[2:]
+# return y
+#
+# def _irregular_adj_nabla(z, k):
+# #applies backwards differences*(-1) and does nonesense at the edge. Thus needs cutting
+# x = z.copy()
+# x[1:-1] /= - k[1:-1] + k[2:]
+# y = -x
+# y[1:] += x[:-1]
+# return y
-def _irregular_nabla(x,k):
- #applies forward differences and does nonesense at the edge. Thus needs cutting
- y = -x
- y[:-1] += x[1:]
- y[1:-1] /= - k[1:-1] + k[2:]
- return y
-
-def _irregular_adj_nabla(z, k):
- #applies backwards differences*(-1) and does nonesense at the edge. Thus needs cutting
- x = z.copy()
- x[1:-1] /= - k[1:-1] + k[2:]
- y = -x
- y[1:] += x[:-1]
- return y
-
-
-class LogLaplaceOperator(EndomorphicOperator):
+class LaplaceOperator(EndomorphicOperator):
def __init__(self, domain,
- default_spaces=None):
- super(LogLaplaceOperator, self).__init__(default_spaces)
+ default_spaces=None, logarithmic = True):
+ super(LaplaceOperator, self).__init__(default_spaces)
if (domain is not None):
if (not isinstance(domain, PowerSpace)):
raise TypeError("The domain has to live over a PowerSpace")
self._domain = self._parse_domain(domain)
+ if logarithmic :
+ self.positions = self.domain[0].kindex.copy()
+ self.positions[1:] = np.log(self.positions[1:])
+ self.positions[0] = -1.
+ else :
+ self.positions = self.domain[0].kindex.copy()
+ self.positions[0] = -1
+
+ self.fwd_dist = self.positions[1:] - self.positions[:-1]
"""
input parameters:
@@ -125,27 +67,35 @@ class LogLaplaceOperator(EndomorphicOperator):
def _times(self, t, spaces):
- l_vec = self.domain[0].kindex.copy()
- l_vec[1:] = np.log(l_vec[1:])
- l_vec[0] = -1.
ret = t.val.copy()
- # ret[2:] *= np.sqrt(np.log(k_vec)-np.log(np.roll(k_vec,1)))[2:]
- ret = _irregular_nabla(ret,l_vec)
- ret = _irregular_adj_nabla(ret,l_vec)
+ ret = self._irregular_laplace(ret)
ret[0:2] = 0
ret[-1] = 0
- ret[2:] *= np.sqrt(l_vec-np.roll(l_vec,1))[2:]
+ ret[2:] *= np.sqrt(self.fwd_dist)[1:]
return Field(self.domain, val=ret).weight(power=-0.5)
def _adjoint_times(self, t, spaces):
ret = t.copy().weight(power=0.5).val
- l_vec = self.domain[0].kindex.copy()
- l_vec[1:] = np.log(l_vec[1:])
- l_vec[0] = -1.
- ret[2:] *= np.sqrt(l_vec-np.roll(l_vec,1))[2:]
+ ret[2:] *= np.sqrt(self.fwd_dist)[1:]
ret[0:2] = 0
ret[-1] = 0
- ret = _irregular_nabla(ret,l_vec)
- ret = _irregular_adj_nabla(ret,l_vec)
- # ret[2:] *= np.sqrt(np.log(k_vec)-np.log(np.roll(k_vec,1)))[2:]
- return Field(self.domain, val=ret).weight(-1)
\ No newline at end of file
+ ret = self._irregular_adj_laplace(ret)
+ return Field(self.domain, val=ret).weight(-1)
+
+ def _irregular_laplace(self, x):
+ ret = np.zeros_like(x)
+ ret[1:-1] = -(x[1:-1] - x[0:-2]) / self.fwd_dist[:-1] \
+ + (x[2:] - x[1:-1]) / self.fwd_dist[1:]
+ ret[1:-1] /= self.positions[2:] - self.positions[:-2]
+ ret *= 2.
+ return ret
+
+ def _irregular_adj_laplace(self, x):
+ ret = np.zeros_like(x)
+ y = x.copy()
+ y[1:-1] /= self.positions[2:] - self.positions[:-2]
+ y *= 2
+ ret[1:-1] = -y[1:-1] / self.fwd_dist[:-1] - y[1:-1] / self.fwd_dist[1:]
+ ret[0:-2] += y[1:-1] / self.fwd_dist[:-1]
+ ret[2:] += y[1:-1] / self.fwd_dist[1:]
+ return ret
diff --git a/nifty/library/operator_library/smoothness_operator.py b/nifty/library/operator_library/smoothness_operator.py
index 24e285a73b582086fbd00108767fe1bf745f44bf..fca9d91e2fa4fd32a5ff604e201ef096fd77eb1d 100644
--- a/nifty/library/operator_library/smoothness_operator.py
+++ b/nifty/library/operator_library/smoothness_operator.py
@@ -1,7 +1,7 @@
from nifty import EndomorphicOperator,\
PowerSpace
-from laplace_operator import LogLaplaceOperator
+from laplace_operator import LaplaceOperator
class SmoothnessOperator(EndomorphicOperator):
@@ -21,7 +21,7 @@ class SmoothnessOperator(EndomorphicOperator):
self._sigma = sigma
- self._Laplace = LogLaplaceOperator(domain=self._domain[0])
+ self._Laplace = LaplaceOperator(domain=self._domain[0])
"""
SmoothnessOperator