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

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

Hi all,
I'm doing a change in the anuga structure, moving the code to

\anuga_core\source\anuga

After you have done an svn update, the PYTHONPATH has to be changed to;
PYTHONPATH = anuga_core/source/

This is part of changes required to make installation of anuga quicker and reducing the size of our sandpits.

If any imports are broken, try fixing them. With adding anuga. to them for example. If this seems to have really broken things, email/phone me.

Cheers
Duncan

File size: 7.3 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
49    theta = acos(p)
50
51    # Correct if v1 is in quadrant 3 or 4 with respect to v2 (as the x-axis)
52    # If v2 was the unit vector [1,0] this would correspond to the test
53    # if v1[1] < 0: theta = 2*pi-theta
54    # In general we use the sign of the projection onto the normal.
55    if c < 0:
57       theta = 2*pi-theta
58
59    return theta
60
61
62def anglediff(v0, v1):
63    """Compute difference between angle of vector v0 (x0, y0) and v1 (x1, y1).
64    This is used for determining the ordering of vertices,
65    e.g. for checking if they are counter clockwise.
66
67    Always return a positive value
68    """
69
70    from math import pi
71
72    a0 = angle(v0)
73    a1 = angle(v1)
74
75    #Ensure that difference will be positive
76    if a0 < a1:
77        a0 += 2*pi
78
79    return a0-a1
80
81def normal_vector(v):
82    """Normal vector to v.
83
84    Returns vector 90 degrees counter clockwise to and of same length as v
85    """
86
87    return array([-v[1], v[0]], Float)
88
89
90#def crossproduct_length(v1, v2):
91#    return v1[0]*v2[1]-v2[0]*v1[1]
92
93
94def mean(x):
95    """Mean value of a vector
96    """
97    return(float(sum(x))/len(x))
98
99
100def cov(x, y=None):
101    """Covariance of vectors x and y.
102
103    If y is None: return cov(x, x)
104    """
105
106    if y is None:
107        y = x
108
109    assert(len(x)==len(y))
110    N = len(x)
111
112    cx = x - mean(x)
113    cy = y - mean(y)
114
115    p = innerproduct(cx,cy)/N
116    return(p)
117
118
119def err(x, y=0, n=2, relative=True):
120    """Relative error of ||x-y|| to ||y||
121       n = 2:    Two norm
122       n = None: Max norm
123
124       If denominator evaluates to zero or
125       if y is omitted or
126       if keyword relative is False,
127       absolute error is returned
128    """
129
130    x = ensure_numeric(x)
131    if y:
132        y = ensure_numeric(y)
133
134    if n == 2:
135        err = norm(x-y)
136        if relative is True:
137            try:
138                err = err/norm(y)
139            except:
140                pass
141
142    else:
143        err = max(abs(x-y))
144        if relative is True:
145            try:
146                err = err/max(abs(y))
147            except:
148                pass
149
150    return err
151
152
153def norm(x):
154    """2-norm of x
155    """
156
157    y = ravel(x)
158    p = sqrt(innerproduct(y,y))
159    return p
160
161
162def corr(x, y=None):
163    """Correlation of x and y
164    If y is None return autocorrelation of x
165    """
166
167    from math import sqrt
168    if y is None:
169        y = x
170
171    varx = cov(x)
172    vary = cov(y)
173
174    if varx == 0 or vary == 0:
175        C = 0
176    else:
177        C = cov(x,y)/sqrt(varx * vary)
178
179    return(C)
180
181
182
183def ensure_numeric(A, typecode = None):
184    """Ensure that sequence is a Numeric array.
185    Inputs:
186        A: Sequence. If A is already a Numeric array it will be returned
187                     unaltered
188                     If not, an attempt is made to convert it to a Numeric
189                     array
190        typecode: Numeric type. If specified, use this in the conversion.
191                                If not, let Numeric decide
192
193    This function is necessary as array(A) can cause memory overflow.
194    """
195
196    if typecode is None:
197        if type(A) == ArrayType:
198            return A
199        else:
200            return array(A)
201    else:
202        if type(A) == ArrayType:
203            if A.typecode == typecode:
204                return array(A)  #FIXME: Shouldn't this just return A?
205            else:
206                return array(A,typecode)
207        else:
208            return array(A,typecode)
209
210
211
212
213def histogram(a, bins, relative=False):
214    """Standard histogram straight from the Numeric manual
215
216    If relative is True, values will be normalised againts the total and
217    thus represent frequencies rather than counts.
218    """
219
220    n = searchsorted(sort(a), bins)
221    n = concatenate( [n, [len(a)]] )
222
223    hist = n[1:]-n[:-1]
224
225    if relative is True:
226        hist = hist/float(sum(hist))
227
228    return hist
229
230def create_bins(data, number_of_bins = None):
231    """Safely create bins for use with histogram
232    If data contains only one point or is constant, one bin will be created.
233    If number_of_bins in omitted 10 bins will be created
234    """
235
236    mx = max(data)
237    mn = min(data)
238
239    if mx == mn:
240        bins = array([mn])
241    else:
242        if number_of_bins is None:
243            number_of_bins = 10
244
245        bins = arange(mn, mx, (mx-mn)/number_of_bins)
246
247    return bins
248
249
250
251####################################################################
252#Python versions of function that are also implemented in numerical_tools_ext.c
253#
254
255def gradient_python(x0, y0, x1, y1, x2, y2, q0, q1, q2):
256    """
257    """
258
259    det = (y2-y0)*(x1-x0) - (y1-y0)*(x2-x0)
260    a = (y2-y0)*(q1-q0) - (y1-y0)*(q2-q0)
261    a /= det
262
263    b = (x1-x0)*(q2-q0) - (x2-x0)*(q1-q0)
264    b /= det
265
266    return a, b
267
268
269def gradient2_python(x0, y0, x1, y1, q0, q1):
271    in the direction orthogonal to (x1-x0), (y1-y0)
272    """
273
274    #Old code
275    #det = x0*y1 - x1*y0
276    #if det != 0.0:
277    #    a = (y1*q0 - y0*q1)/det
278    #    b = (x0*q1 - x1*q0)/det
279
280    #Correct code (ON)
281    det = (x1-x0)**2 + (y1-y0)**2
282    if det != 0.0:
283        a = (x1-x0)*(q1-q0)/det
284        b = (y1-y0)*(q1-q0)/det
285
286    return a, b
287
288
289##############################################
290#Initialise module
291
292from anuga.utilities import compile
293if compile.can_use_C_extension('util_ext.c'):
295else:
298
299
300if __name__ == "__main__":
301    pass
302
303
304
305def angle_obsolete(v):
306    """Compute angle between e1 (the unit vector in the x-direction)
307    and the specified vector v.
308
309    Return a number in [0, 2pi]
310    """
311    from math import acos, pi, sqrt
312
313    # Normalise v
314    v = ensure_numeric(v, Float)
315    v = v/sqrt(sum(v**2))
316
317    # Compute angle
318    theta = acos(v[0])
319
320    if v[1] < 0: