source: inundation-numpy-branch/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: 5.3 KB
Line 
1#!/usr/bin/env python
2"""Auxiliary numerical tools
3
4"""
5
6from numpy import ArrayType, array, sum, innerproduct, ravel, sqrt
7from numpy import searchsorted, sort, concatenate   
8   
9
10def angle(v):
11    """Compute angle between e1 (the unit vector in the x-direction)
12    and the specified vector
13    """
14
15    from math import acos, pi, sqrt
16
17    l = sqrt( sum (array(v)**2))
18    v1 = v[0]/l
19    v2 = v[1]/l
20
21
22    theta = acos(v1)
23   
24    #try:
25    #   theta = acos(v1)
26    #except ValueError, e:
27    #    print 'WARNING (util.py): Angle acos(%s) failed: %s' %(str(v1), e)
28    #
29    #    #FIXME, hack to avoid acos(1.0) Value error
30    #    # why is it happening?
31    #    # is it catching something we should avoid?
32    #    #FIXME (Ole): When did this happen? We need a unit test to
33    #    #reveal this condition or otherwise remove the hack
34    #   
35    #    s = 1e-6
36    #    if (v1+s >  1.0) and (v1-s < 1.0) :
37    #        theta = 0.0
38    #    elif (v1+s >  -1.0) and (v1-s < -1.0):
39    #        theta = 3.1415926535897931
40    #    print 'WARNING (util.py): angle v1 is %f, setting acos to %f '\
41    #          %(v1, theta)
42
43    if v2 < 0:
44        #Quadrant 3 or 4
45        theta = 2*pi-theta
46
47    return theta
48
49
50def anglediff(v0, v1):
51    """Compute difference between angle of vector x0, y0 and x1, y1.
52    This is used for determining the ordering of vertices,
53    e.g. for checking if they are counter clockwise.
54
55    Always return a positive value
56    """
57
58    from math import pi
59
60    a0 = angle(v0)
61    a1 = angle(v1)
62
63    #Ensure that difference will be positive
64    if a0 < a1:
65        a0 += 2*pi
66
67    return a0-a1
68
69
70def mean(x):
71    """Mean value of a vector
72    """
73    return(float(sum(x))/len(x))
74
75
76def cov(x, y=None):
77    """Covariance of vectors x and y.
78
79    If y is None: return cov(x, x)
80    """
81   
82    if y is None:
83        y = x
84
85    assert(len(x)==len(y))
86    N = len(x)
87 
88    cx = x - mean(x) 
89    cy = y - mean(y) 
90
91    p = innerproduct(cx,cy)/N
92    return(p)
93
94
95def err(x, y=0, n=2, relative=True):
96    """Relative error of ||x-y|| to ||y||
97       n = 2:    Two norm
98       n = None: Max norm
99
100       If denominator evaluates to zero or if y is omitted,
101       absolute error is returned
102    """
103
104    x = ensure_numeric(x)
105    if y:
106        y = ensure_numeric(y)       
107
108    if n == 2:
109        err = norm(x-y)
110        if relative is True:
111            try:
112                err = err/norm(y)
113            except:
114                pass
115
116    else:
117        err = max(abs(x-y))
118        if relative is True:
119            try:
120                err = err/max(abs(y))   
121            except:
122                pass
123         
124    return err
125 
126
127def norm(x):
128    """2-norm of x
129    """
130 
131    y = ravel(x)
132    p = sqrt(innerproduct(y,y))
133    return p
134   
135 
136def corr(x, y=None):
137    """Correlation of x and y
138    If y is None return autocorrelation of x
139    """
140
141    from math import sqrt
142    if y is None:
143        y = x
144
145    varx = cov(x)
146    vary = cov(y)
147
148    if varx == 0 or vary == 0:
149        C = 0
150    else: 
151        C = cov(x,y)/sqrt(varx * vary)   
152
153    return(C)
154
155
156       
157def ensure_numeric(A, typecode = None):
158    """Ensure that sequence is a Numeric array.
159    Inputs:
160        A: Sequence. If A is already a Numeric array it will be returned
161                     unaltered
162                     If not, an attempt is made to convert it to a Numeric
163                     array
164        typecode: Numeric type. If specified, use this in the conversion.
165                                If not, let Numeric decide
166
167    This function is necessary as array(A) can cause memory overflow.
168    """
169
170    if typecode is None:
171        if type(A) == ArrayType:
172            return A
173        else:
174            return array(A)
175    else:
176        if type(A) == ArrayType:
177            if A.dtype.char == typecode:
178                return array(A)  #FIXME: Shouldn't this just return A?
179            else:
180                return A.astype(typecode)
181        else:
182            return array(A).astype(typecode)
183
184
185
186
187def histogram(a, bins):
188    """Standard histogram straight from the Numeric manual
189    """
190
191    n = searchsorted(sort(a), bins)
192    n = concatenate( [n, [len(a)]] )
193    return n[1:]-n[:-1]
194
195
196
197####################################################################
198#Python versions of function that are also implemented in numerical_tools_ext.c
199#
200
201def gradient_python(x0, y0, x1, y1, x2, y2, q0, q1, q2):
202    """
203    """
204
205    det = (y2-y0)*(x1-x0) - (y1-y0)*(x2-x0)
206    a = (y2-y0)*(q1-q0) - (y1-y0)*(q2-q0)
207    a /= det
208
209    b = (x1-x0)*(q2-q0) - (x2-x0)*(q1-q0)
210    b /= det
211
212    return a, b
213
214
215def gradient2_python(x0, y0, x1, y1, q0, q1):
216    """Compute radient based on two points and enforce zero gradient
217    in the direction orthogonal to (x1-x0), (y1-y0)
218    """
219
220    #Old code
221    #det = x0*y1 - x1*y0
222    #if det != 0.0:
223    #    a = (y1*q0 - y0*q1)/det
224    #    b = (x0*q1 - x1*q0)/det
225
226    #Correct code (ON)
227    det = (x1-x0)**2 + (y1-y0)**2
228    if det != 0.0:
229        a = (x1-x0)*(q1-q0)/det
230        b = (y1-y0)*(q1-q0)/det
231       
232    return a, b       
233
234
235##############################################
236#Initialise module
237
238#from anuga.utilities import compile
239import compile
240if compile.can_use_C_extension('util_ext.c'):
241    from util_ext import gradient, gradient2
242else:
243    gradient = gradient_python
244    gradient2 = gradient2_python   
245
246
247if __name__ == "__main__":
248    pass
249
Note: See TracBrowser for help on using the repository browser.