# -*- python -*-

'''

recipes - Some famous numerical recipes implemented in Python.

SYNOPSIS

    >>> import recipes

AUTHOR

    Ed Halley (ed@halley.cc) 24 January 2010

SEE ALSO

    Simple Recipes in Python by William Park
    http://www.phys.uu.nl/~haque/computing/WPark_recipes_in_python.html

'''

#----------------------------------------------------------------------------

from math import sqrt

def num(x):
    '''Return an integer or a float, if it can be cast successfully.
    Returns the original value otherwise.  Usually used on strings
    that may or may not parse to a number.
    '''
    try: return int(x)
    except:
        try: return float(x)
        except: pass
    return x

def sign(x):
    '''Return the sign of a value as -1, 0 or +1.'''
    if x < 0: return -1
    if x > 0: return 1
    return 0

def deviation(X, extra=False):
    '''
    Calculate the standard deviation of data.
    If given extra=True, returns a tuple of (mean, stddev).
    '''
    n = float(len(X))
    mean = sum(X) / n
    std = 0
    for a in X:
        std += (a - mean)**2
    std = sqrt(std / (n-1))
    if extra:
        return (mean, std)
    return std

def linefit(X, Y=None, extra=False):
    '''
    Returns a tuple of (a, b) to the regression line "y=ax+b" from X and Y.
    If given only one list, assumes X is a list of incrementing integers.
    If given extra=True, returns additional parameters of the regression,
    as (a, b, R*R, s*s).
    '''
    if Y is None:
        Y = X
        X = map(float, range(len(Y)))
    if len(X) != len(Y):
        raise ValueError, 'unequal length'

    n = len(X)
    Sx = sum(X)
    Sy = sum(Y)
    Sxx = Syy = Sxy = 0.0
    for x, y in zip(X, Y):
        Sxx += x*x
        Syy += y*y
        Sxy += x*y
    det = Sxx*n - Sx*Sx
    a = (Sxy*n - Sy*Sx) / det
    b = (Sxx*Sy - Sx*Sxy) / det

    if extra:
        meanerror = residual = 0.0
        for x, y in zip(X, Y):
            meanerror += (y - Sy/n)**2
            residual += (y - a*x - b)**2
        RR = 1.0 - residual/meanerror
        ss = residual / (n-2)
        #Var_a = ss*n / det
        #Var_b = ss*Sxx / det
        return (a, b, RR, ss) #, Var_a, Var_b)
    
    return a, b

#----------------------------------------------------------------------------

def __test__():
    from testing import __ok__


if __name__ == '__main__':
    from testing import __report__
    __test__()
    __report__()