smbinterp/bin/driver.py

109 lines
3.2 KiB
Python
Raw Normal View History

#!/usr/bin/python
import sys
from optparse import OptionParser
import numpy as np
from baker.tools import rms, exact_func
import grid
def get_mesh(source, destination, use_structured_grid = False):
mesh_source = None
mesh_dest = None
if use_structured_grid:
mesh_source = grid.simple_rect_grid(source, source)
mesh_dest = grid.simple_rect_grid(destination, destination)
else:
mesh_source = grid.simple_random_grid(source)
mesh_dest = grid.simple_random_grid(destination)
if not (mesh_dest and mesh_source):
raise smberror('problem creating mesh objects')
else:
return mesh_source, mesh_dest
if __name__ == '__main__':
parser = OptionParser()
parser.add_option("-o",
"--output-file",
dest="output",
type='str',
default = '/tmp/for_qhull.txt',
help = "qhull output file")
parser.add_option("-e",
"--extra-points",
dest="extra",
type='int',
default = 3,
help = "how many extra points")
parser.add_option("-s",
"--source-total",
dest="source_total",
type='int',
default = 100,
help = "total number of source points for random,\
resolution for structured")
parser.add_option("-d",
"--destination-total",
dest="destination_total",
type='int',
default = 100,
help = "total number of destination points,\
resolution for structured")
parser.add_option("-r",
"--structured",
action = 'store_true',
default = False,
help = "use a structured grid instead of random point cloud")
parser.add_option("-v",
"--verbose",
action = 'store_true',
default = False,
help = "verbosity")
(options, args) = parser.parse_args()
(mesh_source, mesh_dest) = get_mesh(options.source_total, options.destination_total, options.structured)
open(options.output, 'w').write(mesh_source.for_qhull())
if options.verbose:
print >> sys.stderr, "options: %s, args: %s" % (options, args)
print >> sys.stderr, "wrote source mesh output to %s" % options.output
errors = []
success = 0
for X in mesh_dest.points:
answer = mesh_source.run_baker(X)
if answer['a'] == None:
errors.append(0)
continue
exact = exact_func(X[0], X[1])
if np.abs(exact - answer['final']) <= np.abs(exact - answer['qlin']):
success += 1
if options.verbose:
print "current point : %s" % X
print "exact : %0.4f" % exact
print "qlin : %0.4f" % answer['qlin']
print "q_final : %0.4f" % answer['final']
print "qlinerr : %1.4f" % (exact - answer['qlin'],)
print "q_final_err : %0.4f" % (exact - answer['final'],)
cur_error = np.abs(answer['final'] - exact)
errors.append(cur_error)
print rms(errors)
print "%s of %s won" % (success, len(mesh_dest.points))