00001
00002 import sys
00003 import os
00004 import pickle
00005 import math
00006 import numpy
00007 from os.path import join as pjoin
00008 import optparse
00009 from glob import glob
00010
00011 from ufl.algorithms import load_forms
00012
00013 import ufc
00014 import ufc_benchmark
00015
00016
00017
00018
00019
00020
00021
00022 inf = 1e9999
00023 nan = inf*0
00024
00025
00026 from subprocess import Popen, PIPE, STDOUT
00027 def get_status_output(cmd, input=None, cwd=None, env=None):
00028 pipe = Popen(cmd, shell=True, cwd=cwd, env=env, stdout=PIPE, stderr=STDOUT)
00029 (output, errout) = pipe.communicate(input=input)
00030 assert not errout
00031 status = pipe.returncode
00032 return (status, output)
00033
00034 def runcmd(cmd):
00035 get_status_output(cmd)
00036
00037 def write_file(filename, text):
00038 "Write text to a file and close it."
00039 f = open(filename, "w")
00040 f.write(text)
00041 f.close()
00042 print "Wrote file '%s'" % filename
00043
00044
00045
00046 usage = """Compile UFL forms, compute element tensors, and compare with reference values.
00047
00048 Examples:
00049
00050 FIXME: Write usage examples.
00051 """
00052
00053 def opt(long, short, t, default, help):
00054 return optparse.make_option("--%s" % long, "-%s" % short, action="store", type=t, dest=long, default=default, help=help)
00055
00056 option_list = [ \
00057
00058 opt("ufldir", "u", "str", "ufl", "Input directory with .ufl files."),
00059 opt("outputdir", "o", "str", "output", "Output directory to write .ref files to."),
00060 opt("referencedir", "r", "str", "reference", "Reference directory to read .ref files from."),
00061 opt("cachedir", "c", "str", "cache", "Reference directory to read .ref files from."),
00062
00063 opt("skip", "s", "str", "", "Comma-separated list of ufl files to skip."),
00064 opt("write", "w", "int", 0, "Write new reference files."),
00065 opt("debug", "d", "int", 0, "Print debugging info for this script."),
00066
00067 opt("jit_options", "j", "str", "options/default.py", "Python file containing jit options."),
00068
00069 opt("tolerance", "t", "float", 1e-10, "Compare norm of data difference to this tolerance."),
00070 opt("norm", "n", "int", 1, "Compare data with references using the L2 norm of tensor difference."),
00071 opt("eig", "e", "int", 1, "Compare data with references using the eigenvalues."),
00072 opt("random_cell", "a", "int", 1, "Use a (predefined) random cell instead of reference cell."),
00073 opt("benchmark", "b", "int", 1, "Measure the time to call tabulate_tensor."),
00074 ]
00075
00076
00077
00078 def main(args):
00079 "Handle commandline arguments and orchestrate the tests."
00080
00081
00082 parser = optparse.OptionParser(usage=usage, option_list=option_list)
00083 (options, args) = parser.parse_args(args=args)
00084 if args:
00085 print "ERROR: Got additional unknown arguments: ", args
00086 parser.print_usage()
00087 return -1
00088
00089
00090 skip = set(s.strip() for s in options.skip.split(","))
00091 def skipmatch(fn):
00092 if fn in skip: return True
00093 path, name = os.path.split(fn)
00094 if name in skip: return True
00095 basename, ext = os.path.splitext(name)
00096 if basename in skip: return True
00097 return False
00098 uflfiles = glob(pjoin(options.ufldir, "*.ufl"))
00099 uflfiles = [f for f in uflfiles if not skipmatch(f)]
00100 uflfiles = sorted(uflfiles)
00101
00102 if options.debug:
00103 print "."*40
00104 print "Got uflfiles ="
00105 print "\n".join(" " + f for f in uflfiles)
00106
00107
00108 fails = []
00109 passes = []
00110 summaries = []
00111 for filename in uflfiles:
00112 summary, ok = handle_file(filename, options)
00113 summaries.append(summary)
00114 if ok:
00115 passes.append(filename)
00116 else:
00117 fails.append(filename)
00118
00119
00120 print
00121 print "="*80
00122 print "Summaries:",
00123 sep = "\n\n" + "-"*60 + "\n"
00124 print sep + sep.join(summaries)
00125
00126
00127 if passes:
00128 print "="*80
00129 print "The following files passed:"
00130 print "\n".join(" " + f for f in sorted(passes))
00131 if fails:
00132 print "="*80
00133 print "The following files failed:"
00134 print "\n".join(" " + f for f in sorted(fails))
00135
00136 def import_options_iterator(name):
00137 "Import options iterator from a python file."
00138 assert os.path.exists(name)
00139
00140 path, name = os.path.split(name)
00141 basename, dotpy = os.path.splitext(name)
00142 assert dotpy in (".py", "")
00143
00144 sys.path.insert(0, path)
00145
00146
00147 try:
00148 options_module = __import__(basename)
00149 iter_jit_options = options_module.options
00150 finally:
00151 sys.path.pop(0)
00152
00153
00154 return iter_jit_options
00155
00156 def handle_file(filename, options):
00157 "Handle all aspects of testing a single .ufl file."
00158
00159
00160 uflfilename = filename
00161 path, name = os.path.split(uflfilename)
00162 basename, ext = os.path.splitext(name)
00163 if ext != ".ufl":
00164 msg = "Expecting a .ufl file, not %s" % uflfilename
00165 return (msg, False)
00166
00167 if options.debug:
00168 print "."*40
00169 print "In handle_file, filename parsed as:"
00170 print "filename =", filename
00171 print "uflfilename =", uflfilename
00172 print "path =", path
00173 print "name =", name
00174 print "basename =", basename
00175 print "ext =", ext
00176
00177
00178 try:
00179 forms = load_forms(uflfilename)
00180 except:
00181 msg = "Failed to load file, try running\n\n"\
00182 " ufl-analyse %s\n\n"\
00183 "to find the bug in the form file or in UFL." % uflfilename
00184 return (msg, False)
00185
00186 formnames = [form.form_data().name for form in forms]
00187
00188 if options.debug:
00189 print "."*40
00190 print "In handle_file, forms loaded: ", ", ".join(formnames)
00191
00192
00193 iter_jit_options = import_options_iterator(options.jit_options)
00194 if iter_jit_options is None:
00195 msg = "Failed to import options module '%s'." % options.jit_options
00196 return (msg, False)
00197
00198 total_ok = True
00199 summary = ""
00200 for jit, jit_options in iter_jit_options():
00201
00202
00203
00204
00205 try:
00206 jit_result = jit(forms, jit_options)
00207 except:
00208 msg = "Failed to jit-compile forms in file '%s'." % uflfilename
00209 raise
00210
00211
00212 if options.debug:
00213 print ":"*60
00214 print "Jit result:"
00215 print jit_result
00216
00217
00218 data = {}
00219
00220
00221 compiled_forms, module, form_datas = jit_result
00222 for i in range(len(forms)):
00223 form_data = forms[i].form_data()
00224 assert form_data is form_datas[i]
00225 data[form_data.name] = compute_data(compiled_forms[i], form_data, options.random_cell)
00226
00227
00228 benchmark_data = {}
00229 if options.benchmark:
00230 for i in range(len(forms)):
00231 form_data = forms[i].form_data()
00232 assert form_data is form_datas[i]
00233 compiled_form = compiled_forms[i]
00234 result = ufc_benchmark.benchmark_forms([compiled_form], False)
00235 benchmark_data[form_data.name] = result
00236
00237
00238 if options.write:
00239 outputdir = options.referencedir
00240 else:
00241 outputdir = options.outputdir
00242 if options.debug:
00243 print "."*60
00244 print "outputdir =", outputdir
00245
00246 for formname in formnames:
00247 outputfilename = pjoin(outputdir, "%s_%s.ref" % (basename, formname))
00248 if options.debug:
00249 print "Writing to output filename: ", outputfilename
00250 if not data[formname].any():
00251 print "*** Warning: reference tensor only contains zeros!"
00252 write_data(outputfilename, data[formname])
00253
00254 if options.benchmark:
00255 for formname in formnames:
00256 outputfilename = pjoin(outputdir, "%s_%s.timing" % (basename, formname))
00257 if options.debug:
00258 print "Writing to output filename: ", outputfilename
00259 write_data(outputfilename, benchmark_data[formname])
00260
00261
00262 if not options.write:
00263
00264 reference = {}
00265 for formname in formnames:
00266 referencefilename = pjoin(options.referencedir, "%s_%s.ref" % (basename, formname))
00267 if options.debug:
00268 print "Read reference filename: ", referencefilename
00269 reference[formname] = read_data(referencefilename)
00270
00271 if options.benchmark:
00272 benchmark_reference = {}
00273 for formname in formnames:
00274 referencefilename = pjoin(options.referencedir, "%s_%s.timing" % (basename, formname))
00275 if options.debug:
00276 print "Read reference filename: ", referencefilename
00277 benchmark_reference[formname] = read_data(referencefilename)
00278
00279
00280 results = []
00281 for formname in formnames:
00282 msg, ok = compare_data(data[formname], reference[formname], options)
00283
00284 if options.debug:
00285 print "msg =", msg
00286 print "ok = ", ok
00287
00288 if not ok:
00289 total_ok = False
00290 results.append("--- Errors in form %s:\n" % formname)
00291 results.append(msg)
00292 else:
00293 results.append("--- Form %s is ok.\n" % formname)
00294
00295 if options.benchmark:
00296 msg = compare_benchmark(benchmark_data[formname], benchmark_reference[formname], options)
00297 results.extend(msg)
00298
00299 summary += "\n".join(results)
00300
00301
00302 if total_ok:
00303 summary = "%s passed everything." % uflfilename
00304 else:
00305 summary = ("%s has problems:\n" % uflfilename) + summary
00306 if options.write:
00307 summary += "\n\nWrote reference files for %s." % uflfilename
00308 return (summary, total_ok)
00309
00310 def compute_diff_norm(data, reference, options):
00311 if options.debug:
00312 print ":"*40
00313 print "In compute_diff_norm, comparing:"
00314 print "data ="
00315 print data
00316 print "reference ="
00317 print reference
00318
00319
00320 diff = data - reference
00321
00322
00323 d2 = numpy.sum(diff**2)
00324 r2 = numpy.sum(reference**2)
00325 n2 = numpy.sum(data**2)
00326
00327
00328 norm = math.sqrt(d2 / r2)
00329
00330 if options.debug:
00331 print "diff ="
00332 print diff
00333 print "d2, r2, n2="
00334 print d2, r2, n2
00335
00336
00337
00338 return norm
00339
00340 def compute_eigenvalues(data):
00341 sh = data.shape
00342 assert sh[0] == sh[1]
00343 from scipy.linalg import eig
00344 l, v = eig(data)
00345 return numpy.array(sorted(l))
00346
00347 def compare_data(data, reference, options):
00348 norm = nan
00349 eig = nan
00350 msg = ""
00351 if reference is None:
00352 total_ok = False
00353 msg += "No reference to compare to."
00354 else:
00355 total_ok = True
00356
00357 if data.shape != reference.shape:
00358 total_ok = False
00359 msg += "\n ERROR: Data shape is %s, reference shape is %s." % (data.shape, reference.shape)
00360 else:
00361 if options.norm:
00362 norm = compute_diff_norm(data, reference, options)
00363 ok = norm < options.tolerance
00364 if not ok:
00365 total_ok = False
00366 msg += "\n norm = %e >= %e = tol" % (norm, options.tolerance)
00367
00368 if options.eig:
00369 sh = data.shape
00370
00371
00372 if len(sh) == 1:
00373
00374 eig1 = numpy.array(sorted(data))
00375 eig2 = numpy.array(sorted(reference))
00376 eig = sum((eig1-eig2)**2)
00377 ok = eig < options.tolerance
00378 if not ok:
00379 total_ok = False
00380 msg += "\n eig = %e >= %e = tol" % (eig, options.tolerance)
00381
00382 elif sh[0] == sh[1]:
00383
00384 eig1 = compute_eigenvalues(data)
00385 eig2 = compute_eigenvalues(reference)
00386 eig = sum((eig1-eig2)**2)
00387 ok = eig < options.tolerance
00388 if not ok:
00389 total_ok = False
00390 msg += "\n eig = %e >= %e = tol" % (eig, options.tolerance)
00391 else:
00392 if not options.norm:
00393 total_ok = False
00394 msg += "\n WARNING: Not computing eigenvalues of data with shape %s" % repr(sh)
00395
00396
00397 if total_ok:
00398 msg = "Data compared ok."
00399 else:
00400 msg = "Failed because:%s\n\ndata is\n%r\n\nreference is\n%r" % (msg, data, reference)
00401 return (msg, total_ok)
00402
00403 def compare_benchmark(data, reference, options):
00404 msg = []
00405
00406 for (b1, b2) in zip(data, reference):
00407
00408 for (x, y) in zip(b1, b2):
00409
00410 for (r,s) in zip(x,y):
00411 f = r/s
00412 if f > 1:
00413 msg.append("The reference is faster than the new, f = %.2f" % f)
00414 else:
00415 msg.append("The new is faster than the reference, f = %.2f" % f)
00416 return msg
00417
00418 def write_data(fn, data):
00419 try:
00420 f = open(fn, "w")
00421 pickle.dump(data, f)
00422
00423 f.close()
00424 except Exception, what:
00425 print "*** An error occured while trying to write reference file: %s" % fn
00426 raise
00427
00428 def read_data(fn):
00429 try:
00430 f = open(fn)
00431 data = pickle.load(f)
00432
00433 f.close()
00434 except Exception, what:
00435 print "*** An error occured while trying to load reference file: %s" % fn
00436 print "*** Maybe you need to generate the reference? Returning None."
00437 data = None
00438 return data
00439
00440 def make_mesh(cell_shape, random_cell):
00441
00442
00443 if cell_shape == "interval":
00444 mesh = ufc_benchmark.Mesh(1, 1, [2, 1, 0])
00445 if random_cell:
00446 cell = ufc_benchmark.Cell(1, 1, [[-1.445], [0.4713]], [2, 1, 0, 0])
00447 else:
00448 cell = ufc_benchmark.Cell(1, 1, [[0], [1]], [2, 1, 0, 0])
00449 elif cell_shape == "triangle":
00450 mesh = ufc_benchmark.Mesh(2, 2, [3, 3, 1])
00451 if random_cell:
00452 cell = ufc_benchmark.Cell(2, 2, [[-2.2304, -0.88317], [1.3138, -1.0164],\
00453 [0.24622, 1.4431]], [3, 3, 1, 0])
00454 else:
00455 cell = ufc_benchmark.Cell(2, 2, [[0, 0], [1, 0], [1, 1]], [3, 3, 1, 0])
00456 elif cell_shape == "tetrahedron":
00457 mesh = ufc_benchmark.Mesh(3, 3, [4, 6, 4, 1])
00458 if random_cell:
00459 cell = ufc_benchmark.Cell(3, 3, [[-2.2561, -1.6144, -1.7349], [-1.5612, -1.5121, -0.17675],\
00460 [1.6861, -1.1494, 2.4070], [0.52083, 1.1940, 1.8220]], [4, 6, 4, 1])
00461 else:
00462 cell = ufc_benchmark.Cell(3, 3, [[0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 1]], [4, 6, 4, 1])
00463 else:
00464 raise RuntimeError, "Unknown shape " + cell_shape
00465 return mesh, cell
00466
00467 def compute_data(compiled_form, form_data, random_cell):
00468
00469 rank = form_data.rank
00470 num_coefficients = form_data.num_functions
00471 num_arguments = num_coefficients + rank
00472
00473 num_cell_integrals = compiled_form.num_cell_integrals()
00474 num_exterior_facet_integrals = compiled_form.num_exterior_facet_integrals()
00475 num_interior_facet_integrals = compiled_form.num_interior_facet_integrals()
00476
00477
00478 mesh, cell = make_mesh(form_data.cell.domain(), random_cell)
00479 dim = form_data.geometric_dimension
00480 num_facets = cell.num_entities[dim - 1]
00481
00482
00483 dof_maps = [0]*num_arguments
00484 for i in range(num_arguments):
00485 dof_maps[i] = compiled_form.create_dof_map(i)
00486 dof_maps[i].init_mesh(mesh)
00487
00488
00489 w = [0]*num_coefficients
00490 for i in range(num_coefficients):
00491 w[i] = [0]*(dof_maps[rank+i].local_dimension())
00492 for j in range(dof_maps[rank+i].local_dimension()):
00493 w[i][j] = 1.111 + (i + j)/1.111
00494 macro_w = [0]*num_coefficients
00495 for i in range(num_coefficients):
00496 macro_w[i] = [0]*(2*dof_maps[rank+i].local_dimension())
00497 for j in range(2*dof_maps[rank+i].local_dimension()):
00498 macro_w[i][j] = 1.111 + (i + j)/1.111
00499
00500
00501 A = numpy.zeros((1,1))
00502
00503
00504 if num_cell_integrals:
00505 domain = 0
00506
00507 try:
00508 A = ufc_benchmark.tabulate_cell_integral(compiled_form, w, cell, domain)
00509 A = numpy.array(A)
00510 A = numpy.zeros(numpy.shape(A))
00511 for domain in range(num_cell_integrals):
00512 A += ufc_benchmark.tabulate_cell_integral(compiled_form, w, cell, domain)
00513 except Exception, what:
00514 print "*** An error occured while calling tabulate_cell_integral() for domain %d." % domain
00515 raise
00516
00517
00518 if num_exterior_facet_integrals:
00519 domain, facet = (0, 0)
00520 try:
00521 if not numpy.any(A):
00522 A = ufc_benchmark.tabulate_exterior_facet_integral(compiled_form, w, cell, facet, domain)
00523 A = numpy.array(A)
00524 A = numpy.zeros(numpy.shape(A))
00525 for domain in range(num_exterior_facet_integrals):
00526 for facet in range(num_facets):
00527 A += ufc_benchmark.tabulate_exterior_facet_integral(compiled_form, w, cell, facet, domain)
00528 except Exception, what:
00529 print "*** An error occured while calling tabulate_exterior_facet_integral() for domain %d, facet %d." % (domain, facet)
00530 raise
00531
00532
00533
00534
00535 macro_A = numpy.array([0.0])
00536 if num_interior_facet_integrals:
00537 domain, facet0, facet1 = (0,0,0)
00538 try:
00539 macro_A = ufc_benchmark.tabulate_interior_facet_integral(compiled_form, macro_w, cell, cell, facet0, facet1, domain)
00540 macro_A = numpy.array(macro_A)
00541 macro_A = numpy.zeros(numpy.shape(macro_A))
00542 for domain in range(num_interior_facet_integrals):
00543 for facet0 in range(num_facets):
00544 for facet1 in range(num_facets):
00545 macro_A += ufc_benchmark.tabulate_interior_facet_integral(compiled_form, macro_w, cell, cell, facet0, facet1, domain)
00546 except Exception, what:
00547 print "*** An error occured while calling tabulate_interior_facet_integral() for domain %d, facet0 %d, facet1 %d." % (domain, facet0, facet1)
00548 raise
00549
00550
00551
00552 if not macro_A.any():
00553 data = A
00554 elif A.any():
00555 data = macro_A
00556 dims = A.shape
00557 data[:dims[0], :dims[1]] += A
00558
00559 return data
00560
00561
00562
00563 if __name__ == "__main__":
00564 result = main(sys.argv[1:])
00565 sys.exit(result)