Matlab to python conversion result is not equal

Question

Since recently I have been working on a script which I have working in Matlab but need in Python. To give a slightly better explanation: I had a script in Matlab that worked, I wrote the same script but then in Python that also worked as should.

I later found that the matlab script had some problems which made it instable: interpolation over degrees. To solve the problem you get when interpolating between 350 and 10 degrees (the discontinuity issue) I used unwrap in Matlab. This worked fine and stabilized the script.

Currently I have a script in Matlab that works and is stable. I also have a script in Python that works, but doesn't do so well when a discontinuity in angles is found. Now I'm trying to rewrite the Matlab to python.

The old (not always stable) python file does work fine in the situation I'm testing and I can compare values. The new one does not work as expected and I can't find why.

The Matlab results differ slightly from the Python results due to different interpolation alogorithms etc. However the two python files I have should give exactly the same answer, which they don't. They are both in the same folder, use the same functions and the same data.

The Matlab code I'm replicating:

load S_phi.mat
ww=load(strcat(filePrefix,'WaveWindows.txt'));
[sww,nww]=locate(s,x,y,ww(:,1),ww(:,2));

angle1=interp1(sww,unwrap(ww(:,3)*pi/180),sS,'linear','extrap')*180/pi;
angle2=interp1(sww,unwrap(ww(:,4)*pi/180),sS,'linear','extrap')*180/pi;

angle1(angle1<0) = angle1(angle1<0)+360;
angle1(angle1>360) = angle1(angle1>360)-360;

angle2(angle2<0) = angle2(angle2<0)+360;
angle2(angle2>360) = angle2(angle2>360)-360;

%% Compute local S-phi curves
for i=1:ns-1;
% for i=81:121;
    for iphi=1:length(phic)
        if(angle1(i) > angle2(i))
            tangle1 = angle1(i);
            tangle2 = angle2(i)+360;
            wavedir2 = wavedir;
            wavedir2(wavedir <= 180) = wavedir(wavedir <= 180) + 360;
            [wavedir3, I] = sort(wavedir2);
            Splus_temp = Splus(iphi,I);
            Smin_temp = Smin(iphi,I);

            Sploc(i,iphi)=integrate(tangle1,tangle2,wavedir3,Splus_temp);
            Smloc(i,iphi)=integrate(tangle1,tangle2,wavedir3,Smin_temp);
        else
            Sploc(i,iphi)=integrate(angle1(i)*pi/180,angle2(i)*pi/180,unwrap(wavedir*pi/180),Splus(iphi,:));
            Smloc(i,iphi)=integrate(angle1(i)*pi/180,angle2(i)*pi/180,unwrap(wavedir*pi/180),Smin(iphi,:));
            Sploc, (angle1(i)*pi/180), (angle2(i)*pi/180), unwrap(wavedir*pi/180)
        end
    end
end

The Python code of the working (but in some cases unstable) script:

#load S_phi.mat files
phic            = np.load('Sphi_phic_plot.npy')
Snet            = np.load('Sphi_Snet.npy')
Splus           = np.load('Sphi_Splus.npy')
Smin            = np.load('Sphi_Smin.npy')
wavedir         = np.load('Sphi_wavedir.npy')


ww              = pd.read_csv(wave_windows_file, delim_whitespace=True)
locateResult    = coastline.locate(s,x,y,ww['x'],ww['y'])
sww             = locateResult['sp']
nww             = locateResult['np']

fan1            = interpolate.InterpolatedUnivariateSpline(sww, ww['angle1'], k=order) # order = 1: linear
angle1          = fan1(sS)

fan2            = interpolate.InterpolatedUnivariateSpline(sww, ww['angle2'], k=order) # order = 1: linear
angle2          = fan2(sS)

## Compute local S-phi curves
Sploc           = np.zeros((ns-1, phic.shape[0]))
Smloc           = np.zeros((ns-1, phic.shape[0]))
for i in range(ns-1):
    for iphi in range(phic.shape[0]):
        Sploc[i, iphi]      = coastline.integrate(angle1[i], angle2[i], wavedir, Splus[iphi,:])
        Smloc[i, iphi]      = coastline.integrate(angle1[i], angle2[i], wavedir, Smin[iphi,:])

The script I wrote that should give the same results as the python file (and emulate the Matlab script):

# load S_phi.mat files
phic            = np.load('Sphi_phic.npy')
Snet            = np.load('Sphi_Snet.npy')
Splus           = np.load('Sphi_Splus.npy')
Smin            = np.load('Sphi_Smin.npy')
wavedir         = np.load('Sphi_wavedir.npy')

ww              = pd.read_csv(wave_windows_file, delim_whitespace=True)
locateResult    = coastline.locate(s,x,y,ww['x'],ww['y'])
sww             = locateResult['sp']
nww             = locateResult['np']

fan1            = interpolate.InterpolatedUnivariateSpline(sww, np.unwrap(ww['angle1']*m.pi/180), k=1) # order = 1: linear
angle1          = fan1(sS) * 180/m.pi

fan2            = interpolate.InterpolatedUnivariateSpline(sww, np.unwrap(ww['angle2']*m.pi/180), k=1) # order = 1: linear
angle2          = fan2(sS) * 180/m.pi

angle1[angle1<0] = angle1[angle1<0] + 360
angle1[angle1>=360] = angle1[angle1>=360] - 360

angle2[angle2<0] = angle2[angle2<0] + 360
angle2[angle2>=360] = angle2[angle2>=360] - 360

## Compute local S-phi curves
Sploc           = np.zeros((ns-1, phic.shape[0]))
Smloc           = np.zeros((ns-1, phic.shape[0]))
for i in range(ns-1):
    for iphi in range(phic.shape[0]):
        if angle1[i] > angle2[i]:
            tangle1 = deepcopy(angle1[i])
            tangle2 = deepcopy(angle2[i])+360

            wavedir2 = deepcopy(wavedir)
            wavedir2[wavedir <= 180] = wavedir[wavedir <= 180] + 360

            I = np.argsort(wavedir2)
            wavedir3 = np.sort(wavedir2)

            Splus_temp = Splus[iphi,I]
            Smin_temp = Smin[iphi,I]

            Sploc[i, iphi]  = coastline.integrate(tangle1, tangle2, wavedir3, Splus_temp)
            Smloc[i, iphi]  = coastline.integrate(tangle1, tangle2, wavedir3, Smin_temp)
        else:
            Sploc[i, iphi] = coastline.integrate(angle1[i], angle2[i], wavedir, Splus[iphi,:])
            Smloc[i, iphi] = coastline.integrate(angle1[i], angle2[i], wavedir, Smin[iphi,:])

The problem is that the results (Sploc and Smloc) from the second python script are not equal to that of the first script. The input to calculate this is however completely the same. I have the idea that values are overwritten inside the loop.

Does anybody have an idea where things are going wrong?

Help is much appreciated.

Answer 1

I finally found the answer. The problem lied in the comparisson of angle1 with angle2. In the file they're equal but due to the unwrap() and the InterpolatedUnivariateSpline() they differ a little bit (270.0 vs 270.00000000000000006) this caused problems.