SyFi  0.3
verify_tensors/test.py
Go to the documentation of this file.
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)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Defines