Python scipy.sparse.linalg.inv() Examples

def _fit(self):

        # Versions using sparse matrices
        # adj = nx.adjacency_matrix(self._G)
        # ident = sparse.identity(len(self._G.nodes)).tocsc()
        # sim = inv(ident - adj.multiply(self.beta).T) - ident
        # adj = nx.adjacency_matrix(self._G)
        # aux = adj.multiply(-self.beta).T
        # aux.setdiag(1+aux.diagonal(), k=0)
        # sim = inv(aux)
        # sim.setdiag(sim.diagonal()-1)
        # print(sim.nnz)
        # print(adj.nnz)

        # Version using dense matrices
        adj = nx.adjacency_matrix(self._G)
        aux = adj.T.multiply(-self.beta).todense()
        np.fill_diagonal(aux, 1+aux.diagonal())
        sim = np.linalg.inv(aux)
        np.fill_diagonal(sim, sim.diagonal()-1)
        return sim 

def test_inv(self):
        def check(dtype):
            M = array([[1, 0, 2], [0, 0, 3], [-4, 5, 6]], dtype)
            with suppress_warnings() as sup:
                sup.filter(SparseEfficiencyWarning,
                           "spsolve requires A be CSC or CSR matrix format")
                sup.filter(SparseEfficiencyWarning,
                           "spsolve is more efficient when sparse b is in the CSC matrix format")
                sup.filter(SparseEfficiencyWarning,
                           "splu requires CSC matrix format")
                sM = self.spmatrix(M, shape=(3,3), dtype=dtype)
                sMinv = inv(sM)
            assert_array_almost_equal(sMinv.dot(sM).todense(), np.eye(3))
            assert_raises(TypeError, inv, M)
        for dtype in [float]:
            check(dtype) 

def get_from_std(self):
        '''
        Retrieve the matrix that transforms vectors from the standard basis
        to this basis.

        Returns
        -------
        numpy array or scipy sparse matrix
            An array of shape `(size, dim)` where `dim` is the dimension
            of this basis (the length of its vectors) and `size` is the
            size of this basis (its number of vectors).
        '''
        if self.sparse:
            if self.is_complete():
                return _spsl.inv(self.get_to_std().tocsc()).tocsr()
            else:
                assert(self.size < self.dim), "Basis seems to be overcomplete: size > dimension!"
                # we'd need to construct a different pseudo-inverse if the above assert fails

                A = self.get_to_std()  # shape (dim,size) - should have indep *cols*
                Adag = A.getH()        # shape (size, dim)
                invAdagA = _spsl.inv(Adag.tocsr().dot(A.tocsc())).tocsr()
                return invAdagA.dot(Adag.tocsc())
        else:
            if self.is_complete():
                return _inv(self.get_to_std())
            else:
                assert(self.size < self.dim), "Basis seems to be overcomplete: size > dimension!"
                # we'd need to construct a different pseudo-inverse if the above assert fails

                A = self.get_to_std()  # shape (dim,size) - should have indep *cols*
                Adag = A.transpose().conjugate()  # shape (size, dim)
                return _np.dot(_inv(_np.dot(Adag, A)), Adag) 

def inv(self):
        r"""Return the inverse of :math:`\mathbf{A}`.

        Returns
        ----------
        Ainv : :obj:`numpy.ndarray`
            Inverse matrix.

        """
        if isinstance(self.A, np.ndarray):
            Ainv = np.linalg.inv(self.A)
        else:
            Ainv = inv(self.A)

        return Ainv 

def test_inv(self):
        def check(dtype):
            M = array([[1, 0, 2], [0, 0, 3], [-4, 5, 6]], dtype)
            sM = self.spmatrix(M, shape=(3,3), dtype=dtype)
            sMinv = inv(sM)
            assert_array_almost_equal(sMinv.dot(sM).todense(), np.eye(3))
        for dtype in [float]:
            with warnings.catch_warnings():
                warnings.simplefilter("ignore", category=SparseEfficiencyWarning)
                yield check, dtype 

def compute_covariance(self):
        """Compute the covariance matrix after solve has terminated."""
        try:
            # Re-evaluate the precision matrix with the final parameters
            precision, _, _ = self._get_precision_information_and_cost()
            self._covariance_matrix = splinalg.inv(precision.tocsc()).toarray()
        except Exception as e:
            print('Covariance computation failed!\n{}'.format(e)) 

def absorbTime(self):
        P = self.getTransitionMatrix(probabilities=True)
        components,labels = csgraph.connected_components(P, directed=True, connection='strong',return_labels=True)
        if components == 1:
            print("no absorbing states")
            return
            
        transientStates = np.ones(P.shape[0],dtype=bool)

        for component in range(components):
            indices = np.where(labels==component)[0]
            n = len(indices)
            if n==1:
                probSum = P[indices,indices].sum()
            else:            
                probSum = P[np.ix_(indices,indices)].sum()
            if np.isclose(probSum,n):
                transientStates[indices] = False

        indices = np.where(transientStates)[0]
        n = len(indices)
        if n==1:
            Q = P[indices,indices]
        else:
            Q = P[np.ix_(indices,indices)]
        #N will be dense  
        N = inv(eye(n)-Q).A
        N2 = N*(2*N[np.arange(n),np.arange(n)]-np.eye(n))-np.power(N,2)
        t = np.zeros(P.shape[0])
        t[indices] = np.sum(N,axis=1)
        for index in indices:
            print( self.mapping[index],t[index] ) 

def update_Y_closed_form(S_M, Y, Y0, Theta, PiC, gamma, mu):
    # row = []
    # col = []
    # val = []

    for j in range(PiC.shape[1]):
        # for each candidate j, slicing to get submatrix
        mid_list = PiC[:, j].nonzero()[0].tolist()
        Y0_j = Y0[mid_list, :]
        Theta_j = Theta[mid_list, :]
        S_M_j = S_M[mid_list, :][:, mid_list]

        if S_M_j.shape[0] * S_M_j.shape[1] < 2520800000:

            # transform to dense matrix
            tmp = ((1+gamma+mu)*identity(len(mid_list)) - gamma*S_M_j).todense()
            Y_j = npl.inv(tmp) * (Theta_j + mu*Y0_j)
            Y[mid_list, :] = Y_j

            # # sparse 
            # Yc = spsl.inv((1+gamma+mu)*identity(len(mid_list)) - gamma*S_M_j) * (Theta_j + mu*Y0_j)
            # Yc = spsl.spsolve( ((1+gamma+mu)*identity(len(mid_list)) - gamma*S_M_j), (Theta_j + mu*Y0_j) )
            # row_idx, col_idx = Yc.nonzero()
            # for i in range(len(row_idx)):
            #     mid = mid_list[row_idx[i]]
            #     row.append(mid)
            #     col.append(col_idx[i])
            #     val.append(Yc[row_idx[i], col_idx[i]])

        if j % 1000 == 0:
            print 'candidate', j

    # Y = coo_matrix((val, (row, col)), shape = Y0.shape).tocsr()
    return Y 

def test_ptdf(grid):
    """
    Sigma-distances test
    :param grid:
    :return:
    """
    nc = compile_snapshot_circuit(grid)
    islands = split_into_islands(nc)
    inputs = islands[0]  # pick the first island

    PTDF = compute_ptdf(Ybus=inputs.Ybus,
                        Yf=inputs.Yf,
                        Yt=inputs.Yt,
                        Cf=inputs.C_branch_bus_f,
                        Ct=inputs.C_branch_bus_t,
                        V=inputs.Vbus,
                        Ibus=inputs.Ibus,
                        Sbus=inputs.Sbus,
                        pq=inputs.pq,
                        pv=inputs.pv)

    # compose some made up situation
    S2 = inputs.Sbus * 5
    dS = inputs.Sbus - S2
    dSbr = PTDF * dS

    # run a power flow to get the initial branch power and compose the second branch power with the increment
    driver = PowerFlowDriver(grid=grid, options=PowerFlowOptions())
    driver.run()

    Sbr0 = driver.results.Sbranch
    Sbr2 = Sbr0 + dSbr

    PTDFsq = PTDF * PTDF.T
    LODF = PTDFsq * inv(sp.diags(ones(PTDF.shape[0])) - PTDFsq).toarray()

    print('PTDF:')
    print(PTDF.toarray())
    print()
    print('Sbr0')
    print(Sbr0)
    print('Sbr2')
    print(Sbr2) 

def single_short_circuit(self, calculation_inputs: SnapshotCircuit, Vpf, Zf):
        """
        Run a power flow simulation for a single circuit
        @param calculation_inputs:
        @param Vpf: Power flow voltage vector applicable to the island
        @param Zf: Short circuit impedance vector applicable to the island
        @return: short circuit results
        """
        # compute Zbus
        # is dense, so no need to store it as sparse
        if calculation_inputs.Ybus.shape[0] > 1:
            Zbus = inv(calculation_inputs.Ybus).toarray()

            # Compute the short circuit
            V, SCpower = short_circuit_3p(bus_idx=self.options.bus_index,
                                          Zbus=Zbus,
                                          Vbus=Vpf,
                                          Zf=Zf,
                                          baseMVA=calculation_inputs.Sbase)

            # Compute the branches power
            Sbranch, Ibranch, loading, losses = self.compute_branch_results(calculation_inputs=calculation_inputs, V=V)

            # voltage, Sbranch, loading, losses, error, converged, Qpv
            results = ShortCircuitResults(n=calculation_inputs.nbus,
                                          m=calculation_inputs.nbr,
                                          n_tr=calculation_inputs.ntr,
                                          bus_names=calculation_inputs.bus_names,
                                          branch_names=calculation_inputs.branch_names,
                                          transformer_names=calculation_inputs.tr_names,
                                          bus_types=calculation_inputs.bus_types)

            results.Sbus = calculation_inputs.Sbus
            results.voltage = V
            results.Sbranch = Sbranch
            results.Ibranch = Ibranch
            results.losses = losses
            results.SCpower = SCpower


        else:
            nbus = calculation_inputs.Ybus.shape[0]
            nbr = calculation_inputs.nbr

            # voltage, Sbranch, loading, losses, error, converged, Qpv
            results = ShortCircuitResults(n=calculation_inputs.nbus,
                                          m=calculation_inputs.nbr,
                                          n_tr=calculation_inputs.ntr,
                                          bus_names=calculation_inputs.bus_names,
                                          branch_names=calculation_inputs.branch_names,
                                          transformer_names=calculation_inputs.tr_names,
                                          bus_types=calculation_inputs.bus_types)

            results.Sbus = calculation_inputs.Sbus
            results.voltage = np.zeros(nbus, dtype=complex)
            results.Sbranch = np.zeros(nbr, dtype=complex)
            results.Ibranch = np.zeros(nbr, dtype=complex)
            results.losses = np.zeros(nbr, dtype=complex)
            results.SCpower = np.zeros(nbus, dtype=complex)

        return results