source: trunk/anuga_core/source/anuga/utilities/numerical_tools.py @ 7876

Last change on this file since 7876 was 7690, checked in by hudson, 14 years ago

Check for complex polygons raises an exception.

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