text: from __future__ import print_function, division
  
  from itertools import product
  
  from sympy import Tuple, Add, Mul, Matrix, log, expand, Rational
  from sympy.core.trace import Tr
  from sympy.printing.pretty.stringpict import prettyForm
  from sympy.physics.quantum.dagger import Dagger
  from sympy.physics.quantum.operator import HermitianOperator
  from sympy.physics.quantum.represent import represent
  from sympy.physics.quantum.matrixutils import numpy_ndarray, scipy_sparse_matrix, to_numpy
  from sympy.physics.quantum.tensorproduct import TensorProduct, tensor_product_simp
  
  
  class Density(HermitianOperator):
      """Density operator for representing mixed states.
  
      TODO: Density operator support for Qubits
  
      Parameters
      ==========
  
      values : tuples/lists
      Each tuple/list should be of form (state, prob) or [state,prob]
  
      Examples
      ========
  
      Create a density operator with 2 states represented by Kets.
  
      >>> from sympy.physics.quantum.state import Ket
      >>> from sympy.physics.quantum.density import Density
      >>> d = Density([Ket(0), 0.5], [Ket(1),0.5])
      >>> d
      'Density'((|0>, 0.5),(|1>, 0.5))
  
      """
  
      @classmethod
      def _eval_args(cls, args):
          # call this to qsympify the args
          args = super(Density, cls)._eval_args(args)
  
          for arg in args:
              # Check if arg is a tuple
              if not (isinstance(arg, Tuple) and
                      len(arg) == 2):
                  raise ValueError("Each argument should be of form [state,prob]"
                                   " or ( state, prob )")
  
          return args
  
      def states(self):
          """Return list of all states.
  
          Examples
          ========
  
          >>> from sympy.physics.quantum.state import Ket
          >>> from sympy.physics.quantum.density import Density
          >>> d = Density([Ket(0), 0.5], [Ket(1),0.5])
          >>> d.states()
          (|0>, |1>)
  
          """
          return Tuple(*[arg[0] for arg in self.args])
  
      def probs(self):
          """Return list of all probabilities.
  
          Examples
          ========
  
          >>> from sympy.physics.quantum.state import Ket
          >>> from sympy.physics.quantum.density import Density
          >>> d = Density([Ket(0), 0.5], [Ket(1),0.5])
          >>> d.probs()
          (0.5, 0.5)
  
          """
          return Tuple(*[arg[1] for arg in self.args])
  
      def get_state(self, index):
          """Return specific state by index.
  
          Parameters
          ==========
  
          index : index of state to be returned
  
          Examples
          ========
  
          >>> from sympy.physics.quantum.state import Ket
          >>> from sympy.physics.quantum.density import Density
          >>> d = Density([Ket(0), 0.5], [Ket(1),0.5])
          >>> d.states()[1]
          |1>
  
          """
          state = self.args[index][0]
          return state
  
      def get_prob(self, index):
          """Return probability of specific state by index.
  
          Parameters
          ===========
  
          index : index of states whose probability is returned.
  
          Examples
          ========
  
          >>> from sympy.physics.quantum.state import Ket
          >>> from sympy.physics.quantum.density import Density
          >>> d = Density([Ket(0), 0.5], [Ket(1),0.5])
          >>> d.probs()[1]
          0.500000000000000
  
          """
          prob = self.args[index][1]
          return prob
  
      def apply_op(self, op):
          """op will operate on each individual state.
  
          Parameters
          ==========
  
          op : Operator
  
          Examples
          ========
  
          >>> from sympy.physics.quantum.state import Ket
          >>> from sympy.physics.quantum.density import Density
          >>> from sympy.physics.quantum.operator import Operator
          >>> A = Operator('A')
          >>> d = Density([Ket(0), 0.5], [Ket(1),0.5])
          >>> d.apply_op(A)
          'Density'((A*|0>, 0.5),(A*|1>, 0.5))
  
          """
          new_args = [(op*state, prob) for (state, prob) in self.args]
          return Density(*new_args)
  
      def doit(self, **hints):
          """Expand the density operator into an outer product format.
  
          Examples
          ========
  
          >>> from sympy.physics.quantum.state import Ket
          >>> from sympy.physics.quantum.density import Density
          >>> from sympy.physics.quantum.operator import Operator
          >>> A = Operator('A')
          >>> d = Density([Ket(0), 0.5], [Ket(1),0.5])
          >>> d.doit()
          0.5*|0><0| + 0.5*|1><1|
  
          """
  
          terms = []
          for (state, prob) in self.args:
              state = state.expand()  # needed to break up (a+b)*c
              if (isinstance(state, Add)):
                  for arg in product(state.args, repeat=2):
                      terms.append(prob *
                                   self._generate_outer_prod(arg[0], arg[1]))
              else:
                  terms.append(prob *
                               self._generate_outer_prod(state, state))
  
          return Add(*terms)
  
      def _generate_outer_prod(self, arg1, arg2):
          c_part1, nc_part1 = arg1.args_cnc()
          c_part2, nc_part2 = arg2.args_cnc()
  
          if ( len(nc_part1) == 0 or
               len(nc_part2) == 0 ):
              raise ValueError('Atleast one-pair of'
                               ' Non-commutative instance required'
                               ' for outer product.')
  
          # Muls of Tensor Products should be expanded
          # before this function is called
          if (isinstance(nc_part1[0], TensorProduct) and
                  len(nc_part1) == 1 and len(nc_part2) == 1):
              op = tensor_product_simp(nc_part1[0] * Dagger(nc_part2[0]))
          else:
              op = Mul(*nc_part1) * Dagger(Mul(*nc_part2))
  
          return Mul(*c_part1)*Mul(*c_part2)*op
  
      def _represent(self, **options):
          return represent(self.doit(), **options)
  
      def _print_operator_name_latex(self, printer, *args):
          return printer._print(r'\rho', *args)
  
      def _print_operator_name_pretty(self, printer, *args):
          return prettyForm(unichr('\N{GREEK SMALL LETTER RHO}'))
  
      def _eval_trace(self, **kwargs):
          indices = kwargs.get('indices', [])
          return Tr(self.doit(), indices).doit()
  
      def entropy(self):
          """ Compute the entropy of a density matrix.
  
          Refer to density.entropy() method  for examples.
          """
          return entropy(self)
  
  
  def entropy(density):
      """Compute the entropy of a matrix/density object.
  
      This computes -Tr(density*ln(density)) using the eigenvalue decomposition
      of density, which is given as either a Density instance or a matrix
      (numpy.ndarray, sympy.Matrix or scipy.sparse).
  
      Parameters
      ==========
  
      density : density matrix of type Density, sympy matrix,
      scipy.sparse or numpy.ndarray
  
      Examples
      ========
  
      >>> from sympy.physics.quantum.density import Density, entropy
      >>> from sympy.physics.quantum.represent import represent
      >>> from sympy.physics.quantum.matrixutils import scipy_sparse_matrix
      >>> from sympy.physics.quantum.spin import JzKet, Jz
      >>> from sympy import S, log
      >>> up = JzKet(S(1)/2,S(1)/2)
      >>> down = JzKet(S(1)/2,-S(1)/2)
      >>> d = Density((up,0.5),(down,0.5))
      >>> entropy(d)
      log(2)/2
  
      """
      if isinstance(density, Density):
          density = represent(density)  # represent in Matrix
  
      if isinstance(density, scipy_sparse_matrix):
          density = to_numpy(density)
  
      if isinstance(density, Matrix):
          eigvals = density.eigenvals().keys()
          return expand(-sum(e*log(e) for e in eigvals))
      elif isinstance(density, numpy_ndarray):
          import numpy as np
          eigvals = np.linalg.eigvals(density)
          return -np.sum(eigvals*np.log(eigvals))
      else:
          raise ValueError(
              "numpy.ndarray, scipy.sparse or sympy matrix expected")
  
  
  def fidelity(state1, state2):
      """ Computes the fidelity [1]_ between two quantum states
  
      The arguments provided to this function should be a square matrix or a
      Density object. If it is a square matrix, it is assumed to be diagonalizable.
  
      Parameters
      ==========
  
      state1, state2 : a density matrix or Matrix
  
  
      Examples
      ========
  
      >>> from sympy import S, sqrt
      >>> from sympy.physics.quantum.dagger import Dagger
      >>> from sympy.physics.quantum.spin import JzKet
      >>> from sympy.physics.quantum.density import Density, fidelity
      >>> from sympy.physics.quantum.represent import represent
      >>>
      >>> up = JzKet(S(1)/2,S(1)/2)
      >>> down = JzKet(S(1)/2,-S(1)/2)
      >>> amp = 1/sqrt(2)
      >>> updown = (amp * up) + (amp * down)
      >>>
      >>> # represent turns Kets into matrices
      >>> up_dm = represent(up * Dagger(up))
      >>> down_dm = represent(down * Dagger(down))
      >>> updown_dm = represent(updown * Dagger(updown))
      >>>
      >>> fidelity(up_dm, up_dm)
      1
      >>> fidelity(up_dm, down_dm) #orthogonal states
      0
      >>> fidelity(up_dm, updown_dm).evalf().round(3)
      0.707
  
      References
      ==========
  
      .. [1] https://en.wikipedia.org/wiki/Fidelity_of_quantum_states
  
      """
      state1 = represent(state1) if isinstance(state1, Density) else state1
      state2 = represent(state2) if isinstance(state2, Density) else state2
  
      if (not isinstance(state1, Matrix) or
              not isinstance(state2, Matrix)):
          raise ValueError("state1 and state2 must be of type Density or Matrix "
                           "received type=%s for state1 and type=%s for state2" %
                           (type(state1), type(state2)))
  
      if ( state1.shape != state2.shape and state1.is_square):
          raise ValueError("The dimensions of both args should be equal and the "
                           "matrix obtained should be a square matrix")
  
      sqrt_state1 = state1**Rational(1, 2)
      return Tr((sqrt_state1 * state2 * sqrt_state1)**Rational(1, 2)).doit()
  