#pragma TextEncoding	= "UTF-8"
#pragma rtGlobals		= 3
#pragma IgorVersion		= 8.00
#pragma moduleName		= MPFDataFit
#pragma DefaultTab		= {3,20,4}

//________________________________________________________________________________________________
//										MFP Data Fit Function
//	Written by Stephan Thuemer - https://www.wavemetrics.com/user/chozo
//	version 2022-04
//
//	Provides the framework to fit experimental data with Multipeak Fit. After loading MPF, first
//	'Prepare Data Peak Function' needs to be run to set the folder and wave names for the data.
//	The DataPeak fit function is available from within MPF. Here, the last two parameters control
//	the folder and wave number and need to be on Hold at all times. Dial in the desired value
//	to select the data to fit.
//________________________________________________________________________________________________

static StrConstant kMPF_WorkingDir = "root:Packages:MultiPeakFit2:"

Menu "Analysis"
	Submenu MPFDataFit#CreateDataFitMenu()													// add menu entry after loading MPF
		"Prepare Data Peak Function", /Q, MPF_PrepareDataFit()
	End
End

static Function/S CreateDataFitMenu()
	return SelectString(DataFolderExists(kMPF_WorkingDir),"-","Multipeak Fit")
End

// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

Function MPF_PrepareDataFit()																	// setup for data locations
	String setBaseFolder = "", setBaseName = ""
	SVAR/Z baseFolder = $(kMPF_WorkingDir+"MPF_DataFit_baseFolder")								// load previous settings
	SVAR/Z baseName = $(kMPF_WorkingDir+"MPF_DataFit_baseName")
	if (SVAR_Exists(baseFolder))
		setBaseFolder = baseFolder
	endif
	if (SVAR_Exists(baseName))
		setBaseName = baseName
	endif
	
	Prompt setBaseFolder,"Set base folder name without trailing number (e.g., \"root:myFit_\"):"
	Prompt setBaseName,"Set base data name without trailing number (e.g., \"peak_\"):"
	DoPrompt "Setup MPF data fit", setBaseFolder, setBaseName
	if (V_Flag)
		return -1
	endif
	if (!strlen(setBaseFolder) || !strlen(setBaseName))
		Abort "You need to enter both a folder path and base data name."
	endif
	
	int i, j, sourceFound = 0
	if (!DataFolderExists(kMPF_WorkingDir))														// if MPF was not started yet => create basic folder structure
		DFREF saveDF = GetDataFolderDFR()
		SetDataFolder root:
		for (i = 1; i < ItemsInList(kMPF_WorkingDir,":"); i++)
			NewDataFolder/O/S $StringFromList(i,kMPF_WorkingDir,":")
		endfor
		Variable/G currentSetNumber = 0
		SetDataFolder saveDF
	endif
	setBaseFolder = RemoveEnding(setBaseFolder , ":")											// make sure no colon was entered
	String/G $(kMPF_WorkingDir+"MPF_DataFit_baseFolder") = setBaseFolder
	String/G $(kMPF_WorkingDir+"MPF_DataFit_baseName") = setBaseName
	
	for (i = 0; i < 10 && !sourceFound; i++)													// try to look for the data in the first 10 folder / waves
		for (j = 0; j < 10 && !sourceFound; j++)
			Wave/Z source = $(setBaseFolder+num2str(i)+":"+setBaseName+num2str(j))
			sourceFound = WaveExists(source)
		endfor
	endfor
	
	if (!sourceFound)																			// still nothing found? => have a simple check to remind the user
		DoAlert/T="Checking folder and data" 0, "Your first folder and/or data does not seem to exist yet. Make sure to prepare the data in the format: "+setBaseFolder+"[NUM]:"+setBaseName+"[NUM]"
	endif
	return 0
End

//################################################################################################

static Function/Wave getSource(Wave cw)															// get source data from saved base names
	SVAR/Z baseFolder = $(kMPF_WorkingDir+"MPF_DataFit_baseFolder")
	SVAR/Z baseName   = $(kMPF_WorkingDir+"MPF_DataFit_baseName")
	if (SVAR_Exists(baseFolder) && SVAR_Exists(baseName))
		Wave/Z source = $(baseFolder+num2str(cw[3])+":"+baseName+num2str(cw[4]))
		if (WaveExists(source) && WaveType(source,1) == 1 && WaveDims(source) == 1)
			return source
		endif
	endif
	return $""
End

//################################################################################################

static Function MPF_setHoldStateForPeak(Wave cw, String holdStr)								// injects hold state for the current peak
	DFREF DFpath = GetWavesDataFolderDFR(cw)
	int peakNum
	sscanf NameOfWave(cw), "Peak %d Coefs", peakNum
	if (!DataFolderRefStatus(DFpath) || numType(peakNum) != 0 || peakNum < 0)
		return -1
	endif
	SVAR gName = DFPath:GraphName
	String pName = gName+"#MultiPeak2Panel"														// the MPF panel
	if (WinType(pName) == 7)
		Wave/T HoldStrings = DFPath:HoldStrings
		Execute/Z "WMMultiPeakFit#MPF2_RefreshHoldStrings(\""+pName+"\")"						// first, grab the latest hold settings from the panel
		HoldStrings[peakNum+1] = holdStr														// inject hold string
		Execute/Z "WMMultiPeakFit#SetHoldCheckboxesFromWave(WMMultiPeakFit#GetSetNumberFromWinName(\""+pName+"\"))"	// rewrite to panel
	endif
	return 0
End

//################################ MPF peak format functions #####################################

Function/S DataPeak_PeakFuncInfo(Variable which)												// pass names to MPF
	Switch (which)
		case 0:
			return "LocShift;Broaden;Height;Folder;DataNo;"
		case 1:
			return "MPF_DataPeakFunc"
		case 2:
			return "MPF_GaussToDataPeakGuess"
		case 3:
			return "MPF_DataPeakParams"
		case 4:
			return "Location;Height;Area;FWHM;"
	EndSwitch
End

//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

Function MPF_GaussToDataPeakGuess(Wave cw)														// initializes parameters, and set hold state
	// cw values are: "PosShift;GaussWidth;Height;Folder;DataNo;"
	Redimension/N=5 cw
	int i,j, sourceFound = 0
	for (i = 0; i < 10 && !sourceFound; i++)
		for (j = 0; j < 10 && !sourceFound; j++)
			cw[3] = i																			// make sure these are set to the right folder and data before grabbing the source
			cw[4] = j
			Wave/Z source = getSource(cw)														// grabs the first valid source wave, assumes all other peaks behave the same
			sourceFound = WaveExists(source)
		endfor
	endfor
	
	int normInt = 1, normPos = 1
	if (sourceFound)
		Variable wMax  = WaveMax(source)
		Variable left  = leftx(source) < rightx(source) ? leftx(source) : rightx(source)
		Variable right = leftx(source) < rightx(source) ? rightx(source) : leftx(source)
		normInt = wMax > 0.9 && wMax < 1.1 ? 0 : 1
		normPos = left < 0 && right > 0 ? 0 : 1
	endif
	
	cw[0] = normPos ? 0 : cw[0]																	// if data seems centered around zero, peak position can stay
	cw[1] = 0.5																					// initial Gaussian broadening
	cw[2] = normInt ? 1 : cw[2]																	// if data is normalized, peak height can stay
	
	Execute/P/Q "MPFDataFit#MPF_setHoldStateForPeak("+GetWavesDataFolder(cw,2)+", \"00011\")"	// this needs to be executed later or the hold state will be reset
	return 0
End

//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

Function MPF_DataPeakFunc(Wave cw, Wave yw, Wave xw) : FitFunc									// the fit function
	Wave/Z source = getSource(cw)
	if (!WaveExists(source))
		yw = NaN
		if (Exists("MPF2_DisplayStatusMessage"))												// display an error in the MPF graph
			Execute/Z "MPF2_DisplayStatusMessage(\"Base folder setup wrong or data not found.\","+GetWavesDataFolder(cw,1)+")"
		endif
		return -1
	endif
	
	Duplicate/free source, yWave
	SetScale/P x, DimOffset(source,0)+cw[0], DimDelta(source,0), ywave							// prepare shifted convolution wave
	PeakConvolveGaussHelper(yWave, cw[1])														// Gauss convolution
	
	Variable minX = pnt2x(yWave,0), maxX = pnt2x(yWave,numpnts(yWave)-1)						// make sure to not run out-of-bounds (end-point data is repeated)
	yw = xw[p] >= minX ? ( xw[p] <= maxX ? cw[2]*yWave(xw[p]) : cw[2]*yWave(maxX) ) : cw[2]*yWave(minX)
	return 0
End

//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

Function MPF_DataPeakParams(Wave cw, Wave sw, Wave outWave)
	Variable Amplitude = NaN, Location = NaN, FWHM = NaN, PeakArea = NaN						// set to NaN for now
	FindPeakStatsHelper(cw, Amplitude, Location, FWHM, PeakArea)								// approximately calculates all the results from the parameters
	// cw values are: "PosShift;GaussWidth;Height;Folder;DataNo;"
	// output values are: "Location;Height;Area;FWHM;"
	outWave[0][0] = Location		// position
	outWave[0][1] = NaN
	outWave[1][0] = Amplitude		// height
	outWave[1][1] = NaN
	outWave[2][0] = PeakArea		// area
	outWave[2][1] = NaN
	outWave[3][0] = FWHM			// FWHM
	outWave[3][1] = NaN
	return 0
End

//################################################################################################

Static Function PeakConvolveGaussHelper(Wave in, Variable width)								// convolve with a prepared Gaussian
	if (width == 0 || numtype(width) != 0)
		return 0
	endif
	
	Variable pnt, dx  = DimDelta(in,0), size = DimSize(in,0)
	width = abs(width/dx)/sqrt(2)																// width normalized by x scaling and convert to expected Gauss 'width'
	pnt = round(max(abs(10*width),11))															// create a wave 10 times the Gauss width to get the full Gauss in (5 sigma on each side)
	pnt = min(pnt, 2*numpnts(in)+1)																// doesn't have to be larger than two times the full input data
	pnt = pnt+!mod(pnt, 2)																		// force odd size

	Make/FREE/D/N=(pnt) GaussWave	
	GaussWave = Gauss(x, (pnt-1)/2, width)
	Variable A = sum(GaussWave)																	// make sure the Gauss area limited even for a few points
	if (A > 1)
		GaussWave /= A
	endif
	
	Make/FREE/D/N=(2*pnt+size) temp = in[0]														// enlarge the wave to get the full convolution range
	temp[pnt,size+pnt-1] = in[p-pnt]
	temp[size+pnt,size+2*pnt-1] = in[size-1]
	MatrixOP/Free temp = replaceNaNs(temp,0)													// makes sure there are no NaN values in-between
	Convolve/A GaussWave temp
	in[0,size-1] = temp[p+pnt]
End

//################################################################################################

Static Function FindPeakStatsHelper(Wave cw, Variable &Amplitude, Variable &Location, Variable &FWHM, Variable &PeakArea)	// recalculate the (almost) full peak and extract all result values numerically
	Wave/Z sw = getSource(cw)
	if (!WaveExists(sw))
		return -1
	endif
	// cw values are: "PosShift;GaussWidth;Height;Folder;DataNo;"
	Variable left = pnt2x(sw,0)+cw[0], right = pnt2x(sw,numpnts(sw)-1)+cw[0]
	Variable pntSize = DimSize(sw,0)*10															// create a big wave (fine increments) => decides the numeric resolution
	pntSize = pntSize > 100000 ? 100000 : pntSize												// make sure its not too big
	Make/FREE/D/N=(pntSize+1) yw, xw
	Make/FREE/D/N=2 Levels
	SetScale/I x, left, right, xw; xw = x
	//+++++++++++++++++++++++++++++++++++++ find values ++++++++++++++++++++++++++++++++++++++++++
	MPF_DataPeakFunc(cw, yw, xw)																// recalculate the peak with the coef set
	PeakArea = AreaXY(xw, yw, left, right)														// the area is calculated for the full source range
	WaveStats/Q yw
	Amplitude = cw[2] > 0 ? V_max : V_min
	Location  = cw[2] > 0 ? xw[V_maxloc] : xw[V_minloc]
	FindLevels/Q/D=Levels yw, Amplitude/2														// find the half-maximum positions
	if (V_LevelsFound == 2)
		left  = xw[Levels[0]]
		right = xw[Levels[1]]	
		FWHM  = abs(left - right)
	endif
	if (numtype(FWHM) != 0)
		return 0
	endif
	//++++++++++++++++++++++++++ recalculate with higher resolution ++++++++++++++++++++++++++++++
	SetScale/I x, left - 2*FWHM/(pntsize+1), right + 2*FWHM/(pntsize+1), xw; xw = x				// choose a range just a bit bigger than the FWHM
	MPF_DataPeakFunc(cw, yw, xw)
	WaveStats/Q yw
	Amplitude = cw[2] > 0 ? V_max : V_min
	Location  = cw[2] > 0 ? xw[V_maxloc] : xw[V_minloc]
	FindLevels/Q/D=Levels yw, Amplitude/2
	if (V_LevelsFound == 2)
		left  = xw[Levels[0]]
		right = xw[Levels[1]]	
		FWHM  = abs(left - right)
	endif
	return 0
End