Fitting complex functions

Share icon

Sat, 02/11/2012 - 11:35 pm

I have been on a quest to fit complex functions to complex data for quite some time now. While I have had good success with the globalfit dialog included in gencurvefit (a project listed on igorexchange), I wanted some better diagnostic tools, which required abandoning global fit. The idea occurred to me that I could fit my data in the following way, using an all-at-once function and inspired by the global fit method:

Properties that should make this possible:
1. Complex functions take in the same parameters and independent variable for the real and imaginary parts. This simplifies the situation somewhat.
2. The independent variable is the same wave for the real and imaginary parts.

With this in mind, all I need to do is format my input data appropriately and construct my fit function such that

1. It is an all-at-once function. Input waves : parameter wave, frequency wave. Output wave: calculated value wave.
2. The calculated wave is constructed like this: {r0,r1,r2,r3...rn,i0,i1,i2,i3...in} where ri is the real value for point i and ii is the imaginary value for point i.

I have had success with all-at-once functions when I don't have to take apart the input waves or assemble the output wave. Here is an example:

Function ZrealQandles2U(w, realout, frequency) : FitFunc
    Wave w, realout, frequency

    //CurveFitDialog/ Independent Variables 1
    //CurveFitDialog/ frequency
    //CurveFitDialog/ Coefficients 7
    //CurveFitDialog/ w[0] = RP
    //CurveFitDialog/ w[1] = Qdl
    //CurveFitDialog/ w[2] = ndl
    //CurveFitDialog/ w[3] = RS
    //CurveFitDialog/ w[4] = Qbl
    //CurveFitDialog/ w[5] = nbl
    //CurveFitDialog/ w[6] = Rbl   
    Variable RP,Qdl,RS,ndl,Qbl,Rbl,nbl

    RP=w[0]
    Qdl=w[1]
    ndl=w[2]
    RS=w[3]
    Qbl=w[4]
    nbl=w[5]
    Rbl=w[6]
   

    realout = RS+ real(1/(1/RP +Qdl*(2*Pi*10^frequency*sqrt(-1))^ndl)+ 1/(1/Rbl +Qbl*(2*Pi*10^frequency*sqrt(-1))^nbl))
End


here is an example of what I'd like to do, but it doesn't seem to be working:

Function ZrealQandles2U(w, compout, frequency) : FitFunc
    Wave w, compout, frequency
    Wave temp1,temp2,realout,imagout
    variable halfpoint,endpoint
    //CurveFitDialog/ Independent Variables 1
    //CurveFitDialog/ frequency
    //CurveFitDialog/ Coefficients 7
    //CurveFitDialog/ w[0] = RP
    //CurveFitDialog/ w[1] = Qdl
    //CurveFitDialog/ w[2] = ndl
    //CurveFitDialog/ w[3] = RS
    //CurveFitDialog/ w[4] = Qbl
    //CurveFitDialog/ w[5] = nbl
    //CurveFitDialog/ w[6] = Rbl   
    Variable RP,Qdl,RS,ndl,Qbl,Rbl,nbl

    RP=w[0] //assign easy to code variable names
    Qdl=w[1]
    ndl=w[2]
    RS=w[3]
    Qbl=w[4]
    nbl=w[5]
    Rbl=w[6]
   
//calculate some things about the input wave structure
    halfpoint=numpnts(frequency)/2 //this is always an integer, because of the input structure.
    endpoint=numpnts(frequency)

//make two temporary frequency waves       
    duplicate frequency temp1
    duplicate frequency temp2

//make the output wave
    duplicate frequency compout

//trim the temporary frequency waves
    deletepoints halfpoint, halfpoint, temp1
    deletepoints 0, (halfpoint-1), temp2

//compute the real part
    realout = RS+ real(1/(1/RP +Qdl*(2*Pi*10^temp1*sqrt(-1))^ndl)+ 1/(1/Rbl +Qbl*(2*Pi*10^temp1*sqrt(-1))^nbl))
//compute the imagimary part
    imagout = imag(1/(1/RP +Qdl*(2*Pi*10^temp2*sqrt(-1))^ndl)+ 1/(1/Rbl +Qbl*(2*Pi*10^temp2*sqrt(-1))^nbl))
//combine them into the output wave
    concatenate/np/o {realout,imagout},compout
//clean up temporary waves
    killwaves realout, imagout, temp1, temp2
End


Unfortunately, it doesn't work. The fit proceeds as if there are no problems, but the parameters do not vary during the fit so there is obviously something wrong here.

Is there something fundamentally flawed with my method here? Does anyone have any experience fitting complex data without the aid of the global fit panel?

I tried it with a different method,

Function ZrealQandles2Ut(w, compout, frequency) : FitFunc
    Wave w, compout, frequency
    variable halfpoint, temp1, temp2, ii
    //CurveFitDialog/ Independent Variables 1
    //CurveFitDialog/ frequency
    //CurveFitDialog/ Coefficients 7
    //CurveFitDialog/ w[0] = RP
    //CurveFitDialog/ w[1] = Qdl
    //CurveFitDialog/ w[2] = ndl
    //CurveFitDialog/ w[3] = RS
    //CurveFitDialog/ w[4] = Qbl
    //CurveFitDialog/ w[5] = nbl
    //CurveFitDialog/ w[6] = Rbl   
    Variable RP,Qdl,RS,ndl,Qbl,Rbl,nbl

    RP=w[0]
    Qdl=w[1]
    ndl=w[2]
    RS=w[3]
    Qbl=w[4]
    nbl=w[5]
    Rbl=w[6]
   
    halfpoint=numpnts(frequency)/2 //this is always an integer, because of the input structure.
   
    for(ii = 0 ; ii < halfpoint ; ii += 1)
        temp1=frequency[ii]
        compout[ii]=RS+ real(1/(1/RP +Qdl*(2*Pi*10^temp1*sqrt(-1))^ndl)+ 1/(1/Rbl +Qbl*(2*Pi*10^temp1*sqrt(-1))^nbl))
    endfor

    for(ii = 0 ; ii < halfpoint ; ii += 1)
        temp2=frequency[ii+halfpoint]
        compout[ii+halfpoint] = imag(1/(1/RP +Qdl*(2*Pi*10^temp2*sqrt(-1))^ndl)+ 1/(1/Rbl +Qbl*(2*Pi*10^temp2*sqrt(-1))^nbl))  
    endfor
End


But I got the same behavior.

February 13, 2012 at 08:34 am - Permalink

Your second attempt looks reasonable to me. The first is wrong in that you are not allowed to mess with the return wave (the one corresponding to yw in the Igor documentation; compout in your functions).

Your second attempt is basically what you should be doing. You can make it more efficient using wave assignments:
        compout[0, halfpoint-1]=RS+ real(1/(1/RP +Qdl*(2*Pi*10^frequency*sqrt(-1))^ndl)+ 1/(1/Rbl +Qbl*(2*Pi*10^frequency*sqrt(-1))^nbl))
        compout[halfpoint, ] = imag(1/(1/RP +Qdl*(2*Pi*10^ frequency*sqrt(-1))^ndl)+ 1/(1/Rbl +Qbl*(2*Pi*10^ frequency*sqrt(-1))^nbl)) 

But that won't change the result. I find it helpful with things like this to run my fit function outside of FuncFit to make sure that it's producing reasonable output. With an all-at-once function it's a bit tricky: you have to create appropriate wave for compout and frequency in order to run it. But it's easier to debug a fit function on its own than in the context of an actual fit.

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com

February 13, 2012 at 11:14 am - Permalink

Hi John,

Thank you for your assistance. I used the wave assignment method you suggested, which shortened my code considerably. I did some testing and It appears that my function is working as expected outside of any attempts to fit it, so it must be an issue with the function calling the fit. I will look into this and see if I can figure it out.

February 13, 2012 at 01:29 pm - Permalink

Also, in the function you have defined (and John Weeks has cleaned up), my personal taste would be to define and use a complex variable instead of sqrt(-1). For example, in your function define once, and then use
variable/C ki = cmplx(0,1)


This saves frequent recalculation of a complex square-root in several places. It is a pity that Igor's method of Constant definition does not allow complex-number assignments.

February 14, 2012 at 08:47 am - Permalink

OK. All is working well now. Poking around the internet, it seems a lot of people are interested in fitting complex (real+imaginary) data, but are unsure how to do so. Here is my solution:

1. format your data such that the complex part has been seperated into two real waves. concatenate those waves using the /np flag to make a double-length wave containing your data.

2. Do the same for your independent variable, this is usually frequency, making a double-length wave of your frequency spectrum.

3. Format your fit function as below:
Function Zrandles(w,complexcat,logfreq): fitfunc //(takes in w, the wave of coefficients, the output wave (complexcat), and the independent wave (logfreq)).  Don't edit this.
    wave complexcat, w, logfreq //declare wave names.  Don't edit this
    variable halfpoint=numpnts(logfreq)/2 //finding the break between the real and imaginary data.  don't edit this
    variable/c ki = cmplx(0,1) //declare complex constant.  Don't edit this.  Do use "ki" as a replacement for sqrt(-1)
    variable , Rs, Rp, Cdl //declare variable names; you can edit this.  This is optional if you choose to use w[n] instead of easy to read variable names in your function.
   
    Rs=w[0] //assign your variable names to the solution wave.  Edit this.
    Cdl=w[1]
    Rp=w[2]
   
    //put your function in here.  In my case, I've got a simple function, so it just writes complexcat on the next lines.
   
    complexcat[0,halfpoint-1]=Rs+real(1/(1/Rp+1/(ki/-(Cdl*2*Pi*10^logfreq)))) //use real() to get the real part for the first half
    complexcat[halfpoint, ]=imag(1/(1/Rp+1/(-ki/(Cdl*2*Pi*10^logfreq)))) //use imag() to get the imaginary part for the second half
//note that no return command is needed for an all-at-once function.
End


Now use it as you normally would with the global fit dialog, the curvefit dialog, or whatever you use to fit curves these days (I use gencurvefit, a project on igorexchange.)

Thanks to everyone for helping me put this together as I stumbled through the early days of my igor eduction.

I sincerely hope that this saves future researchers from a major headache.

February 16, 2012 at 09:11 am - Permalink

The above works properly, but I realized that I was calculating twice as many times as needed. As such, I decided to use an intermediate complex wave, from which I pull the values to concatenate. All appears to be functioning properly and then at the termination of the fit, I get the following message:

Function Execution Error

While executing a wave read, the following error occured: index out of range.

At andyfaff's suggestion, I tried turning on the procedure debugger, but don't seem to be getting any breaks when running the function, so I started commenting out lines to see where the problem was.

1. the function does not give the error when run independently from gencurvefit
2. if I comment out the line marked below, where the imaginary part is assigned, the error goes away (of course, the curve fit is then bogus)
3. however, if I similarly comment out the line where the real part is assigned, the error stays put.
4. if I replace the last 3 lines with the two commented lines at the end, there is no error.
5. generating the complex wave during the function execution instead of referencing a pre-created wave does not fix the issue.

So I think that I am making a mistake in my wave assignment on that line, specifically "complexer[p-midpointt]", but I can't seem to sort it. I've obtained the same error message with various index shifts and using my best knowledge of igor syntax, it should be correct. The function's code works in its current form, giving a correct result as far as I can see, but kicks out the error message when the fit is complete. All I can figure is something different happens when fitting an allatonce function that I didn't account for. Even more puzzling is the fact that the curve fit gives a perfect result despite the error message. Can anyone shed some light on what may be happening here? Thanks in advance.

Threadsafe Function ZqandlesW(w,complexcat,freqw): fitfunc
    wave complexcat, w, freqw //declare wave names.  Don't edit this
    wave/d/c complexer //pull in temporary complex wave
    variable/g  midpointt //pull in midpoint
    variable/g/c ki, kii //pull in complex constant
    variable Rs, Rp, Qdl,n,sigmaW //declare variable names; you can edit this
   
    Rs=w[0] //assign your variable names to the solution wave.  Edit this.
    Qdl=w[1]
    Rp=w[2]
    n=w[3]
    sigmaW=w[4]

    //computing complex and then assigning
    complexer[0,midpointt-1] = Rs+1/(1/(Rp+sigmaW*kii/sqrt(freqw))+(Qdl*(ki*freqw)^n))
    complexcat[0,midpointt-1]=real(complexer[p])    //if I comment out this line, I still get the error.
    complexcat[midpointt, ]=imag(complexer[p-midpointt]) //if i comment out this line, I no longer get the error.

    //computing seperately
    //complexcat[0,midpointt-1]=real(Rs+1/(1/(Rp+sigmaW*kii/sqrt(freqw))+(Qdl*(ki*freqw)^n)))
    //complexcat[midpointt, ]=imag(Rs+1/(1/(Rp+sigmaW*kii/sqrt(freqw))+(Qdl*(ki*freqw)^n)))
End

April 10, 2012 at 12:31 pm - Permalink

You can't debug a threaded function. For debugging, try removing the Threadsafe keyword and run it on its own.

Be sure that in addition to enabling the debugger, go back to the procedure window right-click menu and enable Debug on Error. It can also help lots of bugs (probably not this one) to enable NVAR, SVAR Wave Checking.

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com

April 10, 2012 at 02:14 pm - Permalink

Ok, with threadsafe off, I am able to get the debugger working and it tells me that the specific wave is "complexer", which is a complex wave. What I realized though, through printing the lengths of various waves (complexcat, complexer and the value for midpoint), I determined that the error is coming up because complexcat is getting re-sized to the value of L (default 200) for the "fit" curve. What happens when you call gencurvefit is the function gets evaluated many times, and then the "fit" values are calculated, based on a wave scaled from the minimum to maximum of the x-values. Because my x-values are actually repeated, instead of always increasing or decreasing ,gencurvefit is having some difficulties. Once the "Fit" values are calculated, the next generation proceeds. I eliminated the error by setting /L=(midpointt*2), such that complexcat's length never changes, and hence, the error message is never triggered.

April 10, 2012 at 04:45 pm - Permalink

mtaylor wrote:
Ok, with threadsafe off, I am able to get the debugger working and it tells me that the specific wave is "complexer", which is a complex wave. What I realized though, through printing the lengths of various waves (complexcat, complexer and the value for midpoint), I determined that the error is coming up because complexcat is getting re-sized to the value of L (default 200) for the "fit" curve. What happens when you call gencurvefit is the function gets evaluated many times, and then the "fit" values are calculated, based on a wave scaled from the minimum to maximum of the x-values. Because my x-values are actually repeated, instead of always increasing or decreasing ,gencurvefit is having some difficulties. Once the "Fit" values are calculated, the next generation proceeds. I eliminated the error by setting /L=(midpointt*2), such that complexcat's length never changes, and hence, the error message is never triggered.


You should never resize/kill/delete/change the type of/etc the waves supplied to the fit function. This will cause an error. The fit waves are not only calculated at the end, they can also be calculated during the fitting process. The fits are (by default) 200 points long and span the range of the input data. The error you were getting is not a gencurvefit problem.

From the help file:

* if /D is not specified, then a fit curve, with the prefix "fit_" to waveName is generated. If the fit is univariate, and the /X flag is not present, then the wave will have destLen points (default 200).
* if the fit is multivariate, or if /X is specified then the fit_wavename wave always has the same number of points as the original data.
*if /D is specified then destWaveName is used as the output for the fit (as described above).

April 10, 2012 at 05:55 pm - Permalink

I think I understand what you are saying, but my function (see above) never changes the length of the wave. It does, however, rely on the length of the wave, which was changed temporarily by gencurvefit to length L, resulting in a non-fatal error message. I experimented with different L values and this was indeed the case as evidenced by having the function print out the length of complexcat every time it was evaluated.

Print "complexcat:"+num2istr(numpnts(complexcat))


during normal function evaluations, the value printed was correct. Periodically (I assume at the end of a generation), it would output as the value of L (default 200). So my code was measuring the length of the wave to be different from what was expected.... and this is what was causing the error message as evidenced by setting L to the same length as complexcat (wy) eliminating the error message, because the wave never changed length. Another way around this would be to calculate the length of the wave within the function, but that wastes a few cpu cycles, so I decided to calculate it beforehand and just reference it.

So yes, you are correct that it is not a bug in gencurvefit, but it is a good thing to know about if you are using a non-standard problem formulation.

April 10, 2012 at 06:37 pm - Permalink

mtaylor wrote:
I think I understand what you are saying, but my function (see above) never changes the length of the wave. It does, however, rely on the length of the wave, which was changed temporarily by gencurvefit to length L, resulting in a non-fatal error message. I experimented with different L values and this was indeed the case as evidenced by having the function print out the length of complexcat every time it was evaluated.

Print "complexcat:"+num2istr(numpnts(complexcat))


during normal function evaluations, the value printed was correct. Periodically (I assume at the end of a generation), it would output as the value of L (default 200). So my code was measuring the length of the wave to be different from what was expected.... and this is what was causing the error message as evidenced by setting L to the same length as complexcat (wy) eliminating the error message, because the wave never changed length. Another way around this would be to calculate the length of the wave within the function, but that wastes a few cpu cycles, so I decided to calculate it beforehand and just reference it.

So yes, you are correct that it is not a bug in gencurvefit, but it is a good thing to know about if you are using a non-standard problem formulation.


The bottom line is that you should not be depending on a given numpnts for the dependent and independent waves, they are variable in size (for gencurvefit and funcfit). You should work out the wave length in the fit function.

April 11, 2012 at 12:10 am - Permalink