Commit d28da9fe authored by Jakob Knollmueller's avatar Jakob Knollmueller

laplace operator adjusted for spaces keyword, critical filter still seems to work

parent 89c555e1
Pipeline #14362 passed with stage
in 5 minutes and 11 seconds
......@@ -3,6 +3,7 @@ setup.cfg
.idea
.DS_Store
*.pyc
*.html
# from https://github.com/github/gitignore/blob/master/Python.gitignore
# Byte-compiled / optimized / DLL files
......
......@@ -62,7 +62,7 @@ if __name__ == "__main__":
# Setting up power space
p_space = PowerSpace(h_space, logarithmic=True,
distribution_strategy=distribution_strategy)#, nbin=5)
distribution_strategy=distribution_strategy)
# Choosing the prior correlation structure and defining correlation operator
p_spec = (lambda k: (.5 / (k + 1) ** 3))
......@@ -140,7 +140,7 @@ if __name__ == "__main__":
m0 = map_energy.position
D0 = map_energy.curvature
# Initializing the power energy with updated parameters
power_energy = CriticalPowerEnergy(position=t0, m=m0, D=D0, sigma=100., samples=3)
power_energy = CriticalPowerEnergy(position=t0, m=m0, D=D0, sigma=10., samples=3)
(power_energy, convergence) = minimizer1(power_energy)
......
......@@ -3,6 +3,7 @@ import numpy as np
from nifty import Field,\
EndomorphicOperator,\
PowerSpace
import nifty.nifty_utilities as utilities
# def _irregular_nabla(x,k):
# #applies forward differences and does nonesense at the edge. Thus needs cutting
......@@ -73,21 +74,57 @@ class LaplaceOperator(EndomorphicOperator):
return False
def _times(self, t, spaces):
ret = t.val.copy()
ret = self._irregular_laplace(ret)
ret[0:2] = 0
ret[-1] = 0
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
ret[2:] *= np.sqrt(self.fwd_dist)[1:]
ret[0:2] = 0
ret[-1] = 0
ret = self._irregular_adj_laplace(ret)
return Field(self.domain, val=ret).weight(-1)
def _times(self, x, spaces):
spaces = utilities.cast_axis_to_tuple(spaces, len(x.domain))
if spaces is None:
# this case means that x lives on only one space, which is
# identical to the space in the domain of `self`. Otherwise the
# input check of LinearOperator would have failed.
axes = x.domain_axes[0]
else:
axes = x.domain_axes[spaces[0]]
axis = axes[0]
prefix = (slice(None),) * axis
fwd_dist = self.fwd_dist.reshape((1,)*axis + self.fwd_dist.shape)
positions = self.positions.reshape((1,)*axis + self.positions.shape)
ret = x.val.copy_empty()
x = x.val
ret[prefix + (slice(1,-1),)] = -(x[prefix + (slice(1,-1),)]
- x[prefix + (slice(0,-2),)]) / fwd_dist[prefix + (slice(0,-1),)] \
+ (x[prefix + (slice(2,None),)] - x[prefix + (slice(1,-1),)]) / fwd_dist[prefix + (slice(1, None),)]
ret[prefix + (slice(1,-1),)] /= positions[prefix + (slice(2,None),)] - positions[prefix + (slice(None,-2),)]
ret *= 2.
ret[prefix + (slice(0,2),)] = 0
ret[prefix + (slice(-1,-1),)] = 0
ret[prefix + (slice(2,None),)] *= np.sqrt(fwd_dist)[prefix + (slice(1,None),)]
return Field(self.domain, val=ret).weight(power=-0.5,spaces=spaces)
def _adjoint_times(self, x, spaces):
spaces = utilities.cast_axis_to_tuple(spaces, len(x.domain))
if spaces is None:
# this case means that x lives on only one space, which is
# identical to the space in the domain of `self`. Otherwise the
# input check of LinearOperator would have failed.
axes = x.domain_axes[0]
else:
axes = x.domain_axes[spaces[0]]
axis = axes[0]
prefix = (slice(None),) * axis
fwd_dist = self.fwd_dist.reshape((1,)*axis + self.fwd_dist.shape)
positions = self.positions.reshape((1,)*axis + self.positions.shape)
y = x.copy().weight(power=0.5).val
y[prefix + (slice(2,None),)] *= np.sqrt(fwd_dist)[prefix + (slice(1,None),)]
y[prefix + (slice(0,2),)] = 0
y[prefix + (slice(-1,-1),)] = 0
ret = y.copy_empty()
y[prefix + (slice(1,-1),)] /= positions[prefix + (slice(2,None),)] - positions[prefix + (slice(None,-2),)]
y *= 2
ret[prefix + (slice(1,-1),)] = -y[prefix + (slice(1,-1),)] / fwd_dist[prefix \
+ (slice(0,-1),)] - y[prefix
+ (slice(1,-1),)] / fwd_dist[prefix + (slice(1, None),)]
ret[prefix + (slice(0,-2),)] += y[prefix + (slice(1,-1),)] / fwd_dist[prefix + (slice(0,-1),)]
ret[prefix + (slice(2,None),)] += y[prefix + (slice(1,-1),)] / fwd_dist[prefix + (slice(1, None),)]
return Field(self.domain, val=ret).weight(-1,spaces=spaces)
def _irregular_laplace(self, x):
ret = np.zeros_like(x)
......
from nifty import EndomorphicOperator,\
PowerSpace
import nifty.nifty_utilities as utilities
from laplace_operator import LaplaceOperator
......@@ -39,7 +40,7 @@ class SmoothnessOperator(EndomorphicOperator):
self._sigma = sigma
self._Laplace = LaplaceOperator(domain=self._domain[0], logarithmic=logarithmic)
self._Laplace = LaplaceOperator(domain=domain, logarithmic=logarithmic)
......@@ -67,6 +68,7 @@ class SmoothnessOperator(EndomorphicOperator):
def self_adjoint(self):
return False
def _times(self, t, spaces):
res = self._Laplace.adjoint_times(self._Laplace.times(t))
return (1/self.sigma)**2*res
def _times(self, x, spaces):
res = self._Laplace.adjoint_times(self._Laplace.times(x,spaces), spaces)
return (1./self.sigma)**2*res
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment