# Standard Python packages
import math
import numbers

# Special dependencies
import numpy

class InfiniteType:
    def __init__(self, multiplier=1.):
        if multiplier == 0.: raise ZeroDivisionError, "Cannot multiply infinity and zero."
        self._multiplier = multiplier

    def __repr__(self):
        if self is Infinity:
            return "Infinity"
        elif self is MinusInfinity:
            return "-Infinity"
        elif self._multiplier > 0.:
            return "Infinity*%g" % self._multiplier
        else:
            return "-Infinity*%g" % abs(self._multiplier)

    def __neg__(self):
        if self is Infinity:
            return MinusInfinity
        elif self is MinusInfinity:
            return Infinity
        else:
            return self * -1.

    def __mul__(self, number):
        if number == 0.: raise ZeroDivisionError, "Cannot multiply infinity and zero."
        return InfiniteType(self._multiplier * number)

    def __div__(self, number):
        if isinstance(number, InfiniteType): raise ZeroDivisionError, "Cannot divide infinity and infinity."
        if number == 0: raise ZeroDivisionError, "Cannot divide infinity and zero."
        return InfiniteType(self._multiplier / number)

    def __truediv__(self, number):
        return self.__div__(number)

#: Symbol representing infinity; can be multiplied by any scalar.
Infinity = InfiniteType()
MinusInfinity = InfiniteType(-1.)

#: A small number (1e-5), used to avoid numerical round-off issues in
#: comparisons.
#:
#: The following can be used to set epsilon (without any
#: multiple-reference issues)::
#: 
#:    import cassius
#:    cassius.epsilon = 1e-10
epsilon = 1e-5

######################################################### Utility functions

def _roundlevel_nsigfigs(num, n):
  if num == 0.: return 1
  return n - int(math.ceil(math.log10(abs(num))))

def str_round(num, n):
  """Round a number to n digits and return the result as a string."""
  num = round(num, n)
  format = "%."+str(max(n, 0))+"f"
  return format % num

def round_sigfigs(num, n):
  "Round a number to n significant figures."
  return round(num, _roundlevel_nsigfigs(num, n))

def str_sigfigs(num, n):
  """Round a number to n significant figures and return the result as
a string."""
  level = _roundlevel_nsigfigs(num, n)
  num = round(num, level)
  format = "%."+str(max(level, 0))+"f"
  return format % num

def round_errpair(num, err, n=2):
  """Round a number and its uncertainty to n significant figures in
  the uncertainty (default is two)."""

  level = _roundlevel_nsigfigs(err, n)
  return round(num, level), round(err, level)

def str_errpair(num, err, n=2):
  """Round a number and its uncertainty to n significant figures in the
uncertainty (default is two) and return the result as a string."""
  level = _roundlevel_nsigfigs(err, n)
  num = round(num, level)
  err = round(err, level)
  format = "%."+str(max(level, 0))+"f"
  return format % num, format % err

def unicode_errpair(num, err, n=2):
  """Round a number and its uncertainty to n significant figures in the
uncertainty (default is two) and return the result joined by a unicode
plus-minus sign."""
  return u"\u00b1".join(str_errpair(num, err, n))

def mean(*values, **kwds):
    """Compute the mean of N values (N > 0).

    Keyword arguments:
       decimals (int or `None`): number of digits after the decimal
       point to found the result, if not `None`

       sigfigs (int or `None`): number of significant digits to round
       the result, if not `None`

       string (bool): return output as a string (forces number of digits)
    """
    decimals = kwds.get("decimals", None)
    sigfigs = kwds.get("sigfigs", None)
    string = kwds.get("string", False)

    if len(values) == 1 and not isinstance(values[0], (numbers.Number, numpy.number)):
        values = values[0]

    sum_1 = 0.
    sum_y = 0.
    for y in values:
        if not isinstance(y, (numbers.Number, numpy.number)):
            raise ValueError, "mean() requires a list of numbers"

        sum_1 += 1.
        sum_y += y
    if sum_1 != 0.:
        output = sum_y / sum_1
        if decimals is not None:
            if string:
                return str_round(output, decimals)
            else:
                return round(output, decimals)
        elif sigfigs is not None:
            if string:
                return str_sigfigs(output, sigfigs)
            else:
                return round_sigfigs(output, sigfigs)
        else:
            if string:
                return str(output)
            else:
                return output
    else:
        raise ValueError, "Cannot take the mean without any values"

def wmean(values, weights, decimals=None, sigfigs=None, string=False):
    """Compute the weighted mean of N values with N weights (N > 0).

    Keyword arguments:
       decimals (int or `None`): number of digits after the decimal
       point to found the result, if not `None`

       sigfigs (int or `None`): number of significant digits to round
       the result, if not `None`

       string (bool): return output as a string (forces number of digits)
    """
    sum_1 = 0.
    sum_y = 0.
    for y, weight in itertools.izip(values, weights):
        if not isinstance(y, (numbers.Number, numpy.number)) or not isinstance(weight, (numbers.Number, numpy.number)):
            raise ValueError, "wmean() requires lists of numbers"

        sum_1 += weight
        sum_y += weight * y
    if sum_1 != 0.:
        outputval, outputerr = sum_y / sum_1, math.sqrt(1. / sum_1)
        if decimals is not None:
            if string:
                return str_round(outputval, decimals), str_round(outputerr, decimals)
            else:
                return round(outputval, decimals), round(outputerr, decimals)
        elif sigfigs is not None:
            if string:
                return str_errpair(outputval, outputerr, sigfigs)
            else:
                return round_errpair(outputval, outputerr, sigfigs)
        else:
            if string:
                return str(outputval), str(outputerr)
            else:
                return outputval, outputerr
    else:
        raise ValueError, "Cannot take the weighted mean without any values"

def linearfit(xvalues, yvalues, weights=None, decimals=None, sigfigs=None, string=False):
    """Compute a linear fit of N x-y pairs with weights (N > 0).

    Keyword arguments:
       weights (list of numbers or `None`): if `None`, weight all
       points equally.

       decimals (int or `None`): number of digits after the decimal
       point to found the result, if not `None`

       sigfigs (int or `None`): number of significant digits to round
       the result, if not `None`

       string (bool): return output as a string (forces number of digits)
    """
    if weights is None:
        weights = numpy.ones(min(len(xvalues), len(yvalues)), dtype=numpy.float)

    sum_1 = 0.
    sum_x = 0.
    sum_xx = 0.
    sum_y = 0.
    sum_xy = 0.
    for x, y, weight in itertools.izip(xvalues, yvalues, weights):
        if not isinstance(x, (numbers.Number, numpy.number)) or not isinstance(y, (numbers.Number, numpy.number)) or not isinstance(weight, (numbers.Number, numpy.number)):
            raise ValueError, "linearfit() requires lists of numbers"

        sum_1 += weight
        sum_x += weight * x
        sum_xx += weight * x**2
        sum_y += weight * y
        sum_xy += weight * x * y
    delta = (sum_1 * sum_xx) - (sum_x * sum_x)
    if delta != 0.:
        intercept = ((sum_xx * sum_y) - (sum_x * sum_xy)) / delta
        intercept_err = math.sqrt(sum_xx / delta)
        slope = ((sum_1 * sum_xy) - (sum_x * sum_y)) / delta
        slope_err = math.sqrt(sum_1 / delta)
        if decimals is not None:
            if string:
                intercept, intercept_err = str_round(intercept, decimals), str_round(intercept_err, decimals)
                slope, slope_err = str_round(slope, decimals), str_round(slope_err, decimals)
            else:
                intercept, intercept_err = round(intercept, decimals), round(intercept_err, decimals)
                slope, slope_err = round(slope, decimals), round(slope_err, decimals)
        elif sigfigs is not None:
            if string:
                intercept, intercept_err = str_errpair(intercept, intercept_err, sigfigs)
                slope, slope_err = str_errpair(slope, slope_err, sigfigs)
            else:
                intercept, intercept_err = round_errpair(intercept, intercept_err, sigfigs)
                slope, slope_err = round_errpair(slope, slope_err, sigfigs)
        elif string:
            intercept, intercept_err = str(intercept), str(intercept_err)
            slope, slope_err = str(slope), str(slope_err)
        return intercept, intercept_err, slope, slope_err
    else:
        raise ValueError, "Cannot take a linear fit without any values"

def rms(*values, **kwds):
    """Compute the root-mean-square of N values (N > 0).

    Keyword arguments:
       decimals (int or `None`): number of digits after the decimal
       point to found the result, if not `None`

       sigfigs (int or `None`): number of significant digits to round
       the result, if not `None`

       string (bool): return output as a string (forces number of digits)
    """
    decimals = kwds.get("decimals", None)
    sigfigs = kwds.get("sigfigs", None)
    string = kwds.get("string", False)

    if len(values) == 1 and not isinstance(values[0], (numbers.Number, numpy.number)):
        values = values[0]

    sum_1 = 0.
    sum_yy = 0.
    for y in values:
        if not isinstance(y, (numbers.Number, numpy.number)):