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

Last change on this file since 3676 was 3676, checked in by duncan, 18 years ago

checked in hack in angle.

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