text: r"""
  Spike Functions
  
  AUTHORS:
  
  - William Stein (2007-07): initial version
  
  - Karl-Dieter Crisman (2009-09): adding documentation and doctests
  """
  #*****************************************************************************
  #       Copyright (C) 2007 William Stein <wstein@gmail.com>
  #       Copyright (C) 2009 Karl-Dieter Crisman <kcrisman@gmail.com>
  #
  #  Distributed under the terms of the GNU General Public License (GPL)
  #  as published by the Free Software Foundation; either version 2 of
  #  the License, or (at your option) any later version.
  #                  http://www.gnu.org/licenses/
  #*****************************************************************************
  import math
  
  from sage.plot.all import line
  from sage.modules.free_module_element import vector
  from sage.rings.all import RDF
  
  class SpikeFunction:
      """
      Base class for spike functions.
  
      INPUT:
  
          -  ``v`` - list of pairs (x, height)
  
          -  ``eps`` - parameter that determines approximation to a true spike
  
      OUTPUT:
  
          a function with spikes at each point ``x`` in ``v`` with the given height.
  
      EXAMPLES::
  
          sage: spike_function([(-3,4),(-1,1),(2,3)],0.001)
          A spike function with spikes at [-3.0, -1.0, 2.0]
  
      Putting the spikes too close together may delete some::
  
          sage: spike_function([(1,1),(1.01,4)],0.1)
          Some overlapping spikes have been deleted.
          You might want to use a smaller value for eps.
          A spike function with spikes at [1.0]
  
      Note this should normally be used indirectly via
      ``spike_function``, but one can use it directly::
  
          sage: from sage.functions.spike_function import SpikeFunction
          sage: S = SpikeFunction([(0,1),(1,2),(pi,-5)])
          sage: S
          A spike function with spikes at [0.0, 1.0, 3.141592653589793]
          sage: S.support
          [0.0, 1.0, 3.141592653589793]
      """
      def __init__(self, v, eps=0.0000001):
          """
          Initializes base class SpikeFunction.
  
          EXAMPLES::
  
              sage: S = spike_function([(-3,4),(-1,1),(2,3)],0.001); S
              A spike function with spikes at [-3.0, -1.0, 2.0]
              sage: S.height
              [4.0, 1.0, 3.0]
              sage: S.eps
              0.00100000000000000
          """
          if len(v) == 0:
             v = [(0,0)]
          v = sorted([(float(x[0]), float(x[1])) for x in v])
          notify = False
  
          for i in reversed(range(len(v)-1)):
              if v[i+1][0] - v[i][0] <= eps:
                  notify = True
                  del v[i+1]
  
          if notify:
              print "Some overlapping spikes have been deleted."
              print "You might want to use a smaller value for eps."
  
          self.v = v
          self.eps = eps
          self.support = [float(x[0]) for x in self.v]
          self.height = [float(x[1]) for x in self.v]
  
      def __repr__(self):
          """
          String representation of a spike function.
  
          EXAMPLES::
  
              sage: spike_function([(-3,4),(-1,1),(2,3)],0.001)
              A spike function with spikes at [-3.0, -1.0, 2.0]
          """
          return "A spike function with spikes at %s"%self.support
  
      def _eval(self, x):
          """
          Evaluates spike function.  Note that when one calls
          the function within the tolerance, the return value
          is the full height at that point.
  
          EXAMPLES::
  
              sage: S = spike_function([(0,5)],eps=.001)
              sage: S(0)
              5.0
              sage: S(.1)
              0.0
              sage: S(.01)
              0.0
              sage: S(.001)
              5.0
          """
          eps = self.eps
          x = float(x)
          for i in range(len(self.support)):
              z = self.support[i]
              if z - eps <= x and x <= z + eps:
                  return self.height[i], i
          return float(0), -1
  
      def __call__(self, x):
          """
          Called when spike function is used as callable function.
  
          EXAMPLES::
  
              sage: S = spike_function([(0,5)],eps=.001)
              sage: S(0)
              5.0
              sage: S(.1)
              0.0
              sage: S(.01)
              0.0
              sage: S(.001)
              5.0
          """
          return self._eval(x)[0]
  
      def plot_fft_abs(self, samples=2**12, xmin=None, xmax=None,  **kwds):
          """
          Plot of (absolute values of) Fast Fourier Transform of
          the spike function with given number of samples.
  
          EXAMPLES::
  
              sage: S = spike_function([(-3,4),(-1,1),(2,3)]); S
              A spike function with spikes at [-3.0, -1.0, 2.0]
              sage: P = S.plot_fft_abs(8)
              sage: p = P[0]; p.ydata
              [5.0, 5.0, 3.367958691924177, 3.367958691924177, 4.123105625617661, 4.123105625617661, 4.759921664218055, 4.759921664218055]
          """
          w = self.vector(samples = samples, xmin=xmin, xmax=xmax)
          xmin, xmax = self._ranges(xmin, xmax)
          z = w.fft()
          k = vector(RDF, [abs(z[i]) for i in range(int(len(z)/2))])
          return k.plot(xmin=0, xmax=1, **kwds)
  
      def plot_fft_arg(self, samples=2**12, xmin=None, xmax=None,  **kwds):
          """
          Plot of (absolute values of) Fast Fourier Transform of
          the spike function with given number of samples.
  
          EXAMPLES::
  
              sage: S = spike_function([(-3,4),(-1,1),(2,3)]); S
              A spike function with spikes at [-3.0, -1.0, 2.0]
              sage: P = S.plot_fft_arg(8)
              sage: p = P[0]; p.ydata
              [0.0, 0.0, -0.211524990023434..., -0.211524990023434..., 0.244978663126864..., 0.244978663126864..., -0.149106180027477..., -0.149106180027477...]
          """
          w = self.vector(samples = samples, xmin=xmin, xmax=xmax)
          xmin, xmax = self._ranges(xmin, xmax)
          z = w.fft()
          k = vector(RDF, [(z[i]).arg() for i in range(int(len(z)/2))])
          return k.plot(xmin=0, xmax=1, **kwds)
  
      def vector(self, samples=2**16, xmin=None, xmax=None):
          """
          Creates a sampling vector of the spike function in question.
  
          EXAMPLES::
  
              sage: S = spike_function([(-3,4),(-1,1),(2,3)],0.001); S
              A spike function with spikes at [-3.0, -1.0, 2.0]
              sage: S.vector(16)
              (4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
          """
          v = vector(RDF, samples) # creates vector of zeros of length 2^16
          xmin, xmax = self._ranges(xmin, xmax)
          delta = (xmax - xmin)/samples
          w = int(math.ceil(self.eps/delta))
          for i in range(len(self.support)):
              x = self.support[i]
              if x > xmax:
                  break
              h = self.height[i]
              j = int((x - xmin)/delta)
              for k in range(j, min(samples, j+w)):
                  v[k] = h
          return v
  
      def _ranges(self, xmin, xmax):
          """
          Quickly find appropriate plotting interval.
  
          EXAMPLES::
  
              sage: S = spike_function([(-1,1),(1,40)])
              sage: S._ranges(None,None)
              (-1.0, 1.0)
          """
          width = (self.support[-1] + self.support[0])/float(2)
          if xmin is None:
             xmin = self.support[0] - width/float(5)
          if xmax is None:
             xmax = self.support[-1] + width/float(5)
          if xmax <= xmin:
             xmax = xmin + 1
          return xmin, xmax
  
      def plot(self, xmin=None, xmax=None, **kwds):
          """
          Special fast plot method for spike functions.
  
          EXAMPLES::
  
              sage: S = spike_function([(-1,1),(1,40)])
              sage: P = plot(S)
              sage: P[0]
              Line defined by 8 points
          """
          v = []
          xmin, xmax = self._ranges(xmin, xmax)
          x = xmin
          eps = self.eps
          while x < xmax:
              y, i = self._eval(x)
              v.append( (x, y) )
              if i != -1:
                  x0 = self.support[i] + eps
                  v.extend([(x0,y), (x0,0)])
                  if i+1 < len(self.support):
                      x = self.support[i+1] - eps
                      v.append( (x, 0) )
                  else:
                      x = xmax
                      v.append( (xmax, 0) )
              else:
                  new_x = None
                  for j in range(len(self.support)):
                      if self.support[j] - eps > x:
                          new_x = self.support[j] - eps
                          break
                  if new_x is None:
                      new_x = xmax
                  v.append( (new_x, 0) )
                  x = new_x
          L = line(v, **kwds)
          L.xmin(xmin-1); L.xmax(xmax)
          return L
  
  
  
  spike_function = SpikeFunction
  