Python sympy.prod() Examples

def _random_coprime_pair(entropy):
  """Returns a pair of random coprime integers."""
  coprime_product = number.integer(entropy, False, min_abs=1)
  factors = sympy.factorint(coprime_product)
  def take():
    prime = random.choice(list(factors.keys()))
    power = factors[prime]
    del factors[prime]
    return prime ** power

  if random.random() < 0.8 and len(factors) >= 2:
    # Disallow trivial factoring where possible.
    count_left = random.randint(1, len(factors) - 1)
    count_right = len(factors) - count_left
  else:
    count_left = random.randint(0, len(factors))
    count_right = len(factors) - count_left

  left = sympy.prod([take() for _ in range(count_left)])
  right = sympy.prod([take() for _ in range(count_right)])
  assert left * right == coprime_product
  return left, right


# @composition.module(number.is_positive_integer) 

def test_boson_states():
    a = BosonOp("a")

    # Fock states
    n = 3
    assert (BosonFockBra(0) * BosonFockKet(1)).doit() == 0
    assert (BosonFockBra(1) * BosonFockKet(1)).doit() == 1
    assert qapply(BosonFockBra(n) * Dagger(a)**n * BosonFockKet(0)) \
        == sqrt(prod(range(1, n+1)))

    # Coherent states
    alpha1, alpha2 = 1.2, 4.3
    assert (BosonCoherentBra(alpha1) * BosonCoherentKet(alpha1)).doit() == 1
    assert (BosonCoherentBra(alpha2) * BosonCoherentKet(alpha2)).doit() == 1
    assert abs((BosonCoherentBra(alpha1) * BosonCoherentKet(alpha2)).doit() -
               exp(-S(1) / 2 * (alpha1 - alpha2) ** 2)) < 1e-12
    assert qapply(a * BosonCoherentKet(alpha1)) == \
        alpha1 * BosonCoherentKet(alpha1) 

def newton_cotes_open(index, **kwargs):
    """
    Open Newton-Cotes formulae.
    <https://math.stackexchange.com/a/1959071/36678>
    """
    points = numpy.linspace(-1.0, 1.0, index + 2)[1:-1]
    degree = index if (index + 1) % 2 == 0 else index - 1
    #
    n = index + 1
    weights = numpy.empty(n - 1)
    t = sympy.Symbol("t")
    for r in range(1, n):
        # Compare with get_weights().
        f = sympy.prod([(t - i) for i in range(1, n) if i != r])
        alpha = (
            2
            * (-1) ** (n - r + 1)
            * sympy.integrate(f, (t, 0, n), **kwargs)
            / (math.factorial(r - 1) * math.factorial(n - 1 - r))
            / n
        )
        weights[r - 1] = alpha
    return C1Scheme("Newton-Cotes (open)", degree, weights, points) 

def probability(self, event):
    # Specializations for optimization.
    if isinstance(event, FiniteProductEvent):
      assert len(self._spaces) == len(event.events)
      return sympy.prod([
          space.probability(event_slice)
          for space, event_slice in zip(self._spaces, event.events)])

    if isinstance(event, CountLevelSetEvent) and self.all_spaces_equal():
      space = self._spaces[0]
      counts = event.counts
      probabilities = {
          value: space.probability(DiscreteEvent({value}))
          for value in six.iterkeys(counts)
      }

      num_events = sum(six.itervalues(counts))
      assert num_events == len(self._spaces)
      # Multinomial coefficient:
      coeff = (
          sympy.factorial(num_events) / sympy.prod(
              [sympy.factorial(i) for i in six.itervalues(counts)]))
      return coeff * sympy.prod([
          pow(probabilities[value], counts[value])
          for value in six.iterkeys(counts)
      ])

    raise ValueError('Unhandled event type {}'.format(type(event))) 

def __init__(self, dims):
    """Constructs a shape whose i-th dim is dims[i].

    Each dim can be one of the following types:
      integer: represents the dimension is a known and fixed.
      string: represents the dimension is an unknown and a sympy dummy symbol is
        used to represent it. Also note that contents of strings only matter for
        logging/printing. Even if the same string is given on multiple
        dimensions, it doesn't mean that they are the same.
      sympy expression: represents a dimension which possibly
        depends on dimensions of other shapes.

    Args:
      dims: A list of either integer, string or sympy.Symbol.
    """
    self._shape = []
    for x in dims:
      assert x is not None, str(dims)
      if isinstance(x, six.string_types):
        # NOTE: Dummy(x) creates a unique symbol. I.e., the value of x has no
        # meaning except for printing, etc.
        self._shape.append(sympy.Dummy(x, integer=True))
      else:
        # Converts x to a sympy type. E.g., int to sympy.Integer.
        self._shape.append(sympy.sympify(x))
    self._size = sympy.prod(self._shape) 

def newton_cotes_closed(index, **kwargs):
    """
    Closed Newton-Cotes formulae.
    <https://en.wikipedia.org/wiki/Newton%E2%80%93Cotes_formulas#Closed_Newton.E2.80.93Cotes_formulae>,
    <http://mathworld.wolfram.com/Newton-CotesFormulas.html>.
    """
    points = numpy.linspace(-1.0, 1.0, index + 1)
    degree = index + 1 if index % 2 == 0 else index

    # Formula (26) from
    # <http://mathworld.wolfram.com/Newton-CotesFormulas.html>.
    # Note that Sympy carries out all operations in rationals, i.e.,
    # _exactly_. Only at the end, the rational is converted into a float.
    n = index
    weights = numpy.empty(n + 1)
    t = sympy.Symbol("t")
    for r in range(n + 1):
        # Compare with get_weights().
        f = sympy.prod([(t - i) for i in range(n + 1) if i != r])
        alpha = (
            2
            * (-1) ** (n - r)
            * sympy.integrate(f, (t, 0, n), **kwargs)
            / (math.factorial(r) * math.factorial(n - r))
            / index
        )
        weights[r] = alpha
    return C1Scheme("Newton-Cotes (closed)", degree, weights, points) 

def symbolic_coeffs_and_diffs(expr,u):
  ''' 
  returns the coefficients for each term containing u or a derivative 
  of u. Also returns the variables that derivatives of u are with 
  respect to
  '''
  # convert expr to a list of terms
  expr = expr.expand()
  expr = expr.as_ordered_terms()
  # throw out terms not containing u
  expr = [i for i in expr if i.has(u)]
  coeffs = []
  diffs = []
  for e in expr:
    # if the expression is a product then expand it into multipliers
    if e.is_Mul:
      e = sp.flatten(e.as_coeff_mul())
    else:
      e = [sp.Integer(1),e]  

    # find multipliers without the queried term
    without_u = [i for i in e if not i.has(u)] 
    coeffs += [without_u]

    # find multipliers with the queried term
    with_u = [i for i in e if i.has(u)]
    if not (len(with_u) == 1):
      raise FormulationError(
        'the term %s has multiple occurrences of %s' % (sp.prod(e),u))

    base,diff = derivative_order(with_u[0])
    if not (base == u):
      raise FormulationError( 
        'cannot express %s as a differential operation of %s' % (base,u))
      
    diffs += diff,

  return coeffs,diffs 

def volume_cell(self):
        """Volume of a single cell e.g  h_x*h_y*h_z in 3D."""
        return prod(d.spacing for d in self.dimensions).subs(self.spacing_map)