Commit 26fb00b8 authored by Philipp Frank's avatar Philipp Frank
Browse files

consistent iwp and wp

parent 54e9bb64
Pipeline #62930 passed with stages
in 6 minutes and 55 seconds
......@@ -216,12 +216,29 @@ def LinearSLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv,
# Go from loglog-space to linear-linear-space
return et @ loglog_ampl
class _Customscale(LinearOperator):
def __init__(self,target,dt):
self._target = DomainTuple.make(target)
self._domain = DomainTuple.scalar_domain()
self._sc = 12/dt**2
self._capability = self.TIMES | self.ADJOINT_TIMES
def apply(self,x,mode):
self._check_input(x,mode)
if mode == self.TIMES:
res = np.zeros(self._target.shape)
res[1] = self._sc * x.to_global_data()
res = from_global_data(self.target,res)
else:
res = (self._sc*x.to_global_data()[1]).sum()
res = Field.scalar(res)
return res
class _TwoLogIntegrations(LinearOperator):
def __init__(self, target):
self._target = makeDomain(target)
self._domain = makeDomain(
UnstructuredDomain(self.target.shape[0] - 2))
UnstructuredDomain((2,self.target.shape[0] - 2)))
self._capability = self.TIMES | self.ADJOINT_TIMES
if not isinstance(self._target[0], PowerSpace):
raise TypeError
......@@ -235,15 +252,21 @@ class _TwoLogIntegrations(LinearOperator):
res = np.empty(self._target.shape)
res[0] = 0
res[1] = 0
res[2:] = np.cumsum(x*self._logvol)
res[2:] = np.cumsum(res[2:]*self._logvol)
res[2:] = np.cumsum(x[0]*np.sqrt(self._logvol))
res[2:] *= self._logvol
res[2:] += (self._logvol*np.sqrt(self._logvol)*
(x[1]/np.sqrt(12) - x[0]/2.))
res[2:] = np.cumsum(res[2:])
return from_global_data(self._target, res)
else:
x = x.to_global_data()
res = np.empty(self._target.shape)
res[2:] = np.cumsum(x[2:][::-1])[::-1]*self._logvol
res[2:] = np.cumsum(res[2:][::-1])[::-1]*self._logvol
return from_global_data(self._domain, res[2:])
x = x.to_global_data_rw()
res = np.empty(self._domain.shape)
x[2:] = np.cumsum(x[2:][::-1])[::-1]
res[0] = self._logvol*np.sqrt(self._logvol)*(-x[2:]/2.)
res[1] = self._logvol*np.sqrt(self._logvol)*(x[2:]/np.sqrt(12))
x[2:] *= self._logvol
res[0] += np.cumsum(x[2:][::-1])[::-1]*np.sqrt(self._logvol)
return from_global_data(self._domain, res)
class _Rest(LinearOperator):
......@@ -269,7 +292,7 @@ class _Rest(LinearOperator):
def LogIntegratedWienerProcess(target, means, stddevs, wienersigmastddev,
wienersigmaprob, keys):
wienersigmaprob, mymean, mystd, keys):
# means and stddevs: zm, slope, yintercept
# keys: rest smooth wienersigma
if not (len(means) == 3 and len(stddevs) == 3 and len(keys) == 3):
......@@ -293,7 +316,17 @@ def LogIntegratedWienerProcess(target, means, stddevs, wienersigmastddev,
sigma = Adder(full(expander.domain, wienermean)) @ (
wienersigmastddev*ducktape(expander.domain, None, keys[2]))
sigma = expander @ sigma.exp()
smooth = _TwoLogIntegrations(target).ducktape(keys[1])*sigma
twolog = _TwoLogIntegrations(target)
ep = ducktape(expander.domain, None, keys[3])
ep = Adder(Field.scalar(mymean)) @ (mystd*ep)
ep =ep.exp()
scale = _Customscale(twolog.domain,twolog._logvol) @ ep
scale = Adder(full(scale.target,1.))@scale
scale = scale**0.5
smooth = sigma*(twolog@(scale*ducktape(scale.target,None,keys[1])))
return rest.ducktape(keys[0]) + smooth
......
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