text: """
  Functions which are common and require SciPy Base and Level 1 SciPy
  (special, linalg)
  """
  
  from __future__ import division, print_function, absolute_import
  
  import numpy
  import numpy as np
  from numpy import (exp, log, asarray, arange, newaxis, hstack, product, array,
                     zeros, eye, poly1d, r_, sum, fromstring, isfinite,
                     squeeze, amax, reshape)
  
  from scipy._lib._version import NumpyVersion
  
  __all__ = ['logsumexp', 'central_diff_weights', 'derivative', 'pade', 'lena',
             'ascent', 'face']
  
  
  _NUMPY_170 = (NumpyVersion(numpy.__version__) >= NumpyVersion('1.7.0'))
  
  
  def logsumexp(a, axis=None, b=None, keepdims=False):
      """Compute the log of the sum of exponentials of input elements.
  
      Parameters
      ----------
      a : array_like
          Input array.
      axis : None or int or tuple of ints, optional
          Axis or axes over which the sum is taken. By default `axis` is None,
          and all elements are summed. Tuple of ints is not accepted if NumPy
          version is lower than 1.7.0.
  
          .. versionadded:: 0.11.0
      keepdims: bool, optional
          If this is set to True, the axes which are reduced are left in the
          result as dimensions with size one. With this option, the result
          will broadcast correctly against the original array.
  
          .. versionadded:: 0.15.0
      b : array-like, optional
          Scaling factor for exp(`a`) must be of the same shape as `a` or
          broadcastable to `a`.
  
          .. versionadded:: 0.12.0
  
      Returns
      -------
      res : ndarray
          The result, ``np.log(np.sum(np.exp(a)))`` calculated in a numerically
          more stable way. If `b` is given then ``np.log(np.sum(b*np.exp(a)))``
          is returned.
  
      See Also
      --------
      numpy.logaddexp, numpy.logaddexp2
  
      Notes
      -----
      Numpy has a logaddexp function which is very similar to `logsumexp`, but
      only handles two arguments. `logaddexp.reduce` is similar to this
      function, but may be less stable.
  
      Examples
      --------
      >>> from scipy.misc import logsumexp
      >>> a = np.arange(10)
      >>> np.log(np.sum(np.exp(a)))
      9.4586297444267107
      >>> logsumexp(a)
      9.4586297444267107
  
      With weights
  
      >>> a = np.arange(10)
      >>> b = np.arange(10, 0, -1)
      >>> logsumexp(a, b=b)
      9.9170178533034665
      >>> np.log(np.sum(b*np.exp(a)))
      9.9170178533034647
      """
      a = asarray(a)
  
      # keepdims is available in numpy.sum and numpy.amax since NumPy 1.7.0
      #
      # Because SciPy supports versions earlier than 1.7.0, we have to handle
      # those old versions differently
  
      if not _NUMPY_170:
          # When support for Numpy < 1.7.0 is dropped, this implementation can be
          # removed. This implementation is a bit hacky. Similarly to old NumPy's
          # sum and amax functions, 'axis' must be an integer or None, tuples and
          # lists are not supported. Although 'keepdims' is not supported by these
          # old NumPy's functions, this function supports it.
  
          # Solve the shape of the reduced array
          if axis is None:
              sh_keepdims = (1,) * a.ndim
          else:
              sh_keepdims = list(a.shape)
              sh_keepdims[axis] = 1
  
          a_max = amax(a, axis=axis)
  
          if a_max.ndim > 0:
              a_max[~isfinite(a_max)] = 0
          elif not isfinite(a_max):
              a_max = 0
  
          if b is not None:
              b = asarray(b)
              tmp = b * exp(a - reshape(a_max, sh_keepdims))
          else:
              tmp = exp(a - reshape(a_max, sh_keepdims))
  
          # suppress warnings about log of zero
          with np.errstate(divide='ignore'):
              out = log(sum(tmp, axis=axis))
  
          out += a_max
  
          if keepdims:
              # Put back the reduced axes with size one
              out = reshape(out, sh_keepdims)
      else:
          # This is a more elegant implementation, requiring NumPy >= 1.7.0
          a_max = amax(a, axis=axis, keepdims=True)
  
          if a_max.ndim > 0:
              a_max[~isfinite(a_max)] = 0
          elif not isfinite(a_max):
              a_max = 0
  
          if b is not None:
              b = asarray(b)
              tmp = b * exp(a - a_max)
          else:
              tmp = exp(a - a_max)
  
          # suppress warnings about log of zero
          with np.errstate(divide='ignore'):
              out = log(sum(tmp, axis=axis, keepdims=keepdims))
  
          if not keepdims:
              a_max = squeeze(a_max, axis=axis)
  
          out += a_max
  
      return out
  
  
  def central_diff_weights(Np, ndiv=1):
      """
      Return weights for an Np-point central derivative.
  
      Assumes equally-spaced function points.
  
      If weights are in the vector w, then
      derivative is w[0] * f(x-ho*dx) + ... + w[-1] * f(x+h0*dx)
  
      Parameters
      ----------
      Np : int
          Number of points for the central derivative.
      ndiv : int, optional
          Number of divisions.  Default is 1.
  
      Notes
      -----
      Can be inaccurate for large number of points.
  
      """
      if Np < ndiv + 1:
          raise ValueError("Number of points must be at least the derivative order + 1.")
      if Np % 2 == 0:
          raise ValueError("The number of points must be odd.")
      from scipy import linalg
      ho = Np >> 1
      x = arange(-ho,ho+1.0)
      x = x[:,newaxis]
      X = x**0.0
      for k in range(1,Np):
          X = hstack([X,x**k])
      w = product(arange(1,ndiv+1),axis=0)*linalg.inv(X)[ndiv]
      return w
  
  
  def derivative(func, x0, dx=1.0, n=1, args=(), order=3):
      """
      Find the n-th derivative of a function at a point.
  
      Given a function, use a central difference formula with spacing `dx` to
      compute the `n`-th derivative at `x0`.
  
      Parameters
      ----------
      func : function
          Input function.
      x0 : float
          The point at which `n`-th derivative is found.
      dx : int, optional
          Spacing.
      n : int, optional
          Order of the derivative. Default is 1.
      args : tuple, optional
          Arguments
      order : int, optional
          Number of points to use, must be odd.
  
      Notes
      -----
      Decreasing the step size too small can result in round-off error.
  
      Examples
      --------
      >>> def f(x):
      ...     return x**3 + x**2
      ...
      >>> derivative(f, 1.0, dx=1e-6)
      4.9999999999217337
  
      """
      if order < n + 1:
          raise ValueError("'order' (the number of points used to compute the derivative), "
                           "must be at least the derivative order 'n' + 1.")
      if order % 2 == 0:
          raise ValueError("'order' (the number of points used to compute the derivative) "
                           "must be odd.")
      # pre-computed for n=1 and 2 and low-order for speed.
      if n == 1:
          if order == 3:
              weights = array([-1,0,1])/2.0
          elif order == 5:
              weights = array([1,-8,0,8,-1])/12.0
          elif order == 7:
              weights = array([-1,9,-45,0,45,-9,1])/60.0
          elif order == 9:
              weights = array([3,-32,168,-672,0,672,-168,32,-3])/840.0
          else:
              weights = central_diff_weights(order,1)
      elif n == 2:
          if order == 3:
              weights = array([1,-2.0,1])
          elif order == 5:
              weights = array([-1,16,-30,16,-1])/12.0
          elif order == 7:
              weights = array([2,-27,270,-490,270,-27,2])/180.0
          elif order == 9:
              weights = array([-9,128,-1008,8064,-14350,8064,-1008,128,-9])/5040.0
          else:
              weights = central_diff_weights(order,2)
      else:
          weights = central_diff_weights(order, n)
      val = 0.0
      ho = order >> 1
      for k in range(order):
          val += weights[k]*func(x0+(k-ho)*dx,*args)
      return val / product((dx,)*n,axis=0)
  
  
  def pade(an, m):
      """
      Return Pade approximation to a polynomial as the ratio of two polynomials.
  
      Parameters
      ----------
      an : (N,) array_like
          Taylor series coefficients.
      m : int
          The order of the returned approximating polynomials.
  
      Returns
      -------
      p, q : Polynomial class
          The pade approximation of the polynomial defined by `an` is
          `p(x)/q(x)`.
  
      Examples
      --------
      >>> from scipy import misc
      >>> e_exp = [1.0, 1.0, 1.0/2.0, 1.0/6.0, 1.0/24.0, 1.0/120.0]
      >>> p, q = misc.pade(e_exp, 2)
  
      >>> e_exp.reverse()
      >>> e_poly = np.poly1d(e_exp)
  
      Compare ``e_poly(x)`` and the pade approximation ``p(x)/q(x)``
  
      >>> e_poly(1)
      2.7166666666666668
  
      >>> p(1)/q(1)
      2.7179487179487181
  
      """
      from scipy import linalg
      an = asarray(an)
      N = len(an) - 1
      n = N - m
      if n < 0:
          raise ValueError("Order of q <m> must be smaller than len(an)-1.")
      Akj = eye(N+1, n+1)
      Bkj = zeros((N+1, m), 'd')
      for row in range(1, m+1):
          Bkj[row,:row] = -(an[:row])[::-1]
      for row in range(m+1, N+1):
          Bkj[row,:] = -(an[row-m:row])[::-1]
      C = hstack((Akj, Bkj))
      pq = linalg.solve(C, an)
      p = pq[:n+1]
      q = r_[1.0, pq[n+1:]]
      return poly1d(p[::-1]), poly1d(q[::-1])
  
  
  def lena():
      """
      Get classic image processing example image, Lena, at 8-bit grayscale
      bit-depth, 512 x 512 size.
  
      Parameters
      ----------
      None
  
      Returns
      -------
      lena : ndarray
          Lena image
  
      Examples
      --------
      >>> import scipy.misc
      >>> lena = scipy.misc.lena()
      >>> lena.shape
      (512, 512)
      >>> lena.max()
      245
      >>> lena.dtype
      dtype('int32')
  
      >>> import matplotlib.pyplot as plt
      >>> plt.gray()
      >>> plt.imshow(lena)
      >>> plt.show()
  
      """
      import pickle
      import os
      fname = os.path.join(os.path.dirname(__file__),'lena.dat')
      f = open(fname,'rb')
      lena = array(pickle.load(f))
      f.close()
      return lena
  
  
  def ascent():
      """
      Get an 8-bit grayscale bit-depth, 512 x 512 derived image for easy use in demos
  
      The image is derived from accent-to-the-top.jpg at
      http://www.public-domain-image.com/people-public-domain-images-pictures/
  
      Parameters
      ----------
      None
  
      Returns
      -------
      ascent : ndarray
         convenient image to use for testing and demonstration
  
      Examples
      --------
      >>> import scipy.misc
      >>> ascent = scipy.misc.ascent()
      >>> ascent.shape
      (512, 512)
      >>> ascent.max()
      255
  
      >>> import matplotlib.pyplot as plt
      >>> plt.gray()
      >>> plt.imshow(ascent)
      >>> plt.show()
  
      """
      import pickle
      import os
      fname = os.path.join(os.path.dirname(__file__),'ascent.dat')
      with open(fname, 'rb') as f:
          ascent = array(pickle.load(f))
      return ascent
  
  
  def face(gray=False):
      """
      Get a 1024 x 768, color image of a raccoon face.
  
      raccoon-procyon-lotor.jpg at http://www.public-domain-image.com
  
      Parameters
      ----------
      gray : bool, optional
          If True then return color image, otherwise return an 8-bit gray-scale
  
      Returns
      -------
      face : ndarray
          image of a racoon face
  
      Examples
      --------
      >>> import scipy.misc
      >>> face = scipy.misc.face()
      >>> face.shape
      (768, 1024, 3)
      >>> face.max()
      230
      >>> face.dtype
      dtype('uint8')
  
      >>> import matplotlib.pyplot as plt
      >>> plt.gray()
      >>> plt.imshow(face)
      >>> plt.show()
  
      """
      import bz2
      import os
      with open(os.path.join(os.path.dirname(__file__), 'face.dat'), 'rb') as f:
          rawdata = f.read()
      data = bz2.decompress(rawdata)
      face = fromstring(data, dtype='uint8')
      face.shape = (768, 1024, 3)
      if gray is True:
          face = (0.21 * face[:,:,0] + 0.71 * face[:,:,1] + 0.07 * face[:,:,2]).astype('uint8')
      return face
  