text: from __future__ import division
  from random import shuffle, random, choice
  import numpy
  
  try:
      from math import factorial
  except ImportError: # python version < 2.6
      from cogent.maths.stats.special import Gamma
      factorial = lambda x: Gamma(x+1)
  
  from cogent.maths.stats.special import igam
  
  __author__ = "Hua Ying, Julien Epps and Gavin Huttley"
  __copyright__ = "Copyright 2007-2012, The Cogent Project"
  __credits__ = ["Julien Epps", "Hua Ying", "Gavin Huttley"]
  __license__ = "GPL"
  __version__ = "1.5.3"
  __maintainer__ = "Gavin Huttley"
  __email__ = "Gavin.Huttley@anu.edu.au"
  __status__ = "Production"
  
  def chi_square(x, p, df=1):
      """returns the chisquare statistic and it's probability"""
      N = len(x)
      end = N
      sim = numpy.logical_not(numpy.logical_xor(x[0:end-p], x[p:end]))*1
      s = ((numpy.ones((N-p,), float)-sim)**2).sum()
      D = s/(N-p)
      p_val = 1 - igam(df/2.0, D/2)
      return D, p_val
  
  def g_statistic(X, p=None, idx=None):
      """
      return g statistic and p value
      arguments: 
          X - the periodicity profile (e.g. DFT magnitudes, autocorrelation etc)
              X needs to contain only those period values being considered, 
              i.e. only periods in the range [llim, ulim]
      """
      # X should be real
      X = abs(numpy.array(X))
      if p is None: 
          power = X.max(0)
          idx = X.argmax(0)
      else:
          assert idx is not None
          power = X[idx]
      g_obs = power/X.sum()
      M = numpy.floor(1/g_obs)
      pmax = len(X)
      
      result = numpy.zeros((int(M+1),), float)
      pmax_fact = factorial(pmax)
      for index in xrange(1, min(pmax, int(M))+1):
          v = (-1)**(index-1)*pmax_fact/factorial(pmax-index)/factorial(index)
          v *= (1-index*g_obs)**(pmax-1)
          result[index] = v
      
      p_val = result.sum()
      return g_obs, p_val
  
  def _seq_to_symbols(seq, motifs, motif_length, result=None):
      """return symbolic represetation of the sequence
      Arguments: 
          - seq: a sequence
          - motifs: a list of sequence motifs
          - motif_length: length of first motif
      """
      if result is None:
          result = numpy.zeros(len(seq), numpy.uint8)
      else:
          result.fill(0)
      
      if motif_length is None:
          motif_length = len(motifs[0])
      
      for i in xrange(len(seq) - motif_length + 1):
          if seq[i: i + motif_length] in motifs:
              result[i] = 1
      
      return result
  
  try:
      from cogent.maths._period import seq_to_symbols
      #raise ImportError
  except ImportError:
      seq_to_symbols = _seq_to_symbols
  
  class SeqToSymbols(object):
      """class for converting all occurrences of motifs in passed sequence
      to 1/0 otherwise"""
      def __init__(self, motifs, length=None, motif_length=None):
          super(SeqToSymbols, self).__init__()
          if type(motifs) == str:
              motifs = [motifs]
          self.motifs = motifs
          self.length = length
          self.motif_length = motif_length or len(motifs[0])
          self.working = None
          if length is not None:
              self.setResultArray(length)
      
      def setResultArray(self, length):
          """sets a result array for length"""
          self.working = numpy.zeros(length, numpy.uint8)
          self.length = length
      
      def __call__(self, seq, result=None):
          if result is None and self.working is None:
              self.setResultArray(len(seq))
          elif self.working is not None:
              if len(seq) != self.working.shape[0]:
                  self.setResultArray(len(seq))
          
          result = self.working
          result.fill(0)
          if type(seq) != str:
              seq = ''.join(seq)
          
          return seq_to_symbols(seq, self.motifs, self.motif_length, result)
      
  
  def circular_indices(vector, start, length, num):
      """docstring for circular_indices"""
      if start > length:
          start = start-length
          
      if start+num < length:
          return vector[start: start+num]
      # get all till end, then from beginning
      return vector[start:] + vector[:start+num-length]
  
  def sampled_places(block_size, length):
      """returns randomly sampled positions with block_size to make a new vector
      with length
      """
      # Main condition is to identify when a draw would run off end, we want to
      # draw from beginning
      num_seg, remainder = divmod(length, block_size)
      vector = range(length)
      result = []
      for seg_num in xrange(num_seg):
          i = choice(vector)
          result += circular_indices(vector, i, length, block_size)
      
      if remainder:
          result += circular_indices(vector, i+block_size, length, remainder)
      
      assert len(result) == length, len(result)
      return result
  
  def blockwise_bootstrap(signal, calc, block_size, num_reps, seq_to_symbols=None, num_stats=None):
      """returns observed statistic and the probability from the bootstrap
      test of observing more `power' by chance than that estimated from the
      observed signal
      
      Arguments:
          - signal: a series, can be a sequence object
          - calc: function to calculate the period power, e.g. ipdft, hybrid,
            auto_corr or any other statistic.
          - block_size: size of contiguous values for resampling
          - num_reps: number of randomly generated permutations
          - seq_to_symbols: function to convert a sequence to 1/0. If not
            provided, the raw data is used.
          - num_stats: the number of statistics being evaluated for each
            interation. Default to 1.
      """
      signal_length = len(signal)
      
      if seq_to_symbols is not None:
          dtype='c'
      else:
          dtype=None # let numpy guess
      
      signal = numpy.array(list(signal), dtype=dtype)
      
      if seq_to_symbols is not None:
          symbolic = seq_to_symbols(signal)
          data = symbolic
      else:
          data = signal
      
      obs_stat = calc(data)
      if seq_to_symbols is not None:
          if sum(symbolic) == 0:
              p = [numpy.array([1.0, 1.0, 1.0]), 1.0][num_stats == 1]
              
              return obs_stat, p
      
      if num_stats is None:
          try:
              num_stats = calc.getNumStats()
          except AttributeError:
              num_stats = 1
      
      if num_stats == 1:
          count = 0
      else:
          count = numpy.zeros(num_stats)
      
      for rep in range(num_reps):
          # get sample positions
          sampled_indices = sampled_places(block_size, signal_length)
          new_signal = signal.take(sampled_indices)
          if seq_to_symbols is not None:
              symbolic = seq_to_symbols(new_signal)
              data = symbolic
          else:
              data = new_signal
          sim_stat = calc(data)
          # count if > than observed
          if num_stats > 1:
              count[sim_stat >= obs_stat] += 1
          elif sim_stat >= obs_stat:
              count += 1
          
      return obs_stat, count / num_reps
  
  
  # def percrb4():
  #     """Return SNR and CRB for periodicity estimates from symbolic signals"""
  #     # TODO: complete the function according to Julien's percrb4.m
  #     pass
  # 
  