Python numpy 模块,broadcast_arrays() 实例源码

我们从Python开源项目中,提取了以下32个代码示例,用于说明如何使用numpy.broadcast_arrays()

项目:uncover-ml    作者:GeoscienceAustralia    | 项目源码 | 文件源码
def __split_apply(self, func_y_lt_f, func_y_gt_f, func_f_lt_0, y, f):

        # Make sure all arrays are of a compatible type
        y, f = np.broadcast_arrays(y, f)
        y = y.astype(float, copy=False)
        f = f.astype(float, copy=False)

        if any(y < 0):
            raise ValueError("y has to be > 0!")

        # get indicators of which likelihoods to apply where
        y_lt_f = np.logical_and(y <= f, f > 0)
        y_gt_f = np.logical_and(y > f, f > 0)
        f_lt_0 = f <= 0

        result = np.zeros_like(y)

        result[y_lt_f] = func_y_lt_f(y[y_lt_f], f[y_lt_f])
        result[y_gt_f] = func_y_gt_f(y[y_gt_f], f[y_gt_f])
        result[f_lt_0] = func_f_lt_0(y[f_lt_0], f[f_lt_0])

        return result
项目:catalyst    作者:enigmampc    | 项目源码 | 文件源码
def compute(self, today, assets, out, dependent, independent):
        alpha = out.alpha
        beta = out.beta
        r_value = out.r_value
        p_value = out.p_value
        stderr = out.stderr

        def regress(y, x):
            regr_results = linregress(y=y, x=x)
            # `linregress` returns its results in the following order:
            # slope, intercept, r-value, p-value, stderr
            alpha[i] = regr_results[1]
            beta[i] = regr_results[0]
            r_value[i] = regr_results[2]
            p_value[i] = regr_results[3]
            stderr[i] = regr_results[4]

        # If `independent` is a Slice or single column of data, broadcast it
        # out to the same shape as `dependent`, then compute column-wise. This
        # is efficient because each column of the broadcasted array only refers
        # to a single memory location.
        independent = broadcast_arrays(independent, dependent)[0]
        for i in range(len(out)):
            regress(y=dependent[:, i], x=independent[:, i])
项目:sparse    作者:mrocklin    | 项目源码 | 文件源码
def _cartesian_product(*arrays):
        """
        Get the cartesian product of a number of arrays.

        Parameters
        ----------
        arrays : Iterable[np.ndarray]
            The arrays to get a cartesian product of. Always sorted with respect
            to the original array.
        Returns
        -------
        out : np.ndarray
            The overall cartesian product of all the input arrays.
        """
        broadcastable = np.ix_(*arrays)
        broadcasted = np.broadcast_arrays(*broadcastable)
        rows, cols = np.prod(broadcasted[0].shape), len(broadcasted)
        dtype = np.result_type(*arrays)
        out = np.empty(rows * cols, dtype=dtype)
        start, end = 0, rows
        for a in broadcasted:
            out[start:end] = a.reshape(-1)
            start, end = end, end + rows
        return out.reshape(cols, rows)
项目:radar    作者:amoose136    | 项目源码 | 文件源码
def test_broadcast_arrays():
    # Test user defined dtypes
    a = np.array([(1, 2, 3)], dtype='u4,u4,u4')
    b = np.array([(1, 2, 3), (4, 5, 6), (7, 8, 9)], dtype='u4,u4,u4')
    result = np.broadcast_arrays(a, b)
    assert_equal(result[0], np.array([(1, 2, 3), (1, 2, 3), (1, 2, 3)], dtype='u4,u4,u4'))
    assert_equal(result[1], np.array([(1, 2, 3), (4, 5, 6), (7, 8, 9)], dtype='u4,u4,u4'))
项目:uncover-ml    作者:GeoscienceAustralia    | 项目源码 | 文件源码
def Ey(self, f, var, z):
        f, z = np.broadcast_arrays(f, z)
        Ey = np.zeros_like(f)
        nz = ~ z

        Ey[z] = self.gaus.Ey(f[z], var)
        Ey[nz] = self.unif.Ey(f[nz])
        return Ey
项目:uncover-ml    作者:GeoscienceAustralia    | 项目源码 | 文件源码
def dp(self, y, f, var, z):
        yz, fz = np.broadcast_arrays(y[z], f[:, z])
        dp = np.zeros_like(f)
        dp[:, z] = self.gaus.dp(yz, fz, var)
        return dp
项目:uncover-ml    作者:GeoscienceAustralia    | 项目源码 | 文件源码
def __split_on_z(self, func_unif, func_gaus, y, f, var, z):

        y, f, z = np.broadcast_arrays(y, f, z)
        val = np.zeros_like(f)
        nz = ~ z
        val[z] = func_gaus(y[z], f[z], var)
        val[nz] = func_unif(y[nz], f[nz])

        return val
项目:catalyst    作者:enigmampc    | 项目源码 | 文件源码
def compute(self, today, assets, out, base_data, target_data):
        # If `target_data` is a Slice or single column of data, broadcast it
        # out to the same shape as `base_data`, then compute column-wise. This
        # is efficient because each column of the broadcasted array only refers
        # to a single memory location.
        target_data = broadcast_arrays(target_data, base_data)[0]
        for i in range(len(out)):
            out[i] = pearsonr(base_data[:, i], target_data[:, i])[0]
项目:catalyst    作者:enigmampc    | 项目源码 | 文件源码
def compute(self, today, assets, out, base_data, target_data):
        # If `target_data` is a Slice or single column of data, broadcast it
        # out to the same shape as `base_data`, then compute column-wise. This
        # is efficient because each column of the broadcasted array only refers
        # to a single memory location.
        target_data = broadcast_arrays(target_data, base_data)[0]
        for i in range(len(out)):
            out[i] = spearmanr(base_data[:, i], target_data[:, i])[0]
项目:svm-street-detector    作者:morris-frank    | 项目源码 | 文件源码
def world2image_uvMat(self, uv_mat):
        '''

        @param XYZ_mat: is a 4 or 3 times n matrix
        '''
        if uv_mat.shape[0] == 2:
            if len(uv_mat.shape)==1:
                uv_mat = uv_mat.reshape(uv_mat.shape + (1,))
            uv_mat = np.vstack((uv_mat, np.ones((1, uv_mat.shape[1]), uv_mat.dtype)))
        result = np.dot(self.Tr33, uv_mat)
        #w0 = -(uv_mat[0]* self.Tr_inv_33[1,0]+ uv_mat[1]* self.Tr_inv[1,1])/self.Tr_inv[1,2]
        resultB = np.broadcast_arrays(result, result[-1,:])
        return resultB[0] / resultB[1]
项目:svm-street-detector    作者:morris-frank    | 项目源码 | 文件源码
def image2world_uvMat(self, uv_mat):
        '''

        @param XYZ_mat: is a 4 or 3 times n matrix
        '''
        #assert self.transformMode == 'kitti' or self.transformMode == 'angles'
        if uv_mat.shape[0] == 2:
            if len(uv_mat.shape)==1:
                uv_mat = uv_mat.reshape(uv_mat.shape + (1,))
            uv_mat = np.vstack((uv_mat, np.ones((1, uv_mat.shape[1]), uv_mat.dtype)))
        result = np.dot(self.Tr33.I, uv_mat)
        #w0 = -(uv_mat[0]* self.Tr_inv_33[1,0]+ uv_mat[1]* self.Tr_inv[1,1])/self.Tr_inv[1,2]
        resultB = np.broadcast_arrays(result, result[-1,:])
        return resultB[0] / resultB[1]
项目:krpcScripts    作者:jwvanderbeck    | 项目源码 | 文件源码
def test_broadcast_arrays():
    # Test user defined dtypes
    a = np.array([(1, 2, 3)], dtype='u4,u4,u4')
    b = np.array([(1, 2, 3), (4, 5, 6), (7, 8, 9)], dtype='u4,u4,u4')
    result = np.broadcast_arrays(a, b)
    assert_equal(result[0], np.array([(1, 2, 3), (1, 2, 3), (1, 2, 3)], dtype='u4,u4,u4'))
    assert_equal(result[1], np.array([(1, 2, 3), (4, 5, 6), (7, 8, 9)], dtype='u4,u4,u4'))
项目:BAG_framework    作者:ucb-art    | 项目源码 | 文件源码
def __call__(self, t, tau, debug=False):
        """Calculate h(t, tau).

        Compute h(t, tau), which is the output at t subject to an impulse
        at time tau. standard numpy broadcasting rules apply.

        Parameters
        ----------
        t : array-like
            the output time.
        tau : array-like
            the input impulse time.
        debug : bool
            True to print debug messages.

        Returns
        -------
        val : :class:`numpy.ndarray`
            the time-varying impulse response evaluated at the given coordinates.
        """
        # broadcast arguments to same shape
        t, tau = np.broadcast_arrays(t, tau)

        # compute impulse using efficient matrix multiply and numpy broadcasting.
        dt = t - tau
        zero_indices = (dt < 0) | (dt > self.k * self.tper)
        t_row = t.reshape((1, -1))
        dt_row = dt.reshape((1, -1))
        tmp = np.dot(self.hmat, np.exp(np.dot(self.n_col, dt_row))) * np.exp(np.dot(self.m_col, t_row))
        result = np.sum(tmp, axis=0).reshape(dt.shape)

        # zero element such that dt < 0 or dt > k * T.
        result[zero_indices] = 0.0

        if debug:
            self._print_debug_msg(result)

        # discard imaginary part
        return np.real(result)
项目:CNN_Mnist    作者:AITTSMD    | 项目源码 | 文件源码
def max_pool_backward_reshape(dout, cache):
  """
  A fast implementation of the backward pass for the max pooling layer that
  uses some clever broadcasting and reshaping.

  This can only be used if the forward pass was computed using
  max_pool_forward_reshape.

  NOTE: If there are multiple argmaxes, this method will assign gradient to
  ALL argmax elements of the input rather than picking one. In this case the
  gradient will actually be incorrect. However this is unlikely to occur in
  practice, so it shouldn't matter much. One possible solution is to split the
  upstream gradient equally among all argmax elements; this should result in a
  valid subgradient. You can make this happen by uncommenting the line below;
  however this results in a significant performance penalty (about 40% slower)
  and is unlikely to matter in practice so we don't do it.
  """
  x, x_reshaped, out = cache

  dx_reshaped = np.zeros_like(x_reshaped)
  out_newaxis = out[:, :, :, np.newaxis, :, np.newaxis]
  mask = (x_reshaped == out_newaxis)
  dout_newaxis = dout[:, :, :, np.newaxis, :, np.newaxis]
  dout_broadcast, _ = np.broadcast_arrays(dout_newaxis, dx_reshaped)
  dx_reshaped[mask] = dout_broadcast[mask]
  dx_reshaped /= np.sum(mask, axis=(3, 5), keepdims=True)
  dx = dx_reshaped.reshape(x.shape)

  return dx
项目:PyDataLondon29-EmbarrassinglyParallelDAWithAWSLambda    作者:SignalMedia    | 项目源码 | 文件源码
def test_broadcast_arrays():
    # Test user defined dtypes
    a = np.array([(1, 2, 3)], dtype='u4,u4,u4')
    b = np.array([(1, 2, 3), (4, 5, 6), (7, 8, 9)], dtype='u4,u4,u4')
    result = np.broadcast_arrays(a, b)
    assert_equal(result[0], np.array([(1, 2, 3), (1, 2, 3), (1, 2, 3)], dtype='u4,u4,u4'))
    assert_equal(result[1], np.array([(1, 2, 3), (4, 5, 6), (7, 8, 9)], dtype='u4,u4,u4'))
项目:aRMSD    作者:armsd    | 项目源码 | 文件源码
def cartesian_product(arrays):
    """ Returns Cartesian product of given arrays (x and y): cartesian_product([x,y]) """

    broadcastable = np.ix_(*arrays)
    broadcasted = np.broadcast_arrays(*broadcastable)
    rows, cols = reduce(np.multiply, broadcasted[0].shape), len(broadcasted)
    out = np.empty(rows * cols, dtype=broadcasted[0].dtype)
    start, end = 0, rows

    for a in broadcasted:
        out[start:end] = a.reshape(-1)
        start, end = end, end + rows

    # Return value(s)
    return out.reshape(cols, rows).T
项目:aws-lambda-numpy    作者:vitolimandibhrata    | 项目源码 | 文件源码
def test_broadcast_arrays():
    # Test user defined dtypes
    a = np.array([(1, 2, 3)], dtype='u4,u4,u4')
    b = np.array([(1, 2, 3), (4, 5, 6), (7, 8, 9)], dtype='u4,u4,u4')
    result = np.broadcast_arrays(a, b)
    assert_equal(result[0], np.array([(1, 2, 3), (1, 2, 3), (1, 2, 3)], dtype='u4,u4,u4'))
    assert_equal(result[1], np.array([(1, 2, 3), (4, 5, 6), (7, 8, 9)], dtype='u4,u4,u4'))
项目:AutoEncoder    作者:situgongyuan    | 项目源码 | 文件源码
def max_pool_backward_reshape(dout, cache):
  """
  A fast implementation of the backward pass for the max pooling layer that
  uses some clever broadcasting and reshaping.

  This can only be used if the forward pass was computed using
  max_pool_forward_reshape.

  NOTE: If there are multiple argmaxes, this method will assign gradient to
  ALL argmax elements of the input rather than picking one. In this case the
  gradient will actually be incorrect. However this is unlikely to occur in
  practice, so it shouldn't matter much. One possible solution is to split the
  upstream gradient equally among all argmax elements; this should result in a
  valid subgradient. You can make this happen by uncommenting the line below;
  however this results in a significant performance penalty (about 40% slower)
  and is unlikely to matter in practice so we don't do it.
  """
  x, x_reshaped, out = cache
  dx_reshaped = np.zeros_like(x_reshaped)
  out_newaxis = out[:, :, :, np.newaxis, :, np.newaxis]
  mask = (x_reshaped == out_newaxis)
  dout_newaxis = dout[:, :, :, np.newaxis, :, np.newaxis]
  dout_broadcast, _ = np.broadcast_arrays(dout_newaxis, dx_reshaped)
  dx_reshaped[mask] = dout_broadcast[mask]
  dx_reshaped /= np.sum(mask, axis=(3, 5), keepdims=True)
  dx = dx_reshaped.reshape(x.shape)

  return dx
项目:optimize-stencil    作者:Ablinne    | 项目源码 | 文件源码
def omega_output(self, fname, parameters):
        """This method writes the dispersion relation to the file fname.

        :param fname: Filename
        :type fname: string
        :param parameters: Array containing the parameters.
        :type parameters: numpy.ndarray
        """

        omega = self.omega(parameters)
        kappa = self.kappamesh
        kmesh = self.kmesh
        k = self.k
        #print(omega, self.dks)
        vgs = np.gradient(omega, *self.dks, edge_order=2)
        vg = np.sqrt(sum(vgc**2 for vgc in vgs))
        if self.dim==2:
            vg[:3,:3] = 1
        elif self.dim==3:
            vg[:3,:3,:3] = 1
        vph = omega/k
        vph.ravel()[0] = 1.
        data = [*np.broadcast_arrays(*kappa), k, omega, vph, vg, *vgs]

        data = np.broadcast_arrays(*data)
        blocks = zip(*map(lambda x: np.vsplit(x, x.shape[0]), data))
        with open(fname, 'wb') as f:
            kappanames = ' '.join(['kx', 'ky', 'kz'][:len(kappa)])
            vgsnames = ' '.join(['vgx', 'vgy', 'vgz'][:len(vgs)])
            f.write(('#' + kappanames + ' k omega vph vg ' + vgsnames + '\n').encode('us-ascii'))

            for block_columns in blocks:
                np.savetxt(f, np.vstack(map(np.ravel, block_columns)).T)
                f.write(b'\n')
项目:lambda-numba    作者:rlhotovy    | 项目源码 | 文件源码
def test_broadcast_arrays():
    # Test user defined dtypes
    a = np.array([(1, 2, 3)], dtype='u4,u4,u4')
    b = np.array([(1, 2, 3), (4, 5, 6), (7, 8, 9)], dtype='u4,u4,u4')
    result = np.broadcast_arrays(a, b)
    assert_equal(result[0], np.array([(1, 2, 3), (1, 2, 3), (1, 2, 3)], dtype='u4,u4,u4'))
    assert_equal(result[1], np.array([(1, 2, 3), (4, 5, 6), (7, 8, 9)], dtype='u4,u4,u4'))
项目:deliver    作者:orchestor    | 项目源码 | 文件源码
def test_broadcast_arrays():
    # Test user defined dtypes
    a = np.array([(1, 2, 3)], dtype='u4,u4,u4')
    b = np.array([(1, 2, 3), (4, 5, 6), (7, 8, 9)], dtype='u4,u4,u4')
    result = np.broadcast_arrays(a, b)
    assert_equal(result[0], np.array([(1, 2, 3), (1, 2, 3), (1, 2, 3)], dtype='u4,u4,u4'))
    assert_equal(result[1], np.array([(1, 2, 3), (4, 5, 6), (7, 8, 9)], dtype='u4,u4,u4'))
项目:Theano-Deep-learning    作者:GeekLiB    | 项目源码 | 文件源码
def _numpy_second(x, y):
    return numpy.broadcast_arrays(x, y)[1]
项目:wavenet    作者:rampage644    | 项目源码 | 文件源码
def __init__(self, *args, mask='B', **kwargs):
        super(MaskedConvolution2D, self).__init__(
            *args, **kwargs
        )

        Cout, Cin, kh, kw = self.W.shape
        pre_mask = self.xp.ones_like(self.W.data).astype('f')
        yc, xc = kh // 2, kw // 2

        # context masking - subsequent pixels won't hav access to next pixels (spatial dim)
        pre_mask[:, :, yc+1:, :] = 0.0
        pre_mask[:, :, yc:, xc+1:] = 0.0

        # same pixel masking - pixel won't access next color (conv filter dim)
        def bmask(i_out, i_in):
            cout_idx = np.expand_dims(np.arange(Cout) % 3 == i_out, 1)
            cin_idx = np.expand_dims(np.arange(Cin) % 3 == i_in, 0)
            a1, a2 = np.broadcast_arrays(cout_idx, cin_idx)
            return a1 * a2

        for j in range(3):
            pre_mask[bmask(j, j), yc, xc] = 0.0 if mask == 'A' else 1.0

        pre_mask[bmask(0, 1), yc, xc] = 0.0
        pre_mask[bmask(0, 2), yc, xc] = 0.0
        pre_mask[bmask(1, 2), yc, xc] = 0.0

        self.mask = pre_mask
项目:wavenet    作者:rampage644    | 项目源码 | 文件源码
def bmask(i_out, i_in):
    cout_idx = np.expand_dims(np.arange(Cout) % 3 == i_out, 1)
    cin_idx = np.expand_dims(np.arange(Cin) % 3 == i_in, 0)
    a1, a2 = np.broadcast_arrays(cout_idx, cin_idx)
    return a1 * a2
项目:DLCourse_AS    作者:TomoyaKung    | 项目源码 | 文件源码
def max_pool_backward_reshape(dout, cache):
  """
  A fast implementation of the backward pass for the max pooling layer that
  uses some clever broadcasting and reshaping.

  This can only be used if the forward pass was computed using
  max_pool_forward_reshape.

  NOTE: If there are multiple argmaxes, this method will assign gradient to
  ALL argmax elements of the input rather than picking one. In this case the
  gradient will actually be incorrect. However this is unlikely to occur in
  practice, so it shouldn't matter much. One possible solution is to split the
  upstream gradient equally among all argmax elements; this should result in a
  valid subgradient. You can make this happen by uncommenting the line below;
  however this results in a significant performance penalty (about 40% slower)
  and is unlikely to matter in practice so we don't do it.
  """
  x, x_reshaped, out = cache

  dx_reshaped = np.zeros_like(x_reshaped)
  out_newaxis = out[:, :, :, np.newaxis, :, np.newaxis]
  mask = (x_reshaped == out_newaxis)
  dout_newaxis = dout[:, :, :, np.newaxis, :, np.newaxis]
  dout_broadcast, _ = np.broadcast_arrays(dout_newaxis, dx_reshaped)
  dx_reshaped[mask] = dout_broadcast[mask]
  dx_reshaped /= np.sum(mask, axis=(3, 5), keepdims=True)
  dx = dx_reshaped.reshape(x.shape)

  return dx
项目:DLCourse_AS    作者:TomoyaKung    | 项目源码 | 文件源码
def max_pool_backward_reshape(dout, cache):
  """
  A fast implementation of the backward pass for the max pooling layer that
  uses some clever broadcasting and reshaping.

  This can only be used if the forward pass was computed using
  max_pool_forward_reshape.

  NOTE: If there are multiple argmaxes, this method will assign gradient to
  ALL argmax elements of the input rather than picking one. In this case the
  gradient will actually be incorrect. However this is unlikely to occur in
  practice, so it shouldn't matter much. One possible solution is to split the
  upstream gradient equally among all argmax elements; this should result in a
  valid subgradient. You can make this happen by uncommenting the line below;
  however this results in a significant performance penalty (about 40% slower)
  and is unlikely to matter in practice so we don't do it.
  """
  x, x_reshaped, out = cache

  dx_reshaped = np.zeros_like(x_reshaped)
  out_newaxis = out[:, :, :, np.newaxis, :, np.newaxis]
  mask = (x_reshaped == out_newaxis)
  dout_newaxis = dout[:, :, :, np.newaxis, :, np.newaxis]
  dout_broadcast, _ = np.broadcast_arrays(dout_newaxis, dx_reshaped)
  dx_reshaped[mask] = dout_broadcast[mask]
  dx_reshaped /= np.sum(mask, axis=(3, 5), keepdims=True)
  dx = dx_reshaped.reshape(x.shape)

  return dx
项目:stream2segment    作者:rizac    | 项目源码 | 文件源码
def min_traveltimes(modelname, source_depths, receiver_depths, distances, phases, callback=None):
    '''computes the minimum travel times using multiprocessing to speed up
    calculations
    min_traveltimes(modelname, ARRAY[N], ARRAY[N], SCALAR, ...) -> [N-length array]
    min_traveltimes(modelname, ARRAY[N], ARRAY[N], ARRAY[M]) -> [NxM] matrix
    min_traveltimes(modelname, SCALAR, SCALAR, ARRAY[M]) -> [M-length array]
    '''
    model = taumodel(modelname)
    source_depths, receiver_depths = np.broadcast_arrays(source_depths, receiver_depths)
    # assert we passed arrays, not matrices:
    if len(source_depths.shape) > 1 or len(distances.shape) > 1:
        raise ValueError("Need to have arrays, not matrices")
    norowdim = source_depths.ndim == 0
    if norowdim:
        source_depths = np.array([source_depths])
        receiver_depths = np.array([receiver_depths])
    nocoldim = distances.ndim == 0
    if nocoldim:
        distances = np.array([distances])
    ttimes = np.full(shape=(source_depths.shape[0], distances.shape[0]), fill_value=np.nan)

    def mp_callback(index, array, _callback=None):
        def _(result):
            array[index] = result
            if _callback is not None:
                _callback()
        return _

    pool = Pool()
    for idx, sd, rd in zip(count(), source_depths, receiver_depths):
        tmp_ttimes = ttimes[idx]
        for i, d in enumerate(distances):
            pool.apply_async(min_traveltime, (model, sd, rd, d, phases),
                             callback=mp_callback(i, tmp_ttimes, callback))
    pool.close()
    pool.join()

    if norowdim or nocoldim:
        ttimes = ttimes.flatten()

    return ttimes
项目:Alfred    作者:jkachhadia    | 项目源码 | 文件源码
def test_broadcast_arrays():
    # Test user defined dtypes
    a = np.array([(1, 2, 3)], dtype='u4,u4,u4')
    b = np.array([(1, 2, 3), (4, 5, 6), (7, 8, 9)], dtype='u4,u4,u4')
    result = np.broadcast_arrays(a, b)
    assert_equal(result[0], np.array([(1, 2, 3), (1, 2, 3), (1, 2, 3)], dtype='u4,u4,u4'))
    assert_equal(result[1], np.array([(1, 2, 3), (4, 5, 6), (7, 8, 9)], dtype='u4,u4,u4'))
项目:lombscargle    作者:jakevdp    | 项目源码 | 文件源码
def periodic_fit(t, y, dy, frequency, t_fit,
                 center_data=True, fit_bias=True, nterms=1):
    """Compute the Lomb-Scargle model fit at a given frequency

    Parameters
    ----------
    t, y, dy : float or array_like
        The times, observations, and uncertainties to fit
    frequency : float
        The frequency at which to compute the model
    t_fit : float or array_like
        The times at which the fit should be computed
    center_data : bool (default=True)
        If True, center the input data before applying the fit
    fit_bias : bool (default=True)
        If True, include the bias as part of the model
    nterms : int (default=1)
        The number of Fourier terms to include in the fit

    Returns
    -------
    y_fit : ndarray
        The model fit evaluated at each value of t_fit
    """
    if dy is None:
        dy = 1

    t, y, dy = np.broadcast_arrays(t, y, dy)
    t_fit = np.asarray(t_fit)
    assert t.ndim == 1
    assert t_fit.ndim == 1
    assert np.isscalar(frequency)

    if center_data:
        w = dy ** -2.0
        y_mean = np.dot(y, w) / w.sum()
        y = (y - y_mean)
    else:
        y_mean = 0

    X = design_matrix(t, frequency, dy=dy, bias=fit_bias, nterms=nterms)
    theta_MLE = np.linalg.solve(np.dot(X.T, X),
                                np.dot(X.T, y / dy))

    X_fit = design_matrix(t_fit, frequency, bias=fit_bias, nterms=nterms)

    return y_mean + np.dot(X_fit, theta_MLE)
项目:lombscargle    作者:jakevdp    | 项目源码 | 文件源码
def lombscargle_scipy(t, y, frequency, normalization='normalized',
                      center_data=True):
    """Lomb-Scargle Periodogram

    This is a wrapper of ``scipy.signal.lombscargle`` for computation of the
    Lomb-Scargle periodogram. This is a relatively fast version of the naive
    O[N^2] algorithm, but cannot handle heteroskedastic errors.

    Parameters
    ----------
    t, y: array_like
        times, values, and errors of the data points. These should be
        broadcastable to the same shape.
    frequency : array_like
        frequencies (not angular frequencies) at which to calculate periodogram
    normalization : string (optional, default='normalized')
        Normalization to use for the periodogram
        TODO: figure out what options to use
    center_data : bool (optional, default=True)
        if True, pre-center the data by subtracting the weighted mean
        of the input data.

    Returns
    -------
    power : array_like
        Lomb-Scargle power associated with each frequency.
        Units of the result depend on the normalization.

    References
    ----------
    .. [1] M. Zechmeister and M. Kurster, A&A 496, 577-584 (2009)
    .. [2] W. Press et al, Numerical Recipies in C (2002)
    .. [3] Scargle, J.D. 1982, ApJ 263:835-853
    """
    if not HAS_SCIPY:
        raise ValueError("scipy must be installed to use lombscargle_scipy")

    t, y = np.broadcast_arrays(t, y)
    frequency = np.asarray(frequency)
    assert t.ndim == 1
    assert frequency.ndim == 1

    if center_data:
        y = y - y.mean()

    # Note: scipy input accepts angular frequencies
    p = signal.lombscargle(t, y, 2 * np.pi * frequency)

    if normalization == 'unnormalized':
        pass
    elif normalization == 'normalized':
        p *= 2 / (t.size * np.mean(y ** 2))
    else:
        raise ValueError("normalization='{0}' "
                         "not recognized".format(normalization))
    return p
项目:DLTK    作者:DLTK    | 项目源码 | 文件源码
def elastic_transform(image, alpha, sigma):
    """
    Elastic deformation of images as described in [1].

    [1] Simard, Steinkraus and Platt, "Best Practices for Convolutional
        Neural Networks applied to Visual Document Analysis", in Proc. of the
        International Conference on Document Analysis and Recognition, 2003.

    Based on gist https://gist.github.com/erniejunior/601cdf56d2b424757de5

    Args:
        image (np.ndarray): image to be deformed
        alpha (list): scale of transformation for each dimension, where larger
            values have more deformation
        sigma (list): Gaussian window of deformation for each dimension, where
            smaller values have more localised deformation

    Returns:
        np.ndarray: deformed image
    """

    assert len(alpha) == len(sigma), \
        "Dimensions of alpha and sigma are different"

    channelbool = image.ndim - len(alpha)
    out = np.zeros((len(alpha) + channelbool, ) + image.shape)

    # Generate a Gaussian filter, leaving channel dimensions zeroes
    for jj in range(len(alpha)):
        array = (np.random.rand(*image.shape) * 2 - 1)
        out[jj] = gaussian_filter(array, sigma[jj],
                                  mode="constant", cval=0) * alpha[jj]

    # Map mask to indices
    shapes = list(map(lambda x: slice(0, x, None), image.shape))
    grid = np.broadcast_arrays(*np.ogrid[shapes])
    indices = list(map((lambda x: np.reshape(x, (-1, 1))), grid + np.array(out)))

    # Transform image based on masked indices
    transformed_image = map_coordinates(image, indices, order=0,
                                        mode='reflect').reshape(image.shape)

    return transformed_image
项目:stream2segment    作者:rizac    | 项目源码 | 文件源码
def min(self, source_depths, receiver_depths, distances, method='linear'):
        '''
        Returns the minimum (minima) travel time(s) for the point(s) identified by
        (source_depths, receiver_depths, distances) by building a grid and
        interpolating on the points identified by each
P[i] = (source_depths[i], receiver_depths[i], distances[i])
    ```
    if the source file has been built with
    receiver depths == [0]. It uses a 2d linear interpolation on a grid. It uses
    scipy griddata
    :param source_depths: numeric or numpy array of length n: the source depth(s), in km
    :param receiver_depths: numeric or numpy array of length n: the receiver depth(s), in km.
    For most applications, this value can be set to zero
    :param distances: numeric or numpy array of length n: the distance(s), in degrees
    :param method: forwarded to `scipy.griddata` function
    :return: a numpy array of length n denoting the minimum travel time(s) of this model for
    each P
    '''
    # Handle the case some arguments are scalars and some arrays:
    source_depths, receiver_depths, distances = \
        np.broadcast_arrays(source_depths, receiver_depths, distances)
    # handle the case all arguments scalars. See
    # https://stackoverflow.com/questions/29318459/python-function-that-handles-scalar-or-arrays
    allscalars = all(_.ndim == 0 for _ in (source_depths, receiver_depths, distances))
    if source_depths.ndim == 0:
        source_depths = source_depths[None]  # Makes x 1D
    if receiver_depths.ndim == 0:
        receiver_depths = receiver_depths[None]  # Makes x 1D
    if distances.ndim == 0:
        distances = distances[None]  # Makes x 1D
    # correct source depths and receiver depths
    source_depths[source_depths < 0] = 0
    receiver_depths[receiver_depths < 0] = 0
    # correct distances to be compatible with obpsy traveltimes calculations:
    distances = distances % 360
    # set values symmetric to 180 degrees if greater than 180:
    _mask = distances > 180
    if _mask.any():  # does this speeds up (allocate mask array once)? FIXME: check
        distances[_mask] = 360 - distances[_mask]
    # set source depths to nan if out of bounds. This prevent method = 'nearest'
    # to return non nan values and be consistent with 'linear' and 'cubic'
    if self._unique_receiver_depth:
        # get values without receiver depth dimension:
        values = np.hstack((self._km2deg(source_depths).reshape(-1, 1),
                            distances.reshape(-1, 1)))
    else:
        values = np.hstack((self._km2deg(source_depths).reshape(-1, 1),
                            self._km2deg(receiver_depths).reshape(-1, 1),
                            distances.reshape(-1, 1)))

    ret = griddata(self._gridpts, self._gridvals, values,
                   method=method, rescale=False, fill_value=np.nan)

    # ret is almost likely a float, so we can set now nans for out of boun values
    # we cannot do it before on any input array cause int arrays do not support nan assignement
    ret[(source_depths > self._sourcedepth_bounds_km[1]) |
        (receiver_depths > self._receiverdepth_bounds_km[1])] = np.nan
    # return scalar if inputs are scalar, array oitherwise
    return np.squeeze(ret) if allscalars else ret

```