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.

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 — Dillion Ecmark
– Dillion Ecmark, Commented Nov 29, 2016 at 8:47

vishes_shell · Accepted Answer · 2016-11-29 09:12:50Z

3

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)))

edited Nov 29, 2016 at 9:12

answered Nov 29, 2016 at 8:45

vishes_shell

23.8k7 gold badges74 silver badges82 bronze badges

Sign up to request clarification or add additional context in comments.

6 Comments

nekomatic Over a year ago

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 Over a year ago

@nekomatic it's good to know.:) And that information is useful for OP.

Dillion Ecmark Over a year ago

@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 Over a year ago

@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 Over a year ago

Also, by coincidence, I read a piece yesterday that cautions against jumping to conclusions about power law behaviour: bactra.org/weblog/491.html

|

Collectives™ on Stack Overflow

Porting a MATLAB code to Python

1 Answer 1

6 Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

6 Comments

Your Answer

Sign up or log in

Post as a guest

Related