source: anuga_work/development/Hinwood_2008/where_close.py @ 5455

Last change on this file since 5455 was 5400, checked in by duncan, 16 years ago

adding a file

File size: 7.0 KB
Line 
1#!/usr/bin/python -tt
2#=======================================================================
3#                        General Documentation
4
5"""Single-function module.
6
7   See function docstring for description.
8"""
9
10#-----------------------------------------------------------------------
11#                       Additional Documentation
12#
13# RCS Revision Code:
14#   $Id: where_close.py,v 1.2 2004/04/28 00:28:15 jlin Exp $
15#
16# Modification History:
17# - 19 Mar 2004:  Original by Johnny Lin, Computation Institute,
18#   University of Chicago.  Passed reasonable tests.
19#
20# Notes:
21# - Written for Python 2.2.2.
22# - Function is based on code from the MA module by Paul F. Dubois.
23#   Some snippets of code in this function are copied directly from
24#   lines in that module.
25# - Module docstrings can be tested using the doctest module.  To
26#   test, execute "python where_close.py".
27# - See import statements throughout for packages/modules required.
28#
29# Copyright (c) 2004 by Johnny Lin.  For licensing, distribution
30# conditions, contact information, and additional documentation see
31# the URL http://www.johnny-lin.com/py_pkgs/gemath/doc/;
32
33# This library is free software; you can redistribute it and/or modify
34# it under the terms of the GNU Lesser General Public License as
35# published by the Free Software Foundation; either version 2.1 of the
36# License, or (at your option) any later version.
37
38# This library is distributed in the hope that it will be useful, but
39# WITHOUT ANY WARRANTY; without even the implied warranty of
40# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
41# Lesser General Public License for more details.
42
43# You should have received a copy of the GNU Lesser General Public
44# License along with this library; if not, write to the Free Software
45# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
46# USA.
47
48# You can contact Johnny Lin at his email address or at the University
49# of Chicago, Department of the Geophysical Sciences, 5734 S. Ellis
50# Ave., Chicago, IL 60637, USA.
51
52#=======================================================================
53
54
55
56#--------------------------- General Function --------------------------
57
58def where_close(x, y, rtol=1.e-5, atol=1.e-8):
59    """Mask of where x and y are element-wise "equal" to each other.
60
61    Returns a long integer array with elements equal to 1 where x and
62    y are "equal", and 0 otherwise.  If x or y are floating point,
63    "equal" means where abs(x-y) <= atol + rtol * abs(y).  This is
64    essentially the same algorithm used in the Numeric function
65    allclose.  If x and y are integer, "equal" means strict equality. 
66    Shape and size of output is the same as x and y; if one is an
67    array and the other is scalar, shape and size of the output is the
68    same as the array.  Output is a Numeric array, unless both inputs
69    are scalar in which the output is a Python integer scalar.
70
71    Positional Input Arguments:
72    * x:  Scalar or Numeric array, Python list/tuple of any size and
73      shape.  Floating or integer type.
74    * y:  Scalar or Numeric array, Python list/tuple of any size and
75      shape.  Floating or integer type.
76
77    Keyword Input Arguments:
78    * rtol:   "Relative" tolerance.  Default is 1.e-5.  Used in the
79              comparison between x and y only if the two are floating
80              point.
81    * atol:   "Absolute" tolerance.  Default is 1.e-8.  Used in the
82              comparison between x and y only if the two are floating
83              point.
84
85    Examples:
86    >>> import Numeric as N
87    >>> from where_close import where_close
88    >>> x = [20.,  -32., -1., 2.            , 5., 29.]
89    >>> y = [20.1, -31., -1., 2.000000000001, 3., 28.99]
90    >>> ind = where_close(x, y)
91    >>> ['%.1g' % ind[i] for i in range(len(ind))]
92    ['0', '0', '1', '1', '0', '0']
93
94    >>> x = N.array([1,  5,  7, -2, 10])
95    >>> y = N.array([1, -5, 17, -2,  0])
96    >>> ind = where_close(x, y)
97    >>> ['%.1g' % ind[i] for i in range(len(ind))]
98    ['1', '0', '0', '1', '0']
99    """
100    import Numeric as N
101    abs = N.absolute
102
103
104    #- Make sure input is Numeric type:
105
106    xN = N.array(x)
107    yN = N.array(y)
108
109
110    #- Safe compare if floating.  Strict compare if integer.  Any other
111    #  type returns an error:
112
113    if (xN.typecode() in N.typecodes['Float']) or \
114       (yN.typecode() in N.typecodes['Float']):
115        return N.less_equal(abs(xN-yN), atol+rtol*abs(yN))
116
117    elif (xN.typecode() in N.typecodes['Integer']) and \
118         (yN.typecode() in N.typecodes['Integer']):
119        return N.equal(xN, yN)
120
121    else:
122        raise ValueError, "where_close:  Inputs must be Float or Integer"
123
124
125
126
127#-------------------------- Main:  Test Module -------------------------
128
129#- Define additional examples for doctest to use:
130
131__test__ = { 'Additional Examples':
132    """
133    >>> from where_close import where_close
134    >>> import Numeric as N
135    >>> x = [20.,  -32., -1., 2.            , 5., 29.]
136    >>> y = [20.1, -31., -1., 2.000000000001, 3., 28.99]
137    >>> x = N.reshape(x, (2,3))
138    >>> y = N.reshape(y, (2,3))
139    >>> ind = where_close(x, y)
140    >>> ['%.1g' % ind[0,i] for i in range(3)]
141    ['0', '0', '1']
142    >>> ['%.1g' % ind[1,i] for i in range(3)]
143    ['1', '0', '0']
144    >>> ind.shape
145    (2, 3)
146    >>> ind.typecode()
147    'l'
148    >>> type(ind)
149    <type 'array'>
150
151    >>> x = [20.,  -32., -1., 2.            , 5., 29.]
152    >>> y = [20.1, -31., -1., 2.000000000001, 3.]
153    >>> ind = where_close(x, y)
154    Traceback (most recent call last):
155        ...
156    ValueError: frames are not aligned
157
158    >>> x = [20, -32, -1]
159    >>> y = [20, -32]
160    >>> ind = where_close(x, y)
161    Traceback (most recent call last):
162        ...
163    ValueError: frames are not aligned
164
165    >>> x = [20,  -32, -1.0, 2, 5, 29]
166    >>> y = 2.
167    >>> ind = where_close(x, y)
168    >>> ['%.1g' % ind[i] for i in range(len(ind))]
169    ['0', '0', '0', '1', '0', '0']
170    >>> x = -32
171    >>> y = N.array([20.,  -32., -1., 2., 5., 29.])
172    >>> ind = where_close(x, y)
173    >>> ['%.1g' % ind[i] for i in range(len(ind))]
174    ['0', '1', '0', '0', '0', '0']
175    >>> x = -32
176    >>> y = -33.
177    >>> ind = where_close(x, y)
178    >>> print ind
179    0
180    >>> type(ind)
181    <type 'int'>
182    >>> x = -33
183    >>> y = -33.
184    >>> ind = where_close(x, y)
185    >>> print ind
186    1
187    >>> x = -33
188    >>> y = -33
189    >>> ind = where_close(x, y)
190    >>> print ind
191    1
192    >>> x = -3.
193    >>> y = -3.
194    >>> ind = where_close(x, y)
195    >>> print ind
196    1
197    """ }
198
199
200#- Execute doctest if module is run from command line:
201
202if __name__ == "__main__":
203    """Test the module.
204
205    Tests the examples in all the module documentation
206    strings, plus __test__.
207
208    Note:  To help ensure that module testing of this file works, the
209    parent directory to the current directory is added to sys.path.
210    """
211    import doctest, sys, os
212    sys.path.append(os.pardir)
213    doctest.testmod(sys.modules[__name__])
214
215
216
217
218# ===== end file =====
219
Note: See TracBrowser for help on using the repository browser.