FitMySpectrum
TRK
#pragma rtGlobals=3 // Use modern global access method and strict wave access.
#include <WMBatchCurveFitIM>
Function FitMySpectrum(ww, xx) : FitFunc
WAVE ww
Variable xx
Variable w1=0,w2=0,w3=0,w4=0,w5=0,w6=0,w7=0,w8=0,w9=0,w10=0,w11=0,w12=0,w13,w14=0,qq=0
Variable T1=0.24, S1=0.047, T2=0.24, S2=0.047, x1=1.91, x2=1.918, S3=0.26, x3=2.20
Variable k, i
Make/O/N=46 A1 //wave to store coefficients A1
Make/O/N=46 A2 //wave to store coefficients A2
Make/O/N=(41,46) fitdata //wavet to store results from fitting
//CurveFitDialog/ These comments were created by the Curve Fitting dialog. Altering them will
//CurveFitDialog/ make the function less convenient to work with in the Curve Fitting dialog.
//CurveFitDialog/ Equation:
//CurveFitDialog/ f(xx) = (A1*(1110*(exp(0.5*(S1/T1)^2-(xx-x1)/T1))*erfc(0.707*(S1/T1-(xx-x1)/S1)))+200*(exp(-((xx-x3)/S3)^2)))+(A2*(373*(exp(0.5*(S2/T2)^2-(xx-x2)/T2)*erfc(0.707*(S2/T2)-(xx-x2)/S2))+200*(exp(-((xx-x3)/S3)^2)))
//CurveFitDialog/ End of Equation
//CurveFitDialog/ Independent Variables 1
//CurveFitDialog/ xx
//CurveFitDialog/ S1 --> Set equal to 0.064
//CurveFitDialog/ T1 --> Set equal to 0.23
//CurveFitDialog/ S2 --> Set equal to 0.12
//CurveFitDialog/ T2 --> Set equal to 0.09
//CurveFitDialog/ x1 --> Set equal to 1.92
//CurveFitDialog/ x2 --> Set equal to 2.17
//CurveFitDialog/ Coefficients 2
//CurveFitDialog/ ww[0] = A1
//CurveFitDialog/ ww[1] = A2
w1 = exp(0.5*(S1/T1)^2-(xx-x1)/T1)
w2 = erfc(0.707*(S1/T1-(xx-x1)/S1))
w3 = 1110
w4 = exp(0.5*(S2/T2)^2-(xx-x2)/T2)
w5 = erfc(0.707*(S2/T2-(xx-x2)/S2))
w6 = 373
w7 = 200
w8 = exp(-((xx-x3)/S3)^2)
w9 = 200
w10 = 560
w11 = ww[0]
w12 = ww[1]
w13 = w3*(w1*w2)+w7*w8
w14 = w6*(w4*w5)+w10*w8
qq = w11*w13+w12*w14
return qq
for(i=0; i<k; i+=1)
String path = "root:ps_2:input" //directory of my 2D wave
Display $path
DFREF dfr = $path
k=DimSize ($path,1)
for(i=0; i<k; i+=1)
Display dfr:f35[][i]
FuncFit/NTHR=0/Q/N FitMySpectrum (f35[][i]) //f35 is the input, 2D wave to which fitting is done
endfor
#include <WMBatchCurveFitIM>
Function FitMySpectrum(ww, xx) : FitFunc
WAVE ww
Variable xx
Variable w1=0,w2=0,w3=0,w4=0,w5=0,w6=0,w7=0,w8=0,w9=0,w10=0,w11=0,w12=0,w13,w14=0,qq=0
Variable T1=0.24, S1=0.047, T2=0.24, S2=0.047, x1=1.91, x2=1.918, S3=0.26, x3=2.20
Variable k, i
Make/O/N=46 A1 //wave to store coefficients A1
Make/O/N=46 A2 //wave to store coefficients A2
Make/O/N=(41,46) fitdata //wavet to store results from fitting
//CurveFitDialog/ These comments were created by the Curve Fitting dialog. Altering them will
//CurveFitDialog/ make the function less convenient to work with in the Curve Fitting dialog.
//CurveFitDialog/ Equation:
//CurveFitDialog/ f(xx) = (A1*(1110*(exp(0.5*(S1/T1)^2-(xx-x1)/T1))*erfc(0.707*(S1/T1-(xx-x1)/S1)))+200*(exp(-((xx-x3)/S3)^2)))+(A2*(373*(exp(0.5*(S2/T2)^2-(xx-x2)/T2)*erfc(0.707*(S2/T2)-(xx-x2)/S2))+200*(exp(-((xx-x3)/S3)^2)))
//CurveFitDialog/ End of Equation
//CurveFitDialog/ Independent Variables 1
//CurveFitDialog/ xx
//CurveFitDialog/ S1 --> Set equal to 0.064
//CurveFitDialog/ T1 --> Set equal to 0.23
//CurveFitDialog/ S2 --> Set equal to 0.12
//CurveFitDialog/ T2 --> Set equal to 0.09
//CurveFitDialog/ x1 --> Set equal to 1.92
//CurveFitDialog/ x2 --> Set equal to 2.17
//CurveFitDialog/ Coefficients 2
//CurveFitDialog/ ww[0] = A1
//CurveFitDialog/ ww[1] = A2
w1 = exp(0.5*(S1/T1)^2-(xx-x1)/T1)
w2 = erfc(0.707*(S1/T1-(xx-x1)/S1))
w3 = 1110
w4 = exp(0.5*(S2/T2)^2-(xx-x2)/T2)
w5 = erfc(0.707*(S2/T2-(xx-x2)/S2))
w6 = 373
w7 = 200
w8 = exp(-((xx-x3)/S3)^2)
w9 = 200
w10 = 560
w11 = ww[0]
w12 = ww[1]
w13 = w3*(w1*w2)+w7*w8
w14 = w6*(w4*w5)+w10*w8
qq = w11*w13+w12*w14
return qq
for(i=0; i<k; i+=1)
String path = "root:ps_2:input" //directory of my 2D wave
Display $path
DFREF dfr = $path
k=DimSize ($path,1)
for(i=0; i<k; i+=1)
Display dfr:f35[][i]
FuncFit/NTHR=0/Q/N FitMySpectrum (f35[][i]) //f35 is the input, 2D wave to which fitting is done
endfor
You also need to declare the wave you want to fit using
Wave MyWave=dfr:xxxxx
December 3, 2017 at 11:59 pm - Permalink
You are missing and "end" after the return command and you need to start a new function afterwards.
You cannot display an entire data folder:
Display $path
vs f35 (Maybe you want to display an empty graph: that would be a singledispay
command.)Try to get single tasks working one by one. Start with the fit function, add the displaying of the raw data, do manual fitting, store the coefficients, put this in loops, append the fits, and list the coefficients.
HJ
December 4, 2017 at 01:08 am - Permalink
You are trying to add code to the fitting function that belongs in a separate function. You need a driver function that looks up the wave with the data in it and loops over the columns, calling FuncFit. The loop in your code needs to be put into a separate function, and then you would call that function.
John Weeks
WaveMetrics, Inc.
support@wavemetrics.com
December 4, 2017 at 09:58 am - Permalink
Thanks Olelytken, HJ, and John for the suggestions. I slightly modified the code to separate fit function and raw spectrum. I can now display raw spectrums in a separate window and I want to call my fit function to fit each of these spectrums. While I am still working on how to call the fit func and append it to the raw spectrum, I would appreciate any insight/comments/help topic to get started.
#include <WMBatchCurveFitIM>
Function FitMySpectrum(ww, xx) : FitFunc
WAVE ww
Variable xx
Variable w1=0,w2=0,w3=0,w4=0,w5=0,w6=0,w7=0,w8=0,w9=0,w10=0,w11=0,w12=0,w13,w14=0,qq=0
Variable T1=0.24, S1=0.047, T2=0.24, S2=0.047, x1=1.91, x2=1.918, S3=0.26, x3=2.20
//Variable k
//Make/O/N=46 A1 //wave to store coefficients A1
//Make/O/N=46 A2 //wave to store coefficients A2
//Make/O/N=(41,46) fitdata //wavet to store results from fitting
//CurveFitDialog/ These comments were created by the Curve Fitting dialog. Altering them will
//CurveFitDialog/ make the function less convenient to work with in the Curve Fitting dialog.
//CurveFitDialog/ Equation:
//CurveFitDialog/ f(xx) = (A1*(1110*(exp(0.5*(S1/T1)^2-(xx-x1)/T1))*erfc(0.707*(S1/T1-(xx-x1)/S1)))+200*(exp(-((xx-x3)/S3)^2)))+(A2*(373*(exp(0.5*(S2/T2)^2-(xx-x2)/T2)*erfc(0.707*(S2/T2)-(xx-x2)/S2))+200*(exp(-((xx-x3)/S3)^2)))
//CurveFitDialog/ End of Equation
//CurveFitDialog/ Independent Variables 1
//CurveFitDialog/ xx
//CurveFitDialog/ S1 --> Set equal to 0.064
//CurveFitDialog/ T1 --> Set equal to 0.23
//CurveFitDialog/ S2 --> Set equal to 0.12
//CurveFitDialog/ T2 --> Set equal to 0.09
//CurveFitDialog/ x1 --> Set equal to 1.92
//CurveFitDialog/ x2 --> Set equal to 2.17
//CurveFitDialog/ Coefficients 2
//CurveFitDialog/ ww[0] = A1
//CurveFitDialog/ ww[1] = A2
w1 = exp(0.5*(S1/T1)^2-(xx-x1)/T1)
w2 = erfc(0.707*(S1/T1-(xx-x1)/S1))
w3 = 1110
w4 = exp(0.5*(S2/T2)^2-(xx-x2)/T2)
w5 = erfc(0.707*(S2/T2-(xx-x2)/S2))
w6 = 373
w7 = 200
w8 = exp(-((xx-x3)/S3)^2)
w9 = 200
w10 = 560
w11 = ww[0]
w12 = ww[1]
w13 = w3*(w1*w2)+w7*w8
w14 = w6*(w4*w5)+w10*w8
qq = w11*w13+w12*w14
return qq
End
Function RawSpectrumDisplay()
variable i
for(i=0;i<46;i+=1)
String path = "root:ps_2"
Display $path
DFREF dfr = $path
Display dfr:f35[][i] //These are the data I want to fit
endfor
End
December 4, 2017 at 10:29 am - Permalink
Thanks Olelytken, HJ, and John for the suggestions. I slightly modified the code to separate fit function and raw spectrum. I can now display raw spectrums in a separate window and want to call my fit function to fit each of these spectrums. While I am still working on how to call the fit func and append it to the raw spectrum, I would appreciate any insight/comments/help topic to get started.
#include <WMBatchCurveFitIM>
Function FitMySpectrum(ww, xx) : FitFunc
WAVE ww
Variable xx
Variable w1=0,w2=0,w3=0,w4=0,w5=0,w6=0,w7=0,w8=0,w9=0,w10=0,w11=0,w12=0,w13,w14=0,qq=0
Variable T1=0.24, S1=0.047, T2=0.24, S2=0.047, x1=1.91, x2=1.918, S3=0.26, x3=2.20
//Variable k
//Make/O/N=46 A1 //wave to store coefficients A1
//Make/O/N=46 A2 //wave to store coefficients A2
//Make/O/N=(41,46) fitdata //wavet to store results from fitting
//CurveFitDialog/ These comments were created by the Curve Fitting dialog. Altering them will
//CurveFitDialog/ make the function less convenient to work with in the Curve Fitting dialog.
//CurveFitDialog/ Equation:
//CurveFitDialog/ f(xx) = (A1*(1110*(exp(0.5*(S1/T1)^2-(xx-x1)/T1))*erfc(0.707*(S1/T1-(xx-x1)/S1)))+200*(exp(-((xx-x3)/S3)^2)))+(A2*(373*(exp(0.5*(S2/T2)^2-(xx-x2)/T2)*erfc(0.707*(S2/T2)-(xx-x2)/S2))+200*(exp(-((xx-x3)/S3)^2)))
//CurveFitDialog/ End of Equation
//CurveFitDialog/ Independent Variables 1
//CurveFitDialog/ xx
//CurveFitDialog/ S1 --> Set equal to 0.064
//CurveFitDialog/ T1 --> Set equal to 0.23
//CurveFitDialog/ S2 --> Set equal to 0.12
//CurveFitDialog/ T2 --> Set equal to 0.09
//CurveFitDialog/ x1 --> Set equal to 1.92
//CurveFitDialog/ x2 --> Set equal to 2.17
//CurveFitDialog/ Coefficients 2
//CurveFitDialog/ ww[0] = A1
//CurveFitDialog/ ww[1] = A2
w1 = exp(0.5*(S1/T1)^2-(xx-x1)/T1)
w2 = erfc(0.707*(S1/T1-(xx-x1)/S1))
w3 = 1110
w4 = exp(0.5*(S2/T2)^2-(xx-x2)/T2)
w5 = erfc(0.707*(S2/T2-(xx-x2)/S2))
w6 = 373
w7 = 200
w8 = exp(-((xx-x3)/S3)^2)
w9 = 200
w10 = 560
w11 = ww[0]
w12 = ww[1]
w13 = w3*(w1*w2)+w7*w8
w14 = w6*(w4*w5)+w10*w8
qq = w11*w13+w12*w14
return qq
End
Function RawSpectrumDisplay()
variable i
for(i=0;i<46;i+=1)
String path = "root:ps_2"
Display $path
DFREF dfr = $path
Display dfr:f35[][i] //These are the data I want to fit
endfor
End
December 4, 2017 at 10:26 am - Permalink
I tried to make the use of driver function within which I looped over the columns of my data. Calling the driver function from the command line just displays the raw data however; I can see the fitting coefficients for the last spectrum being displayed on the command panel. Also I copied the second last line of my code form the help menu and do not understand what ''/D /STRC=" is referring to. I would be thankful for any suggestion on how to get the fits appended over the raw spectrum.
#include <WMBatchCurveFitIM>
Function FitMySpectrum(ww, xx) : FitFunc
WAVE ww
Variable xx
Variable w1=0,w2=0,w3=0,w4=0,w5=0,w6=0,w7=0,w8=0,w9=0,w10=0,w11=0,w12=0,w13,w14=0,qq=0
Variable T1=0.24, S1=0.047, T2=0.24, S2=0.047, x1=1.91, x2=1.918, S3=0.26, x3=2.20
//Variable i, k
//Make/O/N=46 A1 //wave to store coefficients A1
//Make/O/N=46 A2 //wave to store coefficients A2
//Make/O/N=(41,46) fitdata //wavet to store results from fitting
//CurveFitDialog/ These comments were created by the Curve Fitting dialog. Altering them will
//CurveFitDialog/ make the function less convenient to work with in the Curve Fitting dialog.
//CurveFitDialog/ Equation:
//CurveFitDialog/ f(xx) = (A1*(1110*(exp(0.5*(S1/T1)^2-(xx-x1)/T1))*erfc(0.707*(S1/T1-(xx-x1)/S1)))+200*(exp(-((xx-x3)/S3)^2)))+(A2*(373*(exp(0.5*(S2/T2)^2-(xx-x2)/T2)*erfc(0.707*(S2/T2)-(xx-x2)/S2))+200*(exp(-((xx-x3)/S3)^2)))
//CurveFitDialog/ End of Equation
//CurveFitDialog/ Independent Variables 1
//CurveFitDialog/ xx
//CurveFitDialog/ S1 --> Set equal to 0.064
//CurveFitDialog/ T1 --> Set equal to 0.23
//CurveFitDialog/ S2 --> Set equal to 0.12
//CurveFitDialog/ T2 --> Set equal to 0.09
//CurveFitDialog/ x1 --> Set equal to 1.92
//CurveFitDialog/ x2 --> Set equal to 2.17
//CurveFitDialog/ Coefficients 2
//CurveFitDialog/ ww[0] = A1
//CurveFitDialog/ ww[1] = A2
w1 = exp(0.5*(S1/T1)^2-(xx-x1)/T1)
w2 = erfc(0.707*(S1/T1-(xx-x1)/S1))
w3 = 1110
w4 = exp(0.5*(S2/T2)^2-(xx-x2)/T2)
w5 = erfc(0.707*(S2/T2-(xx-x2)/S2))
w6 = 373
w7 = 200
w8 = exp(-((xx-x3)/S3)^2)
w9 = 200
w10 = 560
w11 = ww[0]
w12 = ww[1]
w13 = w3*(w1*w2)+w7*w8
w14 = w6*(w4*w5)+w10*w8
qq = w11*w13+w12*w14
return qq
End
Function RawSpectrumDisplay()
variable i
for(i=0;i<46;i+=1)
String path = "root:ps_2"
Display $path
DFREF dfr = $path
Display dfr:f35[][i] //These are the data I want to fit
endfor
End
// The structure definition
Structure EMGFitStruct
Wave cw // required coefficient wave that should be guessed initially
wave yw // y value which should come from f35[][i]
wave xw // x values already scaled in f35[][i]
variable j
EndStructure
// The driver function that calls Func Fit:
Function EMGFitDriver(cw,yw,xw)
//FUNCREF RawSpectrumDisplay
Wave cw // coefficient waves
Wave yw
Wave xw
variable i
for(i=0;i<46;i+=1)
String path = "root:ps_2"
Display $path
DFREF dfr = $path
Display dfr:f35[][i] //These are the data I want to fit
endfor
STRUCT EMGFitStruct f35
FuncFit FitMySpectrum, cw, yw /X=xw /D /STRC=f35
END
December 4, 2017 at 12:33 pm - Permalink
#include <WMBatchCurveFitIM>
// The structure definition
Structure EMGFitStruct
Wave cw // required coefficient wave that should be guessed initially
wave yw // y value which should come from f35[][i]
wave xw // x values already scaled in f35[][i]
//variable xx, qq
//wave ww
wave inputw
EndStructure
Function FitMySpectrum(ww,xx) : FitFunc
//STRUCT EMGFitStruct &s
//variable xx
//variable qq
WAVE ww
//wave qq
Variable xx
Variable w1=0,w2=0,w3=0,w4=0,w5=0,w6=0,w7=0,w8=0,w9=0,w10=0,w11=0,w12=0,w13,w14=0,qq=0
Variable T1=0.24, S1=0.047, T2=0.24, S2=0.047, x1=1.91, x2=1.918, S3=0.26, x3=2.20
//Variable i, k
//Make/O/N=46 A1 //wave to store coefficients A1
//Make/O/N=46 A2 //wave to store coefficients A2
//Make/O/N=(41,46) qq //wavet to store results from fitting
//CurveFitDialog/ These comments were created by the Curve Fitting dialog. Altering them will
//CurveFitDialog/ make the function less convenient to work with in the Curve Fitting dialog.
//CurveFitDialog/ Equation:
//CurveFitDialog/ f(xx) = (A1*(1110*(exp(0.5*(S1/T1)^2-(xx-x1)/T1))*erfc(0.707*(S1/T1-(xx-x1)/S1)))+200*(exp(-((xx-x3)/S3)^2)))+(A2*(373*(exp(0.5*(S2/T2)^2-(xx-x2)/T2)*erfc(0.707*(S2/T2)-(xx-x2)/S2))+200*(exp(-((xx-x3)/S3)^2)))
//CurveFitDialog/ End of Equation
//CurveFitDialog/ Independent Variables 1
//CurveFitDialog/ xx
//CurveFitDialog/ S1 --> Set equal to 0.064
//CurveFitDialog/ T1 --> Set equal to 0.23
//CurveFitDialog/ S2 --> Set equal to 0.12
//CurveFitDialog/ T2 --> Set equal to 0.09
//CurveFitDialog/ x1 --> Set equal to 1.92
//CurveFitDialog/ x2 --> Set equal to 2.17
//CurveFitDialog/ Coefficients 2
//CurveFitDialog/ ww[0] = A1
//CurveFitDialog/ ww[1] = A2
w1 = exp(0.5*(S1/T1)^2-(xx-x1)/T1)
w2 = erfc(0.707*(S1/T1-(xx-x1)/S1))
w3 = 1110
w4 = exp(0.5*(S2/T2)^2-(xx-x2)/T2)
w5 = erfc(0.707*(S2/T2-(xx-x2)/S2))
w6 = 373
w7 = 200
w8 = exp(-((xx-x3)/S3)^2)
w9 = 200
w10 = 560
w11 = ww[0]
w12 = ww[1]
w13 = w3*(w1*w2)+w7*w8
w14 = w6*(w4*w5)+w10*w8
qq = w11*w13+w12*w14
return qq
//variable i
//for(i=0;i<2;i+=1)
//qq=qq[i]
//endfor
//return qq
End
Function RawSpectrumDisplay (coef, f35)
wave coef
wave f35
variable xx
STRUCT EMGFitStruct fs
variable i
for(i=0;i<2;i+=1)
String path = "root:ps_2"
Display $path
DFREF dfr = $path
Display dfr:f35[][i] //These are the data I want to fit
wave fs.inputw=dfr:f35
endfor
End
// The driver function that calls Func Fit:
Function EMGFitDriver(W_coef,f35,xx)
//FUNCREF RawSpectrumDisplay
Wave W_coef // coefficient waves
Wave f35
Wave xx
variable i
for(i=0;i<3;i+=1)
String path = "root:ps_2:f35[][i]"
Display $path
DFREF dfr = $path
Display dfr:f35[][i] //These are the data I want to fit
endfor
STRUCT EMGFitStruct fs
FuncFit FitMySpectrum, W_coef, f35[][i]/D /STRC=fs //FitMySpectrum is the fitting function
December 5, 2017 at 08:08 am - Permalink
Write a function that only displays the raw data.
Write a fit-function that works with the 'Analysis->Curve Fitting' interface. (I don't think you need the structure feature)
From now on have a close look to the history window:
Fit your data manually.
Manage to store the resulting coefficients manually.
Store the resulting fit wave in a 2d-wave similar to your data.
Add the resulting new fit manually to the display (including coloring etc).
Use the information from the history to modify the displaying function mentioned above. In short words: add these commands within the display loop.
Display the coefficient wave and you should be almost happy.
J
December 5, 2017 at 08:55 am - Permalink
TRK
Thank HJ ! 'Analysis->Curve Fitting' works beautifully and I extract the coefficients from the history/textbox in the graph window. Here I am looking for some automate fitting as I have hundreds of columns in my data (though I am running my loops only for 3 cols in above code for check) and will encounter similar sets of data in future. Initially I tried using simple function but it did not work for me and based on the comments and forum posts, I decided to use structure function. I would be happy and appreciate very much for the help on displaying the fit to every raw spectrums than just the last raw spectum in the loop.
Thanks,
Tika
December 5, 2017 at 09:26 am - Permalink
Once this works, add that command right after the funcfit command in your procedure.
December 5, 2017 at 09:49 am - Permalink
Once again I need help to extract data for each iteration in FuncFit. I checked FuncFit section in help file and it does not explain how to extract results from FuncFit. I just keep on getting the fit only for the last iteration of my loop; it looks like overwriting the fit. Here is what I did; After FuncFit command line I duplicated the wave to get fit data for each iteration (but this does not seem to duplicate the fit data). Then I created a 2D wave to store the fit data. Here is the segment of the code. Thanks
Wave W_coef // coefficient waves in my folder
Wave f35
Wave xx
Wave fit_f35
//Wave outputw
variable i,n
for(i=0;i<10;i+=1)
String path = "root:ps_2:f35[][i]"
Display $path
DFREF dfr = $path
Display dfr:f35[][i] //These are the data I want to fit
endfor
STRUCT EMGFitStruct fs
FuncFit FitMySpectrum, W_coef, f35[][i]/D/STRC=fs //FitMySpectrum is the fitting function
duplicate/o fit_f35, $("fit"+num2istr(n)+"_"+"f35") //I am assuming FucnFit returns the fit data as fit_f35 wave by default because this is what I see in my folder
Make/O/N=(41,46) fitdata
fitdata = fit_f35
SetAxis bottom 1.64, 3.2
//Display outputw[][i]
AppendToGraph f35[][i],fit_f35[][i]
December 5, 2017 at 04:05 pm - Permalink
Wave W_coef // coefficient waves
Wave f35
Wave xx
Make/O/N = (46)fit_f35
//Wave outputw
variable i,n
for(i=0;i<2;i+=1)
String path = "root:ps_2:f35[][i]"
Display $path
DFREF dfr = $path
Display dfr:f35[][i] //These are the data I want to fit
STRUCT EMGFitStruct fs
FuncFit FitMySpectrum, W_coef, f35[][i]/D/STRC=fs //FitMySpectrum is the fitting function
duplicate/o fit_f35, $("fit"+num2istr(n)+"_"+"f35[][i]")
display fit_f35[][i], $("fit"+num2istr(i)+"_"+"f35[][i]")
//wave fitdata=$("fit"+num2istr(n)+"_"+"f35")
Make/O/N=(41,46) fitdata
wave fitdata=$("fit"+num2istr(0)+"_"+"f35[][i]")
SetAxis bottom 1.64, 3.2
endfor
End
December 5, 2017 at 08:41 pm - Permalink
Run your code in the command window and see what happens.
I guess you will have more windows created than you actually need.
I guess your storage wave will be overwritten.
Ill guide you though the development. I assume there are things in the code (e.g., "/TN=") that you don't understand. This is absolutely fine. Consult the manual and if you have additional questions, please feel free to ask.
Let's start.
We do not want to mess with the current data folder (DF): Hence, we save the current one, switch to our data, and at the end of the function we switch back.
I assume you want one window per data set (column). Your initial displaying function should actually look like this:
String DataFolderName, DataName
// This is the full path to your data in two pieces as text. It could also be passed as a wave reference, but this makes the code a little more complicated
Variable i
DFRef OldDF=GetDataFolderDFR() //Store the current DF
SetDataFolder $DataFolderName // Set it to our data
Wave Data=$DataName
For (i=0;i<Dimsize(Data,1);i+=1)
Display /K=1 /W=(10*i,10*i+40, 10*i+440, 10*i+40+210) Data[][i]/TN=$(DataName+"_"+Num2str(i))
ModifyGraph mode(test_2)=2, lsize(test_2)=2
EndFor
SetDataFolder OldDF // Switch DF back
End
Have a look at the manual what the K and W flags are doing.
you should be able to call it like this:
ProcessData("root:ps_2:input:", "f35")
The next step would be to include the fitting function. There is no need for storing in the fit function. All it does is getting one x value and calculating one y value. For each point in your data. In each iteration. You see, it gets called many times. You want the data after all of this is finished. In the meantime, avoid any overhead: Remove all 'Make' command from the fit-function. Variables i & k are also not used - remove this line as well.
Now you want to prepare for storing something. Your fits will share the same data structure as your raw data. Why not use this property?
Duplicate /O Data $("fit_"+DataName")
and show this new wave to the compiler
Wave Fits=$("fit_"+DataName)
Also we need to create the waves for the coefficients; their size is the same as numbers of columns in your data:
Make /O /N=(DimSize(Data,1)) A1, A2
Since this is related to 2d-Data, these commands must be outside the loop:
String DataFolderName, DataName
Variable i
DFRef OldDF=GetDataFolderDFR()
SetDataFolder $DataFolderName
Wave Data=$DataName
Duplicate /O Data $("fit_"+DataName) // Replicate the data structure incl. wavc scaling units etc.
Wave Fits=$("fit_"+DataName)
Make /O /N=(DimSize(Data,1)) A1, A2 // Prepare the coefficient storage
For (i=0;i<Dimsize(Data,1);i+=1)
Display /K=1 /W=(10*i,10*i+40, 10*i+440, 10*i+40+210) Data[][i]/TN=$(DataName+"_"+Num2str(i))
ModifyGraph mode(test_2)=2, lsize(test_2)=2
Endfor
SetDataFolder OldDF
End
Now we can add the fitting. You don't need the data structure. A temporary (/FREE) coefficient wave with initial guesses is sufficient. Also prepare a temporary wave for the fit outside the loop.
Make /FREE /N=(Dimsize(Data,0)) Fit
And add the fitting in the loop:
String DataFolderName, DataName
Variable i
DFRef OldDF=GetDataFolderDFR()
SetDataFolder $DataFolderName
Wave Data=$DataName
Duplicate /O Data $("fit_"+DataName)
Wave Fits=$("fit_"+DataName)
Make /O /N=(DimSize(Data,1)) A1, A2
Make /FREE /D /N=2 w_coef={1,2}
Make /FREE /N=(Dimsize(Data,0)) Fit
For (i=0;i<Dimsize(Data,1);i+=1)
Display /K=1 /W=(10*i,10*i+40, 10*i+440, 10*i+40+210) Data[][i]/TN=$(DataName+"_"+Num2str(i))
ModifyGraph mode($(DataName+"_"+Num2str(i)))=2, lsize($(DataName+"_"+Num2str(i)))=2
FuncFit /Q FitMySpectrum W_coef Data[][i] /D=Fit
Endfor
SetDataFolder OldDF
End
If you test this function, you will realize that there is no fit displayed in the windows. This is ok for now. Let's store our results first.
A1[i]=w_coef[0]
A2[i]=w_coef[1]
Certainly this needs to be done right after the fitting.
Next, we nicely display the resulting fit within the loop
ModifyGraph rgb($(DataName+"_Fit_"+Num2str(i)))=(0,0,0)
Finally we add some output of the coefficeints
Display /K=1 /W=(400,310,800,510) A1
Display /K=1 /W=(400,550,800,750) A2
Altogether it should look like this:
String DataFolderName, DataName
Variable i
DFRef OldDF=GetDataFolderDFR()
SetDataFolder $DataFolderName
Wave Data=$DataName
Duplicate /O Data $("fit_"+DataName)
Wave Fits=$("fit_"+DataName)
Make /O /N=(DimSize(Data,1)) A1, A2
Make /FREE /D /N=2 w_coef={1,2}
Make /FREE /N=(Dimsize(Data,0)) Fit
For (i=0;i<Dimsize(Data,1);i+=1)
Display /K=1 /W=(10*i,10*i+40, 10*i+440, 10*i+40+210) Data[][i]/TN=$(DataName+"_"+Num2str(i))
ModifyGraph mode($(DataName+"_"+Num2str(i)))=2, lsize($(DataName+"_"+Num2str(i)))=2
FuncFit /Q FitMySpectrum W_coef Data[][i] /D=Fit
Fits[][i]=Fit[p]
A1[i]=w_coef[0]
A2[i]=w_coef[1]
AppendToGraph Fits[][i]/TN=$(DataName+"_Fit_"+Num2str(i))
ModifyGraph rgb($(DataName+"_Fit_"+Num2str(i)))=(0,0,0)
Endfor
Edit /K=1 /W=(400,40,800,240) A1 A2
Display /K=1 /W=(400,310,800,510) A1
Display /K=1 /W=(400,550,800,750) A2
SetDataFolder OldDF
End
A final remark: FitMySpectrum is a name that suits for this forum. In Igor, it should get a more meaningful name, e.g. NeuronResponse or DoniachSunjic.
Hope this helps,
HJ
December 6, 2017 at 01:16 am - Permalink
Thank you HJ for step by step development towards the goal. This is what I was looking for exactly. I appreciate your time and effort. It is so good to see how flags after the display command can kill the window (/K=1) with no dialog, determine the location and size of the graph (/w) and provide tracename based on the parent wave(/TN).
At the end, I included a residual wave to visualize how good the fit is. Finally, have a couple of quick questions. I checked that the script "Residuals=Data-Fits" works as it calculates after all the waves are collected. At the same time why "Residuals[][i]=Data[i]-Fits[i] or Residuals[][p]=Data[p]-Fits[p]" does not work inside the loop because I was assuming that the subtraction from the raw data takes place for each iteration and is stored in Residuals.
Next, I know that I can chose not to display the graph window by disabling the display command. I noticed there is /k=3 to hide the window. Other than these is there any quick way to kill the graph window rather than doing it manually?
Once again, thank you very much !
Here is how I tried to calculate Residuals,
String DataFolderName, DataName
Variable i
DFRef OldDF=GetDataFolderDFR()
SetDataFolder $DataFolderName //identifies the path for the current data folder
Wave Data=$DataName // wavename to be executed
Duplicate /O Data $("fit_"+DataName)
Wave Fits=$("fit_"+DataName)
Make /O /N=(DimSize(Data,1)) A1, A2
Make /FREE/D /N=2 w_coef={1,2}
Make /FREE/N=(Dimsize(Data,0)) Fit
Duplicate /O Data $("Residual_"+DataName) //////////////
Wave Residuals=$("Residual_"+DataName) //////////////
variable k
k=Dimsize(Data,1)
For (i=0;i<Dimsize(Data,1);i+=1)
Display /K=1 /W=(10*i,10*i+40, 10*i+440, 10*i+40+210) Data[][i]/TN=$(DataName+"_"+Num2str(i))
ModifyGraph mode($(DataName+"_"+Num2str(i)))=2, lsize($(DataName+"_"+Num2str(i)))=2
FuncFit /Q FitMySpectrum W_coef Data[][i] /D=Fit
Fits[][i]=Fit[p]
Residuals=Data-Fits ////////////////residuals are calculated
//Residuals[][i]=Data[i]-Fits[i]
//Residuals[p][q]=Data[p][q]-Fits[p][q] //this does not work
A1[i]=w_coef[0]
A2[i]=w_coef[1]
AppendToGraph Fits[][i]/TN=$(DataName+"_Fit_"+Num2str(i))
ModifyGraph rgb($(DataName+"_Fit_"+Num2str(i)))=(0,0,0)
Endfor
//Residuals=Data-Fits ////////////////////////
Edit /K=1 /W=(400,100,800,240) A1 A2
Display /K=1 /W=(400,310,800,510) A1
Display /K=1 /W=(400,550,800,750) A2
SetDataFolder OldDF
End
<pre><code class="language-igor">
December 6, 2017 at 05:51 pm - Permalink
Calculating the residuals at the end of the procedure should work that way you wrote it.
Inside the loop you have to pay attention to the dimensions of the waves ( number of [...] pairs ), indices (i), automatic loop indices (p, q, x, y, z etc) and "do-it-for-all-points-in-that-dimension" ( [] )
The right syntax would be
You want to calculate for all rows (which correspond to index p) in column i in Residuals (this is the [][i] part) the difference between row p in column p of data and the p-th row in Fit.
Please read the manual about wave assignment:
displayhelptopic "Multidimensional Wave Assignment"
(Copy these commands to the command window and execute them)
Essentially,
KillWindow
does the job of closing a window. Sadly not a list of windows. So you'd need to write a procedure for this.string WindowList
variable i
for (i=0;i<itemsinlist(WindowList);i+=1)
Killwindow $stringfromlist(i,WindowList)
endfor
end
To kill all graphs run
Killwins(winlist("*",";","WIN:1"))
and to kill all tables
Killwins(winlist("*",";","WIN:2"))
However, it might not be necessary to kill the windows, e.g, if you only want to rerun the fitting. In this case, it might be advisable to separate the processing function into smaller fragments that share the same task. Examples would be creation of data structure, displaying, fitting, summarizing, and cleanup. The same holds, if you want to recreate the graphs later without new fitting. I can give you a step by step example later.
I added a little extra to the processing function:
String DataFolderName, DataName
Variable i
DFRef OldDF=GetDataFolderDFR()
SetDataFolder $DataFolderName
Wave Data=$DataName
Duplicate /O Data $("fit_"+DataName), $("res_"+DataName) // You can duplicate to several destination waves; Keep the names short for automatic extension
Wave Residuals=$("res_"+DataName)
Wave Fits=$("fit_"+DataName)
Make /O /N=(DimSize(Data,1)) A1, A2
Make /FREE /D /N=2 w_coef={1,2}
Make /FREE /N=(Dimsize(Data,0)) Fit
For (i=0;i<Dimsize(Data,1);i+=1)
Display /K=1 /W=(10*i,10*i+40, 10*i+440, 10*i+40+210) Data[][i]/TN=$(DataName+"_"+Num2str(i)) as DataName+"_"+num2str(i)
FuncFit /Q FitMySpectrum W_coef Data[][i] /D=Fit
Fits[][i]=Fit[p]
Residuals[][i]=Data[p][i]-Fit[p]
A1[i]=w_coef[0]
A2[i]=w_coef[1]
AppendToGraph Fits[][i]/TN=$(DataName+"_Fit_"+Num2str(i))
AppendToGraph /L=laxRes Residuals[][i]/TN=$(DataName+"_Res_"+Num2str(i))
ModifyGraph mode($(DataName+"_"+Num2str(i)))=2, lsize($(DataName+"_"+Num2str(i)))=2
ModifyGraph mode($(DataName+"_Res_"+Num2str(i)))=2, lsize($(DataName+"_Res_"+Num2str(i)))=2
ModifyGraph freePos(laxRes)={0,bottom}
ModifyGraph axisEnab(left)={0,0.75}
ModifyGraph axisEnab(laxRes)={0.8,1}
ModifyGraph rgb($(DataName+"_Fit_"+Num2str(i)))=(0,0,0)
Endfor
// Residuals=Data-Fits
Edit /K=1 /W=(400,40,800,240) A1 A2
Display /K=1 /W=(400,310,800,510) A1
Display /K=1 /W=(400,550,800,750) A2
SetDataFolder OldDF
End
HJ
December 7, 2017 at 05:05 am - Permalink
Thank you.
December 7, 2017 at 11:48 pm - Permalink
It is hard going. The first time through it won't make a lot of sense. But some will stick and the next time through more of it will make sense. And it helps to have a project like yours to apply your learning to as you read.
John Weeks
WaveMetrics, Inc.
support@wavemetrics.com
December 8, 2017 at 09:17 am - Permalink
I have two different 2D waves with same number of rows (energy) but different columns (time). These waves overlap in columns, that is first wave is data from 40ps to -5ps with deltay= -1 ps and second wave is data from 300 ps to -50 ps with deltay = -10 ps. I concatenated the two waves to get a new wave destWave( 46*72) which now has data from 300 to 40 ps (col0 to col26) with deltay = -10 and 39 ps to -5 ps (col27 to col72). Since the concatenated wave takes scaling of first wave in the list, ywave scaling mismatches. I tried to run a do while loop but it returns 'nan' in col scaling and the graph is not displayed. Can we scale a 2D wave with two different deltay values?
Here is the script I tried with.
Wave data2, data1
Wave destWave
variable k, dy
//dy =destWave[0] - destWave[1]
Duplicate /O/R=[][1,*] data1, sub1
Duplicate /O/R=[][0,26] data2, sub2
k=Dimsize(destWave,1)
Concatenate/O/NP=1 {sub2,sub1}, destWave
//for (k=0;k<26;k+=1)
do (k<27) //col 0 to col 26; starty 300, deltay= -10, endy 40;
dy = 300-10*k //col 27 to col72, starty39, deltay= -1, endy= -5
while (k>26)
dy = 39-1*k
SetScale/P y, destWave[0], dy, destWave
//SetAxis/A/R left
//endfor
End
<pre><code class="language-igor">
December 9, 2017 at 09:43 am - Permalink
John Weeks
WaveMetrics, Inc.
support@wavemetrics.com
December 11, 2017 at 09:33 am - Permalink
Thanks.
December 11, 2017 at 09:59 pm - Permalink
John Weeks
WaveMetrics, Inc.
support@wavemetrics.com
December 12, 2017 at 09:16 am - Permalink
Please read
displayhelptopic "The XY Model of Data"
December 13, 2017 at 12:52 am - Permalink
First, we need to identify each command with a logical set of actions.
String DataFolderName, DataName
Variable i
DFRef OldDF=GetDataFolderDFR() // Preparation
SetDataFolder $DataFolderName // Preparation
Wave Data=$DataName // Preparation
Duplicate /O Data $("fit_"+DataName), $("res_"+DataName) // Data management
Wave Residuals=$("res_"+DataName) // Preparation
Wave Fits=$("fit_"+DataName) // Preparation
Make /O /N=(DimSize(Data,1)) A1, A2 // Data management
Make /FREE /D /N=2 w_coef={1,2} // Preparation
Make /FREE /N=(Dimsize(Data,0)) Fit // Preparation
For (i=0;i<Dimsize(Data,1);i+=1) // Loop
Display /K=1 /W=(10*i,10*i+40, 10*i+440, 10*i+40+210) Data[][i]/TN=$(DataName+"_"+Num2str(i)) as DataName+"_"+num2str(i) // Display
FuncFit /Q FitMySpectrum W_coef Data[][i] /D=Fit // Fitting
Fits[][i]=Fit[p] // Fitting
Residuals[][i]=Data[p][i]-Fit[p] // Fitting
A1[i]=w_coef[0] // Fitting
A2[i]=w_coef[1] // Fitting
AppendToGraph Fits[][i]/TN=$(DataName+"_Fit_"+Num2str(i)) // Display
AppendToGraph /L=laxRes Residuals[][i]/TN=$(DataName+"_Res_"+Num2str(i)) // Display
ModifyGraph mode($(DataName+"_"+Num2str(i)))=2, lsize($(DataName+"_"+Num2str(i)))=2 // Display
ModifyGraph mode($(DataName+"_Res_"+Num2str(i)))=2, lsize($(DataName+"_Res_"+Num2str(i)))=2 // Display
ModifyGraph freePos(laxRes)={0,bottom} // Display
ModifyGraph axisEnab(left)={0,0.75} // Display
ModifyGraph axisEnab(laxRes)={0.8,1} // Display
ModifyGraph rgb($(DataName+"_Fit_"+Num2str(i)))=(0,0,0) // Display
Endfor // Loop
Edit /K=1 /W=(400,40,800,240) A1 A2 // Summary
Display /K=1 /W=(400,310,800,510) A1 // Summary
Display /K=1 /W=(400,550,800,750) A2 // Summary
SetDataFolder OldDF // Clean up
End
The next steps are to create a copy of the function for each set, rename them appropriately, and delete parts that don't match. Be careful with loops and preparations! They might appear in several places.
(This example is a little overkill, I know. However, it's good for illustration.)
String DataFolderName, DataName
Variable i
DFRef OldDF=GetDataFolderDFR() // Preparation
SetDataFolder $DataFolderName // Preparation
Wave Data=$DataName // Preparation
Duplicate /O Data $("fit_"+DataName), $("res_"+DataName) // Data management
Make /O /N=(DimSize(Data,1)) $(DataName+"_A1"), $(DataName+"_A2") // Data management
SetDataFolder OldDF // Clean up
End
Function FitData(DataFolderName, DataName)
String DataFolderName, DataName
Variable i
DFRef OldDF=GetDataFolderDFR() // Preparation
SetDataFolder $DataFolderName // Preparation
Wave Data=$DataName // Preparation
Wave Residuals=$("res_"+DataName) // Preparation
Wave Fits=$("fit_"+DataName) // Preparation
Wave A1= $(DataName+"_A1"), A2= $(DataName+"_A2") // Preparation <--- This has changed !
Make /FREE /D /N=2 w_coef={1,2} // Preparation
Make /FREE /N=(Dimsize(Data,0)) Fit // Preparation
For (i=0;i<Dimsize(Data,1);i+=1) // Loop
FuncFit /Q FitMySpectrum W_coef Data[][i] /D=Fit // Fitting
Fits[][i]=Fit[p] // Fitting
Residuals[][i]=Data[p][i]-Fit[p] // Fitting
A1[i]=w_coef[0] // Fitting
A2[i]=w_coef[1] // Fitting
Endfor // Loop
SetDataFolder OldDF // Clean up
End
Function DisplayData(DataFolderName, DataName)
String DataFolderName, DataName
Variable i
DFRef OldDF=GetDataFolderDFR() // Preparation
SetDataFolder $DataFolderName // Preparation
Wave Data=$DataName // Preparation
Wave Residuals=$("res_"+DataName) // Preparation
Wave Fits=$("fit_"+DataName) // Preparation
For (i=0;i<Dimsize(Data,1);i+=1) // Loop
Display /K=1 /W=(10*i,10*i+40, 10*i+440, 10*i+40+210) Data[][i]/TN=$(DataName+"_"+Num2str(i)) as DataName+"_"+num2str(i) // Display
AppendToGraph Fits[][i]/TN=$(DataName+"_Fit_"+Num2str(i)) // Display
AppendToGraph /L=laxRes Residuals[][i]/TN=$(DataName+"_Res_"+Num2str(i)) // Display
ModifyGraph mode($(DataName+"_"+Num2str(i)))=2, lsize($(DataName+"_"+Num2str(i)))=2 // Display
ModifyGraph mode($(DataName+"_Res_"+Num2str(i)))=2, lsize($(DataName+"_Res_"+Num2str(i)))=2 // Display
ModifyGraph freePos(laxRes)={0,bottom} // Display
ModifyGraph axisEnab(left)={0,0.75} // Display
ModifyGraph axisEnab(laxRes)={0.8,1} // Display
ModifyGraph rgb($(DataName+"_Fit_"+Num2str(i)))=(0,0,0) // Display
Endfor // Loop
SetDataFolder OldDF // Clean up
End
Function SummaryData(DataFolderName, DataName)
String DataFolderName, DataName
Variable i
DFRef OldDF=GetDataFolderDFR() // Preparation
SetDataFolder $DataFolderName // Preparation
Wave Data=$DataName // Preparation
Duplicate /O Data $("fit_"+DataName), $("res_"+DataName) // Data management
Wave Residuals=$("res_"+DataName) // Preparation
Wave Fits=$("fit_"+DataName) // Preparation
Wave A1= $(DataName+"_A1"), A2= $(DataName+"_A2") // Preparation <--- This has changed from creation to reference
Edit /K=1 /W=(400,40,800,240) A1 A2 // Summary
Display /K=1 /W=(400,310,800,510) A1 // Summary
Display /K=1 /W=(400,550,800,750) A2 // Summary
SetDataFolder OldDF // Clean up
End
Now you have to run the functions one by one
FitData("root:sdf:", "test")
DisplayData("root:sdf:", "test")
SummaryData("root:sdf:", "test")
Note: If you run one of the latter functions on a new set of data, it will fail. It expects the prepared waves to exist, which is not the case for a new data set. Here, some data structure checking would be required to protect the data from the user. Have a look in the manual and read about
if ... else ... endif
-structures and theexists / waveexists
functions.You might have realized that the original function is gone at this point. It was so nice, that it did all of the work with one command. Let's bring it back:
String DataFolderName, DataName
PrepareData(DataFolderName, DataName)
FitData(DataFolderName, DataName)
DisplayData(DataFolderName, DataName)
SummaryData(DataFolderName, DataName)
End
It simply takes your input and calls the other functions. What is it good for? We had the same result using one more complex function.
On the one hand, it gives the programmer a better view on his code. Once the simple functions reach a couple of hundred lines you are happy that it's not in the thousands.
On the other hand, we can re-run single steps!
Assume you did a full run on the data, modified the graphs, added some manual annotation, ..... and you figured out that there was a small typo in the fit function. You don't want to kill all the nice figures. Actually, all you need is to re-run the fit and check your remarks in the graph. This is possible now. The data structure is prepared from the last run. You can call the fitting by itself. Igor updates the graphs as the corresponding waves are altered.
As a next logical step, we want to kill (or hide but this takes some more checking) the fit result windows. In order to do this, we should give the windows names:
and prepare a function that kills the windows related to a data set (based on that name):
String DataFolderName, DataName
Variable i
Wave Data=$(DataFolderName+DataName) // Preparation
For (i=0;i<Dimsize(Data,1);i+=1) // Loop
Killwindow $(DataName+"_"+Num2str(i))
EndFor
End
Be advised not to display the fit windows twice. Killing only works on the first instance. To fix this, you would need some checkup during the displaying routine that skips already displayed fits.
Happy fitting,
HJ
December 13, 2017 at 01:47 am - Permalink
One thing I would like to know is about unequal spaced y wave scaling. As Johnweeks mentioned, I created a y wave scaling (time) called "ywave" to match time for cols in 2D data, "Data". Then I run the following command: "Display; AppendImage Data vs {*,ywave}". Putting cursor in the image now displays the actual time and data however y scale in the displayed image does not match. Where am I doing wrong?
Thank you.
Tika
December 16, 2017 at 05:45 pm - Permalink
You can also send it by private message to me directly. (In the forum, click on my name and go to contact)
HJ
December 17, 2017 at 12:40 am - Permalink
You will recall that my original post about this topic said:
The exception to this statement, that I hinted at but didn't explain, is when you use an auxiliary wave to display an image. For an image, you need the time wave to contain NCOLS+1 rows. This is because each "pixel" in an image display needs a position for the left *and* right edges; the final right edge of the last column requires one extra element.
Please read about this in our help: DisplayHelpTopic "Image X and Y Coordinates"
John Weeks
WaveMetrics, Inc.
support@wavemetrics.com
December 18, 2017 at 10:50 am - Permalink
Thanks,
Tika
December 18, 2017 at 11:36 am - Permalink
Happy Programming,
HJ
December 18, 2017 at 04:47 pm - Permalink