Commit 0aa9e11a authored by Philipp Frank's avatar Philipp Frank
Browse files


parent 00536ff8
......@@ -24,6 +24,7 @@ from ..operators.scaling_operator import ScalingOperator
from ..operators.harmonic_operators import FFTOperator
from ..operators.field_zero_padder import FieldZeroPadder
from ..operators.simple_linear_operators import Realizer,FieldAdapter
from ..operators.diagonal_operator import DiagonalOperator
from ..sugar import makeOp
from ..field import Field
from import UnstructuredDomain
......@@ -38,41 +39,44 @@ def make_dynamic_operator(FFT,harmonic_padding,sm_s0,sm_x0,
Constructs an operator for a dynamic field prior
Constructs an operator for a dynamic field prior.
FFT : FFTOperator
A Fourier Transformation Operator of the space under consideration
harmonic_padding : None, list of float
Amount of central padding in harmonic space in pixels. If None the field is not padded at all.
sm_s0 : float
Cutoff for dynamic smoothness prior
sm_x0 : float, List of float
Scaling of dynamic smoothness along each axis
keys : List of String
keys of input fields of operator.
causal : boolean
Whether or not the reconstructed dynamics should be causal in time
cone : boolean
Whether or not the reconstructed dynamics should be within a light cone
minimum_phase: boolean
Whether or not the reconstructed dynamics should be minimum phase
sigc : float, List of float
variance of light cone parameters.
If cone is False this is ignored
quant : float
Quantization of the light cone in pixels.
If cone is False this is ignored
The Operator encoding the dynamic Greens function in harmonic space when evaluated.
Dictionary of Operator
A collection of sub-chains of Operators which can be used for plotting and evaluation.
Currently only supports RGSpaces.
Note that the first axis of the space is interpreted as the time axis.
ops = {}
if harmonic_padding is None:
......@@ -91,18 +95,16 @@ def make_dynamic_operator(FFT,harmonic_padding,sm_s0,sm_x0,
dists =[0].distances
if isinstance(sm_x0,float):
sm_x0 = list((sm_x0,)*len(dists))
def func(x):
elif len(sm_x0) != len(dists):
raise(ValueError,"Shape mismatch!")
def smoothness_prior_func(x):
res = 1.
for i in range(len(dists)):
res = res + (x[i]/sm_x0[i]/dists[i])**2
return sm_s0/res
Sm = field_from_function(, func)
Sm = makeOp(Sm)
m = Sm(m)
m = FFTB(m)
m = CentralPadd.adjoint(m)
Sm = makeOp(field_from_function(, smoothness_prior_func))
m = CentralPadd.adjoint(FFTB(Sm(m)))
ops[keys[0]+'_k'] = m
m = -m.log()
......@@ -110,14 +112,8 @@ def make_dynamic_operator(FFT,harmonic_padding,sm_s0,sm_x0,
m = m.exp()
ops['Gncc'] = m
if causal:
CRHB = Realizer(
m = FFT.inverse(CRHB.adjoint(m))
def func(x):
res = 1. + np.sign(x[0])
return res
kernel = field_from_function(FFT.domain, func)
kernel = makeOp(kernel)
m = FFT.inverse(Realizer(
kernel = makeOp(field_from_function(FFT.domain, (lambda x: 1.+np.sign(x[0]))))
m = kernel(m)
elif minimum_phase:
raise(ValueError,"minimum phase and not causal not possible!")
......@@ -125,13 +121,17 @@ def make_dynamic_operator(FFT,harmonic_padding,sm_s0,sm_x0,
if cone and len( > 1:
if isinstance(sigc,float):
sigc = list((sigc,)*(len(
cdom = UnstructuredDomain(len(sigc))
c = FieldAdapter(cdom, keys[1])
Sigc = makeOp(Field.from_global_data(, np.array(sigc)))
c = Sigc(c)
c = c.exp()
ops['c'] = c
c = LightConeOperator(,, quant)(c)
elif len(sigc) != len(
raise(ValueError,"Shape mismatch!")
c = FieldAdapter(UnstructuredDomain(len(sigc)), keys[1])
c = makeOp(Field.from_global_data(, np.array(sigc)))(c)
lightspeed = ScalingOperator(-0.5,
scaling =[0].distances[1:]/[0].distances[0]
scaling = DiagonalOperator(Field.from_global_data(,scaling))
ops['lightspeed'] = scaling(lightspeed)
c = LightConeOperator(,, quant)(c.exp())
ops['a'] = c
m = c*m
......@@ -140,5 +140,4 @@ def make_dynamic_operator(FFT,harmonic_padding,sm_s0,sm_x0,
m = FFT(m)
if minimum_phase:
m = m.exp()
ops['G'] = m
return m, ops
......@@ -83,14 +83,11 @@ class LightConeOperator(Operator):
self._sigx = sigx
def apply(self, x):
if not isinstance(x, Linearization):
a , _ = cone_arrays(
x.to_global_data(),, self._sigx,False)
return Field.from_global_data(, a)
a, derivs = cone_arrays(
x.val.to_global_data(),, self._sigx,True)
islin = isinstance(x, Linearization)
val = x.val.to_global_data() if islin else x.to_global_data()
a, derivs = cone_arrays(val,, self._sigx, islin)
res = Field.from_global_data(, a)
if not islin:
return res
jac = LightConeDerivative(,, derivs)(x.jac)
return Linearization(
Field.from_global_data(, a), jac, want_metric=x.want_metric)
return Linearization(res, jac, want_metric=x.want_metric)
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