Skip to content
Navigation Menu
{{ message }}
forked from ipython/ipython
-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathnumutils.py
More file actions
302 lines (224 loc) · 9.81 KB
/
Copy pathnumutils.py
File metadata and controls
302 lines (224 loc) · 9.81 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
# -*- coding: utf-8 -*-
"""
A set of convenient utilities for numerical work.
Most of this module requires Numerical Python or is meant to be used with it.
See http://www.pfdubois.com/numpy for details.
"""
#*****************************************************************************
# Copyright (C) 2001-2005 Fernando Perez <fperez@colorado.edu>
#
# Distributed under the terms of the BSD License. The full license is in
# the file COPYING, distributed as part of this software.
#*****************************************************************************
__all__ = ['sum_flat','mean_flat','rms_flat','base_repr','binary_repr',
'amin','amax','amap','zeros_like','empty_like',
'frange','diagonal_matrix','identity',
'fromfunction_kw','log2','ispower2',
'norm','l1norm','l2norm','exp_safe',
'inf','infty','Infinity',
'Numeric']
#****************************************************************************
# required modules
import __main__
import math
import operator
import sys
import Numeric
from Numeric import *
#*****************************************************************************
# Globals
# useful for testing infinities in results of array divisions (which don't
# raise an exception)
# Python, LaTeX and Mathematica names.
inf = infty = Infinity = (array([1])/0.0)[0]
#****************************************************************************
# function definitions
exp_safe_MIN = math.log(2.2250738585072014e-308)
exp_safe_MAX = 1.7976931348623157e+308
def exp_safe(x):
"""Compute exponentials which safely underflow to zero.
Slow but convenient to use. Note that NumArray will introduce proper
floating point exception handling with access to the underlying
hardware."""
if type(x) is ArrayType:
return exp(clip(x,exp_safe_MIN,exp_safe_MAX))
else:
return math.exp(x)
def amap(fn,*args):
"""amap(function, sequence[, sequence, ...]) -> array.
Works like map(), but it returns an array. This is just a convenient
shorthand for Numeric.array(map(...))"""
return array(map(fn,*args))
def amin(m,axis=0):
"""amin(m,axis=0) returns the minimum of m along dimension axis.
"""
return minimum.reduce(asarray(m),axis)
def amax(m,axis=0):
"""amax(m,axis=0) returns the maximum of m along dimension axis.
"""
return maximum.reduce(asarray(m),axis)
def zeros_like(a):
"""Return an array of zeros of the shape and typecode of a.
If you don't explicitly need the array to be zeroed, you should instead
use empty_like(), which is faster as it only allocates memory."""
return zeros(a.shape,a.typecode())
def empty_like(a):
"""Return an empty (uninitialized) array of the shape and typecode of a.
Note that this does NOT initialize the returned array. If you require
your array to be initialized, you should use zeros_like().
This requires Numeric.empty(), which appeared in Numeric 23.7."""
return empty(a.shape,a.typecode())
def sum_flat(a):
"""Return the sum of all the elements of a, flattened out.
It uses a.flat, and if a is not contiguous, a call to ravel(a) is made."""
if a.iscontiguous():
return Numeric.sum(a.flat)
else:
return Numeric.sum(ravel(a))
def mean_flat(a):
"""Return the mean of all the elements of a, flattened out."""
return sum_flat(a)/float(size(a))
def rms_flat(a):
"""Return the root mean square of all the elements of a, flattened out."""
return math.sqrt(sum_flat(absolute(a)**2)/float(size(a)))
def l1norm(a):
"""Return the l1 norm of a, flattened out.
Implemented as a separate function (not a call to norm() for speed).
Ref: http://mathworld.wolfram.com/L1-Norm.html"""
return sum_flat(absolute(a))
def l2norm(a):
"""Return the l2 norm of a, flattened out.
Implemented as a separate function (not a call to norm() for speed).
Ref: http://mathworld.wolfram.com/L2-Norm.html"""
return math.sqrt(sum_flat(absolute(a)**2))
def norm(a,p=2):
"""norm(a,p=2) -> l-p norm of a.flat
Return the l-p norm of a, considered as a flat array. This is NOT a true
matrix norm, since arrays of arbitrary rank are always flattened.
p can be a number or one of the strings ('inf','Infinity') to get the
L-infinity norm.
Ref: http://mathworld.wolfram.com/VectorNorm.html
http://mathworld.wolfram.com/L-Infinity-Norm.html"""
if p in ('inf','Infinity'):
return max(absolute(a).flat)
else:
return (sum_flat(absolute(a)**p))**(1.0/p)
def frange(xini,xfin=None,delta=None,**kw):
"""frange([start,] stop[, step, keywords]) -> array of floats
Return a Numeric array() containing a progression of floats. Similar to
arange(), but defaults to a closed interval.
frange(x0, x1) returns [x0, x0+1, x0+2, ..., x1]; start defaults to 0, and
the endpoint *is included*. This behavior is different from that of
range() and arange(). This is deliberate, since frange will probably be
more useful for generating lists of points for function evaluation, and
endpoints are often desired in this use. The usual behavior of range() can
be obtained by setting the keyword 'closed=0', in this case frange()
basically becomes arange().
When step is given, it specifies the increment (or decrement). All
arguments can be floating point numbers.
frange(x0,x1,d) returns [x0,x0+d,x0+2d,...,xfin] where xfin<=x1.
frange can also be called with the keyword 'npts'. This sets the number of
points the list should contain (and overrides the value 'step' might have
been given). arange() doesn't offer this option.
Examples:
>>> frange(3)
array([ 0., 1., 2., 3.])
>>> frange(3,closed=0)
array([ 0., 1., 2.])
>>> frange(1,6,2)
array([1, 3, 5])
>>> frange(1,6.5,npts=5)
array([ 1. , 2.375, 3.75 , 5.125, 6.5 ])
"""
#defaults
kw.setdefault('closed',1)
endpoint = kw['closed'] != 0
# funny logic to allow the *first* argument to be optional (like range())
# This was modified with a simpler version from a similar frange() found
# at http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/66472
if xfin == None:
xfin = xini + 0.0
xini = 0.0
if delta == None:
delta = 1.0
# compute # of points, spacing and return final list
try:
npts=kw['npts']
delta=(xfin-xini)/float(npts-endpoint)
except KeyError:
# round() gets npts right even with the vagaries of floating point.
npts=int(round((xfin-xini)/delta+endpoint))
return arange(npts)*delta+xini
def diagonal_matrix(diag):
"""Return square diagonal matrix whose non-zero elements are given by the
input array."""
return diag*identity(len(diag))
def identity(n,rank=2,typecode='l'):
"""identity(n,r) returns the identity matrix of shape (n,n,...,n) (rank r).
For ranks higher than 2, this object is simply a multi-index Kronecker
delta:
/ 1 if i0=i1=...=iR,
id[i0,i1,...,iR] = -|
\ 0 otherwise.
Optionally a typecode may be given (it defaults to 'l').
Since rank defaults to 2, this function behaves in the default case (when
only n is given) like the Numeric identity function."""
iden = zeros((n,)*rank,typecode=typecode)
for i in range(n):
idx = (i,)*rank
iden[idx] = 1
return iden
def base_repr (number, base = 2, padding = 0):
"""Return the representation of a number in any given base."""
chars = '0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ'
if number < base: \
return (padding - 1) * chars [0] + chars [int (number)]
max_exponent = int (math.log (number)/math.log (base))
max_power = long (base) ** max_exponent
lead_digit = int (number/max_power)
return chars [lead_digit] + \
base_repr (number - max_power * lead_digit, base, \
max (padding - 1, max_exponent))
def binary_repr(number, max_length = 1025):
"""Return the binary representation of the input number as a string.
This is more efficient than using base_repr with base 2.
Increase the value of max_length for very large numbers. Note that on
32-bit machines, 2**1023 is the largest integer power of 2 which can be
converted to a Python float."""
assert number < 2L << max_length
shifts = map (operator.rshift, max_length * [number], \
range (max_length - 1, -1, -1))
digits = map (operator.mod, shifts, max_length * [2])
if not digits.count (1): return 0
digits = digits [digits.index (1):]
return ''.join (map (repr, digits)).replace('L','')
def log2(x,ln2 = math.log(2.0)):
"""Return the log(x) in base 2.
This is a _slow_ function but which is guaranteed to return the correct
integer value if the input is an ineger exact power of 2."""
try:
bin_n = binary_repr(x)[1:]
except (AssertionError,TypeError):
return math.log(x)/ln2
else:
if '1' in bin_n:
return math.log(x)/ln2
else:
return len(bin_n)
def ispower2(n):
"""Returns the log base 2 of n if n is a power of 2, zero otherwise.
Note the potential ambiguity if n==1: 2**0==1, interpret accordingly."""
bin_n = binary_repr(n)[1:]
if '1' in bin_n:
return 0
else:
return len(bin_n)
def fromfunction_kw(function, dimensions, **kwargs):
"""Drop-in replacement for fromfunction() from Numerical Python.
Allows passing keyword arguments to the desired function.
Call it as (keywords are optional):
fromfunction_kw(MyFunction, dimensions, keywords)
The function MyFunction() is responsible for handling the dictionary of
keywords it will recieve."""
return function(tuple(indices(dimensions)),**kwargs)
#**************************** end file <numutils.py> ************************
You can’t perform that action at this time.