0

Here is the link to the mat file temp whose sum i want. The sum for python is 1.1230325644146074e-10. The sum in matlab is 1.2189e-010. But when i copy the array from matlab and add them in python terminal i get the same sum as that of matlab. What is the mystery behind this? 

>>> np.sum(temp)

1.1230325644146074e-10

>>> 0 + 1.34369215036371e-13 + 1.44714547828739e-11 + 3.13077747300113e-11 + 0 +         
    7.33156408714693e-10 + 1.07579271945635e-08 + 8.89446393156299e-09 + 0 + 
    8.00303861023109e-08 + 9.10095259764947e-07 + 6.44662571715407e-07 + 0 + 
    6.16697002187199e-06 + 0.000104686113649727 + 0.000240373037717048 + 
    7.07802623602744e-11 + -0.000240372993389479 + -0.000104686106424492 + 
    -6.16697038783319e-06 + 0 + -6.44662640770614e-07 + -9.10095265552302e-07 + 
    -8.00303919930304e-08 + 0 + -8.89446408625544e-09 + -1.07579272042352e-08 +   
    -7.33156409297323e-10 + 0 + -3.13077747365376e-11 + -1.44714547888711e-11 + 
    -1.34369215245131e-13

1.218862069269912e-10

All this is required because num in the code below is derived from temp and then i have a den or denominator component and i am trying to find the division. So even if these are small values but there exists small errors their division creats errors when i am translating matlab code to python.

Here is the matlab code. 

function mth = compute_mean_angle(dEnergy,th)
global NFFT;
sth         =   sin(2*th);
cth         =   cos(2*th);
num         =   sum(sum(dEnergy.*sth));
den         =   sum(sum(dEnergy.*cth));
mth         =   0.5*atan2(num,den);
if(mth <0)
    mth = mth+pi;
end;

%end function compute_mean_angle

Here is the python code. 

def compute_mean_angle(dEnergy, th):
    global NFFT
    sth = np.sin(np.dot(2, th))
    cth = np.cos(np.dot(2, th))
    num = np.sum(np.sum(dEnergy * sth))
    den = np.sum(np.sum(dEnergy * cth))
    mth = np.dot(0.5, math.atan2(num, den))
    if (mth < 0):
        mth = mth + np.pi
    return mth

I am attaching sample files here. Contains all three files temp.mat, dEnergy.mat and th.mat

Roland Smith
  • 42,427
  • 3
  • 64
  • 94
user38751
  • 1

1 Answers1

1

Can't reproduce:

MATLAB

>> load('temp.mat')
>> whos temp
  Name      Size            Bytes  Class     Attributes

  temp      1x32              256  double              

>> sum(temp)
ans =
   1.2189e-10
>> sprintf('%.16e', ans)
ans =
1.2188620688216985e-10

Python

>>> import numpy as np
>>> import scipy.io

>>> t = scipy.io.loadmat("temp.mat")
>>> who(t)
Name            Shape            Bytes            Type
===========================================================

temp            1 x 32           256              float64

Upper bound on total bytes  =       256

>>> np.sum(t["temp"])
1.2188620688216985e-10

EDIT

In response to comments, again the results are pretty much the same:

MALTAB

>> load('dEnergy.mat')
>> load('th.mat')
>> mth = compute_mean_angle(dEnergy,th);
>> sprintf('%.16e', mth)
ans =
1.5707963267948966e+00

Python

>>> m1 = scipy.io.loadmat("dEnergy.mat")
>>> m2 = scipy.io.loadmat("th.mat")
>>> compute_mean_angle(m1["dEnergy"], m2["th"])
1.5707963267948966
Amro
  • 123,847
  • 25
  • 243
  • 454
  • point is: you can't just copy/paste the output and expect to keep all the precision – Amro Jul 13 '14 at 15:13
  • Based on dEnergy.mat and th.mat, you can see that num and den are different for both python and matlab implementation. – user38751 Jul 13 '14 at 17:04
  • @user38751: no I get the exact same output, see my edit... Let me say this again, do not copy the output of variables from MATLAB and paste it in Python, you will loose precision! Save it as MAT-file from MATLAB, then load it in Python as I've done. – Amro Jul 13 '14 at 17:42
  • Its coming same because num is very small and denominator is large. So num/den = 0 and tan inverse of 0 divided by 2 is pi/2 or 1.5707. Check num and den. They are coming different i have another example here you will find the output also different. (http://speedy.sh/eXSQY/matfiles2.rar) I havent copied the results this time. I did as you told. – user38751 Jul 13 '14 at 17:59
  • ah now I understand what you meant! Well the reason you're seeing different results is because with limited floating-point precision, addition is not really commutative (a+b+c != c+b+a). Just in MATLAB alone, you can get two different result depending on how you compute the sum; for instance `sum(sum(dEnergy.*sth))` vs. `sum(dEnergy(:).*sth(:))` gave me different results even though they are mathematically equivalent... This question might provide more info: http://stackoverflow.com/a/10981385/97160 – Amro Jul 13 '14 at 19:48
  • fyi there are specially crafted summation algorithms that minimize numerical errors found in the naive approach: https://en.wikipedia.org/wiki/Kahan_summation_algorithm. Here is an excellent submission on the File Exchange: http://www.mathworks.com/matlabcentral/fileexchange/26800-xsum – Amro Jul 13 '14 at 19:56