nodes: 0, module; 1, function_definition; 2, function_name:spd_eig; 3, parameters; 4, block; 5, identifier:W; 6, default_parameter; 7, default_parameter; 8, default_parameter; 9, expression_statement; 10, comment:# check input; 11, assert_statement; 12, if_statement; 13, expression_statement; 14, comment:# sort them; 15, comment:# determine the cutoff. We know that C0 is an spd matrix,; 16, comment:# so we select the truncation threshold such that everything that is negative vanishes; 17, expression_statement; 18, if_statement; 19, comment:# determine effective rank m and perform low-rank approximations.; 20, expression_statement; 21, expression_statement; 22, expression_statement; 23, if_statement; 24, expression_statement; 25, expression_statement; 26, if_statement; 27, return_statement; 28, identifier:epsilon; 29, float:1e-10; 30, identifier:method; 31, string; 32, identifier:canonical_signs; 33, False; 34, comment:""" Rank-reduced eigenvalue decomposition of symmetric positive definite matrix.
  
      Removes all negligible eigenvalues
  
      Parameters
      ----------
      W : ndarray((n, n), dtype=float)
          Symmetric positive-definite (spd) matrix.
      epsilon : float
          Truncation parameter. Eigenvalues with norms smaller than this cutoff will
          be removed.
      method : str
          Method to perform the decomposition of :math:`W` before inverting. Options are:
  
          * 'QR': QR-based robust eigenvalue decomposition of W
          * 'schur': Schur decomposition of W
      canonical_signs : boolean, default = False
          Fix signs in V, s. t. the largest element of in every row of V is positive.
  
      Returns
      -------
      s : ndarray(k)
          k non-negligible eigenvalues, sorted by descending norms
  
      V : ndarray(n, k)
          k leading eigenvectors
      """; 35, call; 36, string; 37, comparison_operator:method.lower() == 'qr'; 38, block; 39, comment:# compute the Eigenvalues of C0 using Schur factorization; 40, elif_clause; 41, else_clause; 42, assignment; 43, assignment; 44, comparison_operator:evmin < 0; 45, block; 46, assignment; 47, assignment; 48, assignment; 49, comparison_operator:m == 0; 50, block; 51, assignment; 52, assignment; 53, identifier:canonical_signs; 54, comment:# enforce canonical eigenvector signs; 55, block; 56, expression_list; 57, string_content:QR; 58, attribute; 59, argument_list; 60, string_content:W is not a symmetric matrix; 61, call; 62, string; 63, import_from_statement; 64, expression_statement; 65, comparison_operator:method.lower() == 'schur'; 66, block; 67, block; 68, pattern_list; 69, call; 70, identifier:evmin; 71, call; 72, identifier:evmin; 73, integer:0; 74, expression_statement; 75, identifier:evnorms; 76, call; 77, identifier:n; 78, subscript; 79, identifier:m; 80, binary_operator:n - _np.searchsorted(evnorms[::-1], epsilon); 81, identifier:m; 82, integer:0; 83, raise_statement; 84, identifier:Vm; 85, subscript; 86, identifier:sm; 87, subscript; 88, for_statement; 89, identifier:sm; 90, identifier:Vm; 91, identifier:_np; 92, identifier:allclose; 93, attribute; 94, identifier:W; 95, attribute; 96, argument_list; 97, string_content:qr; 98, relative_import; 99, dotted_name; 100, assignment; 101, call; 102, string; 103, import_from_statement; 104, expression_statement; 105, expression_statement; 106, raise_statement; 107, identifier:s; 108, identifier:V; 109, identifier:sort_by_norm; 110, argument_list; 111, attribute; 112, argument_list; 113, assignment; 114, attribute; 115, argument_list; 116, call; 117, integer:0; 118, identifier:n; 119, call; 120, call; 121, identifier:V; 122, slice; 123, slice; 124, identifier:s; 125, slice; 126, identifier:j; 127, call; 128, block; 129, identifier:W; 130, identifier:T; 131, identifier:method; 132, identifier:lower; 133, import_prefix; 134, dotted_name; 135, identifier:eig_qr; 136, pattern_list; 137, call; 138, attribute; 139, argument_list; 140, string_content:schur; 141, dotted_name; 142, dotted_name; 143, assignment; 144, assignment; 145, call; 146, identifier:s; 147, identifier:V; 148, identifier:_np; 149, identifier:min; 150, identifier:s; 151, identifier:epsilon; 152, call; 153, identifier:_np; 154, identifier:abs; 155, identifier:s; 156, attribute; 157, argument_list; 158, attribute; 159, argument_list; 160, identifier:_ZeroRankError; 161, argument_list; 162, integer:0; 163, identifier:m; 164, integer:0; 165, identifier:m; 166, identifier:range; 167, argument_list; 168, expression_statement; 169, expression_statement; 170, identifier:eig_qr; 171, identifier:eig_qr; 172, identifier:s; 173, identifier:V; 174, identifier:eig_qr; 175, argument_list; 176, identifier:method; 177, identifier:lower; 178, identifier:scipy; 179, identifier:linalg; 180, identifier:schur; 181, pattern_list; 182, call; 183, identifier:s; 184, call; 185, identifier:ValueError; 186, argument_list; 187, identifier:max; 188, argument_list; 189, identifier:_np; 190, identifier:shape; 191, identifier:evnorms; 192, identifier:_np; 193, identifier:searchsorted; 194, subscript; 195, identifier:epsilon; 196, binary_operator:'All eigenvalues are smaller than %g, rank reduction would discard all dimensions.'%epsilon; 197, identifier:m; 198, assignment; 199, augmented_assignment; 200, identifier:W; 201, identifier:S; 202, identifier:V; 203, identifier:schur; 204, argument_list; 205, attribute; 206, argument_list; 207, binary_operator:'method not implemented: ' + method; 208, identifier:epsilon; 209, binary_operator:-evmin + 1e-16; 210, identifier:evnorms; 211, slice; 212, string; 213, identifier:epsilon; 214, identifier:jj; 215, call; 216, subscript; 217, call; 218, identifier:W; 219, identifier:_np; 220, identifier:diag; 221, identifier:S; 222, string; 223, identifier:method; 224, unary_operator; 225, float:1e-16; 226, unary_operator; 227, string_content:All eigenvalues are smaller than %g, rank reduction would discard all dimensions.; 228, attribute; 229, argument_list; 230, identifier:Vm; 231, slice; 232, identifier:j; 233, attribute; 234, argument_list; 235, string_content:method not implemented:; 236, identifier:evmin; 237, integer:1; 238, identifier:_np; 239, identifier:argmax; 240, call; 241, identifier:_np; 242, identifier:sign; 243, subscript; 244, attribute; 245, argument_list; 246, identifier:Vm; 247, identifier:jj; 248, identifier:j; 249, identifier:_np; 250, identifier:abs; 251, subscript; 252, identifier:Vm; 253, slice; 254, identifier:j
code: def spd_eig(W, epsilon=1e-10, method='QR', canonical_signs=False):
      """ Rank-reduced eigenvalue decomposition of symmetric positive definite matrix.
  
      Removes all negligible eigenvalues
  
      Parameters
      ----------
      W : ndarray((n, n), dtype=float)
          Symmetric positive-definite (spd) matrix.
      epsilon : float
          Truncation parameter. Eigenvalues with norms smaller than this cutoff will
          be removed.
      method : str
          Method to perform the decomposition of :math:`W` before inverting. Options are:
  
          * 'QR': QR-based robust eigenvalue decomposition of W
          * 'schur': Schur decomposition of W
      canonical_signs : boolean, default = False
          Fix signs in V, s. t. the largest element of in every row of V is positive.
  
      Returns
      -------
      s : ndarray(k)
          k non-negligible eigenvalues, sorted by descending norms
  
      V : ndarray(n, k)
          k leading eigenvectors
      """
      # check input
      assert _np.allclose(W.T, W), 'W is not a symmetric matrix'
  
      if method.lower() == 'qr':
          from .eig_qr.eig_qr import eig_qr
          s, V = eig_qr(W)
      # compute the Eigenvalues of C0 using Schur factorization
      elif method.lower() == 'schur':
          from scipy.linalg import schur
          S, V = schur(W)
          s = _np.diag(S)
      else:
          raise ValueError('method not implemented: ' + method)
  
      s, V = sort_by_norm(s, V) # sort them
  
      # determine the cutoff. We know that C0 is an spd matrix,
      # so we select the truncation threshold such that everything that is negative vanishes
      evmin = _np.min(s)
      if evmin < 0:
          epsilon = max(epsilon, -evmin + 1e-16)
  
      # determine effective rank m and perform low-rank approximations.
      evnorms = _np.abs(s)
      n = _np.shape(evnorms)[0]
      m = n - _np.searchsorted(evnorms[::-1], epsilon)
      if m == 0:
          raise _ZeroRankError('All eigenvalues are smaller than %g, rank reduction would discard all dimensions.'%epsilon)
      Vm = V[:, 0:m]
      sm = s[0:m]
  
      if canonical_signs:
          # enforce canonical eigenvector signs
          for j in range(m):
              jj = _np.argmax(_np.abs(Vm[:, j]))
              Vm[:, j] *= _np.sign(Vm[jj, j])
  
      return sm, Vm