update iwp and scaling

nifty5/library/smooth_linear_amplitude.py

@@ -216,31 +216,11 @@ 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((2, 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
@@ -254,20 +234,18 @@ class _TwoLogIntegrations(LinearOperator):
             res = np.empty(self._target.shape)
             res[0] = 0
             res[1] = 0
-            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(x[1])
+            res[2:] = (res[2:]+res[1:-1])/2*self._logvol + x[0]
             res[2:] = np.cumsum(res[2:])
             return from_global_data(self._target, res)
         else:
             x = x.to_global_data_rw()
-            res = np.empty(self._domain.shape)
+            res = np.zeros(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)
+            res[0] += x[2:]
+            x[2:] *= self._logvol/2.
+            res[1:-1] += res[2:]
+            x[1] += np.cumsum(res[2:][::-1])[::-1]
             return from_global_data(self._domain, res)
 
 
@@ -307,7 +285,7 @@ def LogIntegratedWienerProcess(target, means, stddevs, wienersigmastddev,
     reststddevs = from_global_data(rest.domain, stddevs)
     rest = rest @ Adder(restmeans) @ makeOp(reststddevs)
 
-    expander = VdotOperator(full(target, 1.)).adjoint
+    
     m = means[1]
     L = np.log(target.k_lengths[-1]) - np.log(target.k_lengths[1])
 
@@ -315,20 +293,34 @@ def LogIntegratedWienerProcess(target, means, stddevs, wienersigmastddev,
     wienermean = np.sqrt(3/L)*np.abs(m)/norm.ppf(wienersigmaprob)
     wienermean = np.log(wienermean)
 
+    
+    twolog = _TwoLogIntegrations(target)
+    
+    dt = twolog._logvol
+    sc = np.zeros(twolog.domain.shape)
+    sc[0] = sc[1] = np.sqrt(dt)
+    sc = from_global_data(twolog.domain,sc)
+    expander = VdotOperator(sc).adjoint
+
     sigma = Adder(full(expander.domain, wienermean)) @ (
         wienersigmastddev*ducktape(expander.domain, None, keys[2]))
-    sigma = expander @ sigma.exp()
-    twolog = _TwoLogIntegrations(target)
+    sigmasq = expander @ sigma.exp()
     
     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
+    dist = np.zeros(twolog.domain)
+    dist[0] += 1.
+    dist = from_global_data(twolog.domain,dist)
+    scale = VdotOperator(dist).adjoint @ ep
+    
+    shift = np.ones(scale.target.shape)
+    shift[0] = dt**2/12.
+    shift = from_global_data(scale.target,shift)
+    scale = sigmasq * (Adder(shift)@scale)**0.5
 
-    smooth = sigma*(twolog @ (scale*ducktape(scale.target, None, keys[1])))
+    smooth = twolog@(scale*ducktape(scale.target,None,keys[1]))
     return rest.ducktape(keys[0]) + smooth