source:inundation/utilities/numerical_tools.py@2710

Last change on this file since 2710 was 2710, checked in by ole, 17 years ago

Cleaned up in angle computation

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