source: anuga_core/source/anuga/utilities/numerical_tools.py @ 5729

Last change on this file since 5729 was 5729, checked in by ole, 16 years ago

Refactored get_flow_through_cross_section and added a runtime version to
shallow_water_domain. It still needs some work and testing.

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: 
111       #Quadrant 3 or 4
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   
316
317
318def get_machine_precision():
319    """Calculate the machine precision for Floats
320
321    Depends on static variable machine_precision in this module
322    as this would otherwise require too much computation.
323    """
324
325    global machine_precision
326   
327    if machine_precision is None:
328        epsilon = 1.
329        while epsilon/2 + 1. > 1.:
330            epsilon /= 2
331
332        machine_precision = epsilon
333
334    return machine_precision   
335
336####################################################################
337#Python versions of function that are also implemented in numerical_tools_ext.c
338# FIXME (Ole): Delete these and update tests
339#
340
341def gradient_python(x0, y0, x1, y1, x2, y2, q0, q1, q2):
342    """
343    """
344
345    det = (y2-y0)*(x1-x0) - (y1-y0)*(x2-x0)
346    a = (y2-y0)*(q1-q0) - (y1-y0)*(q2-q0)
347    a /= det
348
349    b = (x1-x0)*(q2-q0) - (x2-x0)*(q1-q0)
350    b /= det
351
352    return a, b
353
354
355def gradient2_python(x0, y0, x1, y1, q0, q1):
356    """Compute radient based on two points and enforce zero gradient
357    in the direction orthogonal to (x1-x0), (y1-y0)
358    """
359
360    #Old code
361    #det = x0*y1 - x1*y0
362    #if det != 0.0:
363    #    a = (y1*q0 - y0*q1)/det
364    #    b = (x0*q1 - x1*q0)/det
365
366    #Correct code (ON)
367    det = (x1-x0)**2 + (y1-y0)**2
368    if det != 0.0:
369        a = (x1-x0)*(q1-q0)/det
370        b = (y1-y0)*(q1-q0)/det
371       
372    return a, b       
373
374
375#-----------------
376#Initialise module
377
378from anuga.utilities import compile
379if compile.can_use_C_extension('util_ext.c'):
380    from util_ext import gradient, gradient2
381else:
382    gradient = gradient_python
383    gradient2 = gradient2_python   
384
385
386if __name__ == "__main__":
387    pass
388
389   
Note: See TracBrowser for help on using the repository browser.