# source:anuga_core/source/anuga/utilities/numerical_tools.py@4249

Last change on this file since 4249 was 4249, checked in by ole, 17 years ago

Made get_machine_precision() fast by using precomputed, global variable.

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