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

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

        -----> G2 -----> G1 ------->

        where G1 is `self`.
        """
        A = numpy.block([[self.A, numpy.zeros((self.A.shape[0], g2.A.shape[1]))],
                         [g2.B1 @ self.C1, g2.A]])

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

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

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

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

        D11 = g2.D11 @ self.D11

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

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

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

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

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

def test_create_block_3_diagonal_matrix():
    np.random.seed(0)
    A = np.empty((4, 3, 3))
    A[:] = np.arange(1, 5)[:, None, None]

    B = np.empty((4, 3, 3))
    B[:] = -np.arange(1, 5)[:, None, None]
    d = 10 * np.arange(10, 15)

    banded = _create_block_3_diagonal_matrix(A, B, d)

    # Convert the banded matrix to the full matrix.
    k, l = list(zip(*product(np.arange(banded.shape[0]),
                             np.arange(banded.shape[1]))))
    k = np.asarray(k)
    l = np.asarray(l)

    i = k - 5 + l
    j = l
    values = banded.ravel()
    mask = (i >= 0) & (i < 15)
    i = i[mask]
    j = j[mask]
    values = values[mask]
    full = np.zeros((15, 15))
    full[i, j] = values

    zero = np.zeros((3, 3))
    eye = np.eye(3)

    # Create the reference full matrix in the most straightforward manner.
    ref = np.block([
        [d[0] * eye, B[0], zero, zero, zero],
        [A[0], d[1] * eye, B[1], zero, zero],
        [zero, A[1], d[2] * eye, B[2], zero],
        [zero, zero, A[2], d[3] * eye, B[3]],
        [zero, zero, zero, A[3], d[4] * eye],
    ])

    assert_allclose(full, ref, atol=1e-19)

 def _compute_jacobian(self, J_eq, J_ineq, s):
     if self.n_ineq == 0:
         return J_eq
     else:
         if sps.issparse(J_eq) or sps.issparse(J_ineq):
             # It is expected that J_eq and J_ineq
             # are already `csr_matrix` because of
             # the way ``BoxConstraint``, ``NonlinearConstraint``
             # and ``LinearConstraint`` are defined.
             J_eq = sps.csr_matrix(J_eq)
             J_ineq = sps.csr_matrix(J_ineq)
             return self._assemble_sparse_jacobian(J_eq, J_ineq, s)
         else:
             S = np.diag(s)
             zeros = np.zeros((self.n_eq, self.n_ineq))
             # Convert to matrix
             if sps.issparse(J_ineq):
                 J_ineq = J_ineq.toarray()
             if sps.issparse(J_eq):
                 J_eq = J_eq.toarray()
             # Concatenate matrices
             return np.block([[J_eq, zeros],
                              [J_ineq, S]])

    def lft(self, other, nu=-1, ny=-1):
        """Return the Linear Fractional Transformation.

        A definition of the LFT operator can be found in Appendix A.7,
        page 512 in the 2nd Edition, Multivariable Feedback Control by
        Sigurd Skogestad.

        An alternative definition can be found here:
        https://www.mathworks.com/help/control/ref/lft.html

        Parameters
        ----------
        other: LTI
            The lower LTI system
        ny: int, optional
            Dimension of (plant) measurement output.
        nu: int, optional
            Dimension of (plant) control input.

        """
        other = _convertToStateSpace(other)
        # maximal values for nu, ny
        if ny == -1:
            ny = min(other.inputs, self.outputs)
        if nu == -1:
            nu = min(other.outputs, self.inputs)
        # dimension check
        # TODO

        # Figure out the sampling time to use
        if (self.dt == None and other.dt != None):
            dt = other.dt       # use dt from second argument
        elif (other.dt == None and self.dt != None) or \
                timebaseEqual(self, other):
            dt = self.dt        # use dt from first argument
        else:
            raise ValueError("Systems have different time bases")

        # submatrices
        A = self.A
        B1 = self.B[:, :self.inputs - nu]
        B2 = self.B[:, self.inputs - nu:]
        C1 = self.C[:self.outputs - ny, :]
        C2 = self.C[self.outputs - ny:, :]
        D11 = self.D[:self.outputs - ny, :self.inputs - nu]
        D12 = self.D[:self.outputs - ny, self.inputs - nu:]
        D21 = self.D[self.outputs - ny:, :self.inputs - nu]
        D22 = self.D[self.outputs - ny:, self.inputs - nu:]

        # submatrices
        Abar = other.A
        Bbar1 = other.B[:, :ny]
        Bbar2 = other.B[:, ny:]
        Cbar1 = other.C[:nu, :]
        Cbar2 = other.C[nu:, :]
        Dbar11 = other.D[:nu, :ny]
        Dbar12 = other.D[:nu, ny:]
        Dbar21 = other.D[nu:, :ny]
        Dbar22 = other.D[nu:, ny:]

        # well-posed check
        F = np.block([[np.eye(ny), -D22], [-Dbar11, np.eye(nu)]])
        if matrix_rank(F) != ny + nu:
            raise ValueError("lft not well-posed to working precision.")
        
        # solve for the resulting ss by solving for [y, u] using [x,
        # xbar] and [w1, w2].
        TH = np.linalg.solve(F, np.block(
            [[C2, np.zeros((ny, other.states)), D21, np.zeros((ny, other.inputs - ny))],
             [np.zeros((nu, self.states)), Cbar1, np.zeros((nu, self.inputs - nu)), Dbar12]]
        ))
        T11 = TH[:ny, :self.states]
        T12 = TH[:ny, self.states: self.states + other.states]
        T21 = TH[ny:, :self.states]
        T22 = TH[ny:, self.states: self.states + other.states]
        H11 = TH[:ny, self.states + other.states: self.states + other.states + self.inputs - nu]
        H12 = TH[:ny, self.states + other.states + self.inputs - nu:]
        H21 = TH[ny:, self.states + other.states: self.states + other.states + self.inputs - nu]
        H22 = TH[ny:, self.states + other.states + self.inputs - nu:]
        
        Ares = np.block([
            [A + B2.dot(T21), B2.dot(T22)],
            [Bbar1.dot(T11), Abar + Bbar1.dot(T12)]
        ])

        Bres = np.block([
            [B1 + B2.dot(H21), B2.dot(H22)],
            [Bbar1.dot(H11), Bbar2 + Bbar1.dot(H12)]
        ])

        Cres = np.block([
            [C1 + D12.dot(T21), D12.dot(T22)],
            [Dbar21.dot(T11), Cbar2 + Dbar21.dot(T12)]
        ])

        Dres = np.block([
            [D11 + D12.dot(H21), D12.dot(H22)],
            [Dbar21.dot(H11), Dbar22 + Dbar21.dot(H12)]
        ])
        return StateSpace(Ares, Bres, Cres, Dres, dt)

    def train_on_texts(self, texts, context_labels=None,
                       batch_size=128,
                       num_epochs=50,
                       verbose=1,
                       new_model=False,
                       gen_epochs=1,
                       train_size=1.0,
                       max_gen_length=300,
                       validation=True,
                       dropout=0.0,
                       via_new_model=False,
                       **kwargs):

        if new_model and not via_new_model:
            self.train_new_model(texts,
                                 context_labels=context_labels,
                                 num_epochs=num_epochs,
                                 gen_epochs=gen_epochs,
                                 batch_size=batch_size,
                                 dropout=dropout,
                                 validation=validation,
                                 **kwargs)
            return

        if context_labels:
            context_labels = LabelBinarizer().fit_transform(context_labels)

        if 'prop_keep' in kwargs:
            train_size = prop_keep

        if self.config['word_level']:
            texts = [text_to_word_sequence(text, filters='') for text in texts]

        # calculate all combinations of text indices + token indices
        indices_list = [np.meshgrid(np.array(i), np.arange(
            len(text) + 1)) for i, text in enumerate(texts)]
        indices_list = np.block(indices_list)

        # If a single text, there will be 2 extra indices, so remove them
        # Also remove first sequences which use padding
        if self.config['single_text']:
            indices_list = indices_list[self.config['max_length']:-2, :]

        indices_mask = np.random.rand(indices_list.shape[0]) < train_size

        gen_val = None
        val_steps = None
        if train_size < 1.0 and validation:
            indices_list_val = indices_list[~indices_mask, :]
            gen_val = generate_sequences_from_texts(
                texts, indices_list_val, self, context_labels, batch_size)
            val_steps = max(
                int(np.floor(indices_list_val.shape[0] / batch_size)), 1)

        indices_list = indices_list[indices_mask, :]

        num_tokens = indices_list.shape[0]
        assert num_tokens >= batch_size, "Fewer tokens than batch_size."

        level = 'word' if self.config['word_level'] else 'character'
        print("Training on {:,} {} sequences.".format(num_tokens, level))

        steps_per_epoch = max(int(np.floor(num_tokens / batch_size)), 1)

        gen = generate_sequences_from_texts(
            texts, indices_list, self, context_labels, batch_size)

        base_lr = 4e-3

        # scheduler function must be defined inline.
        def lr_linear_decay(epoch):
            return (base_lr * (1 - (epoch / num_epochs)))

        if context_labels is not None:
            if new_model:
                weights_path = None
            else:
                weights_path = "{}_weights.hdf5".format(self.config['name'])
                self.save(weights_path)

            self.model = textgenrnn_model(self.num_classes,
                                          dropout=dropout,
                                          cfg=self.config,
                                          context_size=context_labels.shape[1],
                                          weights_path=weights_path)

        self.model.fit_generator(gen, steps_per_epoch=steps_per_epoch,
                                 epochs=num_epochs,
                                 callbacks=[
                                     LearningRateScheduler(
                                         lr_linear_decay),
                                     generate_after_epoch(
                                         self, gen_epochs,
                                         max_gen_length),
                                     save_model_weights(
                                         self.config['name'])],
                                 verbose=verbose,
                                 max_queue_size=2,
                                 validation_data=gen_val,
                                 validation_steps=val_steps
#.........这里部分代码省略.........