# Julia implementation of seqperiod (from MATLAB)

Hi! I’ve been playing with Julia for a while. Suppose you are given a 2000-by10000 Matrix of Float64 type and were told to find the periodicity of each of the 10k columns and store it in some variable. Note that, very close values are assumed to be identical. For that, you are provided with a tolerance value. Here periodicity means the minimum length of cycling…

Ex: 1 2 3 1 2 3 1 2 will give periodicity 3
1.000001 1.999999 1 2 0.999999 will give periodicity 2 for tolerance tol = 0.001.

I have written one code:

``````function seqper(x; tol =0.001)
ind1=findall(x-> x<tol,[abs.(x.-x[1])...]);
period=length(x)
for i=2:length(ind1)
if (maximum(abs.(x[ind1[i]:end].-x[1:end-ind1[i]+1]))<=tol)
period=ind1[i]-1
break
end
end
return period
end
``````

Can anyone help me with an performant implementation?

Why I asked: I am finding reasons for migrating from MATLAB to Julia i.e. translating all my MATLAB codes to Julia and check if it gives me more speed. Recently, I am working on a discrete dynamical system (DDS). I need to solve the system for a range of parameter values (say 100) and (trajectory in DynamicalSystems) and store the last few iterations (say 2000) of one state variable in an Array. This gives me a 2000-by-100 array. I need to find the periodicity of these 100 columns.

`vec(abs.(x .- x[1]))` is much better, you don’t want to splat huge amount of elements

`@views maximum(abs, x[ind1[i]:end] .- x[1:end-ind1[i]+1])` should reduce some allocation for you. You can also turn this into a loop to eliminate allocation.

also I don’t quite understand your algorithm, I think it may be not well defined what you want to find:

``````julia> seqper([1,2,1,2,3, 1,2,1,2,3,4])
ind1 = [1, 3, 6, 8]
11

julia> seqper([1,2,1,2, 1,2,1,2,3])
ind1 = [1, 3, 5, 7]
9
``````

I’d expect periodicity to be 5 in first case and 2 or 4 in second

A quick fix for some allocations:

``````function seqper(x; tol =0.001)
m = period = length(x)
ind1 = findall(abs(x[i]-x[1]) < tol for i = 1:m)
for i = 2:length(ind1)
period = ind1[i] - 1
any(abs(x[j]-x[j-ind1[i]+1]) <= tol for j = ind1[i]:m) && return period
end
return period
end
``````
1 Like

In my use-case, I have consider a time-series data and see whether it repeats periodically or not. The last element ‘3’ in your first example breaks that repeating sequence.

why 12123 as a whole can’t be a pattern? does your data guaranteed not to “cross” itself in the middle of a pattern?

As far as I know, yes. I have never seen any irregularity and the numerical solutions are quite smooth.

Jupyter notebook: public/seqper.ipynb at main · genkuroki/public · GitHub

Generate test data `X` and `Y`:

``````using BenchmarkTools

function testvec(n, per; tol=1e-3)
X = [float(mod1(k, per)) for k in 1:n]
@. X + tol * (rand() - 0.5)
end

X = testvec(10^4, 999);
a = zeros(length(X));
a[end] = 2e-3
Y = X + a;
``````

Original version:

``````function seqper(x; tol=0.001)
ind1 = findall(≤(tol), [abs.(x .- x[1])...])
period = length(x)
for i = 2:length(ind1)
if maximum(abs.(x[ind1[i]:end] .- x[1:end-ind1[i]+1])) ≤ tol
period = ind1[i] - 1
break
end
end
return period
end

@show seqper(X)
@btime seqper(\$X)
@show seqper(Y)
@btime seqper(\$Y);
``````

Results:

``````seqper(X) = 999
228.700 μs (10016 allocations: 602.69 KiB)
seqper(Y) = 10000
314.900 μs (10061 allocations: 1.42 MiB)
``````

My revised simple zero-allocated version:

``````function seqper1(x; tol=1e-3)
@inbounds for k in 2:length(x)
if abs(x[k] - x[1]) ≤ tol
all(j -> abs(x[j] - x[j-k+1]) ≤ tol, k:lastindex(x)) && return k - 1
end
end
return length(x)
end

@show seqper1(X)
@btime seqper1(\$X)
@show seqper1(Y)
@btime seqper1(\$Y);
``````

Results:

``````seqper1(X) = 999
6.640 μs (0 allocations: 0 bytes)
seqper1(Y) = 10000
36.400 μs (0 allocations: 0 bytes)
``````

Case `X`: 230 μs → 6.64 μs
Case `Y`: 315 μs → 36.4 μs

In optimization, aiming for zero-allocation is very important in many cases (Measure performance with @time and pay attention to memory allocation).

4 Likes

Thank you @Seif_Shebl . With some minor changes to your code, I got 5x speed.

``````function seqper(x; tol =0.001)
m = period =length(x)
ind1 = findall([abs(x[i]-x[1]) <= tol for i = 1:m]);
for i=2:length(ind1)
period=ind1[i]-1
all([abs(x[j]-x[j-ind1[i]+1])<=tol for j = ind1[i]:m]) && return period
end
return period
end
``````

Replaced ‘any’ with ‘all’ and added a couple of square brackets inside ‘findall’.

Thanks a lot, @genkuroki . By far this one is the fastest. I will check the links you provided. It’s a huge relief. Here is a part of my code :

``````         NTr=10000;
NIter = 2000;
xpts = ytps = 101
par1range=range(1, 9, length = xpts)
par2range = range(0, 1, length = ytps)
sol_last =zeros(Float64, NIter, xpts*ypts);
# println("i = \$(i)")
for  j = 1:ypts
tr = trajectory(vigi_dis(init; b=par1range[i], v=par2range[j]), NIter; Ttr = NTr);
sol_last[:,(i-1)*xpts+j] = tr[2:end,1];
end
end

@time periods = seqper1.(eachcol(sol_last),tol=0.001)

# ss = mxarray(sol_last) # @time ss = seqper.(eachcol(sol_last),tol=0.001)
# xdata, ydata = mxcall(:meshgrid, 2, par1range, par2range)
# mxcall(:seqperiod, 1, ss)

``````

It’s working faster than ever. Can it improved more?

Edit: Sorry, TAB key published my unfinished writing.

you don’t need this.

use `@view tr[2:end,1]`.

depending how heavy this computation is, you can also multi-thread this

@jling , Using @view gives error:

``````                sol_last[:,(i-1)*xpts+j] = @view tr[2:end,1];
``````
``````ERROR: TaskFailedException

nested task error: ArgumentError: Dataset views only accept indices on one dimension
Stacktrace:
[1] view(::Dataset{2, Float64}, ::UnitRange{Int64}, ::Int64)
@ DelayEmbeddings C:\Users\MHBook\.julia\packages\DelayEmbeddings\DdOZb\src\subdataset.jl:40
[2] macro expansion
@ c:\Users\MHBook\Desktop\ljul\isoperiodic2.jl:49 [inlined]
[4] (::var"#445#threadsfor_fun#138"{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}, Matrix{Float64}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}, UnitRange{Int64}})()
Stacktrace:

top-level scope at c:\Users\MHBook\Desktop\ljul\isoperiodic2.jl

A macro to parallelize a for loop to run with multiple threads. Splits the iteration space among multiple tasks and runs those tasks on threads according to a scheduling policy. A barrier is placed at the end of the loop which waits for all tasks to finish execution.

The schedule argument can be used to request a particular scheduling policy. The only currently supported value is :static, which creates one task per thread and divides the iterations equally among them. Specifying :static is an error if used from inside another @threads loop or from a thread other than 1.

The default schedule (used when no schedule argument is present) is subject to change.

!!! compat "Julia 1.5" The schedule argument is available as of Julia 1.5.