# source:branches/numpy/anuga/utilities/numerical_tools.py@6304

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

Initial commit of numpy changes. Still a long way to go.

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
11NAN = (num.array([1])/0.)[0]
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
76    v1 = ensure_numeric(v1, num.float)
77    v2 = ensure_numeric(v2, num.float)
78
79    # Normalise
80    v1 = v1/num.sqrt(num.sum(v1**2))
81    v2 = v2/num.sqrt(num.sum(v2**2))
82
83    # Compute angle
84    p = num.inner(v1, v2)
85    c = num.inner(v1, normal_vector(v2)) # Projection onto normal
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:
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
127    return num.array([-v[1], v[0]], num.float)
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    """
137    return(float(num.sum(x))/len(x))
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
158    p = num.inner(cx,cy)/N
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
207    y = num.ravel(x)
208    p = num.sqrt(num.inner(y,y))
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
233##
234# @brief Ensure that a sequence is a numeric array of the required type.
235# @param A The sequence object to convert to a numeric array.
236# @param typecode The required numeric type of object A (a numeric dtype).
237# @return A numeric array of the required type.
238def ensure_numeric(A, typecode=None):
239    """Ensure that sequence is a numeric array.
240    Inputs:
241        A: Sequence. If A is already a numeric array it will be returned
242                     unaltered
243                     If not, an attempt is made to convert it to a numeric
244                     array
245        A: Scalar.   Return 0-dimensional array containing that value. Note
246                     that a 0-dim array DOES NOT HAVE A LENGTH UNDER numpy.
247        A: String.   Array of ASCII values (numpy can't handle this)
248
249        typecode:    numeric type. If specified, use this in the conversion.
250                     If not, let numeric package decide.
251                     numpy assumes float64 if no type in A.
252
253    This function is necessary as array(A) can cause memory overflow.
254    """
255
256    if isinstance(A, basestring):
257        msg = 'Sorry, cannot handle string in ensure_numeric()'
258        raise Exception, msg
259
260    if typecode is None:
261        if isinstance(A, num.ndarray):
262            return A
263        else:
264            return num.array(A)
265    else:
266        if isinstance(A, num.ndarray):
267            if A.dtype == typecode:
268#                return num.array(A)  #FIXME: Shouldn't this just return A?
269                return A
270            else:
271                return num.array(A, dtype=typecode)
272        else:
273            return num.array(A, dtype=typecode)
274
275
276
277
278def histogram(a, bins, relative=False):
279    """Standard histogram straight from the numeric manual
280
281    If relative is True, values will be normalised againts the total and
282    thus represent frequencies rather than counts.
283    """
284
285    n = num.searchsorted(num.sort(a), bins)
286    n = num.concatenate( [n, [len(a)]] )
287
288    hist = n[1:]-n[:-1]
289
290    if relative is True:
291        hist = hist/float(num.sum(hist))
292
293    return hist
294
295def create_bins(data, number_of_bins = None):
296    """Safely create bins for use with histogram
297    If data contains only one point or is constant, one bin will be created.
298    If number_of_bins in omitted 10 bins will be created
299    """
300
301    mx = max(data)
302    mn = min(data)
303
304    if mx == mn:
305        bins = num.array([mn])
306    else:
307        if number_of_bins is None:
308            number_of_bins = 10
309
310        bins = num.arange(mn, mx, (mx-mn)/number_of_bins)
311
312    return bins
313
314
315
316def get_machine_precision():
317    """Calculate the machine precision for Floats
318
319    Depends on static variable machine_precision in this module
320    as this would otherwise require too much computation.
321    """
322
323    global machine_precision
324
325    if machine_precision is None:
326        epsilon = 1.
327        while epsilon/2 + 1. > 1.:
328            epsilon /= 2
329
330        machine_precision = epsilon
331
332    return machine_precision
333
334####################################################################
335#Python versions of function that are also implemented in numerical_tools_ext.c
336# FIXME (Ole): Delete these and update tests
337#
338
339def gradient_python(x0, y0, x1, y1, x2, y2, q0, q1, q2):
340    """
341    """
342
343    det = (y2-y0)*(x1-x0) - (y1-y0)*(x2-x0)
344    a = (y2-y0)*(q1-q0) - (y1-y0)*(q2-q0)
345    a /= det
346
347    b = (x1-x0)*(q2-q0) - (x2-x0)*(q1-q0)
348    b /= det
349
350    return a, b
351
352
353def gradient2_python(x0, y0, x1, y1, q0, q1):
355    in the direction orthogonal to (x1-x0), (y1-y0)
356    """
357
358    #Old code
359    #det = x0*y1 - x1*y0
360    #if det != 0.0:
361    #    a = (y1*q0 - y0*q1)/det
362    #    b = (x0*q1 - x1*q0)/det
363
364    #Correct code (ON)
365    det = (x1-x0)**2 + (y1-y0)**2
366    if det != 0.0:
367        a = (x1-x0)*(q1-q0)/det
368        b = (y1-y0)*(q1-q0)/det
369
370    return a, b
371
372################################################################################
373# Decision functions for numeric package objects.
374# It is a little tricky to decide if a numpy thing is of type float.
375# These functions hide numpy-specific details of how we do this.
376################################################################################
377
378##
379# @brief Decide if an object is a numeric package object with datatype of float.
380# @param obj The object to decide on.
381# @return True if 'obj' is a numeric package object, and some sort of float.
382def is_num_float(obj):
383    '''Is an object a numeric package float object?'''
384
385    try:
386        return obj.dtype.char in num.typecodes['Float']
387    except AttributeError:
388        return False
389
390##
391# @brief Decide if an object is a numeric package object with datatype of int.
392# @param obj The object to decide on.
393# @return True if 'obj' is a numeric package object, and some sort of int.
394def is_num_int(obj):
395    '''Is an object a numeric package int object?'''
396
397    try:
398        return obj.dtype.char in num.typecodes['Integer']
399    except AttributeError:
400        return False
401
402
403#-----------------
404#Initialise module
405
406from anuga.utilities import compile
407if compile.can_use_C_extension('util_ext.c'):