Minimization problem: Match parts of two waves via scale and offset

Share icon

I search a solution for the, at first glance probably simple, problem of generating a close-to-zero difference for only parts of two waves, i.e., automatically bring them to match by optimization of some modification parameters. I have attached some demo data for everyone to play with (you can see in the history area how the data was generated). There are two waves, one with a distinct Gaussian feature near the center and one with plain ‘background’ data. Now, the outer ~30% of the two waves can be brought to match by simply adjusting the offset (75 in this demo) and the scale (0.72 here).

I want to find the offset and scale values for the first wave to match programmatically, i.e., a function which calculates these values from the two waves (and optionally how much of the outer values can be used in the calculation).

My ideas are to use a temporary wave to calculate the difference and use some complicated iteration until the difference is very small or abuse functions like 'Optimize' in some way (maybe I’m just too tired today, but I can't find a way ;). Since noise is included, there has to be some criteria for this. Also there may be cases where there is no exact match possible, which also needs to be considered. In Excel there is the solver add-in, where the program fiddles with some values until a minimum/maximum etc. is achieved for some result value. Basically I need this for the sum of the difference for these two waves (in parts). I hope someone can point me to an elegant solution for this. Minimization Test.pxp (29.42 KB) 
Optimize
looks good to me.

A first rough sketch with

#pragma rtGlobals=3		// Use modern global access method and strict wave access.
#include <KBColorizeTraces>

Function func(pw, slope, offset)
	Wave pw
	variable slope, offset
	
	Wave wave1, wave2

	variable pSkipStart = 200
	variable pSkipEnd   = 320
	
	Duplicate/O/FREE/D wave1, sq
	
	sq[] = (p <= pSkipStart || p >= pSkipEnd ? wave1[p] * slope + offset - wave2[p] : 0.0)
	sq = sq[p]^2
	
	variable ret = Sum(sq)
	print ret
	return ret
End

Function doStuff()

	Make/O/N=0/FREE pw
	Optimize/A=0/X={1,0}/T=1e-2 func, pw
	
	Wave wave1, wave2, W_Extremum
	Duplicate/O wave1, result

	variable slope  = W_Extremum[0]
	variable offset = W_Extremum[1]
	result  = wave1[p] * slope + offset
	Display wave1, wave2, result
End


results in following graph.

It is not completely automatized yet as you still need to find pSkipStart and pSkipEnd manually.
Graph4_1.png (110.48 KB) Graph4_1.png (110.48 KB) Graph4_1.png (110.48 KB)

September 3, 2014 at 10:03 am - Permalink
Thanks a lot for the solution. I haven't used 
Optimize
before. It is no problem that you have to give the omitted range manually. I guess there is no way to find this automatically for a wide range of data.
I have modified the code a bit, and giving Optimize a rough parameter size using the /R flag gives a great result. Now I only need to find an elegant way to sneak two arbitrary named waves into the worker function.


Function OptWorker(pw, scale, offset)
	Wave pw
	variable scale, offset
 
	Wave wave1, wave2
	Duplicate/O/D/FREE wave1, sq

 	sq		= wave1 * scale + offset - wave2
 	sq[pw[0],pw[1]]	= 0
	sq		= sq^2
	return  Sum(sq)
End
 
Function doStuff()
	Make/O/N=0/FREE pw={170,350}
	Optimize/Q/A=0/X={1,0}/R = {1, 100} OptWorker, pw
 
	Wave wave1, wave2, W_Extremum
	Duplicate/O wave1, result

	print "Scaling:",W_Extremum[0], "Offset:", W_Extremum[1]
	result = wave1 * W_Extremum[0] + W_Extremum[1]
	KillWaves/Z W_Extremum, W_OptGradient
End
Graph.png (19.92 KB)

September 3, 2014 at 07:42 pm - Permalink
Well for automatically finding the omitted range, you could try 
FindAPeak
and assume some width.

For the sneak in part, unfortunately 
Optimize
does not accept wave reference waves as parameter wave.
So we have to resort to wave notes:

Function OptWorker(pw, scale, offset)
	Wave pw
	variable scale, offset
 
 	string wvNote = note(pw)
 	Wave wave1 = $StringFromList(0,wvNote)
 	Wave wave2 = $StringFromList(1,wvNote)
 	
	Duplicate/D/FREE wave1, sq
 
 	sq		= wave1 * scale + offset - wave2
 	sq[pw[0],pw[1]]	= 0
	sq		= sq^2
	return  Sum(sq)
End
 
Function doStuff()
	Make/FREE pw={170,350}
	Note pw, "root:wave1;root:wave2"
	Optimize/Q/A=0/X={1,0}/R = {1, 100} OptWorker, pw
 
	Wave wave1, wave2, W_Extremum
	Duplicate/O wave1, result
 
	print "Scaling:",W_Extremum[0], "Offset:", W_Extremum[1]
	result = wave1 * W_Extremum[0] + W_Extremum[1]
	KillWaves/Z W_Extremum, W_OptGradient
End

September 4, 2014 at 02:49 am - Permalink
I have a different solution, which is a bit left-field.
I have fitted one wave to the other using an offset and scale, then excluded the points that are most extreme (positive peak). Then repeated the fit and exclusion a few times.
Example code as follows - I hope it is self explanatory:

Function FitDriver()
	wave wave1 // assume wave1 and wave2 already exist
	wave wave2
	variable i,vNum,vNumIter,vIter,vFrac
	vFrac=0.1 // fraction of points to 'mask' out at each iteration
	vNumIter=3 // number of iterations to run
	
	vNum=vFrac*DimSize(wave1,0) // number of pint to mask out each iteration
	
	duplicate/O wave1,waveFit,waveRes
	duplicate/O wave1,wave1iter
	
	Make/D/O/N=(2) W_coef  // make a coeficient wave
	W_coef[0]=0 // offset
	W_coef[1]=1 // scale
	
	for(vIter=0;vIter<vNumIter;vIter+=1) // loop round each iteration
		// run the fit
		FuncFit/Q myFitFunc, W_coef, wave2 /X=wave1iter
		
		waveFit[]=W_coef[0]+W_coef[1]*wave1iter[p]
		waveRes[]=wave2[p]-waveFit[p] // calculate a residuals wave
		
		for(i=0;i<vNum;i+=1) // loop to exclude extreme points from next fit
			waveStats/Q/Z waveRes
			wave1iter[V_maxloc]=NaN
			waveRes[V_maxloc]=NaN
		endfor
	endfor
	// run the fit one last time
	FuncFit/Q myFitFunc,W_coef,wave2 /X=wave1iter
	// calc fit
	waveFit[]=W_coef[0]+W_coef[1]*wave1[p]
	// plot graph
	Display wave1,wave2,waveFit
	ModifyGraph rgb(wave2)=(0,52224,0),rgb(waveFit)=(0,15872,65280)
	TextBox/C/N=text0 "Fit Parameters:\rOffset = "+num2str(W_coef[0])+"\rScale = "+num2str(W_coef[1])
End
Function myFitFunc(pw, yw, xw) : FitFunc
	WAVE pw, yw, xw
	yw = pw[0]+pw[1]*xw
End

Hope this makes sense.
I tried to get this to work using a mask wave, but failed, hence my odd use of 'wave1iter'.
Regards,
Kurt

September 4, 2014 at 07:13 am - Permalink
Wow, that's also a nice solution. Thanks a lot. I never thought that you could take an arbitrary wave like that and just do a fit. Also, this solution doesn't need a 'hint' the Optimize approach (the result magnitude is needed). Maybe we should run a small contest one in while with problems like these to discover new ways of working. Also, "... number of pint ..." made my day. Cheers. XD

September 4, 2014 at 06:14 pm - Permalink
I have realised that the choice to only exclude positive points is not great when the data is noisy as it makes the resultant fit biased to more negative points (you can clearly see this if you run the fit after changing the variables vFrac=0.01 and vNumIter=30 at the start of my old code).
A better solution (at least for this shape of data) is to exclude the most extreme points irrespective of sign.
One could also extend this algorithm to auto stop excluding points when, for example, the largest outlier is 2 x StDev of the residual (say) - this removes another 'human' parameter from the algorithm. (Perhaps this factor should be scaled by number of data points - but my head starts to hurt when I start thinking about that).
New example code as follows:

Function FitDriver()
	wave wave1 // assume wave1 and wave2 already exist
	wave wave2
	variable i,vNum,vMaxNumIter,vIter,vFrac,vMaxSD
	vFrac=0.01 // fraction of points to 'mask' out at each iteration
	vMaxNumIter=50 // maximum number of iterations to run - a safety net
	vMaxSD=2 // Stop excluding points when largest outlier exceeds this number of StdDev's
 
	vNum=vFrac*DimSize(wave1,0) // number of points to mask out each iteration
 
	duplicate/O wave1,waveFit,waveRes
	duplicate/O wave1,wave1iter
 
	Make/D/O/N=(2) W_coef  // make a coefficient wave
	W_coef[0]=0 // offset
	W_coef[1]=1 // scale
 
	vIter=0
	do
		// run the fit
		FuncFit/Q myFitFunc, W_coef, wave2 /X=wave1iter
 
		myFitFunc(W_coef,waveFit,wave1iter)
		waveRes[]=wave2[p]-waveFit[p] // calculate a residuals wave
 
		for(i=0;i<vNum;i+=1) // loop to exclude extreme points from next fit
			waveStats/Q/Z waveRes
			if(abs(V_max)>abs(V_min))
				wave1iter[V_maxloc]=NaN
				waveRes[V_maxloc]=NaN
			else
				wave1iter[V_minloc]=NaN
				waveRes[V_minloc]=NaN
			endif
		endfor
		if(max(abs(V_max),abs(V_min))<vMaxSD*V_sdev) // extreme point not too large...
			break // ... so stop here
		endif
		
		vIter+=1
	while(vIter<vMaxNumIter)
	
	// run the fit one last time
	FuncFit/Q myFitFunc,W_coef,wave2 /X=wave1iter
	// calc fit
	myFitFunc(W_coef,waveFit,wave1)
	// plot graph
	Display wave1,wave2,waveFit
	ModifyGraph rgb(wave2)=(0,52224,0),rgb(waveFit)=(0,15872,65280)
	TextBox/C/N=text0 "Fit Parameters:\rOffset = "+num2str(W_coef[0])+"\rScale = "+num2str(W_coef[1])
End

Function myFitFunc(pw, yw, xw) : FitFunc
	WAVE pw, yw, xw
	yw = pw[0]+pw[1]*xw
End


Apologies for the pint-point typo... perhaps my fingers were hinting at where I should go after work ;)

Regards,
Kurt

September 5, 2014 at 12:33 am - Permalink
Thanks a lot Kurt for the update. Funnily, both versions of your approach give almost the exact same result in my real-world problem unlike the test case. Usually it's the other way around, where you think it works until you discover the solution utterly fails for the real thing. I wonder if there is a more straight forward solution for the excluding-points part. Maybe 'Extract' in some way?

September 8, 2014 at 02:04 am - Permalink