Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Update Katz_FD #36

Open
wants to merge 10 commits into
base: master
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from 9 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 5 additions & 1 deletion antropy/entropy.py
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
"""Entropy functions"""

import numpy as np
from numba import jit, types
from math import factorial, log
Expand Down Expand Up @@ -408,7 +409,10 @@ def _app_samp_entropy(x, order, r, metric="chebyshev", approximate=True):
return phi


@jit((types.Array(types.float64, 1, "C", readonly=True), types.int32, types.float64), nopython=True)
@jit(
(types.Array(types.float64, 1, "C", readonly=True), types.int32, types.float64),
nopython=True,
)
def _numba_sampen(sequence, order, r):
"""
Fast evaluation of the sample entropy using Numba.
Expand Down
36 changes: 30 additions & 6 deletions antropy/fractal.py
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
"""Fractal functions"""

import numpy as np
from numba import jit, types
from math import log, floor
Expand Down Expand Up @@ -183,14 +184,37 @@ def katz_fd(x, axis=-1):
1.0000
"""
x = np.asarray(x)
dists = np.abs(np.diff(x, axis=axis))
ll = dists.sum(axis=axis)
ln = np.log10(ll / dists.mean(axis=axis))
aux_d = x - np.take(x, indices=[0], axis=axis)
d = np.max(np.abs(aux_d), axis=axis)
kfd = np.squeeze(ln / (ln + np.log10(d / ll)))

# euclidian distance calculation
euclidean_distance = np.sqrt(1 + np.square(np.diff(x, axis=axis)))

# total and average path lengths
total_path_length = euclidean_distance.sum(axis=axis)
average_path_length = euclidean_distance.mean(axis=axis)

# max distance from first to all
horizontal_diffs = np.arange(1, x.shape[axis])
vertical_diffs = np.take(x, indices=np.arange(1, x.shape[axis]), axis=axis) - np.take(
x, indices=[0], axis=axis
)

if axis == 1: # reshape if needed
horizontal_diffs = horizontal_diffs.reshape(1, -1)
elif axis == 0:
horizontal_diffs = horizontal_diffs.reshape(-1, 1)

# Euclidean distance and max distance
distances = np.sqrt(np.square(horizontal_diffs) + np.square(vertical_diffs))
max_distance = np.max(distances, axis=axis)

# Katz Fractal Dimension Calculation
full_distance = np.log10(total_path_length / average_path_length)
kfd = np.squeeze(full_distance / (full_distance + np.log10(max_distance / total_path_length)))

# ensure scalar output
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi @PiethonProgram,

I think that this implementation can be simplified, for example by following the proposed new implementation in: #34, or by leveraging the Neurokit2 implementation (which as of present gives the same output as Antropy): https://github.com/neuropsychology/NeuroKit/blob/45c9ad90d863ebf4e9d043b975a10d9f8fdeb06b/neurokit2/complexity/fractal_katz.py#L6

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sounds good, I will make the adjustments.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hello, I have taken a look at the implementations you mentioned. Are you sure it can be simplified? Previous implementations that you mentioned are shorter since they are all single-channel.
If you want to only offer single-channel feature extraction then I can make the changes, but otherwise, unless you want to try and decrease time using Numba, I don't think there is much I can "simplify."

Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good point about the support for ND arrays. I have not yet found the time to do a deep dive into this, but can we just take the existing implementation of Antropy (see below) and replace the distance calculation by the Euclidean distance?

dists = np.abs(np.diff(x, axis=axis))
ll = dists.sum(axis=axis)
ln = np.log10(ll / dists.mean(axis=axis))
aux_d = x - np.take(x, indices=[0], axis=axis)
d = np.max(np.abs(aux_d), axis=axis)
kfd = np.squeeze(ln / (ln + np.log10(d / ll)))

or is there more to your implementation that I'm missing?

Thank you

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think the essence of the code is the same, but some additional "bits" are needed when using Euclidean distance in n-dimensions in order to check for distances from one to the other.
If we were only speaking in 1-dimension, then yes, we can simply just replace the distance calculation line.

if not kfd.ndim:
kfd = kfd.item()

return kfd


Expand Down
1 change: 1 addition & 0 deletions antropy/tests/test_entropy.py
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
"""Test entropy functions."""

import unittest
import numpy as np
from numpy.testing import assert_equal
Expand Down
10 changes: 6 additions & 4 deletions antropy/tests/test_fractal.py
Original file line number Diff line number Diff line change
@@ -1,12 +1,13 @@
"""Test fractal dimension functions."""

import unittest
import numpy as np
from numpy.testing import assert_equal
from numpy import apply_along_axis as aal
from antropy import petrosian_fd, katz_fd, higuchi_fd, detrended_fluctuation


from utils import RANDOM_TS, NORMAL_TS, PURE_SINE, ARANGE, TEST_DTYPES
from utils import RANDOM_TS, NORMAL_TS, PURE_SINE, PURE_COSINE, ARANGE, TEST_DTYPES

PPG_SIGNAL = np.array(
[
Expand Down Expand Up @@ -53,10 +54,11 @@ def test_petrosian_fd(self):

def test_katz_fd(self):
x_k = [0.0, 0.0, 2.0, -2.0, 0.0, -1.0, -1.0, 0.0]
self.assertEqual(np.round(katz_fd(x_k), 3), 5.783)
self.assertEqual(np.round(katz_fd(x_k), 3), 1.503)
# 2D data
assert_equal(aal(katz_fd, axis=1, arr=data), katz_fd(data))
assert_equal(aal(katz_fd, axis=0, arr=data), katz_fd(data, axis=0))
data_kfd = np.vstack((RANDOM_TS, NORMAL_TS, PURE_SINE, PURE_COSINE, ARANGE))
assert_equal(aal(katz_fd, axis=1, arr=data_kfd), katz_fd(data_kfd))
assert_equal(aal(katz_fd, axis=0, arr=data_kfd), katz_fd(data_kfd, axis=0))

def test_higuchi_fd(self):
"""Test for function `higuchi_fd`.
Expand Down
1 change: 1 addition & 0 deletions antropy/tests/utils.py
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@
NORMAL_TS = np.random.normal(size=3000)
RANDOM_TS_LONG = np.random.rand(6000)
PURE_SINE = np.sin(2 * np.pi * 1 * np.arange(3000) / 100)
PURE_COSINE = np.cos(2 * np.pi * 1 * np.arange(3000) / 100)
ARANGE = np.arange(3000)

# Data types for which to test input array compatibility
Expand Down
1 change: 1 addition & 0 deletions antropy/utils.py
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
"""Helper functions"""

import numpy as np
from numba import jit
from math import log, floor
Expand Down
Loading