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

Last change on this file since 5561 was 5561, checked in by duncan, 16 years ago

Clarifying err

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    msg = 'Lengths must be equal: len(x) == %d, len(y) == %d' %(len(x), len(y))
167    assert(len(x)==len(y)), msg
168
169    N = len(x)
170    cx = x - mean(x)
171    cy = y - mean(y)
172
173    p = innerproduct(cx,cy)/N
174    return(p)
175
176
177def err(x, y=0, n=2, relative=True):
178    """Relative error of ||x-y|| to ||y||
179       n = 2:    Two norm
180       n = None: Max norm
181
182       If denominator evaluates to zero or
183       if y is omitted or
184       if keyword relative is False,
185       absolute error is returned
186
187       If there is x and y, n=2 and relative=False, this will calc;
188       sqrt(sum_over_x&y((xi - yi)^2))
189
190       Given this value (err), to calc the root mean square deviation, do
191       err/sqrt(n)
192       where n is the number of elements,(len(x))
193    """
194
195    x = ensure_numeric(x)
196    if y:
197        y = ensure_numeric(y)
198
199    if n == 2:
200        err = norm(x-y)
201        if relative is True:
202            try:
203                err = err/norm(y)
204            except:
205                pass
206
207    else:
208        err = max(abs(x-y))
209        if relative is True:
210            try:
211                err = err/max(abs(y))
212            except:
213                pass
214
215    return err
216
217
218def norm(x):
219    """2-norm of x
220    """
221
222    y = ravel(x)
223    p = sqrt(innerproduct(y,y))
224    return p
225
226
227def corr(x, y=None):
228    """Correlation of x and y
229    If y is None return autocorrelation of x
230    """
231
232    from math import sqrt
233    if y is None:
234        y = x
235
236    varx = cov(x)
237    vary = cov(y)
238
239    if varx == 0 or vary == 0:
240        C = 0
241    else:
242        C = cov(x,y)/sqrt(varx * vary)
243
244    return(C)
245
246
247
248def ensure_numeric(A, typecode = None):
249    """Ensure that sequence is a Numeric array.
250    Inputs:
251        A: Sequence. If A is already a Numeric array it will be returned
252                     unaltered
253                     If not, an attempt is made to convert it to a Numeric
254                     array
255        A: Scalar.   Return 0-dimensional array of length 1, containing that value
256        A: String.   Array of ASCII values
257        typecode: Numeric type. If specified, use this in the conversion.
258                                If not, let Numeric decide
259
260    This function is necessary as array(A) can cause memory overflow.
261    """
262
263    if typecode is None:
264        if type(A) == ArrayType:
265            return A
266        else:
267            return array(A)
268    else:
269        if type(A) == ArrayType:
270            if A.typecode == typecode:
271                return array(A)  #FIXME: Shouldn't this just return A?
272            else:
273                return array(A,typecode)
274        else:
275            return array(A,typecode)
276
277
278
279
280def histogram(a, bins, relative=False):
281    """Standard histogram straight from the Numeric manual
282
283    If relative is True, values will be normalised againts the total and
284    thus represent frequencies rather than counts.
285    """
286
287    n = searchsorted(sort(a), bins)
288    n = concatenate( [n, [len(a)]] )
289
290    hist = n[1:]-n[:-1]
291
292    if relative is True:
293        hist = hist/float(sum(hist))
294
295    return hist
296
297def create_bins(data, number_of_bins = None):
298    """Safely create bins for use with histogram
299    If data contains only one point or is constant, one bin will be created.
300    If number_of_bins in omitted 10 bins will be created
301    """
302
303    mx = max(data)
304    mn = min(data)
305
306    if mx == mn:
307        bins = array([mn])
308    else:
309        if number_of_bins is None:
310            number_of_bins = 10
311
312        bins = arange(mn, mx, (mx-mn)/number_of_bins)
313
314    return bins
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#-----------------
374#Initialise module
375
376from anuga.utilities import compile
377if compile.can_use_C_extension('util_ext.c'):