Python numpy.block函数代码示例

    def test_nested(self):
        one = np.array([1, 1, 1])
        two = np.array([[2, 2, 2], [2, 2, 2], [2, 2, 2]])
        three = np.array([3, 3, 3])
        four = np.array([4, 4, 4])
        five = np.array(5)
        six = np.array([6, 6, 6, 6, 6])
        zero = np.zeros((2, 6))

        result = np.block([
            [
                np.block([
                   [one],
                   [three],
                   [four]
                ]),
                two
            ],
            [five, six],
            [zero]
        ])
        expected = np.array([[1, 1, 1, 2, 2, 2],
                             [3, 3, 3, 2, 2, 2],
                             [4, 4, 4, 2, 2, 2],
                             [5, 6, 6, 6, 6, 6],
                             [0, 0, 0, 0, 0, 0],
                             [0, 0, 0, 0, 0, 0]])

        assert_equal(result, expected)

        def update_basis_values(self, fnEvalBasis):
            if self.basis_values is None:
                self.basis_values = fnEvalBasis(self.basis, self.X)
                return self.basis_values

            prev_basis = self.basis_values.shape[1]
            prev_pts = self.basis_values.shape[0]
            if prev_basis == len(self.basis) and prev_pts == len(self.X):
                return self.basis_values # No change

            X = np.array(self.X)

            A = fnEvalBasis(self.basis, X[prev_pts:, :], 0,
                            prev_basis) if len(X) > prev_pts else None
            B = self.basis_values
            C = fnEvalBasis(self.basis, X, prev_basis) if len(self.basis) > prev_basis else None

            if A is None:
                self.basis_values = np.block([B, C])
            elif C is None:
                self.basis_values = np.block([[B], [A]])
            else:
                self.basis_values = np.block([[B, C[:prev_pts, :]],
                                              [A, C[prev_pts:, :]]])
            return self.basis_values

    def feedback(self, g2=None):
        """
        Calculates the `InternalDelay` object formed when combining two
        `InternalDelay` objects in a feedback loop in the following manner:

        ------>+ o -----> G1 ------->
                 ^-             |
                 |              |
                 |              |
                 |<---- G2 <----|

        where G1 is `self`, and if G2 is not given, it is assumed that
        G2 is an identity matrix.
        """
        if g2 is None:
            g2 = InternalDelay.from_tf_coefficients([1], [1], [0])

        X_inv = numpy.linalg.inv(numpy.eye(g2.D11.shape[0]) + g2.D11 @ self.D11)

        A = numpy.block([
            [self.A - self.B1 @ X_inv @ g2.D11 @ self.C1,
             -self.B1 @ X_inv @ g2.C1],
            [g2.B1 @ self.C1 - g2.B1 @ self.D11 @ X_inv @ g2.D11 @ self.C1,
             g2.A - g2.B1 @ self.D11 @ X_inv @ g2.C1]])

        B1 = numpy.block([[self.B1 - self.B1 @ X_inv @ g2.D11 @ self.D11],
                          [g2.B1 @ self.D11 - g2.B1 @ self.D11 @ X_inv @ g2.D11 @ self.D11]])

        B2 = numpy.block([
            [self.B2 - self.B1 @ X_inv @ g2.D11 @ self.D12,
             -self.B1 @ X_inv @ g2.D12],
            [g2.B1 @ self.D12 - g2.B1 @ self.D11 @ X_inv @ g2.D11 @ self.D12,
             g2.B2 - g2.B1 @ self.D11 @ X_inv @ g2.D12]])

        C1 = numpy.block([self.C1 - self.D11 @ X_inv @ g2.D11 @ self.C1,
                          -self.D11 @ X_inv @ g2.C1])

        C2 = numpy.block([
            [self.C2 - self.D21 @ X_inv @ g2.D11 @ self.C1,
             -self.D21 @ X_inv @ g2.C1],
            [g2.D21 @ self.C1 - g2.D21 @ self.D11 @ X_inv @ g2.D11 @ self.C1,
             g2.C2 - g2.D21 @ self.D11 @ X_inv @ g2.C1]])

        D11 = self.D11 - self.D11 @ X_inv @ g2.D11 @ self.D11

        D12 = numpy.block([self.D12 - self.D11 @ X_inv @ g2.D11 @ self.D12,
                           - self.D11 @ X_inv @ g2.D12])

        D21 = numpy.block([[self.D21 - self.D21 @ X_inv @ g2.D11 @ self.D11],
                           [g2.D21 @ self.D11 - g2.D21 @ self.D11 @ X_inv @ g2.D11 @ self.D11]])

        D22 = numpy.block([
            [self.D22 - self.D21 @ X_inv @ g2.D11 @ self.D12,
             -self.D21 @ X_inv @ g2.D12],
            [g2.D21 @ self.D12 - g2.D21 @ self.D11 @ X_inv @ g2.D11 @ self.D12,
             g2.D22 - g2.D21 @ self.D11 @ X_inv @ g2.D12]])

        delays = numpy.block([self.delays, g2.delays])

        return InternalDelay(A, B1, B2, C1, C2, D11, D12, D21, D22, delays)

def firlp_lowpass1(numtaps, deltap, deltas, cutoff, width, fs, tol=None):
    # Edges of the transition band, expressed as radians per sample.
    wp = np.pi*(cutoff - 0.5*width)/(0.5*fs)
    ws = np.pi*(cutoff + 0.5*width)/(0.5*fs)
    # Grid density.
    density = 16*numtaps/np.pi
    # Number of grid points in the pass band.
    numfreqs_pass = int(np.ceil(wp*density))
    # Number of grid points in the stop band.
    numfreqs_stop = int(np.ceil((np.pi - ws)*density))

    # Grid of frequencies in the pass band.
    wpgrid = np.linspace(0, wp, numfreqs_pass)
    # Remove the first; the inequality associated with this frequency
    # will be replaced by an equality constraint.
    wpgrid = wpgrid[1:]
    # Grid of frequencies in the pass band.
    wsgrid = np.linspace(ws, np.pi, numfreqs_stop)

    # wgrid is the combined array of frequencies.
    wgrid = np.concatenate((wpgrid, wsgrid))

    # The array of weights in the linear programming problem.
    weights = np.concatenate((np.full_like(wpgrid, fill_value=1/deltap),
                              np.full_like(wsgrid, fill_value=1/deltas)))
    # The array of desired frequency responses.
    desired = np.concatenate((np.ones_like(wpgrid),
                              np.zeros_like(wsgrid)))

    R = (numtaps - 1)//2
    C = np.cos(wgrid[:, np.newaxis] * np.arange(R+1))
    V = 1/weights[:, np.newaxis]

    A = np.block([[C, -V], [-C, -V]])
    b = np.block([[desired, -desired]]).T
    c = np.zeros(R+2)
    c[-1] = 1

    # The equality constraint corresponding to H(0) = 1.
    A_eq = np.ones((1, R+2))
    A_eq[:, -1] = 0
    b_eq = np.array([1])

    print("numfreqs_pass =", numfreqs_pass, "  numfreqs_stop =", numfreqs_stop)

    print("R =", R)
    print("c.shape =", c.shape)
    print("A.shape =", A.shape)
    print("b.shape =", b.shape)
    print("A_eq.shape =", A_eq.shape)
    print("b_eq.shape =", b_eq.shape)

    taps_lp = solve_linprog(c, A, b, A_eq=A_eq, b_eq=b_eq, tol=tol)
    return taps_lp

 def time_3d(self, n):
     np.block([
         [
             [self.a000, self.a001],
             [self.a010, self.a011],
         ],
         [
             [self.a100, self.a101],
             [self.a110, self.a111],
         ]
     ])

 def time_nested(self, n):
     np.block([
         [
             np.block([
                [self.one],
                [self.three],
                [self.four]
             ]),
             self.two
         ],
         [self.five, self.six],
         [self.zero]
     ])

def construct_matrix(kernel, verbose=False):
    X = intervals()
    J = abra()
    nondiag_quad = gauss_laguerre()
    collocation_points = nondiag_quad[0]
    if args.skip_generalized:
        diag_quad = [ (nondiag_quad[0], nondiag_quad[1], np.eye(args.degree)) for i in range(args.degree) ]
    else:
        diag_quad = generalized_gauss(collocation_points)
    if verbose:
        print('Mesh with {0}x{0} blocks:'.format(X.size-1))
        H = X[1:] - X[:-1]
        print('Interval lengths: min={}, max={}'.format(minv(H), maxv(H)))

    if verbose:
        print(' -- This is a moment just before constructing matrix A'); start_time = time.time()
    A = np.block([[ make_block(nondiag_quad, diag_quad, kernel, X, i, j)
        for j in range(X.size-1) ] for i in range(X.size-1) ])
    if verbose:
        print(' -- Matrix A is constructed. Elapsed time:', time.time() - start_time)
    if args.plot_matrix:
        splot(XX, A)
        plt.show()
    XX = np.hstack([ make_x(nondiag_quad, X, i) for i in range(X.size-1) ])
    return XX, A

    def laplace_hess(self, rate: Number) -> la.lnarray:
        """Hessian of Laplace transform of SNR memory curve.

        Parameters
        ----------
        rate : float, optional
            Parameter of Laplace transform, ``s``.

        Returns
        -------
        hess : la.lnarray (2n(n-1),2n(n-1))
            Hessian of ``snr_laplace`` at ``s`` with respect to parameters.
        """
        # (p,c), (eta,theta), (Z,Zs,ZQZs)
        rows, cols, mats = self._derivs(rate, True)

        # (n,n,n,n)
        hessww = _dbl_diagsub(_outer3(rows[0], mats[0], cols[1])
                              + _outer3(rows[0], mats[2], cols[0])
                              + _outer3(rows[1], mats[1], cols[0])).sum(0)
        # (2,n,n,n,n)
        hesswq = _dbl_diagsub(_outer3(rows[0], mats[0], cols[0])
                              + _trnsp4(_outer3(rows[0], mats[1], cols[0])))

        # (n(n-1),n(n-1))
        hesspp = tens2mat(hessww + hesswq.sum(0)/self.frac[0])*self.frac[0]**2
        hesspm = tens2mat(hessww - hesswq[0]/self.frac[1]
                          + hesswq[1]/self.frac[0]) * self.frac[0]*self.frac[1]
        hessmm = tens2mat(hessww - hesswq.sum(0)/self.frac[1])*self.frac[1]**2
        # (2n(n-1),2n(n-1))
        return - np.block([[hesspp, hesspm],
                           [hesspm.T, hessmm]])

def genRotationMatrix(alpha=0.0, beta=0.0, gamma=0.0):
    """
    Generate a rotation matrix based on 3 angles

    Implemented based on coordinate system of the rigid body. Imagine a plane
    being pointed at a cardinal direction (NS,EW) by amount `alpha`,
    being tilted up or down to a 'height' by amount `beta`, then rolled
    on its axis by amount `gamma`.

    In order to conserve rotations performed on XYZ across to those performed
    on UVW, the 6x6 rotation matrix is composed of an identical 3x3 rotation
    matrix in the top left and bottom right corners, with zeros everywhere else.

    Parameters
    ----------
    alpha : degree
        angle about Z axis (yaw), angle covers a span of (0,360)
    beta : degree
        angle about Y axis (pitch), angle covers a span of (-90,90)
    gamma : degree
        angle about X axis (roll), angle covers a span of (0,360)
    """
    sub_rotation_matrix = np.dot(r_z(alpha), np.dot(r_y(beta), r_x(gamma)))
    empty_block = np.zeros((3,3))

    rot_mat = np.block([[sub_rotation_matrix, empty_block],
                        [empty_block, sub_rotation_matrix]])

    return rot_mat

def diff(N):
	invD = np.eye(N)
	invD[0][0] = 2
	for i in range(N-2):
		invD[i][i+2] = -1
	upperright = np.dot(la.inv(invD),np.diag(np.arange(1,N+1)*2))
	Mat = np.block([ [np.zeros((N,1)), upperright], [ np.array([[0]]), np.zeros((1,N)) ] ] )
	return Mat

    def test_different_ndims(self):
        a = 1.
        b = 2 * np.ones((1, 2))
        c = 3 * np.ones((1, 1, 3))

        result = np.block([a, b, c])
        expected = np.array([[[1., 2., 2., 3., 3., 3.]]])

        assert_equal(result, expected)

def _build_cauchy_strain_op(bfg):
    dim = bfg.shape[2]
    if dim == 2:
        g1, g2 = bfg[..., 0:1, :], bfg[..., 1:2, :]
        zz = nm.zeros_like(g1)
        out = nm.block([[g1, zz],
                        [zz, g2],
                        [g2, g1]])

    else:
        g1, g2, g3 = bfg[..., 0:1, :], bfg[..., 1:2, :], bfg[..., 2:3, :]
        zz = nm.zeros_like(g1)
        out = nm.block([[g1, zz, zz],
                        [zz, g2, zz],
                        [zz, zz, g3],
                        [g2, g1, zz],
                        [g3, zz, g1],
                        [zz, g3, g2]])

    return out

    def test_3d(self):
        a000 = np.ones((2, 2, 2), int) * 1

        a100 = np.ones((3, 2, 2), int) * 2
        a010 = np.ones((2, 3, 2), int) * 3
        a001 = np.ones((2, 2, 3), int) * 4

        a011 = np.ones((2, 3, 3), int) * 5
        a101 = np.ones((3, 2, 3), int) * 6
        a110 = np.ones((3, 3, 2), int) * 7

        a111 = np.ones((3, 3, 3), int) * 8

        result = np.block([
            [
                [a000, a001],
                [a010, a011],
            ],
            [
                [a100, a101],
                [a110, a111],
            ]
        ])
        expected = array([[[1, 1, 4, 4, 4],
                           [1, 1, 4, 4, 4],
                           [3, 3, 5, 5, 5],
                           [3, 3, 5, 5, 5],
                           [3, 3, 5, 5, 5]],

                          [[1, 1, 4, 4, 4],
                           [1, 1, 4, 4, 4],
                           [3, 3, 5, 5, 5],
                           [3, 3, 5, 5, 5],
                           [3, 3, 5, 5, 5]],

                          [[2, 2, 6, 6, 6],
                           [2, 2, 6, 6, 6],
                           [7, 7, 8, 8, 8],
                           [7, 7, 8, 8, 8],
                           [7, 7, 8, 8, 8]],

                          [[2, 2, 6, 6, 6],
                           [2, 2, 6, 6, 6],
                           [7, 7, 8, 8, 8],
                           [7, 7, 8, 8, 8],
                           [7, 7, 8, 8, 8]],

                          [[2, 2, 6, 6, 6],
                           [2, 2, 6, 6, 6],
                           [7, 7, 8, 8, 8],
                           [7, 7, 8, 8, 8],
                           [7, 7, 8, 8, 8]]])

        assert_array_equal(result, expected)

    def parallel(self, g2):
        """
        Calculates the `InternalDelay` object formed when combining two
        `InternalDelay` objects in parallel in the following manner:

        -----> G1 ------->|
                          |
                          v+
                          o ---->
                          ^+
                          |
        -----> G2 ------->|

        where G1 is `self`.
        """

        A = numpy.block([[self.A, numpy.zeros((self.A.shape[0], g2.A.shape[1]))],
                         [numpy.zeros((g2.A.shape[0], self.A.shape[1])), g2.A]])

        B1 = numpy.block([[self.B1],
                          [g2.B1]])

        B2 = numpy.block([[self.B2, numpy.zeros((self.B2.shape[0], g2.B2.shape[1]))],
                          [numpy.zeros((g2.B2.shape[0], self.B2.shape[1])), g2.B2]])

        C1 = numpy.block([self.C1, g2.C1])

        C2 = numpy.block([[self.C2, numpy.zeros((self.C2.shape[0], g2.C2.shape[1]))],
                          [numpy.zeros((g2.C2.shape[0], self.C2.shape[1])), g2.C2]])

        D11 = self.D11 + g2.D11

        D12 = numpy.block([self.D12, g2.D12])

        D21 = numpy.block([[self.D21],
                           [g2.D21]])

        D22 = numpy.block([[self.D22, numpy.zeros((self.D22.shape[0], g2.D22.shape[1]))],
                           [numpy.zeros((g2.D22.shape[0], self.D22.shape[1])), g2.D22]])

        delays = numpy.block([self.delays, g2.delays])

        return InternalDelay(A, B1, B2, C1, C2, D11, D12, D21, D22, delays)

def _State_freq_resp(mA, mb, sc, f, dt=None):
    """
    This is the low level function to generate the frequency response
    values for a state space representation. The realization must be
    strictly in the observable Hessenberg form.

    Implements the inner loop of Misra, Patel SIMAX 1988 Algo. 3.1 in
    batches of B matrices instead of looping over every column of B.

    Parameters
    ----------

    mA : array_like {n x n}
        The A matrix of the realization in the upper Hessenberg form
    mb : array_like {n x m}
        The B vector of the realization
    sc : float
        The only nonzero coefficient of the o'ble-Hessenberg form
    f : array_like
        The frequency grid
    d : bool, optional
        Evaluate on the imaginary axis or unit circle. Default is imaginary
        axis.

    Returns
    -------
    r  : complex-valued numpy array

    """

    nn, m = mA.shape[0], mb.shape[1]
    r = empty((f.size, m), dtype=complex)
    Ab = block([-mA, mb]).astype(complex)
    U = empty_like(Ab)

    imag_indices = diag_indices(nn)
    # Triangularization of a Hessenberg matrix
    for ind, val in enumerate(f):
        U[:, :] = Ab  # Working copy
        U[imag_indices] += val*1j if dt is None else np.exp(dt*val*1j)
        for x in range(1, nn):
            U[x, x:] -= (U[x, x-1] / U[x-1, x-1]) * U[x-1, x:]

        r[ind, :] = U[-1, -m:] / U[-1, -1-m]

    return r*sc