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 + +