MatrixOp precision?

Share icon

I have found a problem that when using MatrixOp I sometimes get bad values.
When I substitute the equivalent Make command, all is right with the world.
Is this due to the precision of MatrixOp?
Is there some other problem with my code?

(Using 32 bit Igor 6.22A with 64bit Win7)

here is the example code where I find the problem
The MatrixOp command gives some badvals while the commented out Make command generates none:

 
Function test()

Make/D/O/N = 100E6 noise = enoise(0.5)+ 0.5
Extract/INDX/O noise, test0, noise < 0.01
Sort test0, test0 //Make sure the wave is monotonically increasing (should be anyway, but just to make sure I am not dumb)

Make/D/O/N = (numpnts(test0)) test1 = expnoise(100)
Print wavemax(test1) // Make sure that the max val does not exceed 2000

MatrixOp/O results = test0*2000 + test1
//Make/D/O/N=(numpnts(test0)) results = test0*2000 + test1

Differentiate/Meth=2/EP=1 results/D=resultsdiff
Extract/INDX resultsdiff, badvals, resultsdiff < 0

End
The problem in your example is that your test0 is an unsigned integer wave. MatrixOP attempts to maintain unsigned ints as long as possible so the multiplication by 2000 is handled as integer multiplication and this is where the difference between the two approaches comes up. The solution is to redimension test0 as in:
Function test()
 
Make/D/O/N = 100E6 noise = enoise(0.5)+ 0.5
Extract/INDX/O noise, test0, noise < 0.01
Sort test0, test0 //Make sure the wave is monotonically increasing (should be anyway, but just to make sure I am not dumb)
 
Make/D/O/N = (numpnts(test0)) test1 = expnoise(100)
Print wavemax(test1) // Make sure that the max val does not exceed 2000
Redimension/D test0
MatrixOp/O resultsMA = test0*2000 + test1
Make/D/O/N=(numpnts(test0)) resultsOF = test0*2000 + test1

MatrixOP/O aa=sum(abs(resultsMA-resultsOF))
Print aa[0] // check for difference
 
End


I hope this helps,

A.G.
WaveMetrics, Inc.

February 6, 2012 at 12:13 pm - Permalink

I was assuming that MatrixOp would give a double precision wave if any operation involved a double.
This clears things up - thanks.

February 6, 2012 at 02:29 pm - Permalink

BMangum wrote:
I was assuming that MatrixOp would give a double precision wave if any operation involved a double.
This clears things up - thanks.


It is a bit more complicated than that.

MatrixOP needs to support expressions such as:
MatrixOP/O outWave=inImage*factor

Where inImage is /B/U wave and factor is a variable. Now variables in IGOR could be used to represent an integer such as 2 or even a complex number such as cmplx(pi,1) so at runtime the operation looks at the data and decides how to handle this. When it comes to the binary operation inImage*factor it examines factor and finds that it is an integer and also inImage is an integer. The product might fit in an integer wave but to make sure that we do not overflow it is promoted to SP. In your expression there was a following binary operation, i.e., the addition of another wave. At this point, because the wave is DP the whole output is promoted to DP but this is too late as the inaccuracy was already introduced in computing the first product.

FWIW, in IP7 you will be able to explicitly write expressions such as
MatrixOP/O aa=fp64(inImage)

so there would be less ambiguity as to your precise intention.

A.G.
WaveMetrics, Inc.

February 6, 2012 at 03:49 pm - Permalink