source: branches/numpy/anuga/utilities/numerical_tools.py @ 7174

Last change on this file since 7174 was 7174, checked in by ole, 15 years ago

Removed divide-by-zero warning in numpy branch.

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