Python numpy.result_type函数代码示例

    def ShiftConvNumbaFFT(h, N, M, xdtype=np.complex_, powerof2=True):
        # implements Doppler filter:
        # y[n, p] = SUM_k (exp(2*pi*j*n*(k - (L-1))/N) * h[k]) * x[p - k]
        #         = SUM_k (exp(-2*pi*j*n*k/N) * s*[k]) * x[p - (L-1) + k]
        L = len(h)
        outlen = M + L - 1
        nfft = outlen
        if powerof2:
            nfft = pow2(nfft)

        dopplermat = np.exp(2*np.pi*1j*np.arange(N)[:, np.newaxis]*(np.arange(L) - (L - 1))/N)
        dopplermat.astype(np.result_type(h.dtype, np.complex64)) # cast to complex type with precision of h
        hbank = h*dopplermat
        # speed not critical here, just use numpy fft
        hbankpad = zero_pad(hbank, nfft)
        H = np.fft.fft(hbankpad) / nfft # divide by nfft b/c FFTW's ifft does not do this

        xcdtype = np.result_type(xdtype, np.complex64) # cast to complex type with precision of x
        xpad = pyfftw.n_byte_align(np.zeros(nfft, xcdtype), 16)
        X = pyfftw.n_byte_align(np.zeros(nfft, xcdtype), 16)
        xfft = pyfftw.FFTW(xpad, X, threads=_THREADS)

        ydtype = np.result_type(H.dtype, xcdtype)
        Y = pyfftw.n_byte_align_empty(H.shape, 16, ydtype)
        y = pyfftw.n_byte_align_empty(H.shape, 16, ydtype)
        ifft = pyfftw.FFTW(Y, y, direction='FFTW_BACKWARD', threads=_THREADS)

        xtype = numba.__getattribute__(str(np.dtype(xdtype)))

        #htype = numba.__getattribute__(str(H.dtype))
        #xctype = numba.__getattribute__(str(X.dtype))
        #ytype = numba.__getattribute__(str(Y.dtype))
        #@jit(argtypes=[htype[:, ::1], xctype[::1], ytype[:, ::1], xtype[::1]])
        #def fun(H, X, Y, x):
            #xpad[:M] = x
            #xfft.execute() # input is xpad, output is X
            #Y[:, :] = H*X # need expression optimized by numba but that writes into Y
            #ifft.execute() # input is Y, output is y

            #yc = np.array(y)[:, :outlen] # need a copy, which np.array provides
            #return yc

        #@dopplerbank_dec(h, N, M, nfft=nfft, H=H)
        #def shiftconv_numba_fft(x):
            #return fun(H, X, Y, x)

        #@jit(argtypes=[xtype[::1]])
        @jit
        def shiftconv_numba_fft(x):
            xpad[:M] = x
            xfft.execute() # input is xpad, output is X
            Y[:, :] = X*H # need expression optimized by numba but that writes into Y
            ifft.execute() # input is Y, output is y

            yc = np.array(y[:, :outlen]) # need a copy, which np.array provides
            return yc

        shiftconv_numba_fft = dopplerbank_dec(h, N, M, nfft=nfft, H=H)(shiftconv_numba_fft)

        return shiftconv_numba_fft

    def accum (self, key1, key2, value, weight=1):
        index1 = self._key1map[key1]
        index2 = self._key2map[key2]

        if self._m0 is None:
            self._m0 = np.zeros ((self.chunk0size, self.chunk0size), dtype=np.result_type (weight))
            self._m1 = np.zeros ((self.chunk0size, self.chunk0size), dtype=np.result_type (value, weight))
            self._m2 = np.zeros_like (self._m1)

        if index1 >= self._m0.shape[0]:
            self._m0 = np.concatenate ((self._m0, np.zeros_like (self._m0)), axis=0)
            self._m1 = np.concatenate ((self._m1, np.zeros_like (self._m1)), axis=0)
            self._m2 = np.concatenate ((self._m2, np.zeros_like (self._m2)), axis=0)

        if index2 >= self._m0.shape[1]:
            self._m0 = np.concatenate ((self._m0, np.zeros_like (self._m0)), axis=1)
            self._m1 = np.concatenate ((self._m1, np.zeros_like (self._m1)), axis=1)
            self._m2 = np.concatenate ((self._m2, np.zeros_like (self._m2)), axis=1)

        self._m0[index1,index2] += weight
        q = weight * value
        self._m1[index1,index2] += q
        q *= value
        self._m2[index1,index2] += q
        return self

def da_sub(daa, dab):
    """
    subtract 2 DataArrays as cleverly as possible:
      * keep the metadata of the first DA in the result
      * ensures the result has the right type so that no underflows happen
    returns (DataArray): the result of daa - dab
    """
    rt = numpy.result_type(daa, dab) # dtype of result of daa-dab

    dt = None # default is to let numpy decide
    if rt.kind == "f":
        # float should always be fine
        pass
    elif rt.kind in "iub":
        # underflow can happen (especially if unsigned)

        # find the worse case value (could be improved, but would be longer)
        worse_val = int(daa.min()) - int(dab.max())
        dt = numpy.result_type(rt, numpy.min_scalar_type(worse_val))
    else:
        # subtracting such a data is suspicious, but try anyway
        logging.warning("Subtraction on data of type %s unsupported", rt.name)

    res = numpy.subtract(daa, dab, dtype=dt) # metadata is copied from daa
    logging.debug("type = %s, %s", res.dtype.name, daa.dtype.name)
    return res

    def _take_with_fill(self, indices, fill_value=None):
        if fill_value is None:
            fill_value = self.dtype.na_value

        if indices.min() < -1:
            raise ValueError("Invalid value in 'indices'. Must be between -1 "
                             "and the length of the array.")

        if indices.max() >= len(self):
            raise IndexError("out of bounds value in 'indices'.")

        if len(self) == 0:
            # Empty... Allow taking only if all empty
            if (indices == -1).all():
                dtype = np.result_type(self.sp_values, fill_value)
                taken = np.empty_like(indices, dtype=dtype)
                taken.fill(fill_value)
                return taken
            else:
                raise IndexError('cannot do a non-empty take from an empty '
                                 'axes.')

        sp_indexer = self.sp_index.lookup_array(indices)

        if self.sp_index.npoints == 0:
            # Avoid taking from the empty self.sp_values
            taken = np.full(sp_indexer.shape, fill_value=fill_value,
                            dtype=np.result_type(fill_value))
        else:
            taken = self.sp_values.take(sp_indexer)

            # sp_indexer may be -1 for two reasons
            # 1.) we took for an index of -1 (new)
            # 2.) we took a value that was self.fill_value (old)
            new_fill_indices = indices == -1
            old_fill_indices = (sp_indexer == -1) & ~new_fill_indices

            # Fill in two steps.
            # Old fill values
            # New fill values
            # potentially coercing to a new dtype at each stage.

            m0 = sp_indexer[old_fill_indices] < 0
            m1 = sp_indexer[new_fill_indices] < 0

            result_type = taken.dtype

            if m0.any():
                result_type = np.result_type(result_type, self.fill_value)
                taken = taken.astype(result_type)
                taken[old_fill_indices] = self.fill_value

            if m1.any():
                result_type = np.result_type(result_type, fill_value)
                taken = taken.astype(result_type)
                taken[new_fill_indices] = fill_value

        return taken

    def NumbaFFTW(h, M, xdtype=np.complex_, powerof2=True):
        L = len(h)
        outlen = M + L - 1
        nfft = outlen
        if powerof2:
            nfft = pow2(nfft)

        outdtype = np.result_type(h.dtype, xdtype)
        fftdtype = np.result_type(outdtype, np.complex64) # output is always complex, promote using smallest

        # speed not critical here, just use numpy fft
        # cast to outdtype so we use same type of fft as when transforming x
        hpad = zero_pad(h, nfft).astype(outdtype)
        if np.iscomplexobj(hpad):
            H = np.fft.fft(hpad)
        else:
            H = np.fft.rfft(hpad)
        H = (H / nfft).astype(fftdtype) # divide by nfft b/c FFTW's ifft does not do this

        xpad = pyfftw.n_byte_align(np.zeros(nfft, outdtype), 16) # outdtype so same type fft as h->H
        X = pyfftw.n_byte_align(np.zeros(len(H), fftdtype), 16) # len(H) b/c rfft may be used
        xfft = pyfftw.FFTW(xpad, X, threads=_THREADS)

        y = pyfftw.n_byte_align_empty(nfft, 16, outdtype)
        ifft = pyfftw.FFTW(X, y, direction='FFTW_BACKWARD', threads=_THREADS)

        xtype = numba.__getattribute__(str(np.dtype(xdtype)))
        outtype = numba.__getattribute__(str(outdtype))
        ffttype = numba.__getattribute__(str(fftdtype))

        #@jit(restype=outtype[::1],
            #argtypes=[outtype[::1], ffttype[::1], ffttype[::1], outtype[::1], xtype[::1]])
        #def filt(xpad, X, H, y, x):
            #xpad[:M] = x
            #xfft.execute() # input in xpad, result in X
            #X[:] = H*X
            #ifft.execute() # input in X, result in y
            #yc = y[:outlen].copy()
            #return yc

        #@filter_dec(h, M, nfft=nfft, H=H)
        #def numba_fftw(x):
            #return filt(xpad, X, H, y, x)

        #@jit(argtypes=[xtype[::1]])
        @jit
        def numba_fftw(x):
            xpad[:M] = x
            xfft.execute() # input in xpad, result in X
            X[:] = H*X # want expression that is optimized by numba but writes into X
            ifft.execute() # input in X, result in y
            yc = y[:outlen].copy()
            return yc

        numba_fftw = filter_dec(h, M, nfft=nfft, H=H)(numba_fftw)

        return numba_fftw

    def test_result_type(self):
        self.check_promotion_cases(np.result_type)

        f64 = float64(0)
        c64 = complex64(0)
        ## Scalars do not coerce to complex if the value is real
        #assert_equal(np.result_type(c64,array([f64])), np.dtype(float64))
        # But they do if the value is complex
        assert_equal(np.result_type(complex64(3j),array([f64])),
                                                    np.dtype(complex128))

        # Scalars do coerce to complex even if the value is real
        # This is so "a+0j" can be reliably used to make something complex.
        assert_equal(np.result_type(c64,array([f64])), np.dtype(complex128))

    def _convert_list(self, value):
        """Convert a string into a typed numpy array.

        If it is not possible it returns a numpy string.
        """
        try:
            numpy_values = []
            values = value.split(" ")
            types = set([])
            for string_value in values:
                v = self._convert_scalar_value(string_value)
                numpy_values.append(v)
                types.add(v.dtype.type)

            result_type = numpy.result_type(*types)

            if issubclass(result_type.type, (numpy.string_, six.binary_type)):
                # use the raw data to create the result
                return numpy.string_(value)
            elif issubclass(result_type.type, (numpy.unicode_, six.text_type)):
                # use the raw data to create the result
                return numpy.unicode_(value)
            else:
                return numpy.array(numpy_values, dtype=result_type)
        except ValueError:
            return numpy.string_(value)

def test_upcast():
    a0 = csr_matrix([[np.pi, np.pi*1j], [3, 4]], dtype=complex)
    b0 = np.array([256+1j, 2**32], dtype=complex)

    for a_dtype in supported_dtypes:
        for b_dtype in supported_dtypes:
            msg = "(%r, %r)" % (a_dtype, b_dtype)

            if np.issubdtype(a_dtype, np.complexfloating):
                a = a0.copy().astype(a_dtype)
            else:
                a = a0.real.copy().astype(a_dtype)

            if np.issubdtype(b_dtype, np.complexfloating):
                b = b0.copy().astype(b_dtype)
            else:
                b = b0.real.copy().astype(b_dtype)

            if not (a_dtype == np.bool_ and b_dtype == np.bool_):
                c = np.zeros((2,), dtype=np.bool_)
                assert_raises(ValueError, _sparsetools.csr_matvec,
                              2, 2, a.indptr, a.indices, a.data, b, c)

            if ((np.issubdtype(a_dtype, np.complexfloating) and
                 not np.issubdtype(b_dtype, np.complexfloating)) or
                (not np.issubdtype(a_dtype, np.complexfloating) and
                 np.issubdtype(b_dtype, np.complexfloating))):
                c = np.zeros((2,), dtype=np.float64)
                assert_raises(ValueError, _sparsetools.csr_matvec,
                              2, 2, a.indptr, a.indices, a.data, b, c)

            c = np.zeros((2,), dtype=np.result_type(a_dtype, b_dtype))
            _sparsetools.csr_matvec(2, 2, a.indptr, a.indices, a.data, b, c)
            assert_allclose(c, np.dot(a.toarray(), b), err_msg=msg)

 def _normalize_vector_type(self, dtype):
     """Normalize the """
     if self.__at_least_32bits:
         if numpy.issubdtype(dtype, numpy.signedinteger):
             dtype = numpy.result_type(dtype, numpy.uint32)
         if numpy.issubdtype(dtype, numpy.unsignedinteger):
             dtype = numpy.result_type(dtype, numpy.uint32)
         elif numpy.issubdtype(dtype, numpy.floating):
             dtype = numpy.result_type(dtype, numpy.float32)
         elif numpy.issubdtype(dtype, numpy.complexfloating):
             dtype = numpy.result_type(dtype, numpy.complex64)
     if self.__signed_type:
         if numpy.issubdtype(dtype, numpy.unsignedinteger):
             signed = numpy.dtype("%s%i" % ('i', dtype.itemsize))
             dtype = numpy.result_type(dtype, signed)
     return dtype

 def matvec_transp(x):
     if x.shape != (nargout,):
         msg = 'Input has shape ' + str(x.shape)
         msg += ' instead of (%d,)' % self.nargout
         raise ValueError(msg)
     result_type = np.result_type(self.dtype, x.dtype)
     return np.zeros(nargin, dtype=result_type)

def rfft(a):
	n = a.shape[-1]
	b = a.reshape(np.prod(a.shape[:-1]),n)
	fb = np.empty((b.shape[0],b.shape[1]/2+1),dtype=np.result_type(a,0j))
	for i in range(b.shape[0]):
		fb[i] = myfft.rfft(b[i])*n**-0.5
	return np.reshape(fb, list(a.shape[:-1]) + [fb.shape[-1]])

def promote(*operands):
    """
    Take an arbitrary number of graph nodes and produce the promoted
    dtype by discarding all shape information and just looking at the
    measures.

    ::

        5, 5,  | int   |
        2, 5,  | int   |
        1, 3,  | float |

    >>> promote(IntNode(1), Op(IntNode(2))
    int
    >>> promote(FloatNode(1), Op(IntNode(2))
    float
    """

    # Looks something like this...

    # (ArrayNode, IntNode...) -> (dshape('2, int'), dshape('int'))
    # (dshape('2, int', dshape('int')) -> (dshape('int', dshape('int'))
    # (dshape('2, int', dshape('int')) -> (dtype('int', dtype('int'))

    types    = (op.simple_type() for op in operands if op is not None)
    measures = (extract_measure(t) for t in types)
    dtypes   = (to_numpy(m) for m in measures)

    promoted = np.result_type(*dtypes)
    datashape = CType.from_dtype(promoted)
    return datashape

def _contract_plain(mydf, mos, coulG, phase, max_memory):
    cell = mydf.cell
    moiT, mojT, mokT, molT = mos
    nmoi, nmoj, nmok, nmol = [x.shape[0] for x in mos]
    ngrids = moiT.shape[1]
    wcoulG = coulG * (cell.vol/ngrids)
    dtype = numpy.result_type(phase, *mos)
    eri = numpy.empty((nmoi*nmoj,nmok*nmol), dtype=dtype)

    blksize = int(min(max(nmoi,nmok), (max_memory*1e6/16 - eri.size)/2/ngrids/max(nmoj,nmol)+1))
    assert blksize > 0
    buf0 = numpy.empty((blksize,max(nmoj,nmol),ngrids), dtype=dtype)
    buf1 = numpy.ndarray((blksize,nmoj,ngrids), dtype=dtype, buffer=buf0)
    buf2 = numpy.ndarray((blksize,nmol,ngrids), dtype=dtype, buffer=buf0)
    for p0, p1 in lib.prange(0, nmoi, blksize):
        mo_pairs = numpy.einsum('ig,jg->ijg', moiT[p0:p1].conj()*phase,
                                mojT, out=buf1[:p1-p0])
        mo_pairs_G = tools.fft(mo_pairs.reshape(-1,ngrids), mydf.mesh)
        mo_pairs = None
        mo_pairs_G*= wcoulG
        v = tools.ifft(mo_pairs_G, mydf.mesh)
        mo_pairs_G = None
        v *= phase.conj()
        if dtype == numpy.double:
            v = numpy.asarray(v.real, order='C')
        for q0, q1 in lib.prange(0, nmok, blksize):
            mo_pairs = numpy.einsum('ig,jg->ijg', mokT[q0:q1].conj(),
                                    molT, out=buf2[:q1-q0])
            eri[p0*nmoj:p1*nmoj,q0*nmol:q1*nmol] = lib.dot(v, mo_pairs.reshape(-1,ngrids).T)
        v = None
    return eri

    def call(self, args, axis=0, out=None, chunksize=1024 * 1024, **kwargs):
        """ axis is the axis to chop it off.
            if self.altreduce is set, the results will
            be reduced with altreduce and returned
            otherwise will be saved to out, then return out.
        """
        if self.altreduce is not None:
            ret = [None]
        else:
            if out is None :
                if self.outdtype is not None:
                    dtype = self.outdtype
                else:
                    try:
                        dtype = numpy.result_type(*[args[i] for i in self.ins] * 2)
                    except:
                        dtype = None
                out = sharedmem.empty(
                        numpy.broadcast(*[args[i] for i in self.ins] * 2).shape,
                        dtype=dtype)
        if axis != 0:
            for i in self.ins:
                args[i] = numpy.rollaxis(args[i], axis)
            out = numpy.rollaxis(out, axis)
        size = numpy.max([len(args[i]) for i in self.ins])
        with sharedmem.MapReduce() as pool:
            def work(i):
                sl = slice(i, i+chunksize)
                myargs = args[:]
                for j in self.ins:
                    try: 
                        tmp = myargs[j][sl]
                        a, b, c = sl.indices(len(args[j]))
                        myargs[j] = tmp
                    except Exception as e:
                        print tmp
                        print j, e
                        pass
                if b == a: return None
                rt = self.ufunc(*myargs, **kwargs)
                if self.altreduce is not None:
                    return rt
                else:
                    out[sl] = rt
            def reduce(rt):
                if self.altreduce is None:
                    return
                if ret[0] is None:
                    ret[0] = rt
                elif rt is not None:
                    ret[0] = self.altreduce(ret[0], rt)

            pool.map(work, range(0, size, chunksize), reduce=reduce)

        if self.altreduce is None:
            if axis != 0:
                out = numpy.rollaxis(out, 0, axis + 1)
            return out                
        else:
            return ret[0]

    def _inequality(self, other, op, op_name, bad_scalar_msg):
        # Scalar other.
        if isscalarlike(other):
            if 0 == other and op_name in ('_le_', '_ge_'):
                raise NotImplementedError(" >= and <= don't work with 0.")
            elif op(0, other):
                warn(bad_scalar_msg, SparseEfficiencyWarning)
                other_arr = np.empty(self.shape, dtype=np.result_type(other))
                other_arr.fill(other)
                other_arr = self.__class__(other_arr)
                return self._binopt(other_arr, op_name)
            else:
                return self._scalar_binopt(other, op)
        # Dense other.
        elif isdense(other):
            return op(self.todense(), other)
        # Sparse other.
        elif isspmatrix(other):
            #TODO sparse broadcasting
            if self.shape != other.shape:
                raise ValueError("inconsistent shapes")
            elif self.format != other.format:
                other = other.asformat(self.format)
            if op_name not in ('_ge_', '_le_'):
                return self._binopt(other, op_name)

            warn("Comparing sparse matrices using >= and <= is inefficient, "
                 "using <, >, or !=, instead.", SparseEfficiencyWarning)
            all_true = self.__class__(np.ones(self.shape))
            res = self._binopt(other, '_gt_' if op_name == '_le_' else '_lt_')
            return all_true - res
        else:
            raise ValueError("Operands could not be compared.")