SyFi 0.3
|
00001 #!/usr/bin/env python 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 #import sfc 00017 #import logging 00018 #sfc.set_logging_level(logging.DEBUG) 00019 00020 # --- Basic Utilities --- 00021 00022 inf = 1e9999 00023 nan = inf*0 00024 00025 # Taken from http://ivory.idyll.org/blog/mar-07/replacing-commands-with-subprocess 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 # --- Option Parsing Data --- 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 # Directories: 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 # Main behaviour options: 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 # Form compiler options: 00067 opt("jit_options", "j", "str", "options/default.py", "Python file containing jit options."), 00068 # Comparison options: 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 # --- Main routine --- 00077 00078 def main(args): 00079 "Handle commandline arguments and orchestrate the tests." 00080 00081 # Parse commandline arguments 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 # Read input directory and filter filenames 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 # Handle each .ufl file separately 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 # Print summaries 00120 print 00121 print "="*80 00122 print "Summaries:", 00123 sep = "\n\n" + "-"*60 + "\n" 00124 print sep + sep.join(summaries) 00125 00126 # Print files that passed and failed 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 #cwd = os.getcwd() 00146 #os.chdir(path) 00147 try: 00148 options_module = __import__(basename) 00149 iter_jit_options = options_module.options 00150 finally: 00151 sys.path.pop(0) 00152 #os.chdir(cwd) 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 # Split filename 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 # Load forms from this file 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 # Iterate over a given set of jit compilers and options: 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 #jit_options.cache_dir = options.cache_dir # FIXME 00203 00204 # Compile forms with given compiler and options 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 #return (msg, False) 00211 00212 if options.debug: 00213 print ":"*60 00214 print "Jit result:" 00215 print jit_result 00216 00217 # Compute some data for each form from result of jit compilation 00218 data = {} 00219 #compiled_forms = jit_result # Ideally 00220 #compiled_forms, form_datas = jit_result # Previous 00221 compiled_forms, module, form_datas = jit_result # Current 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 # Benchmark the generated tabulate_tensor implementations 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 # Store results for future referencing 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 # Compare to references unless we're writing the references 00262 if not options.write: 00263 # Read reference files 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 # Compare fresh data and reference data 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 # Return final report 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 # Compute difference 00320 diff = data - reference 00321 00322 # Compute sums of squares 00323 d2 = numpy.sum(diff**2) 00324 r2 = numpy.sum(reference**2) 00325 n2 = numpy.sum(data**2) 00326 00327 # Compute normalized norm 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 # Norm from FFC, don't understand the motivation? 00337 #norm = math.sqrt(numpy.sum(diff / (1 + reference))) 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 #assert len(sh) == 2 # code elsewhere ensures data is represented as a matrix 00371 00372 if len(sh) == 1: #sh[0] == 1 or sh[1] == 1: 00373 # Got a vector, compare sorted vectors 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 # Got a matrix, compare matrix 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: # if it has passed the norm test, don't mark it as failed 00393 total_ok = False 00394 msg += "\n WARNING: Not computing eigenvalues of data with shape %s" % repr(sh) 00395 00396 # Make and return summary 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 # For each form 00406 for (b1, b2) in zip(data, reference): 00407 # For each integral type 00408 for (x, y) in zip(b1, b2): 00409 # For each integral 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 #f.write(data) 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 #data = f.read() 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 # Random cells were generated by: 00442 # >>> random.uniform(-2.5,2.5) 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 # FIXME: Define quadrilateral and hexahedral cells 00465 return mesh, cell 00466 00467 def compute_data(compiled_form, form_data, random_cell): 00468 # --- Initialize some variables 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 # --- Initialize geometry variables 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 # --- Initialize dofmaps 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 # --- Generate arbitrary coefficient dofsdofs 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 # --- Target variables 00501 A = numpy.zeros((1,1)) 00502 00503 # --- Add contributions from ALL domains from cell integrals 00504 if num_cell_integrals: 00505 domain = 0 00506 # Get shape of A and reset values 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 # --- Add contributions from ALL domains and facets from exterior integrals 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 # --- Add contributions from ALL domains and facets from interior integrals 00533 # FIXME: this currently makes no sense (integrating interior facets on 1 cell) 00534 # but it should be OK since we just compare numbers. 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 # Add A to the upper left quadrant of macro_A, it makes no sense, 00551 # but the numbers should be OK 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 # --- Execute! --- 00562 00563 if __name__ == "__main__": 00564 result = main(sys.argv[1:]) 00565 sys.exit(result)