source: inundation/utilities/numerical_tools.py @ 2709

Last change on this file since 2709 was 2709, checked in by ole, 18 years ago

Added algorithm for boundary_polygon in the presence of
multiple vertex coordinates. Also tested that new fit_interpolate
can use this new version.

File size: 6.2 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   
11except:
12    #print 'Could not find scipy - using Numeric'
13    from Numeric import ArrayType, array, sum, innerproduct, ravel, sqrt, searchsorted, sort, concatenate, Float   
14
15# Getting an infinite number to use when using Numeric
16INF = (array([1])/0.)[0]
17
18
19
20def angle(v1, v2=None):
21    """Compute angle between 2D vectors v1 and v2.
22   
23    If v2 is not specified it will default
24    to e1 (the unit vector in the x-direction)
25
26    The angle is measured as a number in [0, 2pi] from v2 to v1.
27    """
28    from math import acos, pi, sqrt
29 
30    # Prepare two Numeric vectors
31    if v2 is None:
32        v2 = [1.0, 0.0] # Unit vector along the x-axis
33       
34    v1 = ensure_numeric(v1, Float)
35    v2 = ensure_numeric(v2, Float)   
36   
37    # Normalise
38    v1 = v1/sqrt(sum(v1**2))
39    v2 = v2/sqrt(sum(v2**2))
40   
41    # Compute angle
42    p = innerproduct(v1, v2)
43    c = crossproduct_length(v1, v2)
44   
45    theta = acos(p)
46   
47    # Correct if v1 is in quadrant 3 or 4 with respect to v2 (as the x-axis)
48    # If v2 was the unit vector [1,0] this would correspond to the test
49    # if v1[1] < 0: theta = 2*pi-theta   
50    # In general we use the sign of the cross product length
51    if c > 0: 
52       #Quadrant 3 or 4
53       theta = 2*pi-theta       
54       
55    return theta
56
57   
58def anglediff(v0, v1):
59    """Compute difference between angle of vector v0 (x0, y0) and v1 (x1, y1).
60    This is used for determining the ordering of vertices,
61    e.g. for checking if they are counter clockwise.
62
63    Always return a positive value
64    """
65
66    from math import pi
67
68    a0 = angle(v0)
69    a1 = angle(v1)
70
71    #Ensure that difference will be positive
72    if a0 < a1:
73        a0 += 2*pi
74
75    return a0-a1
76
77def normal_vector(v):
78    """Normal vector to v
79    """
80   
81    return array([-v[1], v[0]], Float)
82
83   
84def crossproduct_length(v1, v2):
85    return v1[0]*v2[1]-v2[0]*v1[1]
86   
87       
88def mean(x):
89    """Mean value of a vector
90    """
91    return(float(sum(x))/len(x))
92
93
94def cov(x, y=None):
95    """Covariance of vectors x and y.
96
97    If y is None: return cov(x, x)
98    """
99   
100    if y is None:
101        y = x
102
103    assert(len(x)==len(y))
104    N = len(x)
105 
106    cx = x - mean(x) 
107    cy = y - mean(y) 
108
109    p = innerproduct(cx,cy)/N
110    return(p)
111
112
113def err(x, y=0, n=2, relative=True):
114    """Relative error of ||x-y|| to ||y||
115       n = 2:    Two norm
116       n = None: Max norm
117
118       If denominator evaluates to zero or if y is omitted,
119       absolute error is returned
120    """
121
122    x = ensure_numeric(x)
123    if y:
124        y = ensure_numeric(y)       
125
126    if n == 2:
127        err = norm(x-y)
128        if relative is True:
129            try:
130                err = err/norm(y)
131            except:
132                pass
133
134    else:
135        err = max(abs(x-y))
136        if relative is True:
137            try:
138                err = err/max(abs(y))   
139            except:
140                pass
141         
142    return err
143 
144
145def norm(x):
146    """2-norm of x
147    """
148 
149    y = ravel(x)
150    p = sqrt(innerproduct(y,y))
151    return p
152   
153 
154def corr(x, y=None):
155    """Correlation of x and y
156    If y is None return autocorrelation of x
157    """
158
159    from math import sqrt
160    if y is None:
161        y = x
162
163    varx = cov(x)
164    vary = cov(y)
165
166    if varx == 0 or vary == 0:
167        C = 0
168    else: 
169        C = cov(x,y)/sqrt(varx * vary)   
170
171    return(C)
172
173
174       
175def ensure_numeric(A, typecode = None):
176    """Ensure that sequence is a Numeric array.
177    Inputs:
178        A: Sequence. If A is already a Numeric array it will be returned
179                     unaltered
180                     If not, an attempt is made to convert it to a Numeric
181                     array
182        typecode: Numeric type. If specified, use this in the conversion.
183                                If not, let Numeric decide
184
185    This function is necessary as array(A) can cause memory overflow.
186    """
187
188    if typecode is None:
189        if type(A) == ArrayType:
190            return A
191        else:
192            return array(A)
193    else:
194        if type(A) == ArrayType:
195            if A.typecode == typecode:
196                return array(A)  #FIXME: Shouldn't this just return A?
197            else:
198                return A.astype(typecode)
199        else:
200            return array(A).astype(typecode)
201
202
203
204
205def histogram(a, bins):
206    """Standard histogram straight from the Numeric manual
207    """
208
209    n = searchsorted(sort(a), bins)
210    n = concatenate( [n, [len(a)]] )
211    return n[1:]-n[:-1]
212
213
214
215####################################################################
216#Python versions of function that are also implemented in numerical_tools_ext.c
217#
218
219def gradient_python(x0, y0, x1, y1, x2, y2, q0, q1, q2):
220    """
221    """
222
223    det = (y2-y0)*(x1-x0) - (y1-y0)*(x2-x0)
224    a = (y2-y0)*(q1-q0) - (y1-y0)*(q2-q0)
225    a /= det
226
227    b = (x1-x0)*(q2-q0) - (x2-x0)*(q1-q0)
228    b /= det
229
230    return a, b
231
232
233def gradient2_python(x0, y0, x1, y1, q0, q1):
234    """Compute radient based on two points and enforce zero gradient
235    in the direction orthogonal to (x1-x0), (y1-y0)
236    """
237
238    #Old code
239    #det = x0*y1 - x1*y0
240    #if det != 0.0:
241    #    a = (y1*q0 - y0*q1)/det
242    #    b = (x0*q1 - x1*q0)/det
243
244    #Correct code (ON)
245    det = (x1-x0)**2 + (y1-y0)**2
246    if det != 0.0:
247        a = (x1-x0)*(q1-q0)/det
248        b = (y1-y0)*(q1-q0)/det
249       
250    return a, b       
251
252
253##############################################
254#Initialise module
255
256from utilities import compile
257if compile.can_use_C_extension('util_ext.c'):
258    from util_ext import gradient, gradient2
259else:
260    gradient = gradient_python
261    gradient2 = gradient2_python   
262
263
264if __name__ == "__main__":
265    pass
266
267
268   
269def angle_obsolete(v):
270    """Compute angle between e1 (the unit vector in the x-direction)
271    and the specified vector v.
272   
273    Return a number in [0, 2pi]
274    """
275    from math import acos, pi, sqrt
276 
277    # Normalise v
278    v = ensure_numeric(v, Float)
279    v = v/sqrt(sum(v**2))
280   
281    # Compute angle
282    theta = acos(v[0])
283     
284    if v[1] < 0: 
285       #Quadrant 3 or 4
286        theta = 2*pi-theta
287   
288    return theta
289   
Note: See TracBrowser for help on using the repository browser.