Newer
Older
import math
import sys
import operator
import types
from pytools.decorator import decorator
Andreas Klöckner
committed
Andreas Klöckner
committed
def delta(x,y):
if x == y:
return 1
else:
return 0
Andreas Klöckner
committed
def factorial(n):
from operator import mul
Andreas Klöckner
committed
assert n == int(n)
return reduce(mul, (i for i in xrange(1,n+1)), 1)
Andreas Klöckner
committed
# Data structures ------------------------------------------------------------
class Reference(object):
def __init__( self, value ):
self.V = value
def get( self ):
return self.V
def set( self, value ):
self.V = value
Andreas Klöckner
committed
def _getattr_(obj, name, default_thunk):
"Similar to .setdefault in dictionaries."
try:
return getattr(obj, name)
except AttributeError:
default = default_thunk()
setattr(obj, name, default)
return default
@decorator
def memoize(func, *args):
# by Michele Simionato
# http://www.phyast.pitt.edu/~micheles/python/
Andreas Klöckner
committed
dic = _getattr_(func, "memoize_dic", dict)
# memoize_dic is created at the first call
if args in dic:
return dic[args]
else:
result = func(*args)
dic[args] = result
return result
FunctionValueCache = memoize
Andreas Klöckner
committed
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
class DictionaryWithDefault(object):
def __init__(self, default_value_generator, start = {}):
self._Dictionary = dict(start)
self._DefaultGenerator = default_value_generator
def __getitem__(self, index):
try:
return self._Dictionary[index]
except KeyError:
value = self._DefaultGenerator(index)
self._Dictionary[index] = value
return value
def __setitem__(self, index, value):
self._Dictionary[index] = value
def __contains__(self, item):
return True
def iterkeys(self):
return self._Dictionary.iterkeys()
def __iter__(self):
return self._Dictionary.__iter__()
def iteritems(self):
return self._Dictionary.iteritems()
class FakeList(object):
def __init__(self, f, length):
self._Length = length
self._Function = f
def __len__(self):
return self._Length
def __getitem__(self, index):
try:
return [self._Function(i)
for i in range(*index.indices(self._Length))]
except AttributeError:
return self._Function(index)
class DependentDictionary(object):
def __init__(self, f, start = {}):
self._Function = f
self._Dictionary = start.copy()
def copy(self):
return DependentDictionary(self._Function, self._Dictionary)
def __contains__(self, key):
try:
self[key]
return True
except KeyError:
return False
def __getitem__(self, key):
try:
return self._Dictionary[key]
except KeyError:
return self._Function(self._Dictionary, key)
def __setitem__(self, key, value):
self._Dictionary[key] = value
def genuineKeys(self):
return self._Dictionary.keys()
def iteritems(self):
return self._Dictionary.iteritems()
def iterkeys(self):
return self._Dictionary.iterkeys()
def itervalues(self):
return self._Dictionary.itervalues()
def add_tuples(t1, t2):
return tuple([t1v + t2v for t1v, t2v in zip(t1, t2)])
def negate_tuple(t1):
return tuple([-t1v for t1v in t1])
def write_1d_gnuplot_graph(f, a, b, steps=100, fname=",,f.data", progress = False):
h = float(b - a)/steps
gnuplot_file = file(fname, "w")
def do_plot(func):
for n in range(steps):
if progress:
sys.stdout.write(".")
sys.stdout.flush()
x = a + h * n
gnuplot_file.write("%f\t%f\n" % (x, func(x)))
do_plot(f)
if progress:
sys.stdout.write("\n")
def write_1d_gnuplot_graphs(f, a, b, steps=100, fnames=None, progress=False):
h = float(b - a)/steps
if not fnames:
result_count = len(f(a))
fnames = [",,f%d.data" % i for i in range(result_count)]
gnuplot_files = [file(fname, "w") for fname in fnames]
for n in range(steps):
if progress:
sys.stdout.write(".")
sys.stdout.flush()
x = a + h * n
for gpfile, y in zip(gnuplot_files, f(x)):
gpfile.write("%f\t%f\n" % (x, y))
if progress:
sys.stdout.write("\n")
Andreas Tester
committed
def write_2d_gnuplot_graph(f, (x0, y0), (x1, y1), (xsteps, ysteps)=(100, 100), fname=",,f.data"):
hx = float(x1 - x0)/xsteps
hy = float(y1 - y0)/ysteps
gnuplot_file = file(fname, "w")
for ny in range(ysteps):
for nx in range(xsteps):
x = x0 + hx * nx
y = y0 + hy * ny
gnuplot_file.write("%g\t%g\t%g\n" % (x, y, f(x, y)))
gnuplot_file.write("\n")
Andreas Klöckner
committed
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
def write_gnuplot_graph(f, a, b, steps = 100, fname = ",,f.data", progress = False):
h = float(b - a)/steps
gnuplot_file = file(fname, "w")
def do_plot(func):
for n in range(steps):
if progress:
sys.stdout.write(".")
sys.stdout.flush()
x = a + h * n
gnuplot_file.write("%f\t%f\n" % (x, func(x)))
if isinstance(f, types.ListType):
for f_index, real_f in enumerate(f):
if progress:
sys.stdout.write("function %d: " % f_index)
do_plot(real_f)
gnuplot_file.write("\n")
if progress:
sys.stdout.write("\n")
else:
do_plot(f)
if progress:
sys.stdout.write("\n")
Andreas Klöckner
committed
# syntactical sugar -----------------------------------------------------------
class InfixOperator:
"""Pseudo-infix operators that allow syntax of the kind `op1 <<operator>> op2'.
Following a recipe from
http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/384122
"""
def __init__(self, function):
self.function = function
def __rlshift__(self, other):
return InfixOperator(lambda x: self.function(other, x))
def __rshift__(self, other):
return self.function(other)
def call(self, a, b):
return self.function(a, b)
Andreas Klöckner
committed
# Generic utilities ----------------------------------------------------------
def flatten(list):
"""For an iterable of sub-iterables, generate each member of each
sub-iterable in turn, i.e. a flattened version of that super-iterable.
Example: Turn [[a,b,c],[d,e,f]] into [a,b,c,d,e,f].
"""
for sublist in list:
for j in sublist:
yield j
Andreas Klöckner
committed
def sum_over(function, arguments):
raise RuntimeError, "Horribly inefficient routine called."
# wherever this is used, it should be replaced by sum() and a generator
# expression.
Andreas Klöckner
committed
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
result = 0
for i in arguments:
result += function(i)
return result
def general_sum(sequence):
return reduce(operator.add, sequence)
def linear_combination(coefficients, vectors):
result = coefficients[0] * vectors[0]
for c,v in zip(coefficients, vectors)[1:]:
result += c*v
return result
def average(sequence):
return general_sum(sequence)/float(len(sequence))
def all_equal(sequence):
item = sequence[0]
for i in sequence[1:]:
if i != item:
return False
return True
def decorate(function, list):
return map(lambda x: (x, function(x)), list)
def partition(criterion, list):
part_true = []
part_false = []
for i in list:
if criterion(i):
part_true.append(i)
else:
part_false.append(i)
return part_true, part_false
def product(list):
return reduce(lambda x,y: x*y, list, 1)
def argmin_f(list, f = lambda x: x):
# deprecated -- the function has become unnecessary because of
# generator expressions
Andreas Klöckner
committed
current_min_index = -1
current_min = f(list[0])
for idx, item in enumerate(list[1:]):
value = f(item)
if value < current_min:
current_min_index = idx
current_min = value
return current_min_index+1
def argmax_f(list, f = lambda x: x):
# deprecated -- the function has become unnecessary because of
# generator expressions
Andreas Klöckner
committed
current_max_index = -1
current_max = f(list[0])
for idx, item in enumerate(list[1:]):
value = f(item)
if value > current_max:
current_max_index = idx
current_max = value
return current_max_index+1
def argmin(list):
current_min_index = -1
it = list.__iter__()
current_min = it.next()
for idx, item in enumerate(it):
value = item
if value < current_min:
current_min_index = idx
current_min = value
return current_min_index+1
def argmax(list):
it = list.__iter__()
current_max = it.next()
for idx, item in enumerate(it):
value = item
if value > current_max:
current_max_index = idx
current_max = value
return current_max_index+1
Andreas Klöckner
committed
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
def cartesian_product(list1, list2):
result = []
for i in list1:
for j in list2:
result.append((i,j))
def cartesian_product_sum(list1, list2):
"""This routine returns a list of sums of each element of
list1 with each element of list2. Also works with lists.
"""
result = []
for i in list1:
for j in list2:
result.append(i+j)
return result
def reverse_dictionary(the_dict):
result = {}
for key, value in the_dict.iteritems():
if value in result:
raise RuntimeError, "non-reversible mapping"
result[value] = key
return result
def wandering_element(length, wanderer=1, landscape=0):
for i in range(length):
yield i*(landscape,) + (wanderer,) + (length-1-i)*(landscape,)
Andreas Klöckner
committed
def generate_nonnegative_integer_tuples_below(n, length, least=0):
Andreas Klöckner
committed
assert length >= 0
if length == 0:
Andreas Klöckner
committed
yield ()
Andreas Klöckner
committed
else:
for i in range(least, n):
Andreas Klöckner
committed
for base in generate_nonnegative_integer_tuples_below(n, length-1, least):
yield (i,) + base
def generate_decreasing_nonnegative_tuples_summing_to(n, length, min=0, max=None):
sig = (n,length,max)
if length == 0:
yield ()
elif length == 1:
if n <= max:
#print "MX", n, max
yield (n,)
else:
return
else:
if max is None or n < max:
max = n
for i in range(min, max+1):
#print "SIG", sig, i
for remainder in generate_decreasing_nonnegative_tuples_summing_to(
n-i, length-1, min, i):
yield (i,) + remainder
Andreas Klöckner
committed
Andreas Klöckner
committed
def generate_nonnegative_integer_tuples_summing_to_at_most(n, length):
"""Enumerate all non-negative integer tuples summing to at most n,
exhausting the search space by varying the first entry fastest,
and the last entry the slowest.
"""
assert length >= 0
if length == 0:
yield ()
else:
for i in range(n+1):
Andreas Klöckner
committed
for remainder in generate_nonnegative_integer_tuples_summing_to_at_most(
n-i, length-1):
yield remainder + (i,)
Andreas Klöckner
committed
def generate_all_nonnegative_integer_tuples(length, least=0):
Andreas Klöckner
committed
assert length >= 0
current_max = least
while True:
for max_pos in range(length):
Andreas Klöckner
committed
for prebase in generate_nonnegative_integer_tuples_below(current_max, max_pos, least):
for postbase in generate_nonnegative_integer_tuples_below(current_max+1, length-max_pos-1, least):
Andreas Klöckner
committed
yield prebase + [current_max] + postbase
current_max += 1
Andreas Klöckner
committed
# backwards compatibility
generate_positive_integer_tuples_below = generate_nonnegative_integer_tuples_below
generate_all_positive_integer_tuples = generate_all_nonnegative_integer_tuples
Andreas Klöckner
committed
def _pos_and_neg_adaptor(tuple_iter):
for tup in tuple_iter:
nonzero_indices = [i for i in range(len(tup)) if tup[i] != 0]
Andreas Klöckner
committed
for do_neg_tup in generate_nonnegative_integer_tuples_below(2, len(nonzero_indices)):
Andreas Klöckner
committed
this_result = list(tup)
for index, do_neg in enumerate(do_neg_tup):
if do_neg:
this_result[nonzero_indices[index]] *= -1
yield tuple(this_result)
def generate_all_integer_tuples_below(n, length, least_abs=0):
Andreas Klöckner
committed
return _pos_and_neg_adaptor(generate_nonnegative_integer_tuples_below(
Andreas Klöckner
committed
n, length, least_abs))
def generate_all_integer_tuples(length, least_abs=0):
Andreas Klöckner
committed
return _pos_and_neg_adaptor(generate_all_nonnegative_integer_tuples(
Andreas Klöckner
committed
length, least_abs))
Andreas Klöckner
committed
def generate_permutations(original):
"""Generate all permutations of the list `original'.
Nicked from http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/252178
"""
if len(original) <=1:
yield original
else:
for perm in generate_permutations(original[1:]):
for i in range(len(perm)+1):
#nb str[0:1] works in both string and list contexts
yield perm[:i] + original[0:1] + perm[i:]
Andreas Klöckner
committed
Andreas Klöckner
committed
def generate_unique_permutations(original):
"""Generate all unique permutations of the list `original'.
"""
had_those = set()
for perm in generate_permutations(original):
if perm not in had_those:
had_those.add(perm)
yield perm
class Table:
"""An ASCII table generator."""
def __init__(self):
self.Rows = []
def add_row(self, row):
self.Rows.append([str(i) for i in row])
def __str__(self):
columns = len(self.Rows[0])
col_widths = [max(len(row[i]) for row in self.Rows)
for i in range(columns)]
lines = [
"|".join([cell.ljust(col_width)
for cell, col_width in zip(row, col_widths)])
for row in self.Rows]
lines[1:1] = ["+".join("-"*col_width
for col_width in col_widths)]
return "\n".join(lines)
Andreas Klöckner
committed
# Obscure stuff --------------------------------------------------------------
def enumerate_basic_directions(dimensions):
coordinate_list = [[0], [1], [-1]]
return reduce(cartesian_product_sum, [coordinate_list] * dimensions)[1:]
Andreas Klöckner
committed
# file system related ---------------------------------------------------------
def assert_not_a_file(name):
import os
if os.access(name, os.F_OK):
raise IOError, "file `%s' already exists" % name