Python scipy.sparse.block_diag() Examples

def batch_of_relational_supports_to_support(batched_relational_supports,
                                            ordering=None):
    if ordering is not None:
        relational_supports_combined = OrderedDict([
            (k, sparse.block_diag([s[k] for s in batched_relational_supports]))
            for k in ordering])
    else:
        first_dict = batched_relational_supports[0]

        if not isinstance(first_dict, OrderedDict):
            ValueError("If you do not provide an ordering, then "
                       "`batched_relational_supports` shguld be list of "
                       "`OrderedDict`. It is a "
                       "list of {}".format(type(first_dict)))

        first_dict_order = list(first_dict.keys())

        relational_supports_combined = OrderedDict([
            (k, sparse.block_diag([s[k] for s in batched_relational_supports]))
            for k in first_dict_order])

    return relational_supports_to_support(relational_supports_combined,
                                          ordering=ordering) 

def compute_dr_wrt(self, wrt):
        if wrt is not self.a:
            return None
    
        Ainv = self.r

        if Ainv.ndim <= 2:
            return -np.kron(Ainv, Ainv.T)
        else:
            Ainv = np.reshape(Ainv,  (-1, Ainv.shape[-2], Ainv.shape[-1]))
            AinvT = np.rollaxis(Ainv, -1, -2)
            AinvT = np.reshape(AinvT, (-1, AinvT.shape[-2], AinvT.shape[-1]))
            result = np.dstack([-np.kron(Ainv[i], AinvT[i]).T for i in range(Ainv.shape[0])]).T
            result = sp.block_diag(result)

        return result 

def aveF2CCV(self):
        """
        averaging operator of x-faces (radial) to cell centered vectors

        Returns
        -------
        scipy.sparse.csr_matrix
            matrix that averages from faces to cell centered vectors
        """
        if getattr(self, '_aveF2CCV', None) is None:
            # n = self.vnC
            if self.isSymmetric is True:
                self._aveF2CCV = sp.block_diag(
                    (self.aveFx2CC, self.aveFz2CC), format="csr"
                )
            else:
                self._aveF2CCV = sp.block_diag(
                    (self.aveFx2CC, self.aveFy2CC, self.aveFz2CC),
                    format="csr"
                )
        return self._aveF2CCV

    ####################################################
    # Deflation Matrices
    #################################################### 

def cellGradBC(self):
        """
        The cell centered Gradient boundary condition matrix
        """
        if getattr(self, '_cellGradBC', None) is None:
            BC = self.setCellGradBC(self._cellGradBC_list)
            n = self.vnC
            if self.dim == 1:
                G = ddxCellGradBC(n[0], BC[0])
            elif self.dim == 2:
                G1 = sp.kron(speye(n[1]), ddxCellGradBC(n[0], BC[0]))
                G2 = sp.kron(ddxCellGradBC(n[1], BC[1]), speye(n[0]))
                G = sp.block_diag((G1, G2), format="csr")
            elif self.dim == 3:
                G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGradBC(n[0], BC[0]))
                G2 = kron3(speye(n[2]), ddxCellGradBC(n[1], BC[1]), speye(n[0]))
                G3 = kron3(ddxCellGradBC(n[2], BC[2]), speye(n[1]), speye(n[0]))
                G = sp.block_diag((G1, G2, G3), format="csr")
            # Compute areas of cell faces & volumes
            S = self.area
            V = self.aveCC2F*self.vol  # Average volume between adjacent cells
            self._cellGradBC = sdiag(S/V)*G
        return self._cellGradBC 

def aveCCV2F(self):
        """
        Construct the averaging operator on cell centers to
        faces as a vector.
        """
        if getattr(self, '_aveCCV2F', None) is None:
            if self.dim == 1:
                self._aveCCV2F = self.aveCC2F
            elif self.dim == 2:
                aveCCV2Fx = sp.kron(speye(self.nCy), av_extrap(self.nCx))
                aveCC2VFy = sp.kron(av_extrap(self.nCy), speye(self.nCx))
                self._aveCCV2F = sp.block_diag((
                    aveCCV2Fx, aveCC2VFy
                ), format="csr")
            elif self.dim == 3:
                aveCCV2Fx = kron3(
                    speye(self.nCz), speye(self.nCy), av_extrap(self.nCx)
                )
                aveCC2VFy = kron3(
                    speye(self.nCz), av_extrap(self.nCy), speye(self.nCx)
                )
                aveCC2BFz = kron3(
                    av_extrap(self.nCz), speye(self.nCy), speye(self.nCx)
                )
                self._aveCCV2F = sp.block_diag((
                        aveCCV2Fx, aveCC2VFy, aveCC2BFz
                ), format="csr")
        return self._aveCCV2F 

def solve(self):
        r"""
        Solve for bound :math:`\\Gamma` using the equation;

        .. math::
                \\mathcal{A}(\\Gamma^n) = u^n

        # ... at constant rotation speed
        ``self.Dfqsdzeta+=scalg.block_diag(*ass.dfqsdzeta_omega(MS.Surfs,MS.Surfs_star))``

        """

        MS = self.MS
        t0 = time.time()

        ### state
        bv = np.dot(self.Ducdu_ext, self.u_ext - self.zeta_dot) + \
             np.dot(self.Ducdzeta, self.zeta)
        self.gamma = np.linalg.solve(self.AIC, -bv)

        ### retrieve gamma over wake
        gamma_star = []
        Ktot = 0
        for ss in range(MS.n_surf):
            Ktot += MS.Surfs[ss].maps.K
            N = MS.Surfs[ss].maps.N
            Mstar = MS.Surfs_star[ss].maps.M
            gamma_star.append(np.concatenate(Mstar * [self.gamma[Ktot - N:Ktot]]))
        gamma_star = np.concatenate(gamma_star)

        ### compute steady force
        self.fqs = np.dot(self.Dfqsdgamma, self.gamma) + \
                   np.dot(self.Dfqsdgamma_star, gamma_star) + \
                   np.dot(self.Dfqsdzeta, self.zeta) + \
                   np.dot(self.Dfqsdu_ext, self.u_ext - self.zeta_dot)

        self.time_sol = time.time() - t0
        print('Solution done in %.2f sec' % self.time_sol) 

def _generate_qp_problem(self):
        '''
        Generate QP problem
        '''

        # Construct the problem
        #       minimize	y' * y + lambda * 1' * t
        #       subject to  y = Ax - b
        #                   -t <= x <= t
        P = spa.block_diag((spa.csc_matrix((self.n, self.n)),
                            2*spa.eye(self.m),
                            spa.csc_matrix((self.n, self.n))), format='csc')
        q = np.append(np.zeros(self.m + self.n),
                      self.lambda_param * np.ones(self.n))
        In = spa.eye(self.n)
        Onm = spa.csc_matrix((self.n, self.m))
        A = spa.vstack([spa.hstack([self.Ad, -spa.eye(self.m),
                                    spa.csc_matrix((self.m, self.n))]),
                        spa.hstack([In, Onm, -In]),
                        spa.hstack([In, Onm, In])]).tocsc()
        l = np.hstack([self.bd, -np.inf * np.ones(self.n), np.zeros(self.n)])
        u = np.hstack([self.bd, np.zeros(self.n), np.inf * np.ones(self.n)])

        problem = {}
        problem['P'] = P
        problem['q'] = q
        problem['A'] = A
        problem['l'] = l
        problem['u'] = u
        problem['m'] = A.shape[0]
        problem['n'] = A.shape[1]

        return problem 

def _generate_qp_problem(self):
        '''
        Generate QP problem
        '''

        # Construct the problem
        #       minimize	y' * y + lambda * 1' * t
        #       subject to  y = Ax - b
        #                   -t <= x <= t
        P = spa.block_diag((spa.csc_matrix((self.n, self.n)),
                            2*spa.eye(self.m),
                            spa.csc_matrix((self.n, self.n))), format='csc')
        q = np.append(np.zeros(self.m + self.n),
                      self.lambda_param * np.ones(self.n))
        In = spa.eye(self.n)
        Onm = spa.csc_matrix((self.n, self.m))
        A = spa.vstack([spa.hstack([self.Ad, -spa.eye(self.m),
                                    spa.csc_matrix((self.m, self.n))]),
                        spa.hstack([In, Onm, -In]),
                        spa.hstack([In, Onm, In])]).tocsc()
        l = np.hstack([self.bd, -np.inf * np.ones(self.n), np.zeros(self.n)])
        u = np.hstack([self.bd, np.zeros(self.n), np.inf * np.ones(self.n)])

        problem = {}
        problem['P'] = P
        problem['q'] = q
        problem['A'] = A
        problem['l'] = l
        problem['u'] = u
        problem['m'] = A.shape[0]
        problem['n'] = A.shape[1]

        return problem 

def _make_sigma(self, TNTs, phiinv):
        return sps.block_diag(TNTs, "csc") + sps.csc_matrix(phiinv) 

def cellGrad(self):
        """
        Cell centered Gradient operator built off of the faceDiv operator.
        Grad =  - (Mf)^{-1} * Div * diag (volume)
        """
        if getattr(self, '_cellGrad', None) is None:

            i_s = self.faceBoundaryInd

            ix = np.ones(self.nFx)
            ix[i_s[0]] = 0.
            ix[i_s[1]] = 0.
            Pafx = sp.diags(ix)

            iy = np.ones(self.nFy)
            iy[i_s[2]] = 0.
            iy[i_s[3]] = 0.
            Pafy = sp.diags(iy)

            MfI = self.getFaceInnerProduct(invMat=True)

            if self.dim == 2:
                Pi = sp.block_diag([Pafx, Pafy])

            elif self.dim == 3:
                iz = np.ones(self.nFz)
                iz[i_s[4]] = 0.
                iz[i_s[5]] = 0.
                Pafz = sp.diags(iz)
                Pi = sp.block_diag([Pafx, Pafy, Pafz])

            self._cellGrad = -Pi * MfI * self.faceDiv.T * sp.diags(self.vol)

        return self._cellGrad 

def aveE2CCV(self):
        "Construct the averaging operator on cell edges to cell centers."
        if getattr(self, '_aveE2CCV', None) is None:
            if self.dim == 1:
                self._aveE2CCV = self.aveEx2CC
            elif self.dim == 2:
                self._aveE2CCV = sp.block_diag(
                    (self.aveEx2CC, self.aveEy2CC), format="csr"
                )
            elif self.dim == 3:
                self._aveE2CCV = sp.block_diag(
                    (self.aveEx2CC, self.aveEy2CC, self.aveEz2CC), format="csr"
                )
        return self._aveE2CCV 

def _lazy_build_elements(self):
        self._elements = []

        for vel in self.vector_elements:
            vstart = 0
            if self.sparse:  # build block-diagonal sparse mx
                diagBlks = []
                for compbasis in self.component_bases:
                    cs = compbasis.elshape
                    comp_vel = vel[vstart:vstart + compbasis.dim]
                    diagBlks.append(comp_vel.reshape(cs))
                    vstart += compbasis.dim
                el = _sps.block_diag(diagBlks, "csr", 'complex')

            else:
                start = [0] * self.elndim
                el = _np.zeros(self.elshape, 'complex')
                for compbasis in self.component_bases:
                    cs = compbasis.elshape
                    comp_vel = vel[vstart:vstart + compbasis.dim]
                    slc = tuple([slice(start[k], start[k] + cs[k]) for k in range(self.elndim)])
                    el[slc] = comp_vel.reshape(cs)
                    vstart += compbasis.dim
                    for k in range(self.elndim): start[k] += cs[k]

            self._elements.append(el)
        if not self.sparse:  # _elements should be an array rather than a list
            self._elements = _np.array(self._elements) 

def aveF2CCV(self):
        "Construct the averaging operator on cell faces to cell centers."
        if getattr(self, '_aveF2CCV', None) is None:
            if self.dim == 1:
                self._aveF2CCV = self.aveFx2CC
            elif self.dim == 2:
                self._aveF2CCV = sp.block_diag((
                    self.aveFx2CC, self.aveFy2CC
                ), format="csr")
            elif self.dim == 3:
                self._aveF2CCV = sp.block_diag((
                    self.aveFx2CC, self.aveFy2CC, self.aveFz2CC
                ), format="csr")
        return self._aveF2CCV 

def symmetry_matrix2D_sign_correction(x_vector: np.array,
                                      y_vector: np.array,
                                      axis: int,
                                      sign: int = 1,
                                      periodicity: bool = False):
    '''
    SymmetryMatrix2D_sign_correction generates a matrixon3 that corrects the sign of every
        parameter component (f, fx, fy, fxy) according to the parametrization being symmetric
        or anti-symmetric.

        The correct symmetry_matrix for [f, fx, fy, fxy] will be of the form:
            sign_matrix@block_diag(4*symmetry matrix)

    Args:
     x_vector: x vector of the fine grid
     y_vector: y vector of the fine grid
     axis: Symmetry axis.
     sign: 1 or -1 depanding whether or not you have symmetry or anti-symmetry.
     periodicity: Indicating whether or not there is periodicity.

    Return:
     sign_correction matrix


    '''
    n_x = len(x_vector)
    n_y = len(y_vector)
    shape = (n_x, n_y)
    return symmetry_matrix_sign_correction(shape, axis, sign, periodicity) 

def duplicate_boundary_data(shape, axis: int) -> sparse.csc.csc_matrix:
    '''
    duplicate_boundary_data creates a sparse matrix that duplicated the vectorized data at
    a boundary along the axis.
    For example: if you have a space A that is 4x5 represented by vector with 20 elements. Then
        duplicate_boundary_data(x_in, y_in, axis: 0 ) will generate a 25x20 matrix that will make
        a vector representation of a 5x5 matrix where the first row is duplicated as last row.
    (used to implement periodicity)
    '''
    if axis == 0:
        block_a = sparse.eye(shape[0] - 1).tocsr()
        block_b = sparse.csr_matrix((1, shape[0] - 1))
        block_b[0, 0] = 1
        block_matrix = sparse.vstack([block_a, block_b])
        per = sparse.block_diag(tuple(shape[1] * [block_matrix]))
        block_c = sparse.csr_matrix((shape[0] - 1, 1))
        block_matrix = sparse.hstack([block_a, block_c])
        per_reverse = sparse.block_diag(tuple(shape[1] * [block_matrix]))
    elif axis == 1:
        block_1 = sparse.eye(int(shape[0] * (shape[1] - 1)))
        block_b = sparse.eye(shape[0]).tocsr()
        block_2 = sparse.hstack([
            block_b,
            sparse.csr_matrix((shape[0], shape[0] * int(shape[1] - 2)))
        ])
        per = sparse.vstack([block_1, block_2])
        block_1 = sparse.eye(int(shape[0] * (shape[1] - 1)))
        block_2 = sparse.csr_matrix((shape[0] * int(shape[1] - 1), shape[0]))
        per_reverse = sparse.hstack([block_1, block_2])
    return per, per_reverse


# Geometry matrix functions
############################ 

def compute_d2(self):
        
        # To stay consistent with numpy, we must upgrade 1D arrays to 2D
        mtx1r = row(self.mtx1.r) if len(self.mtx1.r.shape)<2 else self.mtx1.r
        mtx2r = col(self.mtx2.r) if len(self.mtx2.r.shape)<2 else self.mtx2.r

        if mtx2r.ndim <= 1:
            return self.mtx1r
        elif mtx2r.ndim <= 2:
            return sp.kron(mtx1r, sp.eye(mtx2r.shape[1],mtx2r.shape[1]))
        else:
            mtx1f = mtx1r.reshape((-1, mtx1r.shape[-2], mtx1r.shape[-1]))            
            result = sp.block_diag([np.kron(m1, np.eye(mtx2r.shape[-1],mtx2r.shape[-1])) for m1 in mtx1f])
            assert(result.shape[0] == self.r.size)
            return result 

def compute_d1(self):
        # To stay consistent with numpy, we must upgrade 1D arrays to 2D
        mtx1r = row(self.mtx1.r) if len(self.mtx1.r.shape)<2 else self.mtx1.r
        mtx2r = col(self.mtx2.r) if len(self.mtx2.r.shape)<2 else self.mtx2.r

        if mtx1r.ndim <= 2:
            return sp.kron(sp.eye(mtx1r.shape[0], mtx1r.shape[0]),mtx2r.T)
        else:
            mtx2f = mtx2r.reshape((-1, mtx2r.shape[-2], mtx2r.shape[-1]))
            mtx2f = np.rollaxis(mtx2f, -1, -2) #transpose basically            
            result = sp.block_diag([np.kron(np.eye(mtx1r.shape[-2], mtx1r.shape[-2]),m2) for m2 in mtx2f])
            assert(result.shape[0] == self.r.size)
            return result 

def sparse_factor_solve_kkt(Q_tilde, D_tilde, A_, C_tilde, rx, rs, rz, ry, ns):
    nineq, nz, neq, nBatch = ns

    # H_ = csc_matrix((nz + nineq, nz + nineq))
    # H_[:nz, :nz] = Q_tilde
    # H_[-nineq:, -nineq:] = D_tilde
    H_ = block_diag([Q_tilde, D_tilde], format='csc')
    if neq > 0:
        g_ = torch.cat([rx, rs], 1).squeeze(0).numpy()
        h_ = torch.cat([rz, ry], 1).squeeze(0).numpy()
    else:
        g_ = torch.cat([rx, rs], 1).squeeze(0).numpy()
        h_ = rz.squeeze(0).numpy()

    H_LU = splu(H_)

    invH_A_ = csc_matrix(H_LU.solve(A_.todense().transpose()))
    invH_g_ = H_LU.solve(g_)

    S_ = A_.dot(invH_A_) + C_tilde
    S_LU = splu(S_)
    # t_ = invH_g_[np.newaxis].dot(A_.transpose()).squeeze(0) - h_
    t_ = A_.dot(invH_g_) - h_
    w_ = -S_LU.solve(t_)
    # t_ = -g_ - w_[np.newaxis].dot(A_).squeeze(0)
    t_ = -g_ - A_.transpose().dot(w_)
    v_ = H_LU.solve(t_)

    dx = v_[:nz]
    ds = v_[nz:]
    dz = w_[:nineq]
    dy = w_[nineq:] if neq > 0 else None

    dx = torch.DoubleTensor(dx).unsqueeze(0)
    ds = torch.DoubleTensor(ds).unsqueeze(0)
    dz = torch.DoubleTensor(dz).unsqueeze(0)
    dy = torch.DoubleTensor(dy).unsqueeze(0) if neq > 0 else None

    return dx, ds, dz, dy 

def _test_disjoint_mode(layer, **kwargs):
    A = sp.block_diag([np.ones((N1, N1)), np.ones(
        (N2, N2)), np.ones((N3, N3))]).todense()
    X = np.random.normal(size=(N, F))
    I = np.array([0] * N1 + [1] * N2 + [2] * N3).astype(int)
    sparse = kwargs.pop('sparse', None) is not None

    A_in = Input(shape=(None,), sparse=sparse)
    X_in = Input(shape=(F,))
    I_in = Input(shape=(), dtype=tf.int32)

    layer_instance = layer(**kwargs)
    output = layer_instance([X_in, A_in, I_in])
    model = Model([X_in, A_in, I_in], output)
    output = model([X, A, I])
    X_pool, A_pool, I_pool, mask = output

    N_pool_expected = np.ceil(kwargs['ratio'] * N1) + \
                      np.ceil(kwargs['ratio'] * N2) + \
                      np.ceil(kwargs['ratio'] * N3)
    N_pool_expected = int(N_pool_expected)
    N_pool_true = A_pool.shape[0]

    _check_number_of_nodes(N_pool_expected, N_pool_true)

    assert X_pool.shape == (N_pool_expected, F)
    assert A_pool.shape == (N_pool_expected, N_pool_expected)
    assert I_pool.shape == (N_pool_expected,)

    output_shape = [o.shape for o in output]
    _check_output_and_model_output_shapes(output_shape, model.output_shape) 

def numpy_to_disjoint(X_list, A_list, E_list=None):
    """
    Converts a batch of graphs stored in lists (X, A, and optionally E) to the
    [disjoint mode](https://danielegrattarola.github.io/spektral/data/#disjoint-mode).

    Each entry i of the lists should be associated to the same graph, i.e.,
    `X_list[i].shape[0] == A_list[i].shape[0] == E_list[i].shape[0]`.

    :param X_list: a list of np.arrays of shape `(N, F)`;
    :param A_list: a list of np.arrays or sparse matrices of shape `(N, N)`;
    :param E_list: a list of np.arrays of shape `(N, N, S)`;
    :return:
        -  `X_out`: a rank 2 array of shape `(n_nodes, F)`;
        -  `A_out`: a rank 2 array of shape `(n_nodes, n_nodes)`;
        -  `E_out`: (only if `E_list` is given) a rank 2 array of shape
        `(n_edges, S)`;
    """
    X_out = np.vstack(X_list)
    A_list = [sp.coo_matrix(a) for a in A_list]
    if E_list is not None:
        if E_list[0].ndim == 3:
            E_list = [e[a.row, a.col] for e, a in zip(E_list, A_list)]
        E_out = np.vstack(E_list)
    A_out = sp.block_diag(A_list)
    n_nodes = np.array([x.shape[0] for x in X_list])
    I_out = np.repeat(np.arange(len(n_nodes)), n_nodes)
    if E_list is not None:
        return X_out, A_out, E_out, I_out
    else:
        return X_out, A_out, I_out 

def collate_fn(batch):
        graphs, pmpds, labels = zip(*batch)
        batched_graphs = graph_batch(graphs)
        batched_pmpds = sp.block_diag(pmpds)
        batched_labels = np.concatenate(labels, axis=0)
        return batched_graphs, batched_pmpds, batched_labels 

def periodic_matrix(shape: list, axis: int,
                    periods=int) -> sparse.coo.coo_matrix:
    '''
    periodic_matrix generated a sparse matrix that duplicates a
    parametrization vector according to a give periodicity

    If the periods is 0, it will only duplicate the outer parametrization in
    the axis direction, rather then duplicating the parametrization. (this
    is usefull to describe a structure with periodic boundary conditions)

    Args:
      shape: Shape of the matrix on which you want to make periodic.
      axis: Axis in which you want to have periodicity.
      periods: Amount of periods.

    Returns:
      periodicity matrix
    '''
    n_x = shape[0]
    n_y = shape[1]
    if axis == 0:
        if n_x % periods > 0:
            raise ValueError('x cannot be divided by periods')
        else:
            n_x_periodic = n_x // periods
            block_a = sparse.eye(n_x_periodic).tocsr()
            block_b = sparse.vstack(periods * [block_a])
            geometry_matrix = sparse.block_diag(tuple(n_y * [block_b]))
            block_a = sparse.eye(n_x_periodic).tocsr()
            matrix_list = periods * [sparse.csr_matrix(block_a.shape)]
            matrix_list[0] = block_a
            block_b = sparse.vstack(matrix_list)
            geometry_matrix_reverse = sparse.block_diag(tuple(
                n_y * [block_b])).transpose()
    elif axis == 1:
        if n_y % periods > 0:
            raise ValueError('y cannot be divided by periods')
        else:
            n_y_periodic = n_y // periods
            n = n_y_periodic * n_x
            block_a = sparse.eye(n).tocsr()
            geometry_matrix = sparse.vstack(periods * [block_a])
            matrix_list = periods * [sparse.csr_matrix(block_a.shape)]
            matrix_list[0] = block_a
            geometry_matrix_reverse = sparse.vstack(matrix_list).transpose()

    return geometry_matrix, geometry_matrix_reverse


#deprecated function: use periodic_matrix 

def symmetry_matrix_sign_correction(shape: tuple,
                                    axis: int,
                                    sign: int = 1,
                                    periodicity: bool = False):
    '''
    SymmetryMatrix2D_sign_correction generates a matrix that corrects the sign of every
        parameter component (f, fx, fy, fxy) according to the parametrization being symmetric
        or anti-symmetric.

        The correct symmetry_matrix for [f, fx, fy, fxy] will be of the form:
            sign_matrix@block_diag(4*symmetry matrix)

    Args:
        shape: Shape of the parametrization.
        axis: Symmetry axis.
        sign: 1 or -1 depanding whether or not you have symmetry or anti-symmetry.
        periodicity: Indicating whether or not there is periodicity.

    Returns:
        periodicity matrix
    '''
    n_x = shape[0]
    n_y = shape[1]
    periodic = int(not periodicity)

    if axis == 1:
        periodic = int(not periodicity)
        block_1 = sparse.eye(int(n_x * (n_y // 2)))
        block_c = sign * (sign > 0) * sparse.eye(n_x, n_x)
        block_2 = sign * sparse.eye(int(n_x * (n_y // 2)))
        zero_0 = sparse.diags([periodic] * n_x + [1] * n_x * (n_y - 2) +
                              [periodic] * n_x)
        if bool(n_y % 2):
            sym_corr = sparse.block_diag((block_1, block_c, block_2)) @ zero_0
        else:
            sym_corr = sparse.block_diag((block_1, block_2)) @ zero_0
    elif axis == 0:
        block_1 = sparse.eye(n_x // 2)
        block_c = sign * (sign > 0) * sparse.eye(1, 1)
        block_2 = sign * sparse.eye(n_x // 2)
        zero_0 = sparse.diags([periodic] + [1] * (n_x - 2) + [periodic])
        if bool(n_x % 2):
            sym_corr_row = sparse.block_diag(
                (block_1, block_c, block_2)) @ zero_0
        else:
            sym_corr_row = sparse.block_diag((block_1, block_2)) @ zero_0
        sym_corr = sparse.block_diag(tuple(n_y * [sym_corr_row]))

    return sym_corr


#deprecated function: use symmetry_matrix_sign_correction 

def symmetry_matrix(shape: tuple, axis: int) -> sparse.coo.coo_matrix:
    '''
    SymmetryMatrix generated a symmetry matrix.
        For example, given a parametrization defined by an xy-grid that
        is symmetric in the x axis, this function will return a matrix that
        turn a parametrization that describes the right side of the grid in a
        parametrization over the entire grid. The matrix size will thus be
        (shape[0]*shape[1])x(shape[0]/2*shape[1]).
        Note: you give the symmetry axis if axis is 0 the second coordinate is
            symmetric. If the axis is 1 the first coordinate is symmetric

    Args:
        shape: Shape of the parametrization.
        axis: Symmetry axis.
    '''
    n_x = shape[0]
    n_y = shape[1]

    if axis == 0:
        block_1 = sparse.eye(int(n_x * np.ceil(n_y / 2)))
        block_b = sparse.eye(n_x).tocsr()
        block_b.indices = -1 * block_b.indices + block_b.shape[1] - 1
        block_2 = sparse.block_diag(tuple(int(n_y / 2) * [block_b])).tocsr()
        block_2.indices = -1 * block_2.indices + block_2.shape[1] - 1
        if bool(n_y % 2):
            block_2 = sparse.hstack(
                [block_2, sparse.csr_matrix((n_x * int(n_y / 2), n_x))])
        sym = sparse.vstack([block_1, block_2])
        sym_reverse = sparse.hstack(
            [block_1, sparse.csr_matrix(block_2.shape).T])
    elif axis == 1:
        block_a = sparse.eye(int(np.ceil(n_x / 2))).tocsr()
        block_b = sparse.eye(int(n_x / 2)).tocsr()
        block_b.indices = -1 * block_b.indices + block_b.shape[1] - 1
        if n_x % 2 == 1:
            block_b = sparse.hstack([block_b, sparse.csr_matrix((n_x // 2, 1))])
        block_matrix = sparse.vstack([block_a, block_b])
        sym = sparse.block_diag(tuple(n_y * [block_matrix]))
        block_matrix_reverse = sparse.vstack(
            [block_a, sparse.csr_matrix(block_b.shape)])
        sym_reverse = sparse.block_diag(tuple(
            n_y * [block_matrix_reverse])).transpose()

    return sym, sym_reverse


#deprecated function: use symmetry_matrix 

def get_phi(self, params, cliques=False):
        phis = [signalcollection.get_phi(params) for signalcollection in self._signalcollections]

        # if we found common signals, we'll return a big phivec matrix,
        # otherwise a list of phivec vectors (some of which possibly None)
        if self._commonsignals:
            if np.any([phi.ndim == 2 for phi in phis if phi is not None]):
                # if we have any dense matrices,
                Phi = sl.block_diag(*[np.diag(phi) if phi.ndim == 1 else phi for phi in phis if phi is not None])
            else:
                Phi = np.diag(np.concatenate([phi for phi in phis if phi is not None]))

            # get a dictionary of slices locating each pulsar in Phi matrix
            slices = self._get_slices(phis)

            # self._cliques is a vector of the same size as the Phi matrix
            # for each Phi index i, self._cliques[i] is -1 if row/column
            # belong to no clique, or it gives the clique number otherwise
            if cliques:
                self._resetcliques(Phi.shape[0])
                self._setpulsarcliques(slices, phis)

            # iterate over all common signal classes
            for csclass, csdict in self._commonsignals.items():
                # first figure out which indices are used in this common signal
                # and update the clique index
                if cliques:
                    self._setcliques(slices, csdict)

                # now iterate over all pairs of common signal instances
                pairs = itertools.combinations(csdict.items(), 2)

                for (cs1, csc1), (cs2, csc2) in pairs:
                    crossdiag = csclass.get_phicross(cs1, cs2, params)

                    block1, idx1 = slices[csc1], csc1._idx[cs1]
                    block2, idx2 = slices[csc2], csc2._idx[cs2]

                    if crossdiag.ndim == 1:
                        Phi[block1, block2][idx1, idx2] += crossdiag
                        Phi[block2, block1][idx2, idx1] += crossdiag
                    else:
                        Phi[block1, block2][np.ix_(idx1, idx2)] += crossdiag
                        Phi[block2, block1][np.ix_(idx2, idx1)] += crossdiag

            return Phi
        else:
            return phis 

def _generate_qp_problem(self):
        '''
        Generate QP problem
        '''

        # Construct the problem
        #       minimize	 x.T * x + gamma 1.T * t
        #       subject to  t >= diag(b) A x + 1
        #                   t >= 0

        P = spa.block_diag((spa.eye(self.n),
                            spa.csc_matrix((self.m, self.m))), format='csc')
        q = np.append(np.zeros(self.n), (self.gamma/2) * np.ones(self.m))
        A = spa.vstack([spa.hstack([spa.diags(self.b_svm).dot(self.A_svm),
                                    -spa.eye(self.m)]),
                        spa.hstack([spa.csc_matrix((self.m, self.n)),
                                    spa.eye(self.m)])
                        ]).tocsc()
        l = np.hstack([-np.inf*np.ones(self.m), np.zeros(self.m)])
        u = np.hstack([-np.ones(self.m), np.inf*np.ones(self.m)])

        # Constraints without bounds
        A_nobounds = spa.hstack([spa.diags(self.b_svm).dot(self.A_svm),
                                 -spa.eye(self.m)]).tocsc()
        l_nobounds = -np.inf*np.ones(self.m)
        u_nobounds = -np.ones(self.m)
        bounds_idx = np.arange(self.n, self.n + self.m)

        # Separate bounds
        lx = np.append(-np.inf*np.ones(self.n), np.zeros(self.m))
        ux = np.append(np.inf*np.ones(self.n), np.inf*np.ones(self.m))

        problem = {}
        problem['P'] = P
        problem['q'] = q
        problem['A'] = A
        problem['l'] = l
        problem['u'] = u
        problem['m'] = A.shape[0]
        problem['n'] = A.shape[1]
        problem['A_nobounds'] = A_nobounds
        problem['l_nobounds'] = l_nobounds
        problem['u_nobounds'] = u_nobounds
        problem['bounds_idx'] = bounds_idx
        problem['lx'] = lx
        problem['ux'] = ux

        return problem 

def _generate_qp_problem(self):
        '''
        Generate QP problem
        '''

        # Construct the problem
        #       minimize	x' D x + y' I y - (1/gamma) * mu' x
        #       subject to  1' x = 1
        #                   F' x = y
        #                   0 <= x <= 1
        P = spa.block_diag((2 * self.D, 2 * spa.eye(self.k)), format='csc')
        q = np.append(- self.mu / self.gamma, np.zeros(self.k))
        A = spa.vstack([
                spa.hstack([spa.csc_matrix(np.ones((1, self.n))),
                           spa.csc_matrix((1, self.k))]),
                spa.hstack([self.F.T, -spa.eye(self.k)]),
                spa.hstack((spa.eye(self.n), spa.csc_matrix((self.n, self.k))))
            ]).tocsc()
        l = np.hstack([1., np.zeros(self.k), np.zeros(self.n)])
        u = np.hstack([1., np.zeros(self.k), np.ones(self.n)])

        # Constraints without bounds
        A_nobounds = spa.vstack([
                spa.hstack([spa.csc_matrix(np.ones((1, self.n))),
                            spa.csc_matrix((1, self.k))]),
                spa.hstack([self.F.T, -spa.eye(self.k)]),
                ]).tocsc()
        l_nobounds = np.hstack([1., np.zeros(self.k)])
        u_nobounds = np.hstack([1., np.zeros(self.k)])
        bounds_idx = np.arange(self.n)

        # Separate bounds
        lx = np.hstack([np.zeros(self.n), -np.inf * np.ones(self.k)])
        ux = np.hstack([np.ones(self.n), np.inf * np.ones(self.k)])

        problem = {}
        problem['P'] = P
        problem['q'] = q
        problem['A'] = A
        problem['l'] = l
        problem['u'] = u
        problem['m'] = A.shape[0]
        problem['n'] = A.shape[1]

        problem['A_nobounds'] = A_nobounds
        problem['l_nobounds'] = l_nobounds
        problem['u_nobounds'] = u_nobounds
        problem['bounds_idx'] = bounds_idx
        problem['lx'] = lx
        problem['ux'] = ux

        return problem 

def _generate_qp_problem(self):
        '''
        Generate QP problem
        '''
        # Construct the problem
        #       minimize    1/2 z.T * z + np.ones(m).T * (r + s)
        #       subject to  Ax - b - z = r - s
        #                   r >= 0
        #                   s >= 0
        # The problem reformulation follows from Eq. (24) of the following paper:
        # https://doi.org/10.1109/34.877518
        # x_solver = (x, z, r, s)
        Im = spa.eye(self.m)
        P = spa.block_diag((spa.csc_matrix((self.n, self.n)), Im,
                            spa.csc_matrix((2*self.m, 2*self.m))), format='csc')
        q = np.hstack([np.zeros(self.n + self.m), np.ones(2*self.m)])
        A = spa.bmat([[self.Ad, -Im,   -Im,   Im],
                      [None,     None,  Im,   None],
                      [None,     None,  None, Im]], format='csc')
        l = np.hstack([self.bd, np.zeros(2*self.m)])
        u = np.hstack([self.bd, np.inf*np.ones(2*self.m)])

        # Constraints without bounds
        A_nobounds = spa.hstack([self.Ad, -Im, -Im, Im], format='csc')
        l_nobounds = self.bd
        u_nobounds = self.bd

        # Bounds
        lx = np.hstack([-np.inf * np.ones(self.n + self.m),
                        np.zeros(2*self.m)])
        ux = np.inf*np.ones(self.n + 3*self.m)
        bounds_idx = np.arange(self.n + self.m, self.n + 3*self.m)

        problem = {}
        problem['P'] = P
        problem['q'] = q
        problem['A'] = A
        problem['l'] = l
        problem['u'] = u
        problem['m'] = A.shape[0]
        problem['n'] = A.shape[1]
        problem['A_nobounds'] = A_nobounds
        problem['l_nobounds'] = l_nobounds
        problem['u_nobounds'] = u_nobounds
        problem['bounds_idx'] = bounds_idx
        problem['lx'] = lx
        problem['ux'] = ux

        return problem 

def first_deriv(dts, q, n_segments=1, n_simpson_intervals_per_segment=2, order=4):
    """
    This method differentiates a quantity over time using fourth order finite differencing

    A "segment" is defined as a portion of the quantity vector q with a
    constant delta t (or delta x, etc).
    This routine is designed to be used in the context of Simpson's rule integration
    where segment endpoints are coincident in time and n_points_per_seg = 2*n + 1

    Inputs
    ------
    dts : float
        "Time" step between points for each segment. Length of dts must equal n_segments (vector)
        Note that this is half the corresponding Simpson integration timestep
    q : float
        Quantity to be differentiated (vector)
        Total length of q = n_segments * n_points_per_seg
    n_segments : int
        Number of segments to differentiate. Each segment has constant dt (scalar)
    n_simpson_intervals_per_segment : int
        Number of Simpson's rule intervals per segment. Minimum is 2. (scalar)
        The number of points per segment = 2*n + 1
    order : int
        Order of accuracy (choose 2 or 4)

    Returns
    -------
    dqdt : float
        Derivative of q with respect to time (vector)
    """

    n_int_seg = n_simpson_intervals_per_segment
    n_int_tot = n_segments * n_int_seg
    nn_seg = (n_simpson_intervals_per_segment * 2 + 1)
    nn_tot = n_segments * nn_seg
    if order == 4 and n_int_seg < 2:
        raise ValueError('Must use a minimum of 2 Simpson intervals or 5 points per segment due to fourth-order FD stencil')
    if len(q) != nn_tot:
        raise ValueError('q must be of the correct length')
    if len(dts) != n_segments:
        raise ValueError('must provide same number of dts as segments')

    dqdt = np.zeros(q.shape)
    if order == 4:
        stencil_vec, rowidx, colidx = first_deriv_fourth_order_accurate_stencil(nn_seg)
    elif order == 2:
        stencil_vec, rowidx, colidx = first_deriv_second_order_accurate_stencil(nn_seg)
    else:
        raise ValueError('Must choose second or fourth order accuracy')

    stencil_mat = sp.csr_matrix((stencil_vec, (rowidx, colidx)))
    # now we have a generic stencil for each segment
    # assemble a block diagonal overall stencil for the entire vector q

    block_mat_list = []
    for i in range(n_segments):
        dt_seg = dts[i]
        block_mat_list.append(stencil_mat / dt_seg)
    overall_stencil = sp.block_diag(block_mat_list).toarray()
    dqdt = overall_stencil.dot(q)
    return dqdt