Commit a12138dd authored by Martin Reinecke's avatar Martin Reinecke

demo 3 now working

parent 4356aae3
......@@ -81,7 +81,7 @@ from .minimization.energy_sum import EnergySum
from .sugar import *
from .plotting.plot import plot, plot_finish
from .library.amplitude_model import make_amplitude_model
from .library.amplitude_model import make_amplitude_model, AmplitudeModel
from .library.gaussian_energy import GaussianEnergy
from .library.los_response import LOSResponse
from .library.point_sources import PointSources
......
......@@ -618,6 +618,9 @@ class Field(object):
def unite(self, other):
return self + other
def positive_tanh(self):
return 0.5*(1.+self.tanh())
for op in ["__add__", "__radd__",
"__sub__", "__rsub__",
"__mul__", "__rmul__",
......
......@@ -25,7 +25,9 @@ from ..domains.power_space import PowerSpace
from ..domains.unstructured_domain import UnstructuredDomain
from ..field import Field
from ..multi.multi_field import MultiField
from ..multi.multi_domain import MultiDomain
from ..sugar import makeOp, sqrt
from ..operator import Operator
def _ceps_kernel(dof_space, k, a, k0):
......@@ -144,3 +146,82 @@ def create_cepstrum_amplitude_field(domain, cepstrum):
cepstrum_field[sl2] = np.sum(cepstrum_field[sl], axis=i)
return Field.from_global_data(domain, cepstrum_field)
class AmplitudeModel(Operator):
'''
Computes a smooth power spectrum.
Output lives in PowerSpace.
Parameters
----------
Npixdof : #pix in dof_space
ceps_a, ceps_k0 : Smoothness parameters in ceps_kernel
eg. ceps_kernel(k) = (a/(1+(k/k0)**2))**2
a = ceps_a, k0 = ceps_k0
sm, sv : slope_mean = expected exponent of power law (e.g. -4),
slope_variance (default=1)
im, iv : y-intercept_mean, y-intercept_variance of power_slope
'''
def __init__(self, s_space, Npixdof, ceps_a, ceps_k, sm, sv, im, iv,
keys=['tau', 'phi']):
from ..operators.exp_transform import ExpTransform
from ..operators.qht_operator import QHTOperator
from ..operators.slope_operator import SlopeOperator
from ..operators.symmetrizing_operator import SymmetrizingOperator
from ..models.variable import Variable
from ..models.constant import Constant
from ..models.local_nonlinearity import PointwiseExponential
h_space = s_space.get_default_codomain()
p_space = PowerSpace(h_space)
self._exp_transform = ExpTransform(p_space, Npixdof)
logk_space = self._exp_transform.domain[0]
qht = QHTOperator(target=logk_space)
dof_space = qht.domain[0]
param_space = UnstructuredDomain(2)
sym = SymmetrizingOperator(logk_space)
phi_mean = np.array([sm, im])
phi_sig = np.array([sv, iv])
self._slope = SlopeOperator(param_space, logk_space, phi_sig)
self._norm_phi_mean = Field.from_global_data(param_space, phi_mean/phi_sig)
self._domain = MultiDomain.make({keys[0]: dof_space, keys[1]: param_space})
# fields = {keys[0]: Field.from_random('normal', dof_space),
# keys[1]: Field.from_random('normal', param_space)}
# position = MultiField.from_dict(fields)
# dof_space = position[keys[0]].domain[0]
kern = lambda k: _ceps_kernel(dof_space, k, ceps_a, ceps_k)
cepstrum = create_cepstrum_amplitude_field(dof_space, kern)
ceps = makeOp(sqrt(cepstrum))
self._smooth_op = sym * qht * ceps
self._keys = tuple(keys)
# smooth_spec = smooth_op(Variable(position)[keys[0]])
# phi = Variable(position)[keys[1]] + Constant(position, norm_phi_mean)
# linear_spec = slope(phi)
# loglog_spec = smooth_spec + linear_spec
# xlog_ampl = PointwiseExponential(0.5*loglog_spec)
def __call__(self, x):
smooth_spec = self._smooth_op(x[self._keys[0]])
phi = x[self._keys[1]] + self._norm_phi_mean
linear_spec = self._slope(phi)
loglog_spec = smooth_spec + linear_spec
return self._exp_transform((0.5*loglog_spec).exp())
@property
def domain(self):
return self._domain
@property
def target(self):
return self._exp_transform.target
......@@ -41,6 +41,11 @@ class Linearization(object):
"""Only available if target is a scalar"""
return self._metric
def __getitem__(self, name):
from .operators.field_adapter import FieldAdapter
dom = self._val[name].domain
return Linearization(self._val[name], FieldAdapter(dom, name))
def __neg__(self):
return Linearization(-self._val, self._jac*(-1),
None if self._metric is None else self._metric*(-1))
......@@ -68,11 +73,12 @@ class Linearization(object):
def __mul__(self, other):
from .sugar import makeOp
from .operators.relaxed_sum_operator import RelaxedSumOperator
if isinstance(other, Linearization):
d1 = makeOp(self._val)
d2 = makeOp(other._val)
return Linearization(self._val*other._val,
d2*self._jac + d1*other._jac)
RelaxedSumOperator((d2*self._jac, d1*other._jac)))
if isinstance(other, (int, float, complex)):
# if other == 0:
# return ...
......@@ -104,6 +110,15 @@ class Linearization(object):
tmp = self._val.log()
return Linearization(tmp, makeOp(1./self._val)*self._jac)
def tanh(self):
tmp = self._val.tanh()
return Linearization(tmp, makeOp(1.-tmp**2)*self._jac)
def positive_tanh(self):
tmp = self._val.tanh()
tmp2 = 0.5*(1.+tmp)
return Linearization(tmp2, makeOp(0.5*(1.-tmp**2))*self._jac)
def add_metric(self, metric):
return Linearization(self._val, self._jac, metric)
......
......@@ -7,16 +7,9 @@ from ..multi.multi_field import MultiField
class FieldAdapter(LinearOperator):
def __init__(self, op, name_dom, name_tgt):
if name_dom is None:
self._domain = op.domain
else:
self._domain = MultiDomain.make({name_dom: op.domain})
if name_tgt is None:
self._target = op.target
else:
self._target = MultiDomain.make({name_tgt: op.target})
self._op = op
def __init__(self, dom, name_dom):
self._domain = MultiDomain.make({name_dom: dom})
self._target = dom
@property
def domain(self):
......@@ -28,16 +21,11 @@ class FieldAdapter(LinearOperator):
@property
def capability(self):
return self._op.capability
return self._all_ops
def apply(self, x, mode):
self._check_input(x, mode)
dom = self._dom(mode)
if x.domain is not dom:
x = dom.values()[0]
res = self._op.apply(x. mode)
tgt = self._tgt(mode)
if res.domain is not tgt:
res = MultiField(tgt, (res,))
return res
if mode == self.TIMES:
return x.values()[0]
return MultiField(self._domain, (x,))
......@@ -63,3 +63,13 @@ class RelaxedSumOperator(LinearOperator):
tmp = op.apply(x.extract(op._dom(mode)), mode)
res = tmp if res is None else res.unite(tmp)
return res
def draw_sample(self, from_inverse=False, dtype=np.float64):
if from_inverse:
raise NotImplementedError(
"cannot draw from inverse of this operator")
res = None
for op in self._ops:
tmp = op.draw_sample(from_inverse, dtype)
res = tmp if res is None else res.unite(tmp)
return 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