#pragma rtGlobals=1		// Use modern global access method.

//---------------------------------------------------------------------------------
//Killgraphs
Function KillGraphs(numstart,numend)
variable numstart,numend
variable i
String str1

i = numstart
Do
	str1 = "Graph"+num2str(i)
	DoWindow/K $str1
	i = i+1
while (i < (numend+1))
end

Function KillTables(numstart,numend)
variable numstart,numend
variable i
String str1

i = numstart
Do
	str1 = "Table"+num2str(i)
	DoWindow/K $str1
	i = i+1
while (i < (numend+1))
End

//Medianfilter the original data
Function MedianSmooth(inWave,filterwidth)
Wave inWave
variable filterwidth
	Duplicate/O inWave, tempwave
	Wildpoint/F tempwave,filterwidth,0
	Duplicate/O tempwave, $NameofWave(inwave)+"_S"
	killwaves/Z tempwave
End

//Median filter and histogram a wave
Function SmoothAndBinWave(inwave,d,binsize)
Wave inwave
Variable d,binsize        
String Histname
variable start,nbins
String str1
Wave temp, temp_hist
Duplicate/O inwave, temp
MedianSmooth(temp,d)								//Median filter of wave
Make/O/n=0 temp_hist
WaveStats/Q temp_S									//Calculate histogram limits
start=floor(V_min)
nbins=(ceil(V_max)-start)/binsize
Histogram/b={start,binsize,nbins} temp_S temp_hist		//Bin the data		
//Display temp_S
//Display inwave
End

//histogram fitting to double gaussian to estimate the step location and stepsize
Function DoubleGaus(w,x) : FitFunc	
	Wave w
	Variable x
	 return w[0]*exp(-((x-w[1])^2/(2*(w[2]^2))))+w[3]*exp(-((x-w[4])^2/(2*(w[5]^2))))
End

Function Heaviside(w,x) : FitFunc  //heaviside function fitting
	wave w 
	variable x
	return w[3] + w[2]*(1/(1+exp(-2*w[1]*(x-w[0]))))
   	//w[0]=steplocn, w[1] controls fitting slope, w[2]=stepsize, w[3]=offset    
End

//-----------------------------------------------------------------------------------------------------------------------------------
//Program finds step locations for 2-state folding by fitting to Heaviside function, It also works for more steps
//winsize is size of window for calculating t-valules of data from constant extension value.
//d is Median smoothing parameter, binsize is histogram parameter. 
Function TwoStateStepFit(inwave,winsize,d,binsize,smthhist,levelthreshold)
Wave inwave
Variable winsize,d,binsize,smthhist,levelthreshold
Variable fitslope		//For setting the slope of Heaviside fit. 0.5 is reasonable number
Variable dt			//time scaling in data wave
Wave temp,temp_hist,temp_S,temp1,temp1_hist,smooth_hist,temp1_S
Variable i,j,k,m, n
Variable wavelength,dx
Variable level1,min1,level2,peak1,peak2

//Filter and bin the wave 
Duplicate/O inwave,temp1
//Display temp1
dt = deltax(inwave)
SmoothAndBinWave(temp1,d,binsize)
Duplicate/O temp_S temp1_S
Duplicate/O temp_hist temp1_hist
wavelength = numpnts(temp1)
//Calculate the extensions of the two states by fitting the extension histogram
Duplicate/O  temp1_hist smooth_hist
Smooth smthhist, smooth_hist
FindPeak/Q/M=25 smooth_hist					//Find first maximum in histogram
level1 = V_PeakLoc
peak1 = V_PeakVal
FindPeak/Q/N/R=(level1,inf) smooth_hist			//Find minimum
min1 = V_PeakLoc
FindPeak/Q/M=25/R=(min1,inf) smooth_hist		//Find second maximum
level2 = V_PeakLoc
peak2 = V_PeakVal
//Print level1,peak1,level2,peak2
Make/D/O/N=6 W_coef
W_coef = {peak1*1.3,level1,1,peak2*1.3,level2,1}
FuncFit/Q DoubleGaus W_coef temp1_hist /D
level1=W_coef[1]
level2=W_coef[4]
dx = level2-level1	
//Display smooth_hist
//print dx							//Distance between states
//Finding steps
Wave win,temp2
Wave FStepLocn,UStepLocn,FTime,UTime,FStepSize, UStepsize		
Variable startavg,stepFlag,steplocn,stepstartlevel,varthresh
Variable winext,fitthresh,stepsize

Make /O/N=(numpnts(inwave)) tValues=nan
Make /O/N=(winsize) wBefore,wAfter
Make/O/N=(winsize) win,temp2
temp2=temp1_S         //temp2 is the sliding average window wave, but it doesn't move forward in dowhile loop, don't know why
WaveStats/Q temp2    // how to change the data contents of temp2?
startavg = V_Avg	
//Find out if starts in state 1 or 2, set constant value for win to that state
If (abs(startavg-level1)<abs(startavg-level2))
	stepstartlevel = level1
    Else    
	stepstartlevel = level2
EndIf
//This part ends up truncating the wave, we lose the last bit. Find a better way in a subsequent version!
//Make/O/N=(numpnts(temp1)-windowSize) var 
Make/O/N=0 FStepLocn,FStepSize,Ftime,UStepLocn,UStepSize,Utime
//i=0				//counter for point positions in the wave
//j=0			//counter for number of steps
variable averagediff, sqrtvalues
Variable Fcntr,Ucntr	//counters for setting up waves to report results
Variable pvariance       //calculate pooled variance
Fcntr=0        //number of folding events
Ucntr=0		//number of unfolding events
	
//Finding putative steps
//Function StepFinder(inWave,threshold,windowSize)
Make /O/N=(numpnts(inwave)) tValues=nan
Make /O/N=(winsize) wBefore
Make/O/N=(winsize) wAfter
 i=winsize 	// 'i' will be the location of the putative step.
Do 						// Slide window until a step is found in the source wave 'w'.
	wBefore = temp1_S[p+i-winsize]
	wAfter = temp1_S[p+i]
	averagediff=mean(wBefore)-mean(wAfter)
	sqrtvalues=sqrt( variance(wBefore)/winsize+ variance(wAfter)/winsize)
	tValues[i] = averagediff/sqrtvalues
		i+=1
while (i<(numpnts(inwave)-winsize))	 //Stop after the population changed.
Display tValues
Make /O/N=0 Usteps,Fsteps
Findlevels /q/d=Usteps/edge=1 tValues,-levelthreshold 		// Find crossings down through t=-threshold (Folding steps).
Findlevels /q/d=Fsteps/edge=2 tValues,levelthreshold 			// Find crossings up through t=threshold (Unfolding steps).
//Edit Usteps,Fsteps
	 	
//Fit to heaviside function around step 
    	Make/O/N=(3*(winsize)) temp3,temp4			//Unscaled data!!!
Do
  	i=Usteps[Ucntr]
  	temp3 = temp1[p + i - 1.5*winsize]		//May want to change the start position of this wave for fitting!
//	Display temp3
	Wavestats/Q temp3
//	winext = V_max-V_min				//Find the distance from highest to lowest extension values in the fitting window
//	fitthresh = V_min + (winext/2)
	fitthresh =(V_min + V_max)/2
	FindLevel/Q/P temp3, fitthresh
 	Print "level for fitting: ",V_levelx
	fitslope = 0.5
	Make/D/N=4/O W_coef
	W_coef = {V_levelx, fitslope,dx,level1}		//changing win[0] to stepstartlevel, then changing stepstartlevel=level1 for unfolding steps
  	Make/O/T/N=2 T_Constraints; 				//Constrain w[1] to be positive. then sign of w[2] tells direction of fold/unfold
  	T_Constraints[0] = {"K0>0","K1 > 0"}; 			
  	FuncFit/Q/NTHR=0/X Heaviside W_coef temp3 /D/C=T_Constraints
	Print "Fit successful, W_coef: ", W_coef
//	If (W_coef[2] > 0)							//extension increases: unfold
		InsertPoints inf,1,UStepLocn,UStepSize,Ftime
		UStepLocn[Ucntr] = i-1.5*winsize+W_coef[0]		//point location of step in original wave
		UstepSize[Ucntr] = W_coef[2]
		If (Fsteps[0]<Usteps[0])
		Ftime[Ucntr]= (Usteps[Ucntr] - FSteps[Ucntr])*dt
		Print Usteps[Ucntr],Fsteps[Fcntr]
		ElseIf (Ucntr < 0.5)
			Ftime[Ucntr] = Usteps[Ucntr]*dt
			Else
			Ftime[Ucntr] = (Usteps[Ucntr] - FSteps[Ucntr-1])*dt
				Else
				Print "Problem with unfolding time at step", Fcntr
		EndIf
			Print "Location: ", UStepLocn[Ucntr], "Size: ",UstepSize[Ucntr], "Lifetime: ", Ftime[Ucntr],"U counter:", Ucntr,"F Counter:", Fcntr
//	Endif
	Ucntr+=1
While (Ucntr<numpnts(Usteps))

Do
  	i=Fsteps[Fcntr]
  	temp4 = temp1[p + i - 1.5*winsize]
//	Display temp4
	Wavestats/Q temp4
//	winext = V_max-V_min				//Find the distance from highest to lowest extension values in the fitting window
//	fitthresh = V_min + (winext/2)
	fitthresh =(V_min + V_max)/2
	FindLevel/Q/P temp4, fitthresh
 	Print "level for fitting: ",V_levelx
	fitslope = 0.5
	Make/D/N=4/O W_coef
	W_coef = {V_levelx, fitslope,dx,level2}		//changing win[0] to stepstartlevel, then changing stepstartlevel to level2 for folding steps
  	Make/O/T/N=1 T_Constraints; 				//Constrain w[1] to be positive. then sign of w[2] tells direction of fold/unfold
  	T_Constraints[0] = {"K1 > 0"}; 			
  	FuncFit/Q/NTHR=0/X Heaviside W_coef temp3 /D/C=T_Constraints
	Print "Fit successful, W_coef: ", W_coef
//	If (W_coef[2] < 0)								//extension decreases: fold
		InsertPoints inf,1,FStepLocn,FStepSize,Utime
		FStepLocn[Fcntr] = i-1.5*winsize+W_coef[0]	//point location of step in original wave
		FstepSize[Fcntr] = W_coef[2]
		If (Fsteps[0]>Usteps[0])
		Utime[Fcntr]= (FstepLocn[Fcntr] - UStepLocn[Fcntr])*dt
		ElseIf (Fcntr < 0.5)
			Utime[Fcntr] = FStepLocn[Fcntr]*dt
			ElseIf (Fcntr > 0.5)
			Utime[Fcntr] = (FstepLocn[Fcntr] - UStepLocn[Fcntr-1])*dt
				Else
				Print "Problem with unfolding time at step", Fcntr
		EndIf
		Print "Location: ", FStepLocn[Fcntr], "Size: ",FstepSize[Fcntr], "Lifetime: ", Utime[Fcntr],"F Counter:", Fcntr,"U counter:", Ucntr
//	Endif
		Fcntr+=1
While (Fcntr<numpnts(Fsteps))
Edit FStepLocn,FstepSize,Utime,UStepLocn,UstepSize,Ftime
End