Average data according to a category (or key or index) wave

// get the desired summary statistics

#pragma rtGlobals=3		// Use modern global access method and strict wave access.
function /WAVE averageByKey(waveToAvg, dataKeys[, keyTol])
// Given a wave waveToAvg and a wave dataKeys,
// determines the number of unique entries (categories/keys/...) in dataKeys
// and averages waveToAvg according to those entries.
// The standard error of the mean, # of points averaged, and corresponding
// categories are also saved.
// e.g. 
//	averageByKey(calibrationData, calibrationSolution)	// Key by entry in calibrationSolution
//	make /n=(numpnts(signal)) signalInt=round(signal); averageByKey(signal,signalInt) // Average to nearest integer
//
// OUTPUT:
	// averageByKey(data, keys), produces:
	// wave data_Keys -- a list of the unique data keys found in keyWave
	// wave data_Kavg -- the mean of the data corresponding to the keys
	// wave data_KSE  -- the standard error of the mean
	// wave data_Knum -- the # of points used for calculating the Avg & SD
	
wave waveToAvg
wave dataKeys
variable keyTol // Optional: the tolerance of keys, e.g. is 1.001 a "1" or not?
if (ParamIsDefault(keyTol))
	keyTol= 1e-07
endif
	
	// get list of unique keys
	wave keys = getKeys(dataKeys, tol=keyTol, wname=nameOfWave(waveToAvg) + "_Keys")
	variable nkeys = numpnts(keys)
	
	// try to predict memory overflow
	variable np=numpnts(waveToAvg) 
	if (np*nkeys > 2147000000)
		abort "averageByKey() was given too many keys / data points. Make a smaller selection."
	endif
	
	// prepare waves for saving summary stats
	make /d/n=(nkeys)/o $(nameOfWave(waveToAvg) + "_KAvg")		/WAVE=out_avg
	make /d/n=(nkeys)/o $(nameOfWave(waveToAvg) + "_KSE")		/WAVE=out_SE
	make /d/n=(nkeys)/o $(nameOfWave(waveToAvg) + "_Knum")		/WAVE=out_n
	out_avg=nan
	out_SE=nan
	out_n=nan
	
	// prepare output matrix
	make /FREE /d/n=(np,nkeys) outmat=nan // each column represents one key & contains only that key's data
	variable i,colkey
	for (i=0; i<nkeys; i+=1)
		SetDimLabel 1, i, $num2str(keys[i]), outmat  // label for later sorting
	endfor
	
	// go through the data, insert appropriate results into outmat
	string Sthiskey
	for (i=0; i<np; i+=1)
		Sthiskey = num2str(dataKeys[i])
		outmat[i][% $Sthiskey]=waveToAvg[i]
	endfor
	
	// get the desired summary statistics
	for (i=0; i<nkeys; i+=1)
		ImageStats /G={0, np-1, i,i} outmat
		if (V_flag < 0) // probably unnecessary: avoid saving the previous loop's stats
			abort "averageByKey(): unknown error" 
		endif
		out_n[i]=V_npnts
		if (V_npnts>0) // for V_npnts==1, V_sdev is set to NaN
			out_avg[i]=V_avg
			out_SE[i]=V_sdev / sqrt(V_npnts)
		endif
	endfor
	return keys
end
 
function /WAVE getKeys(keys[, tol, wname])
// Given a list of keys (e.g. laser on/off = 0/1),
// return a wave of the unique keys only (in e.g. above, {0,1})
 
wave keys
variable tol // optional, allow 1.0001==1 to be true.
string wname // optional name of output wave
if (ParamIsDefault(tol))
	tol= 1e-07
endif
if (ParamIsDefault(wname))
	wname= nameOfWave(keys) + "_uniqKeys"
endif
 
	// create output wave
	make /d/o/n=0 $wname /WAVE=outw
	outw=nan
 
	// special case of an empty key wave, or one with only NaN/inf
	wavestats /Q/M=1 keys
	if (V_npnts==0)
		return outw
	endif
 
	// go through input wave & save any new entries
	variable i=0, np=numpnts(keys)
	do
		if (getIndex(keys[i], outw, tol=tol) == -1)
			appendWave(outw, {keys[i]})
			SetDimLabel 0, (numpnts(outw)-1), $num2str(keys[i]), outw // label for ease of use (but tol!=0 must be considered)
		endif
		i += 1
	while (i<np)
 
	return outw
end
 
function getIndex(v, w, [tol, rev, start])
// Finds the index of the given variable var in the given wave w.
// Optional variable tolerance gives the ± acceptance range. Default 1.0E-7.
// Optional variable rev searches in reverse if set.
// Capable of finding NaNs.
// Returns -1 if not found.
// see also FindValue
 
	variable v
	wave w
	variable tol
	variable rev
	variable start // default == 0
 
 
	if (paramIsDefault(tol)) // TODO: change this to an autodetected spacing from w
		tol = 1.0E-7
	elseif (tol < 0)
		tol = abs(tol)
	endif
 
	variable i, n = numpnts(w)
 
	if ((rev))
		if (ParamIsDefault(start))
			start = n-1
		endif
		if (start < 0)
			abort "Cannot start at a negative index." 
			// Strictly speaking getIndex() is not responsible for this, but it helps with debugging.
		endif
 
		for (i =start; i >= 0; i -= 1)
			if (abs(w[i] - v) < tol)
				return i
			endif
		endfor
 
	else
 
		for (i =start; i < n; i += 1)
			if (abs(w[i] - v) < tol)
				return i
			endif
		endfor
 
	endif
 
	return -1
end
 
function appendWave(taker, appendThis)
// Concatenate alias. Given a wave "appendThis", appends its data to "taker". 1D waves only.
	wave taker, appendThis
	concatenate /NP=0 {appendThis}, taker
end