The following are code examples for showing how to use . They are extracted from open source Python projects. You can vote up the examples you like or vote down the exmaples you don’t like. You can also save this page to your account.

Example 1

def correlate(self, imgfft):
        #Very much related to the convolution theorem, the cross-correlation
        #theorem states that the Fourier transform of the cross-correlation of
        #two functions is equal to the product of the individual Fourier
        #transforms, where one of them has been complex conjugated:  
        
        
        if self.imgfft is not 0 or imgfft.imgfft is not 0:
            imgcj = np.conjugate(self.imgfft)
            imgft = imgfft.imgfft
            
            prod = deepcopy(imgcj)
            for x in range(imgcj.shape[0]):
                for y in range(imgcj.shape[0]):
                    prod[x][y] = imgcj[x][y] * imgft[x][y]
            
            cc = Corr( np.real(fft.ifft2(fft.fftshift(prod)))) # real image of the correlation
            
            # adjust to center
            cc.data = np.roll(cc.data, int(cc.data.shape[0] / 2), axis = 0)
            cc.data = np.roll(cc.data, int(cc.data.shape[1] / 2), axis = 1)
        else:
            raise FFTnotInit()
        return cc 

Example 2

def genSpectra(time,dipole,signal):

    fw, frequency = pade(time,dipole)
    fw_sig, frequency = pade(time,signal,alternate=True)

    fw_re = np.real(fw)
    fw_im = np.imag(fw)
    fw_abs = fw_re**2 + fw_im**2

    #spectra = (fw_re*17.32)/(np.pi*field*damp_const)
    #spectra = (fw_re*17.32*514.220652)/(np.pi*field*damp_const)
    #numerator = np.imag((fw*np.conjugate(fw_sig)))
    numerator = np.imag(fw_abs*np.conjugate(fw_sig))
    #numerator = np.abs((fw*np.conjugate(fw_sig)))
    #numerator = np.abs(fw)
    denominator = np.real(np.conjugate(fw_sig)*fw_sig)
    #denominator = 1.0 
    spectra = ((4.0*27.21138602*2*frequency*np.pi*(numerator))/(3.0*137.036*denominator))
    spectra *= 1.0/100.0
    #plt.plot(frequency*27.2114,fourier)
    #plt.show()
    return frequency, spectra 

Example 3

def test_generic_methods(self):
        # Tests some MaskedArray methods.
        a = array([1, 3, 2])
        assert_equal(a.any(), a._data.any())
        assert_equal(a.all(), a._data.all())
        assert_equal(a.argmax(), a._data.argmax())
        assert_equal(a.argmin(), a._data.argmin())
        assert_equal(a.choose(0, 1, 2, 3, 4), a._data.choose(0, 1, 2, 3, 4))
        assert_equal(a.compress([1, 0, 1]), a._data.compress([1, 0, 1]))
        assert_equal(a.conj(), a._data.conj())
        assert_equal(a.conjugate(), a._data.conjugate())

        m = array([[1, 2], [3, 4]])
        assert_equal(m.diagonal(), m._data.diagonal())
        assert_equal(a.sum(), a._data.sum())
        assert_equal(a.take([1, 2]), a._data.take([1, 2]))
        assert_equal(m.transpose(), m._data.transpose()) 

Example 4

def test_testArrayMethods(self):
        a = array([1, 3, 2])
        self.assertTrue(eq(a.any(), a._data.any()))
        self.assertTrue(eq(a.all(), a._data.all()))
        self.assertTrue(eq(a.argmax(), a._data.argmax()))
        self.assertTrue(eq(a.argmin(), a._data.argmin()))
        self.assertTrue(eq(a.choose(0, 1, 2, 3, 4),
                           a._data.choose(0, 1, 2, 3, 4)))
        self.assertTrue(eq(a.compress([1, 0, 1]), a._data.compress([1, 0, 1])))
        self.assertTrue(eq(a.conj(), a._data.conj()))
        self.assertTrue(eq(a.conjugate(), a._data.conjugate()))
        m = array([[1, 2], [3, 4]])
        self.assertTrue(eq(m.diagonal(), m._data.diagonal()))
        self.assertTrue(eq(a.sum(), a._data.sum()))
        self.assertTrue(eq(a.take([1, 2]), a._data.take([1, 2])))
        self.assertTrue(eq(m.transpose(), m._data.transpose())) 

Example 5

def test_spherical_harmonic(scheme):
    '''Assert the norm of the spherical harmonic

    Y_1^1(phi, theta) = -1/2 sqrt(3/2/pi) * exp(i*phi) * sin(theta)

    is indeed 1, i.e.,

    int_0^2pi int_0^pi
        Y_1^1(phi, theta) * conj(Y_1^1(phi, theta)) * sin(theta)
        dphi dtheta = 1.
    '''
    def spherical_harmonic_11(azimuthal, polar):
        # y00 = 1.0 / numpy.sqrt(4*numpy.pi)
        y11 = -0.5 * numpy.sqrt(3.0/2.0/numpy.pi) \
            * numpy.exp(1j*azimuthal) * numpy.sin(polar)
        return y11 * numpy.conjugate(y11)

    val = quadpy.sphere.integrate_spherical(
        spherical_harmonic_11,
        rule=scheme
        )

    assert abs(val - 1.0) < 1.0e-15
    return 

Example 6

def __init__(self, n):
        self.degree = n

        self.points = -numpy.cos(numpy.pi * (numpy.arange(n) + 0.5) / n)

        # n -= 1
        N = numpy.arange(1, n, 2)
        length = len(N)
        m = n - length
        K = numpy.arange(m)

        v0 = numpy.concatenate([
            2 * numpy.exp(1j*numpy.pi*K/n) / (1 - 4*K**2),
            numpy.zeros(length+1)
            ])
        v1 = v0[:-1] + numpy.conjugate(v0[:0:-1])

        w = numpy.fft.ifft(v1)
        assert max(w.imag) < 1.0e-15
        self.weights = w.real

        return 

Example 7

def CoreShellScatteringFunction(mCore,mShell,wavelength,dCore,dShell,minAngle=0, maxAngle=180, angularResolution=0.5, normed=False):
#  http://pymiescatt.readthedocs.io/en/latest/forwardCS.html#CoreShellScatteringFunction
  xCore = np.pi*dCore/wavelength
  xShell = np.pi*dShell/wavelength
  theta = np.linspace(minAngle,maxAngle,int((maxAngle-minAngle)/angularResolution))*np.pi/180
  thetaSteps = len(theta)
  SL = np.zeros(thetaSteps)
  SR = np.zeros(thetaSteps)
  SU = np.zeros(thetaSteps)
  for j in range(thetaSteps):
    u = np.cos(theta[j])
    S1,S2 = CoreShellS1S2(mCore,mShell,xCore,xShell,u)
    SL[j] = (np.sum((np.conjugate(S1)*S1))).real
    SR[j] = (np.sum((np.conjugate(S2)*S2))).real
    SU[j] = (SR[j]+SL[j])/2
  if normed:
    SL /= np.max(SL)
    SR /= np.max(SR)
    SU /= np.max(SU)
  return theta,SL,SR,SU 

Example 8

def CoreShellScatteringFunction(mCore,mShell,wavelength,dCore,dShell,minAngle=0, maxAngle=180, angularResolution=0.5, normed=False):
#  http://pymiescatt.readthedocs.io/en/latest/forwardCS.html#CoreShellScatteringFunction
  xCore = np.pi*dCore/wavelength
  xShell = np.pi*dShell/wavelength
  theta = np.linspace(minAngle,maxAngle,int((maxAngle-minAngle)/angularResolution))*np.pi/180
  thetaSteps = len(theta)
  SL = np.zeros(thetaSteps)
  SR = np.zeros(thetaSteps)
  SU = np.zeros(thetaSteps)
  for j in range(thetaSteps):
    u = np.cos(theta[j])
    S1,S2 = CoreShellS1S2(mCore,mShell,xCore,xShell,u)
    SL[j] = (np.sum((np.conjugate(S1)*S1))).real
    SR[j] = (np.sum((np.conjugate(S2)*S2))).real
    SU[j] = (SR[j]+SL[j])/2
  if normed:
    SL /= np.max(SL)
    SR /= np.max(SR)
    SU /= np.max(SU)
  return theta,SL,SR,SU 

Example 9

def gain_substitution_scalar(gain, x, xwt):
    nants, nchan, nrec, _ = gain.shape
    newgain = numpy.ones_like(gain, dtype='complex')
    gwt = numpy.zeros_like(gain, dtype='float')
    
    # We are going to work with Jones 2x2 matrix formalism so everything has to be
    # converted to that format
    x = x.reshape(nants, nants, nchan, nrec, nrec)
    xwt = xwt.reshape(nants, nants, nchan, nrec, nrec)
    
    for ant1 in range(nants):
        for chan in range(nchan):
            # Loop over e.g. 'RR', 'LL, or 'xx', 'YY' ignoring cross terms
            top = numpy.sum(x[:, ant1, chan, 0, 0] *
                            gain[:, chan, 0, 0] * xwt[:, ant1, chan, 0, 0], axis=0)
            bot = numpy.sum((gain[:, chan, 0, 0] * numpy.conjugate(gain[:, chan, 0, 0]) *
                             xwt[:, ant1, chan, 0, 0]).real, axis=0)
            
            if bot > 0.0:
                newgain[ant1, chan, 0, 0] = top / bot
                gwt[ant1, chan, 0, 0] = bot
            else:
                newgain[ant1, chan, 0, 0] = 0.0
                gwt[ant1, chan, 0, 0] = 0.0
    return newgain, gwt 

Example 10

def convolve_scalestack(scalestack, img):
    """Convolve img by the specified scalestack, returning the resulting stack

    :param scalestack: stack containing the scales
    :param img: Image to be convolved
    :return: stack
    """

    convolved = numpy.zeros(scalestack.shape)
    ximg = numpy.fft.fftshift(numpy.fft.fft2(numpy.fft.fftshift(img)))

    nscales = scalestack.shape[0]
    for iscale in range(nscales):
        xscale = numpy.fft.fftshift(numpy.fft.fft2(numpy.fft.fftshift(scalestack[iscale, :, :])))
        xmult = ximg * numpy.conjugate(xscale)
        convolved[iscale, :, :] = numpy.real(numpy.fft.ifftshift(numpy.fft.ifft2(numpy.fft.ifftshift(xmult))))
    return convolved 

Example 11

def test_generic_methods(self):
        # Tests some MaskedArray methods.
        a = array([1, 3, 2])
        assert_equal(a.any(), a._data.any())
        assert_equal(a.all(), a._data.all())
        assert_equal(a.argmax(), a._data.argmax())
        assert_equal(a.argmin(), a._data.argmin())
        assert_equal(a.choose(0, 1, 2, 3, 4), a._data.choose(0, 1, 2, 3, 4))
        assert_equal(a.compress([1, 0, 1]), a._data.compress([1, 0, 1]))
        assert_equal(a.conj(), a._data.conj())
        assert_equal(a.conjugate(), a._data.conjugate())

        m = array([[1, 2], [3, 4]])
        assert_equal(m.diagonal(), m._data.diagonal())
        assert_equal(a.sum(), a._data.sum())
        assert_equal(a.take([1, 2]), a._data.take([1, 2]))
        assert_equal(m.transpose(), m._data.transpose()) 

Example 12

def test_testArrayMethods(self):
        a = array([1, 3, 2])
        self.assertTrue(eq(a.any(), a._data.any()))
        self.assertTrue(eq(a.all(), a._data.all()))
        self.assertTrue(eq(a.argmax(), a._data.argmax()))
        self.assertTrue(eq(a.argmin(), a._data.argmin()))
        self.assertTrue(eq(a.choose(0, 1, 2, 3, 4),
                           a._data.choose(0, 1, 2, 3, 4)))
        self.assertTrue(eq(a.compress([1, 0, 1]), a._data.compress([1, 0, 1])))
        self.assertTrue(eq(a.conj(), a._data.conj()))
        self.assertTrue(eq(a.conjugate(), a._data.conjugate()))
        m = array([[1, 2], [3, 4]])
        self.assertTrue(eq(m.diagonal(), m._data.diagonal()))
        self.assertTrue(eq(a.sum(), a._data.sum()))
        self.assertTrue(eq(a.take([1, 2]), a._data.take([1, 2])))
        self.assertTrue(eq(m.transpose(), m._data.transpose())) 

Example 13

def test_generic_methods(self):
        # Tests some MaskedArray methods.
        a = array([1, 3, 2])
        assert_equal(a.any(), a._data.any())
        assert_equal(a.all(), a._data.all())
        assert_equal(a.argmax(), a._data.argmax())
        assert_equal(a.argmin(), a._data.argmin())
        assert_equal(a.choose(0, 1, 2, 3, 4), a._data.choose(0, 1, 2, 3, 4))
        assert_equal(a.compress([1, 0, 1]), a._data.compress([1, 0, 1]))
        assert_equal(a.conj(), a._data.conj())
        assert_equal(a.conjugate(), a._data.conjugate())

        m = array([[1, 2], [3, 4]])
        assert_equal(m.diagonal(), m._data.diagonal())
        assert_equal(a.sum(), a._data.sum())
        assert_equal(a.take([1, 2]), a._data.take([1, 2]))
        assert_equal(m.transpose(), m._data.transpose()) 

Example 14

def test_testArrayMethods(self):
        a = array([1, 3, 2])
        self.assertTrue(eq(a.any(), a._data.any()))
        self.assertTrue(eq(a.all(), a._data.all()))
        self.assertTrue(eq(a.argmax(), a._data.argmax()))
        self.assertTrue(eq(a.argmin(), a._data.argmin()))
        self.assertTrue(eq(a.choose(0, 1, 2, 3, 4),
                           a._data.choose(0, 1, 2, 3, 4)))
        self.assertTrue(eq(a.compress([1, 0, 1]), a._data.compress([1, 0, 1])))
        self.assertTrue(eq(a.conj(), a._data.conj()))
        self.assertTrue(eq(a.conjugate(), a._data.conjugate()))
        m = array([[1, 2], [3, 4]])
        self.assertTrue(eq(m.diagonal(), m._data.diagonal()))
        self.assertTrue(eq(a.sum(), a._data.sum()))
        self.assertTrue(eq(a.take([1, 2]), a._data.take([1, 2])))
        self.assertTrue(eq(m.transpose(), m._data.transpose())) 

Example 15

def cumsum_snr(freqs, detector, data, hpf, hxf, theta, phi, psi, zeroFreq=False, normalizeTemplate=True ):
    """
    returns the cumulative sum of the snr as a function of frequency
    does NOT maximize over the phase at coalescence
    """
    template = detector.project(freqs, hpf, hxf, theta, phi, psi, zeroFreq=zeroFreq)
    PSD = detector.PSD(freqs)
    deltaF = freqs[1]-freqs[0]

    ans = 2*np.cumsum(deltaF*np.conjugate(data)*template/PSD).real 
    if normalizeTemplate:
        ans /= np.sum(deltaF*np.conjugate(template)*template/PSD).real**0.5

    return ans

#------------------------ 

Example 16

def __newlambdagammaC(self, theta, l):
        """ Apply Eqns. 25-27 in Vidal 2003 to update lambda^C and gamma^C
            (lambda and gamma of this qbit).
        """
        gamma_ket = self.coefs[l+1].lam
        gamma_bra = np.conjugate(gamma_ket)
        Gamma_star = np.conjugate(self.coefs[l+1].gamma)
        inputs = [Gamma_star, theta, gamma_bra, gamma_ket]
        Gamma_star_idx = [1, -3, -2]
        theta_idx = [-1, 1, -4, -5]
        gamma_bra_idx = [-6]
        gamma_ket_idx = [-7]
        idx = [Gamma_star_idx, theta_idx, gamma_bra_idx, gamma_ket_idx]
        contract_me = scon(inputs, idx)
        svd_me = np.einsum('agibggg', contract_me)
        evals, evecs = la.eigh(svd_me)
        return evals, evecs 

Example 17

def paulidouble(i, j, tensor=True):
    pauli_i = pauli(i)
    pauli_j = pauli(j)
    outer = np.zeros((4, 4), dtype=np.complex)
    outer[:2, :2] = pauli_i
    outer[2:, 2:] = pauli_j
    if tensor:
        outer.shape = (2, 2, 2, 2)
    return outer

# def paulitwo_left(i):
    # return np.kron(pauli(i), pauli(0)) 

# def paulitwo_right(i):
    # return np.kron(pauli(0), pauli(i))

# def newrho_DK(Lket, theta_ij):
    # Lbra = np.conjugate(L_before)
    # theta_star = np.conjugate(theta_ij)
    # in_bracket = scon(Lbra, Lket, theta_ij, theta_star,
                # [1], [1], [2, 3, 1, 

Example 18

def test_generic_methods(self):
        # Tests some MaskedArray methods.
        a = array([1, 3, 2])
        assert_equal(a.any(), a._data.any())
        assert_equal(a.all(), a._data.all())
        assert_equal(a.argmax(), a._data.argmax())
        assert_equal(a.argmin(), a._data.argmin())
        assert_equal(a.choose(0, 1, 2, 3, 4), a._data.choose(0, 1, 2, 3, 4))
        assert_equal(a.compress([1, 0, 1]), a._data.compress([1, 0, 1]))
        assert_equal(a.conj(), a._data.conj())
        assert_equal(a.conjugate(), a._data.conjugate())

        m = array([[1, 2], [3, 4]])
        assert_equal(m.diagonal(), m._data.diagonal())
        assert_equal(a.sum(), a._data.sum())
        assert_equal(a.take([1, 2]), a._data.take([1, 2]))
        assert_equal(m.transpose(), m._data.transpose()) 

Example 19

def test_testArrayMethods(self):
        a = array([1, 3, 2])
        self.assertTrue(eq(a.any(), a._data.any()))
        self.assertTrue(eq(a.all(), a._data.all()))
        self.assertTrue(eq(a.argmax(), a._data.argmax()))
        self.assertTrue(eq(a.argmin(), a._data.argmin()))
        self.assertTrue(eq(a.choose(0, 1, 2, 3, 4),
                           a._data.choose(0, 1, 2, 3, 4)))
        self.assertTrue(eq(a.compress([1, 0, 1]), a._data.compress([1, 0, 1])))
        self.assertTrue(eq(a.conj(), a._data.conj()))
        self.assertTrue(eq(a.conjugate(), a._data.conjugate()))
        m = array([[1, 2], [3, 4]])
        self.assertTrue(eq(m.diagonal(), m._data.diagonal()))
        self.assertTrue(eq(a.sum(), a._data.sum()))
        self.assertTrue(eq(a.take([1, 2]), a._data.take([1, 2])))
        self.assertTrue(eq(m.transpose(), m._data.transpose())) 

Example 20

def qisunitary(self,op):
		mat = op[1]
		(r,c) = mat.shape
		if r != c:
			return False
		invmat = np.conjugate(np.transpose(mat))
		pmat = np.asarray(mat * invmat)
		for i in range(r):
			for j in range(c):
				if i != j:
					if np.absolute(pmat[i][j]) > self.maxerr:
						return False
				else:
					if np.absolute(pmat[i][j]-1.0) > self.maxerr:
						return False
		return True 

Example 21

def __init__(self, gain):
        self.gain = gain * np.pi  # pi term puts angle output to pi range

        # components
        self.conjugate = Conjugate()
        self.complex_mult = ComplexMultiply()
        self.angle = Angle()
        self.out = Sfix()

        # specify component delay
        self._delay = self.conjugate._delay + \
                      self.complex_mult._delay + \
                      self.angle._delay + 1

        # constants
        self.gain_sfix = Const(Sfix(self.gain, 3, -14)) 

Example 22

def test_generic_methods(self):
        # Tests some MaskedArray methods.
        a = array([1, 3, 2])
        assert_equal(a.any(), a._data.any())
        assert_equal(a.all(), a._data.all())
        assert_equal(a.argmax(), a._data.argmax())
        assert_equal(a.argmin(), a._data.argmin())
        assert_equal(a.choose(0, 1, 2, 3, 4), a._data.choose(0, 1, 2, 3, 4))
        assert_equal(a.compress([1, 0, 1]), a._data.compress([1, 0, 1]))
        assert_equal(a.conj(), a._data.conj())
        assert_equal(a.conjugate(), a._data.conjugate())

        m = array([[1, 2], [3, 4]])
        assert_equal(m.diagonal(), m._data.diagonal())
        assert_equal(a.sum(), a._data.sum())
        assert_equal(a.take([1, 2]), a._data.take([1, 2]))
        assert_equal(m.transpose(), m._data.transpose()) 

Example 23

def test_testArrayMethods(self):
        a = array([1, 3, 2])
        self.assertTrue(eq(a.any(), a._data.any()))
        self.assertTrue(eq(a.all(), a._data.all()))
        self.assertTrue(eq(a.argmax(), a._data.argmax()))
        self.assertTrue(eq(a.argmin(), a._data.argmin()))
        self.assertTrue(eq(a.choose(0, 1, 2, 3, 4),
                           a._data.choose(0, 1, 2, 3, 4)))
        self.assertTrue(eq(a.compress([1, 0, 1]), a._data.compress([1, 0, 1])))
        self.assertTrue(eq(a.conj(), a._data.conj()))
        self.assertTrue(eq(a.conjugate(), a._data.conjugate()))
        m = array([[1, 2], [3, 4]])
        self.assertTrue(eq(m.diagonal(), m._data.diagonal()))
        self.assertTrue(eq(a.sum(), a._data.sum()))
        self.assertTrue(eq(a.take([1, 2]), a._data.take([1, 2])))
        self.assertTrue(eq(m.transpose(), m._data.transpose())) 

Example 24

def test_poly(self):
        assert_array_almost_equal(np.poly([3, -np.sqrt(2), np.sqrt(2)]),
                                  [1, -3, -2, 6])

        # From matlab docs
        A = [[1, 2, 3], [4, 5, 6], [7, 8, 0]]
        assert_array_almost_equal(np.poly(A), [1, -6, -72, -27])

        # Should produce real output for perfect conjugates
        assert_(np.isrealobj(np.poly([+1.082j, +2.613j, -2.613j, -1.082j])))
        assert_(np.isrealobj(np.poly([0+1j, -0+-1j, 1+2j,
                                      1-2j, 1.+3.5j, 1-3.5j])))
        assert_(np.isrealobj(np.poly([1j, -1j, 1+2j, 1-2j, 1+3j, 1-3.j])))
        assert_(np.isrealobj(np.poly([1j, -1j, 1+2j, 1-2j])))
        assert_(np.isrealobj(np.poly([1j, -1j, 2j, -2j])))
        assert_(np.isrealobj(np.poly([1j, -1j])))
        assert_(np.isrealobj(np.poly([1, -1])))

        assert_(np.iscomplexobj(np.poly([1j, -1.0000001j])))

        np.random.seed(42)
        a = np.random.randn(100) + 1j*np.random.randn(100)
        assert_(np.isrealobj(np.poly(np.concatenate((a, np.conjugate(a)))))) 

Example 25

def test_generic_methods(self):
        # Tests some MaskedArray methods.
        a = array([1, 3, 2])
        assert_equal(a.any(), a._data.any())
        assert_equal(a.all(), a._data.all())
        assert_equal(a.argmax(), a._data.argmax())
        assert_equal(a.argmin(), a._data.argmin())
        assert_equal(a.choose(0, 1, 2, 3, 4), a._data.choose(0, 1, 2, 3, 4))
        assert_equal(a.compress([1, 0, 1]), a._data.compress([1, 0, 1]))
        assert_equal(a.conj(), a._data.conj())
        assert_equal(a.conjugate(), a._data.conjugate())

        m = array([[1, 2], [3, 4]])
        assert_equal(m.diagonal(), m._data.diagonal())
        assert_equal(a.sum(), a._data.sum())
        assert_equal(a.take([1, 2]), a._data.take([1, 2]))
        assert_equal(m.transpose(), m._data.transpose()) 

Example 26

def test_testArrayMethods(self):
        a = array([1, 3, 2])
        self.assertTrue(eq(a.any(), a._data.any()))
        self.assertTrue(eq(a.all(), a._data.all()))
        self.assertTrue(eq(a.argmax(), a._data.argmax()))
        self.assertTrue(eq(a.argmin(), a._data.argmin()))
        self.assertTrue(eq(a.choose(0, 1, 2, 3, 4),
                           a._data.choose(0, 1, 2, 3, 4)))
        self.assertTrue(eq(a.compress([1, 0, 1]), a._data.compress([1, 0, 1])))
        self.assertTrue(eq(a.conj(), a._data.conj()))
        self.assertTrue(eq(a.conjugate(), a._data.conjugate()))
        m = array([[1, 2], [3, 4]])
        self.assertTrue(eq(m.diagonal(), m._data.diagonal()))
        self.assertTrue(eq(a.sum(), a._data.sum()))
        self.assertTrue(eq(a.take([1, 2]), a._data.take([1, 2])))
        self.assertTrue(eq(m.transpose(), m._data.transpose())) 

Example 27

def get_pmf_xi(self):
        """Return the values of the variable ``xi``.

        The components ``xi`` make up the probability mass function, i.e.
        :math:`\\xi(k) = pmf(k) = Pr(X = k)`.
        """
        chi = np.empty(self.number_trials + 1, dtype=complex)
        chi[0] = 1
        half_number_trials = int(
            self.number_trials / 2 + self.number_trials % 2)
        # set first half of chis:
        chi[1:half_number_trials + 1] = self.get_chi(
            np.arange(1, half_number_trials + 1))
        # set second half of chis:
        chi[half_number_trials + 1:self.number_trials + 1] = np.conjugate(
            chi[1:self.number_trials - half_number_trials + 1] [::-1])
        chi /= self.number_trials + 1
        xi = np.fft.fft(chi)
        if self.check_xi_are_real(xi):
            xi = xi.real
        else:
            raise TypeError("pmf / xi values have to be real.")
        xi += np.finfo(type(xi[0])).eps
        return xi 

Example 28

def test_generic_methods(self):
        # Tests some MaskedArray methods.
        a = array([1, 3, 2])
        assert_equal(a.any(), a._data.any())
        assert_equal(a.all(), a._data.all())
        assert_equal(a.argmax(), a._data.argmax())
        assert_equal(a.argmin(), a._data.argmin())
        assert_equal(a.choose(0, 1, 2, 3, 4), a._data.choose(0, 1, 2, 3, 4))
        assert_equal(a.compress([1, 0, 1]), a._data.compress([1, 0, 1]))
        assert_equal(a.conj(), a._data.conj())
        assert_equal(a.conjugate(), a._data.conjugate())

        m = array([[1, 2], [3, 4]])
        assert_equal(m.diagonal(), m._data.diagonal())
        assert_equal(a.sum(), a._data.sum())
        assert_equal(a.take([1, 2]), a._data.take([1, 2]))
        assert_equal(m.transpose(), m._data.transpose()) 

Example 29

def test_testArrayMethods(self):
        a = array([1, 3, 2])
        self.assertTrue(eq(a.any(), a._data.any()))
        self.assertTrue(eq(a.all(), a._data.all()))
        self.assertTrue(eq(a.argmax(), a._data.argmax()))
        self.assertTrue(eq(a.argmin(), a._data.argmin()))
        self.assertTrue(eq(a.choose(0, 1, 2, 3, 4),
                           a._data.choose(0, 1, 2, 3, 4)))
        self.assertTrue(eq(a.compress([1, 0, 1]), a._data.compress([1, 0, 1])))
        self.assertTrue(eq(a.conj(), a._data.conj()))
        self.assertTrue(eq(a.conjugate(), a._data.conjugate()))
        m = array([[1, 2], [3, 4]])
        self.assertTrue(eq(m.diagonal(), m._data.diagonal()))
        self.assertTrue(eq(a.sum(), a._data.sum()))
        self.assertTrue(eq(a.take([1, 2]), a._data.take([1, 2])))
        self.assertTrue(eq(m.transpose(), m._data.transpose())) 

Example 30

def _solve_lens_equation(self, complex_value):
        """
        Solve the lens equation for the given point (in complex coordinates).
        """
        complex_conjugate = np.conjugate(complex_value)
        return complex_value - (1. / (1. + self.q)) * (
            (1./complex_conjugate) + (self.q / (complex_conjugate - self.s))) 

Example 31

def _jacobian_determinant_ok_WM95(self, source_x, source_y):
        """determinants of lens equation Jacobian for verified roots"""
        roots_ok_bar = np.conjugate(self._polynomial_roots_ok_WM95(
                                   source_x=source_x, source_y=source_y))
        # Variable X_bar is conjugate of variable X.
        denominator_1 = self._position_z1_WM95 - roots_ok_bar
        add_1 = self.mass_1 / denominator_1**2
        denominator_2 = self._position_z2_WM95 - roots_ok_bar
        add_2 = self.mass_2 / denominator_2**2
        derivative = add_1 + add_2
        # Can we use Utils.complex_fsum()-like function here?

        return 1.-derivative*np.conjugate(derivative)
        # Can we use Utils.complex_fsum()-like function here? 

Example 32

def adj(self,x):
        """Returns Hermitian adjoint of a matrix"""
        assert x.shape[0] == x.shape[1]
        return np.conjugate(x).T 

Example 33

def test_conjugate(self):
        a = np.array([1-1j, 1+1j, 23+23.0j])
        ac = a.conj()
        assert_equal(a.real, ac.real)
        assert_equal(a.imag, -ac.imag)
        assert_equal(ac, a.conjugate())
        assert_equal(ac, np.conjugate(a))

        a = np.array([1-1j, 1+1j, 23+23.0j], 'F')
        ac = a.conj()
        assert_equal(a.real, ac.real)
        assert_equal(a.imag, -ac.imag)
        assert_equal(ac, a.conjugate())
        assert_equal(ac, np.conjugate(a))

        a = np.array([1, 2, 3])
        ac = a.conj()
        assert_equal(a, ac)
        assert_equal(ac, a.conjugate())
        assert_equal(ac, np.conjugate(a))

        a = np.array([1.0, 2.0, 3.0])
        ac = a.conj()
        assert_equal(a, ac)
        assert_equal(ac, a.conjugate())
        assert_equal(ac, np.conjugate(a))

        a = np.array([1-1j, 1+1j, 1, 2.0], object)
        ac = a.conj()
        assert_equal(ac, [k.conjugate() for k in a])
        assert_equal(ac, a.conjugate())
        assert_equal(ac, np.conjugate(a))

        a = np.array([1-1j, 1, 2.0, 'f'], object)
        assert_raises(AttributeError, lambda: a.conj())
        assert_raises(AttributeError, lambda: a.conjugate()) 

Example 34

def test_var_values(self):
        for mat in [self.rmat, self.cmat, self.omat]:
            for axis in [0, 1, None]:
                msqr = _mean(mat * mat.conj(), axis=axis)
                mean = _mean(mat, axis=axis)
                tgt = msqr - mean * mean.conjugate()
                res = _var(mat, axis=axis)
                assert_almost_equal(res, tgt) 

Example 35

def test_oddfeatures_1(self):
        # Test of other odd features
        x = arange(20)
        x = x.reshape(4, 5)
        x.flat[5] = 12
        assert_(x[1, 0] == 12)
        z = x + 10j * x
        assert_equal(z.real, x)
        assert_equal(z.imag, 10 * x)
        assert_equal((z * conjugate(z)).real, 101 * x * x)
        z.imag[...] = 0.0

        x = arange(10)
        x[3] = masked
        assert_(str(x[3]) == str(masked))
        c = x >= 8
        assert_(count(where(c, masked, masked)) == 0)
        assert_(shape(where(c, masked, masked)) == c.shape)

        z = masked_where(c, x)
        assert_(z.dtype is x.dtype)
        assert_(z[3] is masked)
        assert_(z[4] is not masked)
        assert_(z[7] is not masked)
        assert_(z[8] is masked)
        assert_(z[9] is masked)
        assert_equal(x, z) 

Example 36

def test_basic_ufuncs(self):
        # Test various functions such as sin, cos.
        (x, y, a10, m1, m2, xm, ym, z, zm, xf) = self.d
        assert_equal(np.cos(x), cos(xm))
        assert_equal(np.cosh(x), cosh(xm))
        assert_equal(np.sin(x), sin(xm))
        assert_equal(np.sinh(x), sinh(xm))
        assert_equal(np.tan(x), tan(xm))
        assert_equal(np.tanh(x), tanh(xm))
        assert_equal(np.sqrt(abs(x)), sqrt(xm))
        assert_equal(np.log(abs(x)), log(xm))
        assert_equal(np.log10(abs(x)), log10(xm))
        assert_equal(np.exp(x), exp(xm))
        assert_equal(np.arcsin(z), arcsin(zm))
        assert_equal(np.arccos(z), arccos(zm))
        assert_equal(np.arctan(z), arctan(zm))
        assert_equal(np.arctan2(x, y), arctan2(xm, ym))
        assert_equal(np.absolute(x), absolute(xm))
        assert_equal(np.angle(x + 1j*y), angle(xm + 1j*ym))
        assert_equal(np.angle(x + 1j*y, deg=True), angle(xm + 1j*ym, deg=True))
        assert_equal(np.equal(x, y), equal(xm, ym))
        assert_equal(np.not_equal(x, y), not_equal(xm, ym))
        assert_equal(np.less(x, y), less(xm, ym))
        assert_equal(np.greater(x, y), greater(xm, ym))
        assert_equal(np.less_equal(x, y), less_equal(xm, ym))
        assert_equal(np.greater_equal(x, y), greater_equal(xm, ym))
        assert_equal(np.conjugate(x), conjugate(xm)) 

Example 37

def test_testUfuncs1(self):
        # Test various functions such as sin, cos.
        (x, y, a10, m1, m2, xm, ym, z, zm, xf, s) = self.d
        self.assertTrue(eq(np.cos(x), cos(xm)))
        self.assertTrue(eq(np.cosh(x), cosh(xm)))
        self.assertTrue(eq(np.sin(x), sin(xm)))
        self.assertTrue(eq(np.sinh(x), sinh(xm)))
        self.assertTrue(eq(np.tan(x), tan(xm)))
        self.assertTrue(eq(np.tanh(x), tanh(xm)))
        with np.errstate(divide='ignore', invalid='ignore'):
            self.assertTrue(eq(np.sqrt(abs(x)), sqrt(xm)))
            self.assertTrue(eq(np.log(abs(x)), log(xm)))
            self.assertTrue(eq(np.log10(abs(x)), log10(xm)))
        self.assertTrue(eq(np.exp(x), exp(xm)))
        self.assertTrue(eq(np.arcsin(z), arcsin(zm)))
        self.assertTrue(eq(np.arccos(z), arccos(zm)))
        self.assertTrue(eq(np.arctan(z), arctan(zm)))
        self.assertTrue(eq(np.arctan2(x, y), arctan2(xm, ym)))
        self.assertTrue(eq(np.absolute(x), absolute(xm)))
        self.assertTrue(eq(np.equal(x, y), equal(xm, ym)))
        self.assertTrue(eq(np.not_equal(x, y), not_equal(xm, ym)))
        self.assertTrue(eq(np.less(x, y), less(xm, ym)))
        self.assertTrue(eq(np.greater(x, y), greater(xm, ym)))
        self.assertTrue(eq(np.less_equal(x, y), less_equal(xm, ym)))
        self.assertTrue(eq(np.greater_equal(x, y), greater_equal(xm, ym)))
        self.assertTrue(eq(np.conjugate(x), conjugate(xm)))
        self.assertTrue(eq(np.concatenate((x, y)), concatenate((xm, ym))))
        self.assertTrue(eq(np.concatenate((x, y)), concatenate((x, y))))
        self.assertTrue(eq(np.concatenate((x, y)), concatenate((xm, y))))
        self.assertTrue(eq(np.concatenate((x, y, x)), concatenate((x, ym, x)))) 

Example 38

def test_testUfuncRegression(self):
        f_invalid_ignore = [
            'sqrt', 'arctanh', 'arcsin', 'arccos',
            'arccosh', 'arctanh', 'log', 'log10', 'divide',
            'true_divide', 'floor_divide', 'remainder', 'fmod']
        for f in ['sqrt', 'log', 'log10', 'exp', 'conjugate',
                  'sin', 'cos', 'tan',
                  'arcsin', 'arccos', 'arctan',
                  'sinh', 'cosh', 'tanh',
                  'arcsinh',
                  'arccosh',
                  'arctanh',
                  'absolute', 'fabs', 'negative',
                  'floor', 'ceil',
                  'logical_not',
                  'add', 'subtract', 'multiply',
                  'divide', 'true_divide', 'floor_divide',
                  'remainder', 'fmod', 'hypot', 'arctan2',
                  'equal', 'not_equal', 'less_equal', 'greater_equal',
                  'less', 'greater',
                  'logical_and', 'logical_or', 'logical_xor']:
            try:
                uf = getattr(umath, f)
            except AttributeError:
                uf = getattr(fromnumeric, f)
            mf = getattr(np.ma, f)
            args = self.d[:uf.nin]
            with np.errstate():
                if f in f_invalid_ignore:
                    np.seterr(invalid='ignore')
                if f in ['arctanh', 'log', 'log10']:
                    np.seterr(divide='ignore')
                ur = uf(*args)
                mr = mf(*args)
            self.assertTrue(eq(ur.filled(0), mr.filled(0), f))
            self.assertTrue(eqmask(ur.mask, mr.mask)) 

Example 39

def power_process(data, sfreq, toffset, modulus, integration, log_scale, zscale, title):
    """ Break power by modulus and display each block. Integration here acts
    as a pure average on the power level data.
    """
    if modulus:
        block = 0
        block_size = integration * modulus
        block_toffset = toffset
        while block < len(data) / block_size:

            vblock = data[block * block_size:block * block_size + modulus]
            pblock = (vblock * numpy.conjugate(vblock)).real

            # complete integration
            for idx in range(1, integration):

                vblock = data[block * block_size + idx *
                              modulus:block * block_size + idx * modulus + modulus]
                pblock += (vblock * numpy.conjugate(vblock)).real

            pblock /= integration

            yield power_plot(pblock, sfreq, block_toffset, log_scale, zscale, title)

            block += 1
            block_toffset += block_size / sfreq

    else:
        pdata = (data * numpy.conjugate(data)).real
        yield power_plot(pdata, sfreq, toffset, log_scale, zscale, title) 

Example 40

def MatrixElements(m,wavelength,diameter,mu):
#  http://pymiescatt.readthedocs.io/en/latest/forward.html#MatrixElements
  x = np.pi*diameter/wavelength
  # B&H eqs. 4.77, where mu=cos(theta)
  S1,S2 = MieS1S2(m,x,mu)
  S11 = 0.5*(np.abs(S2)**2+np.abs(S1)**2)
  S12 = 0.5*(np.abs(S2)**2-np.abs(S1)**2)
  S33 = 0.5*(np.conjugate(S2)*S1+S2*np.conjugate(S1))
  S34 = 0.5j*(S1*np.conjugate(S2)-S2*np.conjugate(S1))
  return S11, S12, S33, S34 

Example 41

def CoreShellS1S2(mCore,mShell,xCore,xShell,mu):
#  http://pymiescatt.readthedocs.io/en/latest/forwardCS.html#CoreShellS1S2
  nmax = np.round(2+xShell+4*(xShell**(1/3)))
  an,bn = CoreShell_ab(mCore,mShell,xCore,xShell)
  pin,taun = MiePiTau(mu,nmax)
  n = np.arange(1,int(nmax)+1)
  n2 = (2*n+1)/(n*(n+1))
  pin *= n2
  taun *= n2
  S1=np.sum(an*np.conjugate(pin))+np.sum(bn*np.conjugate(taun))
  S2=np.sum(an*np.conjugate(taun))+np.sum(bn*np.conjugate(pin))
  return S1,S2 

Example 42

def CoreShellMatrixElements(mCore,mShell,xCore,xShell,mu):
#  http://pymiescatt.readthedocs.io/en/latest/forwardCS.html#CoreShellMatrixElements
  S1,S2 = CoreShellS1S2(mCore,mShell,xCore,xShell,mu)
  S11 = 0.5*(np.abs(S2)**2+np.abs(S1)**2)
  S12 = 0.5*(np.abs(S2)**2-np.abs(S1)**2)
  S33 = 0.5*(np.conjugate(S2)*S1+S2*np.conjugate(S1))
  S34 = 0.5j*(S1*np.conjugate(S2)-S2*np.conjugate(S1))
  return S11, S12, S33, S34 

Example 43

def MatrixElements(m,wavelength,diameter,mu):
#  http://pymiescatt.readthedocs.io/en/latest/forward.html#MatrixElements
  x = np.pi*diameter/wavelength
  # B&H eqs. 4.77, where mu=cos(theta)
  S1,S2 = MieS1S2(m,x,mu)
  S11 = 0.5*(np.abs(S2)**2+np.abs(S1)**2)
  S12 = 0.5*(np.abs(S2)**2-np.abs(S1)**2)
  S33 = 0.5*(np.conjugate(S2)*S1+S2*np.conjugate(S1))
  S34 = 0.5j*(S1*np.conjugate(S2)-S2*np.conjugate(S1))
  return S11, S12, S33, S34 

Example 44

def CoreShellS1S2(mCore,mShell,xCore,xShell,mu):
#  http://pymiescatt.readthedocs.io/en/latest/forwardCS.html#CoreShellS1S2
  nmax = np.round(2+xShell+4*(xShell**(1/3)))
  an,bn = CoreShell_ab(mCore,mShell,xCore,xShell)
  pin,taun = MiePiTau(mu,nmax)
  n = np.arange(1,int(nmax)+1)
  n2 = (2*n+1)/(n*(n+1))
  pin *= n2
  taun *= n2
  S1=np.sum(an*np.conjugate(pin))+np.sum(bn*np.conjugate(taun))
  S2=np.sum(an*np.conjugate(taun))+np.sum(bn*np.conjugate(pin))
  return S1,S2 

Example 45

def CoreShellMatrixElements(mCore,mShell,xCore,xShell,mu):
#  http://pymiescatt.readthedocs.io/en/latest/forwardCS.html#CoreShellMatrixElements
  S1,S2 = CoreShellS1S2(mCore,mShell,xCore,xShell,mu)
  S11 = 0.5*(np.abs(S2)**2+np.abs(S1)**2)
  S12 = 0.5*(np.abs(S2)**2-np.abs(S1)**2)
  S33 = 0.5*(np.conjugate(S2)*S1+S2*np.conjugate(S1))
  S34 = 0.5j*(S1*np.conjugate(S2)-S2*np.conjugate(S1))
  return S11, S12, S33, S34 

Example 46

def sum_visibility(vis: Visibility, direction: SkyCoord) -> numpy.array:
    """ Direct Fourier summation in a given direction

    :param vis: Visibility to be summed
    :param direction: Direction of summation
    :return: flux[nch,npol], weight[nch,pol]
    """
    # TODO: Convert to Visibility or remove?
    
    assert isinstance(vis, Visibility) or isinstance(vis, BlockVisibility), vis
    
    svis = copy_visibility(vis)
    
    l, m, n = skycoord_to_lmn(direction, svis.phasecentre)
    phasor = numpy.conjugate(simulate_point(svis.uvw, l, m))
    
    # Need to put correct mapping here
    _, frequency = get_frequency_map(svis, None)
    
    frequency = list(frequency)
    
    nchan = max(frequency) + 1
    npol = svis.polarisation_frame.npol
    
    flux = numpy.zeros([nchan, npol])
    weight = numpy.zeros([nchan, npol])
    
    coords = svis.vis, svis.weight, phasor, list(frequency)
    for v, wt, p, ic in zip(*coords):
        for pol in range(npol):
            flux[ic, pol] += numpy.real(wt[pol] * v[pol] * p)
            weight[ic, pol] += wt[pol]
    
    flux[weight > 0.0] = flux[weight > 0.0] / weight[weight > 0.0]
    flux[weight <= 0.0] = 0.0
    return flux, weight 

Example 47

def gain_substitution_vector(gain, x, xwt):
    nants, nchan, nrec, _ = gain.shape
    newgain = numpy.ones_like(gain, dtype='complex')
    if nrec > 0:
        newgain[..., 0, 1] = 0.0
        newgain[..., 1, 0] = 0.0
    
    gwt = numpy.zeros_like(gain, dtype='float')
    
    # We are going to work with Jones 2x2 matrix formalism so everything has to be
    # converted to that format
    x = x.reshape(nants, nants, nchan, nrec, nrec)
    xwt = xwt.reshape(nants, nants, nchan, nrec, nrec)
    
    if nrec > 0:
        gain[..., 0, 1] = 0.0
        gain[..., 1, 0] = 0.0
    
    for ant1 in range(nants):
        for chan in range(nchan):
            # Loop over e.g. 'RR', 'LL, or 'xx', 'YY' ignoring cross terms
            for rec in range(nrec):
                top = numpy.sum(x[:, ant1, chan, rec, rec] *
                                gain[:, chan, rec, rec] * xwt[:, ant1, chan, rec, rec], axis=0)
                bot = numpy.sum((gain[:, chan, rec, rec] * numpy.conjugate(gain[:, chan, rec, rec]) *
                                 xwt[:, ant1, chan, rec, rec]).real, axis=0)
                
                if bot > 0.0:
                    newgain[ant1, chan, rec, rec] = top / bot
                    gwt[ant1, chan, rec, rec] = bot
                else:
                    newgain[ant1, chan, rec, rec] = 0.0
                    gwt[ant1, chan, rec, rec] = 0.0
    
    return newgain, gwt 

Example 48

def gain_substitution_matrix(gain, x, xwt):
    nants, nchan, nrec, _ = gain.shape
    newgain = numpy.ones_like(gain, dtype='complex')
    gwt = numpy.zeros_like(gain, dtype='float')
    
    # We are going to work with Jones 2x2 matrix formalism so everything has to be
    # converted to that format
    x = x.reshape(nants, nants, nchan, nrec, nrec)
    xwt = xwt.reshape(nants, nants, nchan, nrec, nrec)
    
    # Write these loops out explicitly. Derivation of these vector equations is tedious but they are
    # structurally identical to the scalar case with the following changes
    # Vis -> 2x2 coherency vector, g-> 2x2 Jones matrix, *-> matmul, conjugate->Hermitean transpose (.H)
    for ant1 in range(nants):
        for chan in range(nchan):
            top = 0.0
            bot = 0.0
            for ant2 in range(nants):
                if ant1 != ant2:
                    xmat = x[ant2, ant1, chan]
                    xwtmat = xwt[ant2, ant1, chan]
                    g2 = gain[ant2, chan]
                    top += xmat * xwtmat * g2
                    bot += numpy.conjugate(g2) * xwtmat * g2
            newgain[ant1, chan][bot > 0.0] = top[bot > 0.0] / bot[bot > 0.0]
            newgain[ant1, chan][bot <= 0.0] = 0.0
            gwt[ant1, chan] = bot.real
    return newgain, gwt 

Example 49

def solution_residual_scalar(gain, x, xwt):
    """Calculate residual across all baselines of gain for point source equivalent visibilities
    
    :param gain: gain [nant, ...]
    :param x: Point source equivalent visibility [nant, ...]
    :param xwt: Point source equivalent weight [nant, ...]
    :return: residual[...]
    """
    
    nants, nchan, nrec, _ = gain.shape
    x = x.reshape(nants, nants, nchan, nrec, nrec)
    
    xwt = xwt.reshape(nants, nants, nchan, nrec, nrec)
    
    residual = numpy.zeros([nchan, nrec, nrec])
    sumwt = numpy.zeros([nchan, nrec, nrec])
    
    for ant1 in range(nants):
        for ant2 in range(nants):
            for chan in range(nchan):
                error = x[ant2, ant1, chan, 0, 0] - \
                    gain[ant1, chan, 0, 0] * numpy.conjugate(gain[ant2, chan, 0, 0])
                residual += (error * xwt[ant2, ant1, chan, 0, 0] * numpy.conjugate(error)).real
                sumwt += xwt[ant2, ant1, chan, 0, 0]
    
    residual[sumwt > 0.0] = numpy.sqrt(residual[sumwt > 0.0] / sumwt[sumwt > 0.0])
    residual[sumwt <= 0.0] = 0.0
        
    return residual 

Example 50

def solution_residual_vector(gain, x, xwt):
    """Calculate residual across all baselines of gain for point source equivalent visibilities
    
    Vector case i.e. off-diagonals of gains are zero

    :param gain: gain [nant, ...]
    :param x: Point source equivalent visibility [nant, ...]
    :param xwt: Point source equivalent weight [nant, ...]
    :return: residual[...]
    """
    
    nants, nchan, nrec, _ = gain.shape
    x = x.reshape(nants, nants, nchan, nrec, nrec)
    x[..., 1, 0] = 0.0
    x[..., 0, 1] = 0.0
    
    xwt = xwt.reshape(nants, nants, nchan, nrec, nrec)
    xwt[..., 1, 0] = 0.0
    xwt[..., 0, 1] = 0.0
    
    residual = numpy.zeros([nchan, nrec, nrec])
    sumwt = numpy.zeros([nchan, nrec, nrec])
    
    for ant1 in range(nants):
        for ant2 in range(nants):
            for chan in range(nchan):
                for rec in range(nrec):
                    error = x[ant2, ant1, chan, rec, rec] - \
                        gain[ant1, chan, rec, rec] * numpy.conjugate(gain[ant2, chan, rec, rec])
                    residual += (error * xwt[ant2, ant1, chan, rec, rec] * numpy.conjugate(error)).real
                    sumwt += xwt[ant2, ant1, chan, rec, rec]

    residual[sumwt > 0.0] = numpy.sqrt(residual[sumwt > 0.0] / sumwt[sumwt > 0.0])
    residual[sumwt <= 0.0] = 0.0

    return residual