W_sigma not being updated

Share icon

I have a procedure that involves multiple curve fittings. I call the fitting function with the command

FuncFit/H="01"/NTHR=0 cloudExpansion W_coefX averagesX /X=averagesXTimes /W=averagesXSigmas

and the sigma for each coefficient associated with the fitting is saved in W_sigma. However, when I call the function again on another wave

FuncFit/NTHR=0 cloudExpansion W_coefY averagesY /Y=averagesYTimes /W=averagesYSigmas

W_sigma still contains exactly the same value as it did after the first fit (even though different fit values are produced, and the data stored in averagesXTimes, averagesYTimes, averagesYSigmas, and averagesXSigmas are all different). The coefficients saved in W_coefY and W_coefX are different, but how can W_sigma not change? If I kill the wave W_sigma before calling the second FuncFit, it is not recreated.
Is the current data folder by any chance different when you perform the two fits? W_Sigma will be created in the current data folder. If you view W_sigma in a table, then change the current data folder, the table won't change because that's not the same W_sigma.

If that's not the explanation, please post a copy of your experiment file so I can investigate.

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com

January 20, 2014 at 09:32 am - Permalink

I don't believe the data folder is ever being changed. Everything is being handled in the root folder.

Here is my code:
#pragma rtGlobals=3     // Use modern global access method and strict wave access.
Function openTemperatureData(pathname) 
    string pathname
    variable binariesIndex = 0
    variable stringIndex
    string currentFilename
   
   
    variable uselessNumber, timeInMs, iterationNumber
    // Get a semicolon-separated list of all files in the folder
    String fileList = IndexedFile($pathName, -1, ".dat")
   

    make/N=500 times//A wave storing the times associated with data points in data waves
    make/N=500 iterations//A wave storing the iteration number associated with data points in data waves

   
    make/N=500 timesToUse//A wave storing the iteration number associated with data points in data waves
    make/N=500 FWHMsXToUse//A wave storing the FWHMs associated with each data wave
    make/N=500 FWHMsYToUse
    make/N=500 FWHMsXSigmaToUse
    make/N=500 FWHMsYSigmaToUse
    Make/D/O/N=7 waveSizeCoefficients
    make averagesX//The average FWHM at each timestep
    make averagesY
    // Sort using combined alpha and numeric sort
    fileList = SortList(fileList, ";", 16)
    variable filledPoint = 0
    variable filledPointToUse = 0
    make W_sigma
    variable cropSizeX
    variable cropSizeY
    Make/N=(200,200)/W cropped
    variable lastTime = 0////The last time that holds data that can be fit
    variable firstTime = 9999 //The first time that holds data that can be fit 
    variable V_FitError//Holds error status of fit routines.  
    make temp
    make temp3
   
   
    for(binariesIndex = 0; binariesIndex < itemsInList(fileList)-1; binariesIndex += 1)//itemsInList(fileList)-1
        currentFilename = StringFromList(binariesIndex,fileList)
        stringIndex=0
        //FindListItem(currentFilename, fileList)
        sscanf currentFilename, "%d_%dms_a_%d", uselessNumber, timeInMs, iterationNumber
       
        GBLoadWave/Q/B=0/A=binary/T={16,16}/W=1/P=$pathName currentFilename
        string wavename = "binary" + num2str(filledPoint)  
        times[filledPoint] = timeInMs
        Redimension/N=(640,480) $wavename
        wave temp=$wavename
        cropped = temp[x+350][y+100]
       
               
        iterations[filledPoint]= iterationNumber
        CurveFit/NTHR=0/Q/N Gauss2D kwCWave=waveSizeCoefficients,  cropped
        if(V_FitError > 0.0002)
            print "Error encountered when attempting to fit to " + currentFilename+".  "  
            print "This is usually caused by a picture that does not contain a coherent cloud.  This file will not be used in analysis.  "
            V_FitError = 0
        endif

           
        if((waveSizeCoefficients[3] < 100) && (waveSizeCoefficients[3] > 0.0001))
            FWHMsXToUse[filledPointToUse]=waveSizeCoefficients[3]*6.7*3*10^-6
            timesToUse[filledPointToUse]=times[filledPoint]
            if (times[filledPoint] > lastTime)
                lastTime = times[filledPoint]
            endif
            if (times[filledPoint] < firstTime)
                firstTime = times[filledPoint]
            endif                  
            FWHMsYToUse[filledPointToUse]=waveSizeCoefficients[5]*6.7*3*10^-6
            FWHMsXSigmaToUse[filledPointToUse]=W_sigma[3]*6.7*3*10^-6
            FWHMsYSigmaToUse[filledPointToUse]=W_sigma[5]*6.7*3*10^-6
            filledPointToUse = filledPointToUse + 1
        endif
           
                //Determine time
            filledPoint = filledPoint + 1
        //KillWaves temp3
    endfor
   
   
   
    variable timeXIndex = firstTime
    make averagesXTimes
    make averagesXSigmas
    variable averagesXIndex = 0//Index for array of averaged X data
    variable populationXSize
    variable averageXIndex//The index for finding the original data
   
   
    do
        averageXIndex = 0
        populationXSize = 0
            do
                if( (timeXIndex < timesToUse[averageXIndex] + 0.001) && (timeXIndex > timesToUse[averageXIndex] - 0.001) )
                    averagesX[averagesXIndex] = averagesX[averagesXIndex] + (FWHMsXToUse[averageXIndex] * (1/(FWHMsXSigmaToUse[averageXIndex])^2))
                    populationXSize = populationXSize + 1/((FWHMsXSigmaToUse[averageXIndex])^2)
                endif
                averageXIndex = averageXIndex + 1
            while(averageXIndex < filledPointToUse - 0.001)//Now we have sum of members for a particular time
            if(populationXSize > 0.0001)//To handle float rounding errors (Igor doesn't have int variables)
                averagesX[averagesXIndex] = averagesX[averagesXIndex]/populationXSize
                averagesXTimes[averagesXIndex] = timeXIndex//Will place the time into the space with the same index number as the x and y data corresponding to that time.
                averagesXTimes = averagesXTimes - averagesXTimes[0]
                averagesXSigmas[averagesXIndex] = sqrt(1/populationXSize)
                if(averagesXIndex > 4)
                    averagesXSigmas[averagesXIndex] = averagesXSigmas[averagesXIndex]*2
                endif
                averagesXIndex = averagesXIndex + 1
            endif
        timeXIndex = timeXIndex + 1
    while(timeXIndex <= lastTime)//handles out of order times

   
    Make/D/N=2/O W_coefX
    Duplicate/D timesToUse,velocity_fit
    W_coefX = {2,0.000322167}
    FuncFit/H="01"/NTHR=0 cloudExpansion W_coefX  averagesX /X=averagesXTimes /W=averagesXSigmas /I=1
   

   
   
   
    variable velocityX = W_coefX[0] * 1000
    print "velocityX: " + num2str(velocityX)
    variable temperatureX = 2*(85*1.67*10^-27*0.5*velocityX^2)/(1.38*10^-23)//2 * (85*1.67*10^-27*0.5*(W_coefX[0]*1000)^2) / (1.38 * 10^-23)
    print "Temperature X: " + num2str(temperatureX) + " +/- " + num2str(W_sigma[0]) + " Kelvins"
   

    variable timeYIndex = firstTime
    make averagesYTimes
    make averagesYSigmas
    variable averagesYIndex = 0//Index for array of averaged Y data
    variable populationYSize
    variable averageYIndex//The index for finding the original data
   
   
    do
        averageYIndex = 0
        populationYSize = 0
            do
                if( (timeYIndex < timesToUse[averageYIndex] + 0.001) && (timeYIndex > timesToUse[averageYIndex] - 0.001) )
                    averagesY[averagesYIndex] = averagesY[averagesYIndex] + (FWHMsYToUse[averageYIndex] * (1/(FWHMsYSigmaToUse[averageYIndex])^2))
                    populationYSize = populationYSize + 1/((FWHMsYSigmaToUse[averageYIndex])^2)
                endif
                averageYIndex = averageYIndex + 1
            while(averageYIndex < filledPointToUse - 0.001)//Now we have sum of members for a particular time
            if(populationYSize > 0.0001)//To handle float rounding errors (Igor doesn't have int variables)
                averagesY[averagesYIndex] = averagesY[averagesYIndex]/populationYSize
                averagesYTimes[averagesYIndex] = timeYIndex//Will place the time into the space with the same index number as the x and y data corresponding to that time.
                averagesYTimes = averagesYTimes - averagesYTimes[0]
                averagesYSigmas[averagesYIndex] = sqrt(1/populationYSize)
                averagesYIndex = averagesYIndex + 1
            endif
        timeYIndex = timeYIndex + 1
    while(timeYIndex <= lastTime)//handles out of order times
    Make/D/N=2/O W_coefY
    W_coefY = {2,0.000248354}
    //killwaves W_sigma
    FuncFit/NTHR=0 cloudExpansion W_coefY averagesY /Y=averagesYTimes /W=averagesYSigmas /I=1
   


   
    variable velocityY = W_coefX[0] * 1000
    variable temperatureY = 2*(85*1.67*10^-27*0.5*velocityY^2)/(1.38*10^-23)
    print "Temperature Y: " + num2str(temperatureY) + " +/- " + num2str(W_sigma[0]) + " Kelvins"

    print "Average: " + num2str((temperatureY + temperatureX) / 2)
End

Make/D/N=2/O W_coef
W_coef[0] = {50,0.000322167}
FuncFit/NTHR=0 cloudExpansion W_coef  averagesY /X=averagesYTimes /W=averagesXSigmas /I=1 /F={0.950000, 4}


--Now let's all stay positive, and do some science.

January 20, 2014 at 07:59 pm - Permalink

I think it might fix your problem if you remove Make W_sigma (Igor will make that wave whether you want it or not) and ADD WAVE W_sigma AFTER the call to FuncFit.

It is important that the WAVE statement be after FuncFit- it must be after the wave has been created by FuncFit. To learn why, read this in the Igor help:
DisplayHelpTopic "Wave References"
The explanation starts with the sentence "A Wave statement has both a compile-time and a runtime effect."

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com

January 21, 2014 at 09:44 am - Permalink

Are you saying to add "WAVE W_sigma" after each call to FuncFit? I did this, and removed the make W_sigma command at the beginning of the function, and the problem still persists.

I read the topic "Wave References", which doesn't explain why I would need to use "WAVE W_sigma" after calls to FuncFit, since I already know the name of the wave that is going to be created. I can just reference it as W_sigma even if I haven't used "WAVE W_sigma", right?

--Now let's all stay positive, and do some science.

January 21, 2014 at 11:01 am - Permalink

DimensionalDifficulties wrote:
Are you saying to add "WAVE W_sigma" after each call to FuncFit? I did this, and removed the make W_sigma command at the beginning of the function, and the problem still persists.

Yes, that's what I was recommending.
Quote:
I read the topic "Wave References", which doesn't explain why I would need to use "WAVE W_sigma" after calls to FuncFit,

It does, actually.
Quote:
since I already know the name of the wave that is going to be created. I can just reference it as W_sigma even if I haven't used "WAVE W_sigma", right?

No. The WAVE statement has two functions, one at compile time and one at runtime. That's what is described in that help topic.

At compile time, the compiler encounters a name like W_Sigma and it has to know what type of object it is (wave, variable, string) to compile the right kind of code.
At runtime, the WAVE statement actually causes Igor to go look up the wave and store a pointer to the wave in the wave reference. This is why it has to be after FuncFit- the wave doesn't exist (or may be a different wave) before FuncFit runs, so if the WAVE statement is before, then it can't look it up.

Just to confuse the whole issue, there are certain cases where the compiler knows what kind of object is required and will generate the right kind of code even without the WAVE statement. But let's not worry about that detail just yet.

Try to boil down the code to something I can run easily (I don't want to have to load the data set from an external file), make sure the code is copied into the experiment file, and post the file here. I'll try to take a look at it. Guidance is always appreciated in terms of instructions on how to invoke it, initial guesses if it's required, etc.

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com

January 21, 2014 at 12:51 pm - Permalink

Alright, I've found a syntax error in my code that seems to have been preventing the FuncFit command from working properly. Fixing that, and utilizing your advice concerning the make command, I was able to get the function to work properly and update the W_sigma wave as expected. Thank you for your help. Here is my function in a form that works:

#pragma rtGlobals=3     // Use modern global access method and strict wave access.
Function openTemperatureData(pathname) 
    string pathname
    variable binariesIndex = 0
    variable stringIndex
    string currentFilename
   
   
    variable uselessNumber, timeInMs, iterationNumber
    // Get a semicolon-separated list of all files in the folder
    String fileList = IndexedFile($pathName, -1, ".dat")
   

    make/N=500 times//A wave storing the times associated with data points in data waves
    make/N=500 iterations//A wave storing the iteration number associated with data points in data waves

   
    make/N=500 timesToUse//A wave storing the iteration number associated with data points in data waves
    make/N=500 FWHMsXToUse//A wave storing the FWHMs associated with each data wave
    make/N=500 FWHMsYToUse
    make/N=500 FWHMsXSigmaToUse
    make/N=500 FWHMsYSigmaToUse
    Make/D/O/N=7 waveSizeCoefficients
    make averagesX//The average FWHM at each timestep
    make averagesY
    // Sort using combined alpha and numeric sort
    fileList = SortList(fileList, ";", 16)
    variable filledPoint = 0
    variable filledPointToUse = 0
   
    variable cropSizeX
    variable cropSizeY
    Make/N=(200,200)/W cropped
    variable lastTime = 0////The last time that holds data that can be fit
    variable firstTime = 9999 //The first time that holds data that can be fit 
    variable V_FitError//Holds error status of fit routines.  
    make temp
    make temp3
   
   
    for(binariesIndex = 0; binariesIndex < itemsInList(fileList)-1; binariesIndex += 1)//itemsInList(fileList)-1
        currentFilename = StringFromList(binariesIndex,fileList)
        stringIndex=0
        //FindListItem(currentFilename, fileList)
        sscanf currentFilename, "%d_%dms_a_%d", uselessNumber, timeInMs, iterationNumber
       
        GBLoadWave/Q/B=0/A=binary/T={16,16}/W=1/P=$pathName currentFilename
        string wavename = "binary" + num2str(filledPoint)  
        times[filledPoint] = timeInMs
        Redimension/N=(640,480) $wavename
        wave temp=$wavename
        cropped = temp[x+350][y+100]
       
               
        iterations[filledPoint]= iterationNumber
        CurveFit/NTHR=0/Q/N Gauss2D kwCWave=waveSizeCoefficients,  cropped
        WAVE W_sigma
        if(V_FitError > 0.0002)
            print "Error encountered when attempting to fit to " + currentFilename+".  "  
            print "This is usually caused by a picture that does not contain a coherent cloud.  This file will not be used in analysis.  "
            V_FitError = 0
        endif

           
        if((waveSizeCoefficients[3] < 100) && (waveSizeCoefficients[3] > 0.0001))
            FWHMsXToUse[filledPointToUse]=waveSizeCoefficients[3]*6.7*3*10^-6
            timesToUse[filledPointToUse]=times[filledPoint]
            if (times[filledPoint] > lastTime)
                lastTime = times[filledPoint]
            endif
            if (times[filledPoint] < firstTime)
                firstTime = times[filledPoint]
            endif                  
            FWHMsYToUse[filledPointToUse]=waveSizeCoefficients[5]*6.7*3*10^-6
            FWHMsXSigmaToUse[filledPointToUse]=W_sigma[3]*6.7*3*10^-6
            FWHMsYSigmaToUse[filledPointToUse]=W_sigma[5]*6.7*3*10^-6
            filledPointToUse = filledPointToUse + 1
        endif
           
                //Determine time
            filledPoint = filledPoint + 1
        //KillWaves temp3
    endfor
   
   
   
    variable timeXIndex = firstTime
    make averagesXTimes
    make averagesXSigmas
    variable averagesXIndex = 0//Index for array of averaged X data
    variable populationXSize
    variable averageXIndex//The index for finding the original data
   
   
    do
        averageXIndex = 0
        populationXSize = 0
            do
                if( (timeXIndex < timesToUse[averageXIndex] + 0.001) && (timeXIndex > timesToUse[averageXIndex] - 0.001) )
                    averagesX[averagesXIndex] = averagesX[averagesXIndex] + (FWHMsXToUse[averageXIndex] * (1/(FWHMsXSigmaToUse[averageXIndex])^2))
                    populationXSize = populationXSize + 1/((FWHMsXSigmaToUse[averageXIndex])^2)
                endif
                averageXIndex = averageXIndex + 1
            while(averageXIndex < filledPointToUse - 0.001)//Now we have sum of members for a particular time
            if(populationXSize > 0.0001)//To handle float rounding errors (Igor doesn't have int variables)
                averagesX[averagesXIndex] = averagesX[averagesXIndex]/populationXSize
                averagesXTimes[averagesXIndex] = timeXIndex//Will place the time into the space with the same index number as the x and y data corresponding to that time.
                averagesXTimes = averagesXTimes - averagesXTimes[0]
                averagesXSigmas[averagesXIndex] = sqrt(1/populationXSize)
                if(averagesXIndex > 4)
                    averagesXSigmas[averagesXIndex] = averagesXSigmas[averagesXIndex]*2
                endif
                averagesXIndex = averagesXIndex + 1
            endif
        timeXIndex = timeXIndex + 1
    while(timeXIndex <= lastTime)//handles out of order times

   
    Make/D/N=2/O W_coefX
    Duplicate/D timesToUse,velocity_fit
    W_coefX = {2,0.000322167}
   

    KillWaves W_sigma
    FuncFit/Q/H="01"/NTHR=0 cloudExpansion W_coefX  averagesX[0,9] /X=averagesXTimes[0,9] /W=averagesXSigmas[0,9] /I=1
    print "X"
    print W_sigma
    WAVE W_sigma
    print W_sigma
   
    variable velocityX = W_coefX[0] * 1000
    print "velocityX: " + num2str(velocityX)
    variable temperatureX = 2*(85*1.67*10^-27*0.5*velocityX^2)/(1.38*10^-23)//2 * (85*1.67*10^-27*0.5*(W_coefX[0]*1000)^2) / (1.38 * 10^-23)
    print "Temperature X: " + num2str(temperatureX) + " +/- " + num2str(W_sigma[0]) + " Kelvins"
   

    variable timeYIndex = firstTime
    make averagesYTimes
    make averagesYSigmas
    variable averagesYIndex = 0//Index for array of averaged Y data
    variable populationYSize
    variable averageYIndex//The index for finding the original data
   
   
    do
        averageYIndex = 0
        populationYSize = 0
            do
                if( (timeYIndex < timesToUse[averageYIndex] + 0.001) && (timeYIndex > timesToUse[averageYIndex] - 0.001) )
                    averagesY[averagesYIndex] = averagesY[averagesYIndex] + (FWHMsYToUse[averageYIndex] * (1/(FWHMsYSigmaToUse[averageYIndex])^2))
                    populationYSize = populationYSize + 1/((FWHMsYSigmaToUse[averageYIndex])^2)
                endif
                averageYIndex = averageYIndex + 1
            while(averageYIndex < filledPointToUse - 0.001)//Now we have sum of members for a particular time
            if(populationYSize > 0.0001)//To handle float rounding errors (Igor doesn't have int variables)
                averagesY[averagesYIndex] = averagesY[averagesYIndex]/populationYSize
                averagesYTimes[averagesYIndex] = timeYIndex//Will place the time into the space with the same index number as the x and y data corresponding to that time.
                averagesYTimes = averagesYTimes - averagesYTimes[0]
                averagesYSigmas[averagesYIndex] = sqrt(1/populationYSize)
                averagesYIndex = averagesYIndex + 1
            endif
        timeYIndex = timeYIndex + 1
    while(timeYIndex <= lastTime)//handles out of order times
    Make/D/N=2/O W_coefY
   

   
    W_coefY = {2,0.000248354}
    KillWaves W_sigma
    FuncFit/H="01"/Q/NTHR=0 cloudExpansion W_coefY averagesY[0,9] /X=averagesYTimes[0,9] /W=averagesYSigmas[0,9] /I=1


    print "Y"
    print W_sigma
    WAVE W_sigma
    print W_sigma
   

   
    variable velocityY = W_coefX[0] * 1000
    variable temperatureY = 2*(85*1.67*10^-27*0.5*velocityY^2)/(1.38*10^-23)
    print "Temperature Y: " + num2str(temperatureY) + " +/- " + num2str(W_sigma[0]) + " Kelvins"

    print "Average: " + num2str((temperatureY + temperatureX) / 2)
End


--Now let's all stay positive, and do some science.

January 21, 2014 at 04:32 pm - Permalink