code: def generate_pseudo(strain_states, order=3):
      
      s = sp.Symbol('s')
      nstates = len(strain_states)
      ni = np.array(strain_states)*s
      mis, absent_syms = [], []
      for degree in range(2, order + 1):
          cvec, carr = get_symbol_list(degree)
          sarr = np.zeros((nstates, 6), dtype=object)
          for n, strain_v in enumerate(ni):
              # Get expressions
              exps = carr.copy()
              for i in range(degree - 1):
                  exps = np.dot(exps, strain_v)
              exps /= np.math.factorial(degree - 1)
              sarr[n] = [sp.diff(exp, s, degree - 1) for exp in exps]
          svec = sarr.ravel()
          present_syms = set.union(*[exp.atoms(sp.Symbol) for exp in svec])
          absent_syms += [set(cvec) - present_syms]
          m = np.zeros((6*nstates, len(cvec)))
          for n, c in enumerate(cvec):
              m[:, n] = v_diff(svec, c)
          mis.append(np.linalg.pinv(m))
      return mis, absent_syms
docs: Generates the pseudoinverse for a given set of strains.
  
      Args:
          strain_states (6xN array like): a list of voigt-notation
              "strain-states", i. e. perturbed indices of the strain
              as a function of the smallest strain e. g. (0, 1, 0, 0, 1, 0)
          order (int): order of pseudoinverse to calculate
  
      Returns:
          mis: pseudo inverses for each order tensor, these can
              be multiplied by the central difference derivative
              of the stress with respect to the strain state
          absent_syms: symbols of the tensor absent from the PI
              expression