Porting a MATLAB code to Python

Question

UPDATED TO REFLECT COMMENTS

I am required to port a MATLAB code into python. Unfortunately, despite lot os sleepless nights of debugging and googling, I cannot get my code to run. Here is the MATLAB code in question:

clc;
serie =0:50;
serie = serie - mean(serie);
y = cumsum(serie);
L = length(y);
%Calculate the detrended fluctuation short term coefficient
npuntos = 10;
f1=zeros(1,npuntos);
for n1=4:16
%Segmentation
seg =(1:n1:L);
%disp(length(seg))
yn = zeros(1,L);
 for k = 1:length(seg)-1
 yaux = y(seg(k):seg(k+1)-1);
 x = 1:length(yaux);
 A=[sum(x.^2),sum(x); sum(x),length(x)];
 C=[sum(x.*yaux);sum(yaux)];
 v=inv(A)*C;
 a=v(1); b=v(2);
 pen=a;
 ord=b;
 ytrend = pen*x + ord;
 yn(seg(k):seg(k+1)-1) = ytrend';
 end
f1(n1) = sqrt((1/seg(end)).*sum((y(1:seg(end)-1)-yn(1:seg(end)-1)).^2));
end
n1=4:16;
f1=f1(4:end);
p1 = polyfit(log10(n1(1:end-2)),log10(f1(1:end-2)),1);
alpha1 = p1(1);
disp(alpha1)

My attempt to translate the code into python goes as follows:

import numpy as np
data = np.arange(start=0, stop=51, step=1)
data = data.transpose()
data = data - np.mean(data)
y = np.cumsum(data)
L = len(y)
# Estimate the value of alpha1
npuntos = 12
f1 = [0] * npuntos
for i, n1 in enumerate(np.arange(start=4, stop=16, step=1)):
 seg = np.arange(start=0, stop=L, step=n1) # Potential error
 yn = [0] * L
 for k in np.arange(start=0, stop=len(seg)-1, step=1): # Potential Error
 yaux = y[seg[k]:seg[k + 1]-1] # Potential Error
 x = np.arange(start=1, stop=len(yaux) + 1, step=1)
 A = [[sum(x ** 2), sum(x)], [sum(x), len(x)]]
 C = [[sum(x * yaux)], [sum(yaux)]]
 v = (np.linalg.inv(A)).dot(C)
 pen = v[0]
 ord = v[1]
 ytrend = pen * x + ord
 yn[seg[k]: seg[k + 1] - 1] = ytrend
 f1[i] = np.sqrt((1 / seg[-1]) * sum((y[1:seg[-1] - 1] - yn[1:seg[-1] - 1]) ** 2))
n1 =np.arange(start=4, stop=16, step=1)
f1 = f1[4:]
xx =np.log10(n1[1: - 2])
yy=np.log10(f1[1: - 2])
print(len(xx))
print(len(yy))
p1 = np.polyfit(xx, yy, 1)
alpha1 = p1[1]
print(alpha1)

Unfortunately, I get TypeError when the program is executing this line

p1 = np.polyfit(xx, yy, 1)

This is indeed expected as xx is of length 9 while xx is just 5. By using try/catch block as suggested in the comment,

try:
 f1[n1] = np.sqrt((1 / seg[-1]) * sum((y[1:seg[-1] - 1] - yn[1:seg[-1] - 1]) ** 2))
except IndexError:
 f1.append(np.sqrt((1 / seg[-1]) * sum((y[1:seg[-1] - 1] - yn[1:seg[-1] - 1]) ** 2)))

the error is fixed by the output is completely wrong.

I have gone through the debugger but I cannot quite how the error. Can someone please help? P.S -The above snippet is supposed to calculate the Detrended Fluctuation Analysis (DFA) if anyone is interested.

• Dillion Ecmark: That I have tried but it gives me the same error. Do you mean something like this f1[n1] = np.sqrt((1 / seg[-1]) * sum((y[0:seg[-1] - 1] - yn[0:seg[-1] - 1]) ** 2))? This did not work either. I get the same IndexError: list assignment index out of range
• TheBlackCat: Which version of python are you using?

Answer 1

Thats because you have npuntos = 10 so that f1 = [0] * npuntos makes your f1 list's size equal 10. And then you iterate over

for n1 in np.arange(start=4, stop=16, step=1):

And accessing f1[n1] which from 10 to 15 will give you an IndexError

UPDATE

First of all you can use np.zeros((5,), dtype=np.int) since you already using np module.

Secondly to figure out the thing with IndexError. Personally i don't want to fall into mathematical problem that you are solving, so the solution won't be the best. Just a slight change. I believe that you know that Python indexing is zero based. Since that you will start populating your at 5th element. I'm not sure that you want it. You can do enumerate(np.arange(start=4, stop=16, step=1)) to create zero based indexing to your list:

for i, n1 in enumerate(np.arange(start=4, stop=16, step=1)):
 ...
 f1[i] = np.sqrt((1 / seg[-1]) * sum((y[1:seg[-1] - 1] - yn[1:seg[-1] - 1]) ** 2))

But len(np.arange(start=4, stop=16, step=1)) is 12 not the size you creating your f1(10). So from that point you can create 12 element list.

npuntos = 12
f1 = [0] * npuntos # or np.zeros((5,), dtype=np.int)

Or you can do appending values as it done in MATLAB(as @nekomatic noted) if it's required.

So you need to wrap

f1[n1] = np.sqrt((1 / seg[-1]) * sum((y[1:seg[-1] - 1] - yn[1:seg[-1] - 1]) ** 2))

in try / except:

try:
 f1[n1] = np.sqrt((1 / seg[-1]) * sum((y[1:seg[-1] - 1] - yn[1:seg[-1] - 1]) ** 2))
except IndexError:
 f1.append(np.sqrt((1 / seg[-1]) * sum((y[1:seg[-1] - 1] - yn[1:seg[-1] - 1]) ** 2)))

• nekomatic: Note that the MATLAB code does the same thing, but in MATLAB a matrix or vector automatically expands when you assign to an index outside its current bounds; the equivalent structure in Python (and most other languages) doesn't do this.
• vishes_shell: @nekomatic it's good to know.:) And that information is useful for OP.
• Dillion Ecmark: @nekomatic, how to you suggest I go about solving this? I have tried to create an empty list and appending the values but I also get an error
• nekomatic: @JohnP.Smith if the MATLAB code gives the correct result then your immediate error should be fixed by initialising f1 with 17 elements, not 10 (17 because Python addresses them as 0...16). You should check whether there are any other changes you need to make for zero-based vs one-based indexing though. Also the fact that the MATLAB code contains this error would make me want to check carefully that it was actually implementing the algorithm correctly. I assume 4...16 are the intended values of n (the number of samples in each 'box' as described on your link)?
• nekomatic: Also, by coincidence, I read a piece yesterday that cautions against jumping to conclusions about power law behaviour: bactra.org/weblog/491.html