Commit f468ec84 authored by Martin Reinecke's avatar Martin Reinecke
Browse files

cosmetics

parent 76c0839f
from ...energies.energy import Energy
from ...operators.smoothness_operator import SmoothnessOperator
from ...operators.diagonal_operator import DiagonalOperator
from ...operators.linear_operator import LinearOperator
from ...operators.power_projection_operator import PowerProjection
from . import CriticalPowerCurvature
from ...memoization import memo
......@@ -72,7 +71,8 @@ class CriticalPowerEnergy(Energy):
strength=smoothness_prior,
logarithmic=logarithmic)
self.rho = self.position.domain[0].rho
self.P = PowerProjection(domain=self.m.domain,target=self.position.domain)
self.P = PowerProjection(domain=self.m.domain,
target=self.position.domain)
self._w = w if w is not None else None
if inverter is None:
preconditioner = DiagonalOperator(self._theta.domain,
......@@ -80,7 +80,7 @@ class CriticalPowerEnergy(Energy):
copy=False)
inverter = ConjugateGradient(preconditioner=preconditioner)
self._inverter = inverter
self.one = Field(self.position.domain,val=1.)
self.one = Field(self.position.domain, val=1.)
self.constants = self.one/2. + self.alpha - 1
@property
......@@ -152,7 +152,6 @@ class CriticalPowerEnergy(Energy):
def _theta(self):
return exp(-self.position) * (self.q + self.w / 2.)
@property
@memo
def _Tt(self):
......
......@@ -39,3 +39,5 @@ from .response_operator import ResponseOperator
from .laplace_operator import LaplaceOperator
from .smoothness_operator import SmoothnessOperator
from .power_projection_operator import PowerProjection
from .power_projection_operator import PowerProjection
\ No newline at end of file
from .power_projection_operator import PowerProjection
from ... import Field,\
FieldArray
from ... import Field
from ..linear_operator import LinearOperator
from ...spaces.power_space import PowerSpace
class PowerProjection(LinearOperator):
def __init__(self, domain, target, default_spaces=None):
self._domain = self._parse_domain(domain)
self._target = self._parse_domain(target)
if len(self._domain)!=1 or len(self._target)!=1:
if len(self._domain) != 1 or len(self._target) != 1:
raise ValueError("Operator only works over one space")
if not self._domain[0].harmonic:
raise ValueError("domain must be a harmonic space")
......@@ -16,18 +16,18 @@ class PowerProjection(LinearOperator):
self.pindex = self.target[0].pindex
super(PowerProjection, self).__init__(default_spaces)
def _times(self,x,spaces):
def _times(self, x, spaces):
if spaces is None:
spaces = 0
projected_x = self.pindex.bincount(
weights=x.weight(1,spaces=spaces).val.real,
weights=x.weight(1, spaces=spaces).val.real,
axis=x.domain_axes[spaces])
tgt_domain = list(x.domain)
tgt_domain[spaces] = self._target[0]
y = Field(tgt_domain, val=projected_x).weight(-1, spaces=spaces)
return y
def _adjoint_times(self,x,spaces):
def _adjoint_times(self, x, spaces):
if spaces is None:
spaces = 0
tgt_domain = list(x.domain)
......@@ -38,11 +38,11 @@ class PowerProjection(LinearOperator):
spec = self._spec_to_rescaler(spec, spaces, axes)
y.val.apply_scalar_function(lambda x: x * spec.real,
inplace=True)
inplace=True)
return y
def _spec_to_rescaler(self, spec, power_space_index,axes):
def _spec_to_rescaler(self, spec, power_space_index, axes):
# weight the random fields with the power spectrum
# therefore get the pindex from the power space
......@@ -67,6 +67,7 @@ class PowerProjection(LinearOperator):
# here, the power_spectrum is distributed into the new shape
local_rescaler = spec[local_blow_up]
return local_rescaler
@property
def domain(self):
return self._domain
......
......@@ -78,7 +78,7 @@ def create_power_operator(domain, power_spectrum, dtype=None,
f = f.real
# f **= 2
diag = Field(domain,f).weight()
diag = Field(domain, f).weight()
return DiagonalOperator(domain, diag)
......@@ -111,7 +111,7 @@ def generate_posterior_sample(mean, covariance):
power = sqrt(S.diagonal().weight(-1))
mock_signal = Field.from_random(random_type="normal", domain=S.domain,
dtype=power.dtype)
dtype=power.dtype)
mock_signal = power * mock_signal
noise = N.diagonal().weight(-1)
......
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