#pragma TextEncoding = "UTF-8"
#pragma rtGlobals=3				// Use modern global access method and strict wave access
#pragma DefaultTab={3,20,4}		// Set default tab width in Igor Pro 9 and later


function MTspeedtest(string func, variable order, variable datalength)
	
	make/O/N=(50,50,datalength) w3D, w3Dfits
	
	// create some data
	int i, j
	make/O/N=3 coef
	for (i=0;i<50;i++)
		for(j=0;j<50;j++)
			coef = enoise(10)
			w3D[i][j][] = poly(coef, r) + enoise(5000)
		endfor
	endfor
	
	Make/free/N=(DimSize(w3D, 0),DimSize(w3D,1))/wave wwFits
	wave/Z mask = $""
	
	variable timeST, timeMT
	tic()
	multithread wwFits = FF3DWorker(func, w3D, p, q, mask, polyorder=order)
	timeMT = toc()
	
	tic()
	wwFits = FF3DWorker(func, w3D, p, q, mask, polyorder=order)
	timeST = toc()
	
	printf "Unmasked %s fit. ", func
	printf "poly order: %g, MT time: %g, ST time: %g, Speedup: %g\r", order, timeMT, timeST, timeST/timeMT
		
	make/free/N=(datalength) mask
	mask[0,datalength/2] = 1
	
	tic()
	multithread wwFits = FF3DWorker(func, w3D, p, q, mask, polyorder=order)
	timeMT = toc()
	
	tic()
	wwFits = FF3DWorker(func, w3D, p, q, mask, polyorder=order)
	timeST = toc()
	
	
	printf "Masked %s fit. ", func
	printf "poly order: %g, MT time: %g, ST time: %g, Speedup: %g\r", order, timeMT, timeST, timeST/timeMT

	// populate the output wave
	for(i=DimSize(w3Dfits, 0)-1;i>=0;i--)
		for(j=DimSize(w3Dfits, 1)-1;j>=0;j--)
			wave/Z w = wwFits[i][j]
			if (WaveExists(w))
				w3Dfits[i][j][] = w[r]
			else
				w3Dfits[i][j][] = NaN
			endif
		endfor
	endfor
	
end

threadsafe function/wave FF3DWorker(string strFunction, wave w_3D, variable pp, variable qq, wave/Z w_mask, [variable polyorder])
	
	DFREF dfSav= GetDataFolderDFR()
	DFREF dfFree = NewFreeDataFolder()
	SetDataFolder dfFree
	
	Make/N=(DimSize(w_3D, 2)) w_data
	w_data = w_3D[pp][qq][p]
	
	Duplicate w_data dfFree:$NameOfWave(w_data) + num2str(pp) + num2str(qq) /wave=wFit
	
	try
		strswitch (strFunction)
			case "constant" :
				Make/O/D w_coef = {0,0}
				CurveFit/Q/N/H="01" line, kwCWave=w_coef, w_data /M=w_mask/NWOK; AbortOnRTE
				// the wave assignment
				FastOp wFit = (w_coef[0])
				w_coef = {w_coef[0]}
				break
			case "line" :
				CurveFit/Q/N line, w_data /M=w_mask/NWOK; AbortOnRTE
				wave/D w_coef
				wFit = w_coef[0] + w_coef[1] * x
				break
			case "poly" :
				CurveFit/Q/N poly polyorder, w_data /M=w_mask/NWOK; AbortOnRTE
				wave/D w_coef
				wFit = poly(w_coef, x)
				break
			case "gauss" :
				CurveFit/Q/N Gauss, w_data /M=w_mask/NWOK; AbortOnRTE
				wave/D w_coef
				wFit = w_coef[0] + w_coef[1] * exp(-((x - w_coef[2])/w_coef[3])^2)
				break
			case "lor" :
				CurveFit/Q/N lor, w_data /M=w_mask/NWOK; AbortOnRTE
				wave/D w_coef
				wFit = w_coef[0] + w_coef[1]/((x - w_coef[2])^2 + w_coef[3])
				break
			case "voigt" :
				CurveFit/Q/N voigt, w_data /M=w_mask/NWOK; AbortOnRTE
				wave/D w_coef
				wFit = VoigtPeak(w_coef, x)
				break
			case "sin" :
				CurveFit/Q/N sin, w_data /M=w_mask/NWOK; AbortOnRTE
				wave/D w_coef
				wFit = w_coef[0] + w_coef[1]*sin(w_coef[2]*x + w_coef[3])
				break
			case "sigmoid" :
				CurveFit/Q/N sigmoid, w_data /M=w_mask/NWOK; AbortOnRTE
				wave/D w_coef
				wFit = w_coef[0] + w_coef[1]/(1+exp(-(x - w_coef[2])/w_coef[3]))
				break
			case "exp" :
				CurveFit/Q/N exp, w_data /M=w_mask/NWOK; AbortOnRTE
				wave/D w_coef
				wFit = w_coef[0] + w_coef[1]*exp(-w_coef[2]*x)
				break
			case "dblexp" :
				CurveFit/Q/N dblexp, w_data /M=w_mask/NWOK; AbortOnRTE
				wave/D w_coef
				wFit = w_coef[0]+w_coef[1]*exp(-w_coef[2]*x)+w_coef[3]*exp(-w_coef[4]*x)
				break
			case "dblexp_peak" :
				CurveFit/Q/N dblexp_peak, w_data /M=w_mask/NWOK; AbortOnRTE
				wave/D w_coef
				wFit = W_coef[0]+W_coef[1]*(-exp(-(x-W_coef[4])/W_coef[2])+exp(-(x-W_coef[4])/W_coef[3]))
				break
			case "hillequation" :
				CurveFit/Q/N hillequation, w_data /M=w_mask/NWOK; AbortOnRTE
				wave/D w_coef
				wFit = w_coef[0]+(w_coef[1]-w_coef[0])*(x^w_coef[2]/(1+(x^w_coef[2]+w_coef[3]^w_coef[2])))
				break
			case "power" :
				CurveFit/Q/N power, w_data /M=w_mask/NWOK; AbortOnRTE
				wave/D w_coef
				wFit = w_coef[0] + w_coef[1]*x^w_coef[2]
				break
			case "log" :
				CurveFit/Q/N log, w_data /M=w_mask/NWOK; AbortOnRTE
				wave/D w_coef
				wFit = w_coef[0] + w_coef[1]*log(x)
				break
			case "lognormal" :
				CurveFit/Q/N lognormal, w_data /M=w_mask/NWOK; AbortOnRTE
				wave/D w_coef
				wFit = w_coef[0] + w_coef[1]*exp(-(ln(x/w_coef[2])/w_coef[3])^2)
				break
			default:
				SetDataFolder dfSav
				return $""
		endswitch

	catch
		if (V_AbortCode == -4)
			// Clear the error silently.
			variable CFerror = GetRTError(1)	// 1 to clear the error
			FastOp wFit = (NaN)
		endif
	endtry
		
	SetDataFolder dfSav
	return wFit // wFit becomes free wave when free data folder goes out of scope
end

function tic()
	variable/G tictoc = StartMSTimer
end

function toc()
	NVAR/Z tictoc
	variable ttTime = StopMSTimer(tictoc)
	KillVariables/Z tictoc
	return ttTime
end