Commit ccc06713 authored by Philipp Frank's avatar Philipp Frank

added dynamic and light cone prior

parent 19fbc7e0
......@@ -78,6 +78,8 @@ from .plot import Plot
from .library.amplitude_model import AmplitudeModel
from .library.inverse_gamma_model import InverseGammaModel
from .library.los_response import LOSResponse
from .library.dynamic_operator import make_dynamic_operator
from .library.light_cone_operator import LightConeOperator
from .library.wiener_filter_curvature import WienerFilterCurvature
from .library.correlated_fields import CorrelatedField, MfCorrelatedField
......
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
# Copyright(C) 2013-2018 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.
from __future__ import absolute_import, division, print_function
import numpy as np
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 ..sugar import makeOp
from ..field import Field
from ..domains.unstructured_domain import UnstructuredDomain
from .light_cone_operator import LightConeOperator,field_from_function
def make_dynamic_operator(FFT,harmonic_padding,
keys=['f', 'c'],
causal=True,
cone=True,
minimum_phase=False,
sm_s0=15.,
sm_x0=[0.18, 0.18],
sigc=[3.],
quant=5.):
ops = {}
if harmonic_padding is None:
CentralPadd = ScalingOperator(1.,FFT.target)
else:
shp = ()
for i in range(len(FFT.target.shape)):
shp += (FFT.target.shape[i] + harmonic_padding[i],)
CentralPadd = FieldZeroPadder(FFT.target,shp,central=True)
ops['CentralPadd'] = CentralPadd
sdom = CentralPadd.target[0].get_default_codomain()
FFTB = FFTOperator(sdom)(Realizer(sdom))
m = FieldAdapter(sdom, keys[0])
dists = m.target[0].distances
def func(x):
res = 1.
for i in range(len(sm_x0)):
res = res + (x[i]/sm_x0[i]/dists[i])**2
return sm_s0/res
Sm = field_from_function(m.target, func)
Sm = makeOp(Sm)
m = Sm(m)
m = FFTB(m)
m = CentralPadd.adjoint(m)
ops[keys[0]+'_k'] = m
m = -m.log()
if not minimum_phase:
m = m.exp()
ops['Gncc'] = m
if causal:
CRHB = Realizer(FFT.target)
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 = kernel(m)
elif minimum_phase:
raise(ValueError,"minimum phase and not causal not possible!")
if cone:
if len(m.target.shape) < 2:
raise(ValueError,"Light cone requires dimensionality >= 2")
cdom = UnstructuredDomain(len(sigc))
c = FieldAdapter(cdom, keys[1])
Sigc = makeOp(Field(c.target, np.array(sigc)))
c = Sigc(c)
c = c.exp()
ops['c'] = c
c = LightConeOperator(c.target, m.target, quant)(c)
ops['a'] = c
m = c*m
ops['Gx'] = m
if causal:
m = FFT(m)
if minimum_phase:
m = m.exp()
ops['G'] = m
return m, ops
from __future__ import absolute_import, division, print_function
from ..field import Field
from ..domain_tuple import DomainTuple
from ..linearization import Linearization
from ..operators.linear_operator import LinearOperator
from ..operators.operator import Operator
import numpy as np
def make_coords(domain, absolute=False):
domain = DomainTuple.make(domain)
dim = len(domain.shape)
dist = domain[0].distances
shape = domain.shape
k_array = np.zeros((dim,) + shape)
for i in range(dim):
ks = np.minimum(shape[i] - np.arange(shape[i]), np.arange(
shape[i]))*dist[i]
if not absolute:
ks[int(shape[i]/2) + 1:] *= -1
fst_dims = (1,)*i
lst_dims = (1,)*(dim - i - 1)
k_array[i] += ks.reshape(fst_dims + (shape[i],) + lst_dims)
return k_array
def field_from_function(domain, func, absolute=False):
domain = DomainTuple.make(domain)
k_array = make_coords(domain, absolute=absolute)
return Field(domain, val=func(k_array))
class LightConeDerivative(LinearOperator):
def __init__(self, domain, target, derivatives):
super(LightConeDerivative, self).__init__()
self._domain = domain
self._target = target
self._derivatives = derivatives
self._capability = self.TIMES | self.ADJOINT_TIMES
def apply(self, x, mode):
self._check_input(x, mode)
x = x.val
res = np.zeros(self._tgt(mode).shape, dtype=self._derivatives.dtype)
for i in range(self.domain.shape[0]):
if mode == self.TIMES:
res += self._derivatives[i]*x[i]
else:
res[i] = np.sum(self._derivatives[i]*x)
return Field(self._tgt(mode), val=res)
def cone_arrays(c, domain, sigx,want_gradient):
x = make_coords(domain)
a = np.zeros(domain.shape, dtype=np.complex)
if want_gradient:
derivs = np.zeros((c.size,) + domain.shape, dtype=np.complex)
else:
derivs = None
a -= (x[0]/(sigx*domain[0].distances[0]))**2
for i in range(c.size):
res = (x[i + 1]/(sigx*domain[0].distances[i + 1]))**2
a += c[i]*res
if want_gradient:
derivs[i] = res
a = np.sqrt(a)
if want_gradient:
derivs *= -0.5
for i in range(c.size):
derivs[i][a == 0] = 0.
derivs[i][a != 0] /= a[a != 0]
a = a.real
if want_gradient:
derivs *= a
a = np.exp(-0.5*a**2)
if want_gradient:
derivs = a*derivs.real
return a, derivs
class LightConeOperator(Operator):
def __init__(self, domain, target, sigx):
self._domain = domain
self._target = target
self._sigx = sigx
def apply(self, x):
if not isinstance(x, Linearization):
a , _ = cone_arrays(
x.to_global_data(), self.target, self._sigx,False)
return Field(self.target, a)
a, derivs = cone_arrays(
x.val.to_global_data(), self.target, self._sigx,True)
jac = LightConeDerivative(x.jac.target, self.target, derivs)(x.jac)
return Linearization(
Field(self.target, a), 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