code: #!/usr/bin/python
  # -*- coding: utf-8 -*-
  ##
  # utils.py : some auxiliary functions
  ##
  # © 2017, Chris Ferrie (csferrie@gmail.com) and
  #         Christopher Granade (cgranade@cgranade.com).
  #
  # Redistribution and use in source and binary forms, with or without
  # modification, are permitted provided that the following conditions are met:
  #
  #     1. Redistributions of source code must retain the above copyright
  #        notice, this list of conditions and the following disclaimer.
  #
  #     2. Redistributions in binary form must reproduce the above copyright
  #        notice, this list of conditions and the following disclaimer in the
  #        documentation and/or other materials provided with the distribution.
  #
  #     3. Neither the name of the copyright holder nor the names of its
  #        contributors may be used to endorse or promote products derived from
  #        this software without specific prior written permission.
  #
  # THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
  # AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  # IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  # ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
  # LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
  # CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
  # SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
  # INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
  # CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
  # ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
  # POSSIBILITY OF SUCH DAMAGE.
  ##
  
  ## FEATURES ###################################################################
  
  from __future__ import absolute_import
  from __future__ import print_function
  from __future__ import division
  
  ## IMPORTS ####################################################################
  
  from builtins import range
  
  import warnings
  
  import numpy as np
  import numpy.linalg as la
  
  from scipy.linalg import eigh
  
  from scipy.stats import logistic, binom
  from scipy.special import gammaln, gamma, expit, logit
  from scipy.linalg import sqrtm
  
  from numpy.testing import assert_almost_equal
  
  from qinfer._due import due, Doi
  from qinfer._exceptions import ApproximationWarning
  
  ## FUNCTIONS ##################################################################
  
  def get_qutip_module(required_version='3.2'):
      """
      Attempts to return the qutip module, but
      silently returns ``None`` if it can't be
      imported, or doesn't have version at
      least ``required_version``.
  
      :param str required_version: Valid input to
          ``distutils.version.LooseVersion``.
      :return: The qutip module or ``None``.
      :rtype: ``module`` or ``NoneType``
      """
      try:
          import qutip as qt
          from distutils.version import LooseVersion
          _qt_version = LooseVersion(qt.version.version)
          if _qt_version < LooseVersion(required_version):
              return None
      except ImportError:
          return None
  
      return qt
  
  def check_qutip_version(required_version='3.2'):
      """
      Returns ``true`` iff the imported qutip
      version exists and has ``LooseVersion``
      of at least ``required_version``.
  
      :param str required_version: Valid input to
          ``distutils.version.LooseVersion``.
      :rtype: ``bool``
      """
      try:
          qt = get_qutip_module(required_version)
          return qt is not None
      except:
          # In any other case (including something other
          # than ImportError) we say it's not good enough
          return False
  
  
  def binomial_pdf(N,n,p):
      r"""
      Returns the PDF of the binomial distribution
      :math:`\operatorname{Bin}(N, p)` evaluated at :math:`n`.
      """
      return binom(N, p).pmf(n)
  
  def multinomial_pdf(n,p):
      r"""
      Returns the PDF of the multinomial distribution
      :math:`\operatorname{Multinomial}(N, n, p)=
          \frac{N!}{n_1!\cdots n_k!}p_1^{n_1}\cdots p_k^{n_k}`
  
      :param np.ndarray n : Array of outcome integers
          of shape ``(sides, ...)`` where sides is the number of
          sides on the dice and summing over this index indicates
          the number of rolls for the given experiment.
      :param np.ndarray p : Array of (assumed) probabilities
          of shape ``(sides, ...)`` or ``(sides-1,...)``
          with the rest of the dimensions the same as ``n``.
          If ``sides-1``, the last probability is chosen so that the
          probabilities of all sides sums to 1. If ``sides``
          is the last index, these probabilities are assumed
          to sum to 1.
  
      Note that the numbers of experiments don't need to be given because
      they are implicit in the sum over the 0 index of ``n``.
      """
  
      # work in log space to avoid overflow
      log_N_fac = gammaln(np.sum(n, axis=0) + 1)[np.newaxis,...]
      log_n_fac_sum = np.sum(gammaln(n + 1), axis=0)
  
      # since working in log space, we need special
      # consideration at p=0. deal with p=0, n>0 later.
      def nlogp(n,p):
          result = np.zeros(p.shape)
          mask = p!=0
          result[mask] = n[mask] * np.log(p[mask])
          return result
  
      if p.shape[0] == n.shape[0] - 1:
          ep = np.empty(n.shape)
          ep[:p.shape[0],...] = p
          ep[-1,...] = 1-np.sum(p,axis=0)
      else:
          ep = p
      log_p_sum = np.sum(nlogp(n, ep), axis=0)
  
      probs = np.exp(log_N_fac - log_n_fac_sum + log_p_sum)
  
      # if n_k>0 but p_k=0, the whole probability must be 0
      mask = np.sum(np.logical_and(n!=0, ep==0), axis=0) == 0
      probs = mask * probs
  
      return probs[0,...]
  
  def sample_multinomial(N, p, size=None):
      r"""
      Draws fixed number of samples N from different
      multinomial distributions (with the same number dice sides).
  
      :param int N: How many samples to draw from each distribution.
      :param np.ndarray p: Probabilities specifying each distribution.
          Sum along axis 0 should be 1.
      :param size: Output shape. ``int`` or tuple of
          ``int``s. If the given shape is,
          e.g., ``(m, n, k)``, then m * n * k samples are drawn
          for each distribution.
          Default is None, in which case a single value
          is returned for each distribution.
  
      :rtype: np.ndarray
      :return: Array of shape ``(p.shape, size)`` or p.shape if
          size is ``None``.
      """
      # ensure s is array
      s = np.array([1]) if size is None else np.array([size]).flatten()
  
      def take_samples(ps):
          # we have to flatten to make apply_along_axis work.
          return np.random.multinomial(N, ps, np.prod(s)).flatten()
  
      # should have shape (prod(size)*ps.shape[0], ps.shape[1:])
      samples = np.apply_along_axis(take_samples, 0, p)
      # should have shape (size, p.shape)
      samples = samples.reshape(np.concatenate([s, p.shape]))
      # should have shape (p.shape, size)
      samples = samples.transpose(np.concatenate(
          [np.arange(s.ndim, p.ndim+s.ndim), np.arange(s.ndim)]
      ))
  
      if size is None:
          # get rid of trailing singleton dimension.
          samples = samples[...,0]
  
      return samples
  
  
  def outer_product(vec):
      r"""
      Returns the outer product of a vector :math:`v`
      with itself, :math:`v v^\T`.
      """
      return (
          np.dot(vec[:, np.newaxis], vec[np.newaxis, :])
          if len(vec.shape) == 1 else
          np.dot(vec, vec.T)
          )
  
  def particle_meanfn(weights, locations, fn=None):
      r"""
      Returns the mean of a function :math:`f` over model
      parameters.
  
      :param numpy.ndarray weights: Weights of each particle.
      :param numpy.ndarray locations: Locations of each
          particle.
      :param callable fn: Function of model parameters to
          take the mean of. If `None`, the identity function
          is assumed.
      """
      warnings.warn('particle_meanfn is deprecated, please use distributions.ParticleDistribution',
                    DeprecationWarning)
      fn_vals = fn(locations) if fn is not None else locations
      return np.sum(weights * fn_vals.transpose([1, 0]),
          axis=1)
  
  
  def particle_covariance_mtx(weights,locations):
      """
      Returns an estimate of the covariance of a distribution
      represented by a given set of SMC particle.
  
      :param weights: An array containing the weights of each
          particle.
      :param location: An array containing the locations of
          each particle.
      :rtype: :class:`numpy.ndarray`, shape
          ``(n_modelparams, n_modelparams)``.
      :returns: An array containing the estimated covariance matrix.
      """
      # TODO: add shapes to docstring.
      warnings.warn('particle_covariance_mtx is deprecated, please use distributions.ParticleDistribution',
                    DeprecationWarning)
  
      # Find the mean model vector, shape (n_modelparams, ).
      mu = particle_meanfn(weights, locations)
  
      # Transpose the particle locations to have shape
      # (n_modelparams, n_particles).
      xs = locations.transpose([1, 0])
      # Give a shorter name to the particle weights, shape (n_particles, ).
      ws = weights
  
      cov = (
          # This sum is a reduction over the particle index, chosen to be
          # axis=2. Thus, the sum represents an expectation value over the
          # outer product $x . x^T$.
          #
          # All three factors have the particle index as the rightmost
          # index, axis=2. Using the Einstein summation convention (ESC),
          # we can reduce over the particle index easily while leaving
          # the model parameter index to vary between the two factors
          # of xs.
          #
          # This corresponds to evaluating A_{m,n} = w_{i} x_{m,i} x_{n,i}
          # using the ESC, where A_{m,n} is the temporary array created.
          np.einsum('i,mi,ni', ws, xs, xs)
          # We finish by subracting from the above expectation value
          # the outer product $mu . mu^T$.
          - np.dot(mu[..., np.newaxis], mu[np.newaxis, ...])
      )
  
      # The SMC approximation is not guaranteed to produce a
      # positive-semidefinite covariance matrix. If a negative eigenvalue
      # is produced, we should warn the caller of this.
      assert np.all(np.isfinite(cov))
      if not np.all(la.eig(cov)[0] >= 0):
          warnings.warn('Numerical error in covariance estimation causing positive semidefinite violation.', ApproximationWarning)
  
      return cov
  
  
  def ellipsoid_volume(A=None, invA=None):
      """
      Returns the volume of an ellipsoid given either its
      matrix or the inverse of its matrix.
      """
  
      if invA is None and A is None:
          raise ValueError("Must pass either inverse(A) or A.")
  
      if invA is None and A is not None:
          invA = la.inv(A)
  
      # Find the unit sphere volume.
      # http://en.wikipedia.org/wiki/Unit_sphere#General_area_and_volume_formulas
      n  = invA.shape[0]
      Vn = (np.pi ** (n/2)) / gamma(1 + (n/2))
  
      return Vn * la.det(sqrtm(invA))
  
  @due.dcite(
      Doi("10.1016/j.dam.2007.02.013"),
      description="Khachiyan algorithm",
      tags=["implementation"]
  )
  def mvee(points, tol=0.001):
      """
      Returns the minimum-volume enclosing ellipse (MVEE)
      of a set of points, using the Khachiyan algorithm.
      """
  
      # This function is a port of the matlab function by
      # Nima Moshtagh found here:
      # https://www.mathworks.com/matlabcentral/fileexchange/9542-minimum-volume-enclosing-ellipsoid
      # with accompanying writup here:
      # https://www.researchgate.net/profile/Nima_Moshtagh/publication/254980367_MINIMUM_VOLUME_ENCLOSING_ELLIPSOIDS/links/54aab5260cf25c4c472f487a.pdf
  
      N, d = points.shape
  
      Q = np.zeros([N,d+1])
      Q[:,0:d] = points[0:N,0:d]
      Q[:,d] = np.ones([1,N])
  
      Q = np.transpose(Q)
      points = np.transpose(points)
      count = 1
      err = 1
      u = (1/N) * np.ones(shape = (N,))
  
      while err > tol:
  
          X = np.dot(np.dot(Q, np.diag(u)), np.transpose(Q))
          M = np.diag( np.dot(np.dot(np.transpose(Q), la.inv(X)),Q))
          jdx = np.argmax(M)
          step_size = (M[jdx] - d - 1)/((d+1)*(M[jdx] - 1))
          new_u = (1 - step_size)*u
          new_u[jdx] = new_u[jdx] + step_size
          count = count + 1
          err = la.norm(new_u - u)
          u = new_u
  
      U = np.diag(u)
      c = np.dot(points,u)
      A = (1/d) * la.inv(np.dot(np.dot(points,U), np.transpose(points)) - np.outer(c,c) )
      return A, np.transpose(c)
  
  def in_ellipsoid(x, A, c):
      """
      Determines which of the points ``x`` are in the
      closed ellipsoid with shape matrix ``A`` centered at ``c``.
      For a single point ``x``, this is computed as
  
          .. math::
              (c-x)^T\cdot A^{-1}\cdot (c-x) \leq 1
  
      :param np.ndarray x: Shape ``(n_points, dim)`` or ``n_points``.
      :param np.ndarray A: Shape ``(dim, dim)``, positive definite
      :param np.ndarray c: Shape ``(dim)``
      :return: `bool` or array of bools of length ``n_points``
      """
      if x.ndim == 1:
          y = c - x
          return np.einsum('j,jl,l', y, np.linalg.inv(A), y) <= 1
      else:
          y = c[np.newaxis,:] - x
          return np.einsum('ij,jl,il->i', y, np.linalg.inv(A), y) <= 1
  
  def uniquify(seq):
      """
      Returns the unique elements of a sequence ``seq``.
      """
      #from http://stackoverflow.com/a/480227/1205799
      seen = set()
      seen_add = seen.add
      return [ x for x in seq if x not in seen and not seen_add(x)]
  
  def assert_sigfigs_equal(x, y, sigfigs=3):
      """
      Tests if all elements in x and y
      agree up to a certain number of
      significant figures.
  
      :param np.ndarray x: Array of numbers.
      :param np.ndarray y: Array of numbers you want to
          be equal to ``x``.
      :param int sigfigs: How many significant
          figures you demand that they share.
          Default is 3.
      """
      # determine which power of 10 best describes x
      xpow =  np.floor(np.log10(x))
      # now rescale 1 \leq x < 9
      x = x * 10**(- xpow)
      # scale y by the same amount
      y = y * 10**(- xpow)
  
      # now test if abs(x-y) < 0.5 * 10**(-sigfigs)
      assert_almost_equal(x, y, sigfigs)
  
  def format_uncertainty(value, uncertianty, scinotn_break=4):
      """
      Given a value and its uncertianty, format as a LaTeX string
      for pretty-printing.
  
      :param int scinotn_break: How many decimal points to print
          before breaking into scientific notation.
      """
      if uncertianty == 0:
          # Return the exact number, without the ± annotation as a fixed point
          # number, since all digits matter.
          # FIXME: this assumes a precision of 6; need to select that dynamically.
          return "{0:f}".format(value)
      else:
          # Return a string of the form "0.00 \pm 0.01".
          mag_unc = int(np.log10(np.abs(uncertianty)))
          # Zero should be printed as a single digit; that is, as wide as str "1".
          mag_val = int(np.log10(np.abs(value))) if value != 0 else 0
          n_digits = max(mag_val - mag_unc, 0)
  
  
          if abs(mag_val) < abs(mag_unc) and abs(mag_unc) > scinotn_break:
              # We're formatting something close to zero, so recale uncertianty
              # accordingly.
              scale = 10**mag_unc
              return r"({{0:0.{0}f}} \pm {{1:0.{0}f}}) \times 10^{{2}}".format(
                  n_digits
              ).format(
                  value / scale,
                  uncertianty / scale,
                  mag_unc
             )
          if abs(mag_val) <= scinotn_break:
              return r"{{0:0.{n_digits}f}} \pm {{1:0.{n_digits}f}}".format(n_digits=n_digits).format(value, uncertianty)
          else:
              scale = 10**mag_val
              return r"({{0:0.{0}f}} \pm {{1:0.{0}f}}) \times 10^{{2}}".format(
                  n_digits
              ).format(
                  value / scale,
                  uncertianty / scale,
                  mag_val
             )
  
  def compactspace(scale, n):
      r"""
      Returns points :math:`x` spaced in the open interval
      :math:`(-\infty, \infty)`  by linearly spacing in the compactified
      coordinate :math:`s(x) = e^{-\alpha x} / (1 + e^{-\alpha x})^2`,
      where :math:`\alpha` is a scale factor.
      """
      logit = logistic(scale=scale).ppf
      compact_xs = np.linspace(0, 1, n + 2)[1:-1]
      return logit(compact_xs)
  
  def to_simplex(y):
      r"""
      Interprets the last index of ``y`` as stick breaking fractions 
      in logit space and returns a non-negative array of 
      the same shape where the last dimension always sums to unity.
      
      A unit simplex is a list of non-negative numbers :math:`(x_1,...,x_K)`
      that sum to one, :math:`\sum_{k=1}^K x_k=1`, for example, the 
      probabilities of an K-sided die.
      It is sometimes desireable to parameterize this object with variables 
      that are unconstrained and "decorrelated".
      To this end, we imagine :math:`\vec{x}` as a partition of the unit 
      stick :math:`[0,1]` with :math:`K-1` break points between 
      :math:`K` successive intervals of respective lengths :math:`(x_1,...,x_K)`.
      Instead of storing the interval lengths, we start from the left-most break 
      point and iteratively store the breaking fractions, :math:`z_k`, 
      of the remaining stick.
      This gives the formula 
      :math:`z_k=x_k / (1-\sum_{k'=1}^{k-1}x_k)` with the convention 
      :math:`x_0:=0`, 
      which has an inverse formula :math:`x_k = z_k(1-z_{k-1})\cdots(1-z_1)`.
      Note that :math:`z_K=1` since the last stick is not broken; this is the 
      result of the redundant information imposed by :math:`\sum_{k=1}^K x_k=1`.
      To unbound the parameters :math:`z_k` into the real line, 
      we pass through the logit function, 
      :math:`\operatorname{logit}(p)=\log\frac{p}{1-p}`, 
      to end up with the parameterization 
      :math:`y_k=\operatorname{logit}(z_k)+\log(K-k)`, with the convention 
      :math:`y_K=0`.
      The shift by :math:`\log(K-k)` is largely asthetic and causes the 
      uniform simplex :math:`\vec{x}=(1/K,1/K,...,1/K)` to be mapped to 
      :math:`\vec{x}=(0,0,...,0)`.
  
      Inverse to :func:`from_simplex`.
  
      :param np.ndarray: Array of logit space stick breaking 
          fractions along the last index.
  
      :rtype: ``np.ndarray``
      """
      n = y.shape[-1]
      # z are the stick breaking fractions in [0,1]
      z = expit(y - np.log(n - np.arange(1, n+1)))
      x = np.empty(y.shape)
      x[..., 0] = z[..., 0]
      x[..., 1:] = z[..., 1:] * (1 - z[..., :-1]).cumprod(axis=-1)
      return x
  
  def from_simplex(x):
      r"""
      Inteprets the last index of x as unit simplices and returns a
      real array of the sampe shape in logit space.
  
      Inverse to :func:`to_simplex` ; see that function for more details.
  
      :param np.ndarray: Array of unit simplices along the last index.
      
      :rtype: ``np.ndarray``
      """
      n = x.shape[-1]
      # z are the stick breaking fractions in [0,1]
      # the last one is always 1, so don't worry about it
      z = np.empty(shape=x.shape)
      z[..., 0] = x[..., 0]
      z[..., 1:-1] = x[..., 1:-1] / (1 - x[..., :-2].cumsum(axis=-1))
  
      # now z are the logit-transformed breaking fractions
      z[..., :-1] = logit(z[..., :-1]) - logit(1 / (n - np.arange(n-1, dtype=np.float)))
      # set this to 0 manually to avoid subtracting inf-inf
      z[..., -1] = 0
      return z
  
  def pretty_time(secs, force_h=False, force_m=False):
      if secs > 86400:
          return "{d} days, ".format(d=int(secs//86400)) + pretty_time(secs % 86400, force_h=True)
      elif force_h or secs > 3600:
          return "{h}:".format(h=int(secs//3600)) + pretty_time(secs % 3600, force_m=True)
      elif force_m or secs > 60:
          return (
              "{m:0>2}:{s:0>2}" if force_m else "{m}:{s:0>2}"
          ).format(m=int(secs//60), s=int(secs%60))
      else:
          return "{0:0.2f} seconds".format(secs)
  
  def safe_shape(arr, idx=0, default=1):
      shape = np.shape(arr)
      return shape[idx] if idx < len(shape) else default
  
  def join_struct_arrays(arrays):
      """
      Takes a list of possibly structured arrays, concatenates their
      dtypes, and returns one big array with that dtype. Does the
      inverse of ``separate_struct_array``.
  
      :param list arrays: List of ``np.ndarray``s
      """
      # taken from http://stackoverflow.com/questions/5355744/numpy-joining-structured-arrays
      sizes = np.array([a.itemsize for a in arrays])
      offsets = np.r_[0, sizes.cumsum()]
      shape = arrays[0].shape
      joint = np.empty(shape + (offsets[-1],), dtype=np.uint8)
      for a, size, offset in zip(arrays, sizes, offsets):
          joint[...,offset:offset+size] = np.atleast_1d(a).view(np.uint8).reshape(shape + (size,))
      dtype = sum((a.dtype.descr for a in arrays), [])
      return joint.ravel().view(dtype)
  
  def separate_struct_array(array, dtypes):
      """
      Takes an array with a structured dtype, and separates it out into
      a list of arrays with dtypes coming from the input ``dtypes``.
      Does the inverse of ``join_struct_arrays``.
  
      :param np.ndarray array: Structured array.
      :param dtypes: List of ``np.dtype``, or just a ``np.dtype`` and the number of
          them is figured out automatically by counting bytes.
      """
      try:
          offsets = np.cumsum([np.dtype(dtype).itemsize for dtype in dtypes])
      except TypeError:
          dtype_size = np.dtype(dtypes).itemsize
          num_fields = int(array.nbytes / (array.size * dtype_size))
          offsets = np.cumsum([dtype_size] * num_fields)
          dtypes = [dtypes] * num_fields
      offsets = np.concatenate([[0], offsets]).astype(int)
      uint_array = array.view(np.uint8).reshape(array.shape + (-1,))
      return [
          uint_array[..., offsets[idx]:offsets[idx+1]].flatten().view(dtype)
          for idx, dtype in enumerate(dtypes)
      ]
  
  def sqrtm_psd(A, est_error=True, check_finite=True):
      """
      Returns the matrix square root of a positive semidefinite matrix,
      truncating negative eigenvalues.
      """
      w, v = eigh(A, check_finite=check_finite)
      mask = w <= 0
      w[mask] = 0
      np.sqrt(w, out=w)
      A_sqrt = (v * w).dot(v.conj().T)
  
      if est_error:
          return A_sqrt, np.linalg.norm(np.dot(A_sqrt, A_sqrt) - A, 'fro')
      else:
          return A_sqrt
  
  def decorate_init(init_decorator):
      """
      Given a class definition and a decorator that acts on methods,
      applies that decorator to the class' __init__ method.
      Useful for decorating __init__ while still allowing __init__ to be
      inherited.
      """
  
      def class_decorator(cls):
          cls.__init__ = init_decorator(cls.__init__)
          return cls
  
      return class_decorator
  
  #==============================================================================
  #Test Code
  if __name__ == "__main__":
  
      from mpl_toolkits.mplot3d import Axes3D
      from mpl_toolkits.mplot3d.art3d import Poly3DCollection
      import matplotlib.pyplot as plt
      from scipy.spatial import Delaunay
  
      #some random points
      points = np.array([[ 0.53135758, -0.25818091, -0.32382715],
      [ 0.58368177, -0.3286576,  -0.23854156,],
      [ 0.18741533,  0.03066228, -0.94294771],
      [ 0.65685862, -0.09220681, -0.60347573],
      [ 0.63137604, -0.22978685, -0.27479238],
      [ 0.59683195, -0.15111101, -0.40536606],
      [ 0.68646128,  0.0046802,  -0.68407367],
      [ 0.62311759,  0.0101013,  -0.75863324]])
  
      # compute mvee
      A, centroid = mvee(points)
      print(A)
  
      # point it and some other stuff
      U, D, V = la.svd(A)
  
      rx, ry, rz = [1/np.sqrt(d) for d in D]
      u, v = np.mgrid[0:2*np.pi:20j,-np.pi/2:np.pi/2:10j]
  
      x=rx*np.cos(u)*np.cos(v)
      y=ry*np.sin(u)*np.cos(v)
      z=rz*np.sin(v)
  
      for idx in range(x.shape[0]):
          for idy in range(y.shape[1]):
              x[idx,idy],y[idx,idy],z[idx,idy] = np.dot(np.transpose(V),np.array([x[idx,idy],y[idx,idy],z[idx,idy]])) + centroid
  
  
      fig = plt.figure()
      ax = fig.add_subplot(111, projection='3d')
      ax.scatter(points[:,0],points[:,1],points[:,2])
      ax.plot_surface(x, y, z, cstride = 1, rstride = 1, alpha = 0.1)
      plt.show()
  
  def binom_est_p(n, N, hedge=float(0)):
      r"""
      Given a number of successes :math:`n` and a number of trials :math:`N`,
      estimates the binomial distribution parameter :math:`p` using the
      hedged maximum likelihood estimator of [FB12]_.
  
      :param n: Number of successes.
      :type n: `numpy.ndarray` or `int`
      :param int N: Number of trials.
      :param float hedge: Hedging parameter :math:`\beta`.
      :rtype: `float` or `numpy.ndarray`.
      :return: The estimated binomial distribution parameter :math:`p` for each
          value of :math:`n`.
      """
      return (n + hedge) / (N + 2 * hedge)
  
  def binom_est_error(p, N, hedge = float(0)):
      r"""
      """
  
      # asymptotic np.sqrt(p * (1 - p) / N)
      return np.sqrt(p*(1-p)/(N+2*hedge+1))
  
apis: numpy.ones numpy.sum numpy.diag scipy.special.logit numpy.dtype scipy.special.gamma numpy.log numpy.apply_along_axis numpy.isfinite numpy.testing.assert_almost_equal numpy.transpose scipy.stats.binom matplotlib.pyplot.figure numpy.logical_and numpy.abs numpy.cos scipy.linalg.eigh numpy.log10 numpy.linspace numpy.linalg.eig scipy.special.gammaln numpy.zeros numpy.dot numpy.argmax numpy.arange numpy.prod numpy.einsum numpy.array numpy.linalg.norm numpy.cumsum numpy.empty scipy.stats.logistic numpy.linalg.inv numpy.exp numpy.linalg.svd numpy.atleast_1d matplotlib.pyplot.show numpy.shape numpy.sqrt numpy.sin numpy.concatenate numpy.outer scipy.linalg.sqrtm