source: anuga_core/source/anuga/utilities/numerical_tools.py @ 6174

Last change on this file since 6174 was 6174, checked in by rwilson, 15 years ago

Changed .array(A) to .array(A, num.Int) where appropriate. Helps when going to numpy.

File size: 8.6 KB
RevLine 
[5897]1#!/usr/bin/env python
2"""Auxiliary numerical tools
3
4"""
5
6from math import acos, pi, sqrt
7from warnings import warn
8
[6158]9import Numeric as num
[5897]10
[6158]11NAN = (num.array([1])/0.)[0]
[5897]12# if we use a package that has NAN, this should be updated to use NAN.
13
14# Static variable used by get_machine_precision
15machine_precision = None
16
17
18def safe_acos(x):
19    """Safely compute acos
20
21    Protect against cases where input argument x is outside the allowed
22    interval [-1.0, 1.0] by no more than machine precision
23    """
24
25    error_msg = 'Input to acos is outside allowed domain [-1.0, 1.0].'+\
26                'I got %.12f' %x
27    warning_msg = 'Changing argument to acos from %.18f to %.1f' %(x, sign(x))
28
29    eps = get_machine_precision() # Machine precision
30    if x < -1.0:
31        if x < -1.0 - eps:
32            raise ValueError, errmsg
33        else:
34            warn(warning_msg)
35            x = -1.0
36
37    if x > 1.0:
38        if x > 1.0 + eps:
39            raise ValueError, errmsg
40        else:
41            print 'NOTE: changing argument to acos from %.18f to 1.0' %x
42            x = 1.0
43
44    return acos(x)
45
46
47def sign(x):
48    if x > 0: return 1
49    if x < 0: return -1
50    if x == 0: return 0   
51
52
53def is_scalar(x):
54    """True if x is a scalar (constant numeric value)
55    """
56
57    from types import IntType, FloatType
58    if type(x) in [IntType, FloatType]:
59        return True
60    else:
61        return False
62
63def angle(v1, v2=None):
64    """Compute angle between 2D vectors v1 and v2.
65   
66    If v2 is not specified it will default
67    to e1 (the unit vector in the x-direction)
68
69    The angle is measured as a number in [0, 2pi] from v2 to v1.
70    """
71 
72    # Prepare two Numeric vectors
73    if v2 is None:
74        v2 = [1.0, 0.0] # Unit vector along the x-axis
75       
[6158]76    v1 = ensure_numeric(v1, num.Float)
77    v2 = ensure_numeric(v2, num.Float)   
[5897]78   
79    # Normalise
[6158]80    v1 = v1/num.sqrt(num.sum(v1**2))
81    v2 = v2/num.sqrt(num.sum(v2**2))
[5897]82
83    # Compute angle
[6158]84    p = num.innerproduct(v1, v2)
85    c = num.innerproduct(v1, normal_vector(v2)) # Projection onto normal
[5897]86                                            # (negative cross product)
87       
88    theta = safe_acos(p)
89           
90   
91    # Correct if v1 is in quadrant 3 or 4 with respect to v2 (as the x-axis)
92    # If v2 was the unit vector [1,0] this would correspond to the test
93    # if v1[1] < 0: theta = 2*pi-theta   
94    # In general we use the sign of the projection onto the normal.
95    if c < 0: 
96       #Quadrant 3 or 4
97       theta = 2*pi-theta       
98       
99    return theta
100
101   
102def anglediff(v0, v1):
103    """Compute difference between angle of vector v0 (x0, y0) and v1 (x1, y1).
104    This is used for determining the ordering of vertices,
105    e.g. for checking if they are counter clockwise.
106
107    Always return a positive value
108    """
109
110    from math import pi
111
112    a0 = angle(v0)
113    a1 = angle(v1)
114
115    #Ensure that difference will be positive
116    if a0 < a1:
117        a0 += 2*pi
118
119    return a0-a1
120
121def normal_vector(v):
122    """Normal vector to v.
123
124    Returns vector 90 degrees counter clockwise to and of same length as v
125    """
126   
[6158]127    return num.array([-v[1], v[0]], num.Float)
[5897]128
129   
130#def crossproduct_length(v1, v2):
131#    return v1[0]*v2[1]-v2[0]*v1[1]
132   
133       
134def mean(x):
135    """Mean value of a vector
136    """
[6158]137    return(float(num.sum(x))/len(x))
[5897]138
139
140def cov(x, y=None):
141    """Covariance of vectors x and y.
142
143    If y is None: return cov(x, x)
144    """
145   
146    if y is None:
147        y = x
148
149    x = ensure_numeric(x)
150    y = ensure_numeric(y)
151    msg = 'Lengths must be equal: len(x) == %d, len(y) == %d' %(len(x), len(y))
152    assert(len(x)==len(y)), msg
153
154    N = len(x) 
155    cx = x - mean(x) 
156    cy = y - mean(y) 
157
[6158]158    p = num.innerproduct(cx,cy)/N
[5897]159    return(p)
160
161
162def err(x, y=0, n=2, relative=True):
163    """Relative error of ||x-y|| to ||y||
164       n = 2:    Two norm
165       n = None: Max norm
166
167       If denominator evaluates to zero or
168       if y is omitted or
169       if keyword relative is False,
170       absolute error is returned
171
172       If there is x and y, n=2 and relative=False, this will calc;
173       sqrt(sum_over_x&y((xi - yi)^2))
174
175       Given this value (err), to calc the root mean square deviation, do
176       err/sqrt(n)
177       where n is the number of elements,(len(x))
178    """
179
180    x = ensure_numeric(x)
181    if y:
182        y = ensure_numeric(y)       
183
184    if n == 2:
185        err = norm(x-y)
186        if relative is True:
187            try:
188                err = err/norm(y)
189            except:
190                pass
191
192    else:
193        err = max(abs(x-y))
194        if relative is True:
195            try:
196                err = err/max(abs(y))   
197            except:
198                pass
199         
200    return err
201 
202
203def norm(x):
204    """2-norm of x
205    """
206 
[6158]207    y = num.ravel(x)
208    p = num.sqrt(num.innerproduct(y,y))
[5897]209    return p
210   
211 
212def corr(x, y=None):
213    """Correlation of x and y
214    If y is None return autocorrelation of x
215    """
216
217    from math import sqrt
218    if y is None:
219        y = x
220
221    varx = cov(x)
222    vary = cov(y)
223
224    if varx == 0 or vary == 0:
225        C = 0
226    else: 
227        C = cov(x,y)/sqrt(varx * vary)   
228
229    return(C)
230
231
232       
233def ensure_numeric(A, typecode = None):
234    """Ensure that sequence is a Numeric array.
235    Inputs:
236        A: Sequence. If A is already a Numeric array it will be returned
237                     unaltered
238                     If not, an attempt is made to convert it to a Numeric
239                     array
240        A: Scalar.   Return 0-dimensional array of length 1, containing that value
241        A: String.   Array of ASCII values
242        typecode: Numeric type. If specified, use this in the conversion.
243                                If not, let Numeric decide
244
245    This function is necessary as array(A) can cause memory overflow.
246    """
247
248    if typecode is None:
[6158]249        if type(A) == num.ArrayType:
[5897]250            return A
251        else:
[6158]252            return num.array(A)
[5897]253    else:
[6158]254        if type(A) == num.ArrayType:
[5897]255            if A.typecode == typecode:
[6174]256#                return num.array(A)  #FIXME: Shouldn't this just return A?
257                return A
[5897]258            else:
[6158]259                return num.array(A,typecode)
[5897]260        else:
[6158]261            return num.array(A,typecode)
[5897]262
263
264
265
266def histogram(a, bins, relative=False):
267    """Standard histogram straight from the Numeric manual
268
269    If relative is True, values will be normalised againts the total and
270    thus represent frequencies rather than counts.
271    """
272
[6158]273    n = num.searchsorted(num.sort(a), bins)
274    n = num.concatenate( [n, [len(a)]] )
[5897]275
276    hist = n[1:]-n[:-1]
277
278    if relative is True:
[6158]279        hist = hist/float(num.sum(hist))
[5897]280       
281    return hist
282
283def create_bins(data, number_of_bins = None):
284    """Safely create bins for use with histogram
285    If data contains only one point or is constant, one bin will be created.
286    If number_of_bins in omitted 10 bins will be created
287    """
288
289    mx = max(data)
290    mn = min(data)
291
292    if mx == mn:
[6158]293        bins = num.array([mn])
[5897]294    else:
295        if number_of_bins is None:
296            number_of_bins = 10
297           
[6158]298        bins = num.arange(mn, mx, (mx-mn)/number_of_bins)
[5897]299
300    return bins
301   
302
303
304def get_machine_precision():
305    """Calculate the machine precision for Floats
306
307    Depends on static variable machine_precision in this module
308    as this would otherwise require too much computation.
309    """
310
311    global machine_precision
312   
313    if machine_precision is None:
314        epsilon = 1.
315        while epsilon/2 + 1. > 1.:
316            epsilon /= 2
317
318        machine_precision = epsilon
319
320    return machine_precision   
321
322####################################################################
323#Python versions of function that are also implemented in numerical_tools_ext.c
324# FIXME (Ole): Delete these and update tests
325#
326
327def gradient_python(x0, y0, x1, y1, x2, y2, q0, q1, q2):
328    """
329    """
330
331    det = (y2-y0)*(x1-x0) - (y1-y0)*(x2-x0)
332    a = (y2-y0)*(q1-q0) - (y1-y0)*(q2-q0)
333    a /= det
334
335    b = (x1-x0)*(q2-q0) - (x2-x0)*(q1-q0)
336    b /= det
337
338    return a, b
339
340
341def gradient2_python(x0, y0, x1, y1, q0, q1):
342    """Compute radient based on two points and enforce zero gradient
343    in the direction orthogonal to (x1-x0), (y1-y0)
344    """
345
346    #Old code
347    #det = x0*y1 - x1*y0
348    #if det != 0.0:
349    #    a = (y1*q0 - y0*q1)/det
350    #    b = (x0*q1 - x1*q0)/det
351
352    #Correct code (ON)
353    det = (x1-x0)**2 + (y1-y0)**2
354    if det != 0.0:
355        a = (x1-x0)*(q1-q0)/det
356        b = (y1-y0)*(q1-q0)/det
357       
358    return a, b       
359
360
361#-----------------
362#Initialise module
363
364from anuga.utilities import compile
365if compile.can_use_C_extension('util_ext.c'):
366    from util_ext import gradient, gradient2
367else:
368    gradient = gradient_python
369    gradient2 = gradient2_python   
370
371
[6119]372if __name__ == '__main__':
[5897]373    pass
374
375   
Note: See TracBrowser for help on using the repository browser.