# source:inundation/utilities/numerical_tools.py@2709

Last change on this file since 2709 was 2709, checked in by ole, 17 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:
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):
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'):
259else:
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: