source: inundation-numpy-branch/utilities/numerical_tools.py @ 3236

#!/usr/bin/env python
"""Auxiliary numerical tools

"""

from numpy import ArrayType, array, sum, innerproduct, ravel, sqrt
from numpy import searchsorted, sort, concatenate    
    

def angle(v):
    """Compute angle between e1 (the unit vector in the x-direction)
    and the specified vector
    """

    from math import acos, pi, sqrt

    l = sqrt( sum (array(v)**2))
    v1 = v[0]/l
    v2 = v[1]/l


    theta = acos(v1)
    
    #try:
    #   theta = acos(v1)
    #except ValueError, e:
    #    print 'WARNING (util.py): Angle acos(%s) failed: %s' %(str(v1), e)
    #
    #    #FIXME, hack to avoid acos(1.0) Value error
    #    # why is it happening?
    #    # is it catching something we should avoid?
    #    #FIXME (Ole): When did this happen? We need a unit test to
    #    #reveal this condition or otherwise remove the hack
    #    
    #    s = 1e-6
    #    if (v1+s >  1.0) and (v1-s < 1.0) :
    #        theta = 0.0
    #    elif (v1+s >  -1.0) and (v1-s < -1.0):
    #        theta = 3.1415926535897931
    #    print 'WARNING (util.py): angle v1 is %f, setting acos to %f '\
    #          %(v1, theta)

    if v2 < 0:
        #Quadrant 3 or 4
        theta = 2*pi-theta

    return theta


def anglediff(v0, v1):
    """Compute difference between angle of vector x0, y0 and x1, y1.
    This is used for determining the ordering of vertices,
    e.g. for checking if they are counter clockwise.

    Always return a positive value
    """

    from math import pi

    a0 = angle(v0)
    a1 = angle(v1)

    #Ensure that difference will be positive
    if a0 < a1:
        a0 += 2*pi

    return a0-a1


def mean(x):
    """Mean value of a vector
    """
    return(float(sum(x))/len(x))


def cov(x, y=None):
    """Covariance of vectors x and y.

    If y is None: return cov(x, x) 
    """
    
    if y is None:
        y = x 

    assert(len(x)==len(y))
    N = len(x)
 
    cx = x - mean(x)  
    cy = y - mean(y)  

    p = innerproduct(cx,cy)/N
    return(p)


def err(x, y=0, n=2, relative=True):
    """Relative error of ||x-y|| to ||y||
       n = 2:    Two norm
       n = None: Max norm

       If denominator evaluates to zero or if y is omitted,
       absolute error is returned
    """

    x = ensure_numeric(x)
    if y:
        y = ensure_numeric(y)        

    if n == 2:
        err = norm(x-y)
        if relative is True:
            try:
                err = err/norm(y)
            except:
                pass

    else:
        err = max(abs(x-y))
        if relative is True:
            try:
                err = err/max(abs(y))    
            except:
                pass
          
    return err
  

def norm(x):
    """2-norm of x
    """
  
    y = ravel(x)
    p = sqrt(innerproduct(y,y))
    return p
    
  
def corr(x, y=None):
    """Correlation of x and y
    If y is None return autocorrelation of x
    """

    from math import sqrt
    if y is None:
        y = x 

    varx = cov(x)
    vary = cov(y)

    if varx == 0 or vary == 0:
        C = 0
    else:  
        C = cov(x,y)/sqrt(varx * vary)   

    return(C)


        
def ensure_numeric(A, typecode = None):
    """Ensure that sequence is a Numeric array.
    Inputs:
        A: Sequence. If A is already a Numeric array it will be returned
                     unaltered
                     If not, an attempt is made to convert it to a Numeric
                     array
        typecode: Numeric type. If specified, use this in the conversion.
                                If not, let Numeric decide

    This function is necessary as array(A) can cause memory overflow.
    """

    if typecode is None:
        if type(A) == ArrayType:
            return A
        else:
            return array(A)
    else:
        if type(A) == ArrayType:
            if A.dtype.char == typecode:
                return array(A)  #FIXME: Shouldn't this just return A?
            else:
                return A.astype(typecode)
        else:
            return array(A).astype(typecode)




def histogram(a, bins):
    """Standard histogram straight from the Numeric manual
    """

    n = searchsorted(sort(a), bins)
    n = concatenate( [n, [len(a)]] )
    return n[1:]-n[:-1]



####################################################################
#Python versions of function that are also implemented in numerical_tools_ext.c
#

def gradient_python(x0, y0, x1, y1, x2, y2, q0, q1, q2):
    """
    """

    det = (y2-y0)*(x1-x0) - (y1-y0)*(x2-x0)
    a = (y2-y0)*(q1-q0) - (y1-y0)*(q2-q0)
    a /= det

    b = (x1-x0)*(q2-q0) - (x2-x0)*(q1-q0)
    b /= det

    return a, b


def gradient2_python(x0, y0, x1, y1, q0, q1):
    """Compute radient based on two points and enforce zero gradient
    in the direction orthogonal to (x1-x0), (y1-y0) 
    """

    #Old code
    #det = x0*y1 - x1*y0
    #if det != 0.0:
    #    a = (y1*q0 - y0*q1)/det
    #    b = (x0*q1 - x1*q0)/det

    #Correct code (ON)
    det = (x1-x0)**2 + (y1-y0)**2
    if det != 0.0:
        a = (x1-x0)*(q1-q0)/det
        b = (y1-y0)*(q1-q0)/det
        
    return a, b        


##############################################
#Initialise module

#from utilities import compile
import compile
if compile.can_use_C_extension('util_ext.c'):
    from util_ext import gradient, gradient2
else:
    gradient = gradient_python
    gradient2 = gradient2_python    


if __name__ == "__main__":
    pass