[Commits] [svn:einsteintoolkit] pyGWAnalysis/trunk/DiscreteFunction/ (Rev. 7)

[reisswig at tapir.caltech.edu]

User: reisswig
Date: 2012/02/01 07:08 PM

Modified:
 /trunk/DiscreteFunction/
  DiscreteFunction.py

Log:
 Added function to align two DiscreteFunction objects
 by minimizing the L2 norm of the difference over a spcified range.

File Changes:

Directory: /trunk/DiscreteFunction/
===================================

File [modified]: DiscreteFunction.py
Delta lines: +74 -0
===================================================================
--- trunk/DiscreteFunction/DiscreteFunction.py	2012-02-02 01:00:03 UTC (rev 6)
+++ trunk/DiscreteFunction/DiscreteFunction.py	2012-02-02 01:08:28 UTC (rev 7)
@@ -692,3 +692,77 @@
     return DF2.__class__(x, f)
 
 
+
+def AlignOverRange(DF1, DF2, x0, x1, start_xoff, delta, dx, accuracy=1e-2, monotonic=False):
+    """Aligns two functions by minimizing their difference in the L2 norm and 
+       returns the x-shift that needs to be applied.
+       [x0,x1]: The minimization interval.
+       start_xoff: An initial guess to speed up computation.
+       delta: the shift will be applied in the range [start_xoff-delta, start_xoff+delta]. 
+              The smaller delta, the faster (but keep in mind that you need to be sure that the minimum is located in that range!
+       dx: The delta spacing of interpolated function values on which the pointwise difference is taken.
+       accuracy: The accuracy in finding the shift that minimizes the L2 norm of the error ('the delta spacing of the shift').
+       monotonic: Set to 'True' if function is monotonic (this will greatly speed up computation)"""
+
+    # Keep DF1 fixed. Find offset x_off that needs to be applied to DF2.x
+    # in order to make the L2-norm of the difference over a range [x0,x1]
+    # as small as possible.
+
+    x_off = start_xoff-delta
+    minErr = 10e10
+    minX = x_off
+    trueMinX = x_off
+
+    # If function is monotonic, we can use a faster algorithm
+    if (monotonic==True):
+	factor = 2.0*delta/accuracy*0.1
+	maxX = start_xoff+delta
+	while (factor >= 1.0):
+	    x_off = minX
+	    minErr = 10e10
+	    while (x_off < maxX):
+		DF2shift = DiscreteFunction(DF2.x+x_off, DF2.f)
+		subDF1, subDF2 = MakeConsistent(DF1, DF2shift, x0, x1, dx)
+		err = 0
+		for ii in range (0, subDF1.Length()):
+		    err += (subDF1.f[ii]-subDF2.f[ii])**2
+		err = sqrt(err)
+		#print minErr, err
+		if (err < minErr):
+		    minErr = err
+		    minX = x_off-accuracy*factor
+		    trueMinX = x_off
+		else:
+		    maxX = x_off
+		    break
+		x_off += accuracy*factor
+	    factor *= 0.1
+	return trueMinX
+	
+	print "ERROR! ...using brute force method!"
+    
+    # Proceed with brute force algorithm...
+
+    x_off = start_xoff-delta
+
+    minErr = 10e10
+    minX = x_off
+
+
+    # We shift DF2 and compute the error over range [x0,x1]
+    while (x_off < start_xoff+delta):
+	DF2shift = DiscreteFunction(DF2.x+x_off, DF2.f)
+	subDF1, subDF2 = MakeConsistent(DF1, DF2shift, x0, x1, dx)
+	err = 0
+	for ii in range (0, subDF1.Length()):
+	    err += (subDF1.f[ii]-subDF2.f[ii])**2
+	err = sqrt(err)
+	#print err
+	if (err < minErr):
+	    minErr = err
+	    minX = x_off
+	x_off += accuracy
+
+    return minX
+
+