Compute fixed point residual in Fourier domain.
434220
def rsdl(self):
        """"""

        diff = self.Xf - self.Yfprv
        return sl.rfl2norm2(diff, self.X.shape, axis=self.cri.axisN)
155622
def _compute_fixed(self):
        ''''''
        try:
            lon, lat = np.meshgrid(self.lon, self.lat)
        except:
            lat = self.lat
        phi = np.deg2rad(lat)
        try:
            albedo = self.a0 + self.a2 * P2(np.sin(phi))
        except:
            albedo = np.zeros_like(phi)
        # make sure that the diagnostic has the correct field dimensions.
        #dom = self.domains['default']
        #  this is a more robust way to get the single value from dictionary:
        dom = next(iter(self.domains.values()))
        self.albedo = Field(albedo, domain=dom)
1137340
def _call(self, x):
        """"""
        return scipy.signal.fftconvolve(self.kernel, x, mode='same')
925481
def inverse(self):
        """"""
        sign = '+' if self.sign == '-' else '-'
        return FourierTransform(
            domain=self.range, range=self.domain, impl=self.impl,
            axes=self.axes, halfcomplex=self.halfcomplex, shift=self.shifts,
            sign=sign, tmp_r=self._tmp_r, tmp_f=self._tmp_f)
56400
def _compute_dlt(self):
        """"""
        res = super()._compute_dlt()
        for rec in self:
            ltaf_to_apply = self.env['ddmrp.adjustment'].search(
                rec._ltaf_to_apply_domain())
            if ltaf_to_apply:
                ltaf = 1
                values = ltaf_to_apply.mapped('value')
                for val in values:
                    ltaf *= val
                prev = rec.dlt
                rec.dlt *= ltaf
                _logger.debug(
                    "LTAF=%s applied to %s. DLT: %s -> %s" %
                    (ltaf, rec.name, prev, rec.dlt))
        return res
923658
def _calc_fourier_spectrum(self, fa_length=None):
        """"""

        if fa_length is None:
            # Use the next power of 2 for the length
            n = 1
            while n < self.accels.size:
                n <<= 1
        else:
            n = fa_length

        self._fourier_amps = np.fft.rfft(self._accels, n)

        freq_step = 1. / (2 * self._time_step * (n / 2))
        self._freqs = freq_step * np.arange(1 + n / 2)
1091446
def tilt_residual(params, data, mask):
    """"""
    bg = tilt_model(params, shape=data.shape)
    res = (data - bg)[mask]
    return res.flatten()
589656
def fft(ts):
    
    t_step = ts.index[1] - ts.index[0]
    oc = np.abs(np.fft.fftshift(np.fft.fft(ts.values))) / len(ts.values)
    t = np.fft.fftshift(np.fft.fftfreq(len(oc), d=t_step))
    return Trace(oc, t)
882826
def _call(self, x):
        """"""
        if isinstance(self.domain, ProductSpace):
            norm = PointwiseNorm(self.domain, 2)(x)
        else:
            norm = x.ufuncs.absolute()

        if self.gamma > 0:
            tmp = norm.ufuncs.square()
            tmp *= 1 / (2 * self.gamma)

            index = norm.ufuncs.greater_equal(self.gamma)
            tmp[index] = norm[index] - self.gamma / 2
        else:
            tmp = norm

        return tmp.inner(tmp.space.one())
755064
def spa_tmplt_precondition(length, delta_f, kmin=0):
    
    global _prec
    if _prec is None or _prec.delta_f != delta_f or len(_prec) < length:
        v = numpy.arange(0, (kmin+length*2), 1.0) * delta_f
        v = numpy.power(v[1:len(v)], -7.0/6.0)
        _prec = FrequencySeries(v, delta_f=delta_f, dtype=float32)
    return _prec[kmin:kmin + length]
108202
def obfn_dfd(self):
        
        """

        Ef = sl.inner(self.Zf, self.obfn_fvarf(), axis=self.cri.axisM) \
          - self.Sf
        return (np.linalg.norm(self.W * sl.irfftn(Ef, self.cri.Nv,
                                                  self.cri.axisN))**2) / 2.0
211524
def fft_image(shape, sd=None, decay_power=1):
    """"""

    sd = sd or 0.01
    batch, h, w, ch = shape
    freqs = rfft2d_freqs(h, w)
    init_val_size = (2, ch) + freqs.shape

    images = []
    for _ in range(batch):
        # Create a random variable holding the actual 2D fourier coefficients
        init_val = np.random.normal(size=init_val_size, scale=sd).astype(np.float32)
        spectrum_real_imag_t = tf.Variable(init_val)
        spectrum_t = tf.complex(spectrum_real_imag_t[0], spectrum_real_imag_t[1])

        # Scale the spectrum. First normalize energy, then scale by the square-root
        # of the number of pixels to get a unitary transformation.
        # This allows to use similar leanring rates to pixel-wise optimisation.
        scale = 1.0 / np.maximum(freqs, 1.0 / max(w, h)) ** decay_power
        scale *= np.sqrt(w * h)
        scaled_spectrum_t = scale * spectrum_t

        # convert complex scaled spectrum to shape (h, w, ch) image tensor
        # needs to transpose because irfft2d returns channels first
        image_t = tf.transpose(tf.spectral.irfft2d(scaled_spectrum_t), (1, 2, 0))

        # in case of odd spatial input dimensions we need to crop
        image_t = image_t[:h, :w, :ch]

        images.append(image_t)

    batched_image_t = tf.stack(images) / 4.0  # TODO: is that a magic constant?
    return batched_image_t
594429
def poly2o_residual(params, data, mask):
    """"""
    bg = poly2o_model(params, shape=data.shape)
    res = (data - bg)[mask]
    return res.flatten()
1295670
def _reshape_like(F, x, y):
    """"""
    return x.reshape(y.shape) if F is ndarray else F.reshape_like(x, y)
759036
def damping_kraus_map(p=0.10):
    
    damping_op = np.sqrt(p) * np.array([[0, 1],
                                        [0, 0]])

    residual_kraus = np.diag([1, np.sqrt(1 - p)])
    return [residual_kraus, damping_op]
61098
def rsdl_rn(self, AX, Y):
        """"""

        # The primal residual normalisation term is
        # max( ||A x^(k)||_2, ||B y^(k)||_2 ) and B = -(I I I ...)^T.
        # The scaling by sqrt(Nb) of the l2 norm of Y accounts for the
        # block replication introduced by multiplication by B
        return max((np.linalg.norm(AX), np.sqrt(self.Nb) * np.linalg.norm(Y)))
537614
def normalize_residuals(self, points):
        """"""
        residuals = self.evaluate_residual(points)
        solutions = self.evaluate_solution(points)
        return [resid / soln for resid, soln in zip(residuals, solutions)]