Code in R to generate the vector of entropy of a signal

This is a document that contains all the functions involved to compute the vector of entropy of a given signal together with explanations and examples.

Geometric Pattern Transformation

Given a signal \(Y=\{y_1,y_2, \dots ,y_n\}\), the Geometric Pattern Transformation (GPT) \(T_{MS}:\mathbb{R}^3 \to \Delta M \times \Delta S\) is defined as: \[\begin{equation} \label{eq:GPTMS} T_{MS}(y_{n-2},y_{n-1},y_n) = (2y_n-y_{n-1}-y_{n-2}, y_n-y_{n-2}), \end{equation}\] and can be computed using the following code.

# INPUT

# signal: vector with the values of the signal

# OUTPUT

# GPT points of the signal.

GPT <- function(signal){
  
  n <- length(signal) - 2 # amount of points after applying the GPT
  
  points <- matrix(0, nrow=n, ncol=2) # initialization
  
  for (i in 1:n){
    points[i,1] <- 2*signal[i+2]-signal[i+1]-signal[i] # first component in the plane
    points[i,2] <- signal[i+2]-signal[i] # second component in the plane
  }
  
  return(points) 
}

Consider the following signal \(Y\).

# Example

y <- c(1,3,2,5,6,8,2) # signal
GPT(y)
##      [,1] [,2]
## [1,]    0    1
## [2,]    5    2
## [3,]    5    4
## [4,]    5    3
## [5,]  -10   -4
# INPUT

# x: x coordinate of a point
# y: y coordinate of a point

# OUTPUT

# Polar coordinates (r, t) from cartesian coordinates (x, y), where r is the radius and t the angle.

cartesian2polar <- function(x, y){
  
  # Radius
  r <- sqrt(x^2 + y^2)
  
  # Origin
  if (x == 0 & y == 0){
    t <- 0
  }
  
  # Angle of a point over the posive y axis
  if (x == 0 & y > 0){
    t <- pi / 2
  }
  
  # Angle of a point over the negative y axis
  if (x == 0 & y < 0){
    t <- 3 * pi / 2
  }
  
  # Angle of a point over the posive x axis
  if (y == 0 & x > 0){
    t <- 0
  }
  
  # Angle of a point over the negative x axis
  if (y == 0 & x < 0){
    t <- pi
  }
  
  # Angle of a point in the first quadrant
  if (x > 0 & y > 0){
    t <- atan(y / x)
  }
  
  # Angle of a point in the second quadrant
  if (x < 0 & y > 0){
    t <- pi + atan(y / x)
  }
  
  # Point in the third quadrant
  if (x < 0 & y < 0){
    t <- pi + atan(y / x)
  }
  
  # Angle of a point in the fourth quadrant
  if (x > 0 & y < 0){
    t <- 2*pi + atan(y / x)
  }
  
  return(c(r,t))
}
# Example

pol_coord_1 <- cartesian2polar(0, 0)
pol_coord_2 <- cartesian2polar(0, 1)
pol_coord_3 <- cartesian2polar(1, 1)
pol_coord_4 <- cartesian2polar(-2, 0)

cat(paste0("The polar coordinates of the point (0,0) are ", as.data.frame(pol_coord_1), ". \nThe polar coordinates of the point (0,1) are ", as.data.frame(pol_coord_2), ". \nThe polar coordinates of the point (1,1) are ", as.data.frame(pol_coord_3), ". \nThe polar coordinates of the point (-2,0) are ", as.data.frame(pol_coord_4[1]), "."))
## The polar coordinates of the point (0,0) are c(0, 0). 
## The polar coordinates of the point (0,1) are c(1, 1.5707963267949). 
## The polar coordinates of the point (1,1) are c(1.4142135623731, 0.785398163397448). 
## The polar coordinates of the point (-2,0) are 2.
# INPUT

# v: vector in the plane

# OUTPUT

# Norm of the vector v.

vec_norm <- function(v){
  
  return(sqrt(v[1] * v[1] + v[2] * v[2])) 
  
}
# Example

v <- c(0 ,1)
w <- c(1 ,1)

vnorm <- vec_norm(v)
wnorm <- vec_norm(w)

cat(paste0("The norm of the vector (0,1) is ", vnorm, ". \nThe norm of the vector (1,1) is ", wnorm, "."))
## The norm of the vector (0,1) is 1. 
## The norm of the vector (1,1) is 1.4142135623731.
# INPUT

# v: vector in the plane
# w: vector in the plane

# OUTPUT

# Angle between vectors v and w.

angle <- function(v, w){
  
  cdot <- (v[1] * w[1] + v[2] * w[2]) # inner product
  vnorm <- vec_norm(v) # norm of v
  wnorm <- vec_norm(w) # norm of w

  theta <- acos(cdot / (vnorm * wnorm)) # angle bewteen v and w
  
  return(theta)
}
# Example

v <- c(0 ,1)
w <- c(1 ,0)
u <- c(1 ,1)

ang_vw <- angle(v, w)
ang_vu <- angle(v, u)

cat(paste0("The angle between the vectors v and w is ", ang_vw, ". \nThe angle between the vectors v and u is ", ang_vu, "."))
## The angle between the vectors v and w is 1.5707963267949. 
## The angle between the vectors v and u is 0.785398163397448.

We proved that the points of the GPT are located in region of the \(\Delta M \times \Delta S\) plane defined by the ordinal patterns of the signal as it is shown in the two figures below.

# INPUT

# ang: angle from the polar coordinates of a point

# OUTPUT

# Corresponding ordinal pattern: 123, 132, 213, 232, 312 and 321.

# Angles between lines that delimite the regions occupied by each ordinal pattern. 
theta0 <- 0
theta1 <- angle(c(1, 0), c(2, 1))
theta2 <- angle(c(2, 1), c(1, 1))
theta3 <- angle(c(1, 1), c(-1, 0))
theta4 <- angle(c(-1, 0), c(-2, -1))
theta5 <- angle(c(-2, -1), c(-1, -1))
theta6 <- angle(c(-1, -1), c(1, 0))


# Delimiters angles
ang0 <- 0
ang1 <- theta1
ang2 <- ang1 + theta2
ang3 <- ang2 + theta3
ang4 <- ang3 + theta4
ang5 <- ang4 + theta5
ang6 <- ang5 + theta6

# Function to obtain ordinal pattern
pattern <- function(ang){
  
  if (ang0 <= ang & ang < ang1){
    pat <- 213
  }
  
  if (ang1 <= ang & ang < ang2){
    pat <- 123
  }
  
  if (ang2 <= ang & ang < ang3){
    pat <- 132
  }
  
  if (ang3 <= ang & ang < ang4){
    pat <- 231
  }
  
  if (ang4 <= ang & ang < ang5){
    pat <- 321
  }
  
  if (ang5 <= ang & ang < ang6){
    pat <- 312
  }
  
  return(pat)
}
# Example

pat_1 <- pattern(0)
pat_2 <- pattern(pi / 2)
pat_3 <- pattern(3 * pi / 2)

cat(paste0("The angle 0 corresponds to the pattern ", pat_1, ". \nThe angle 90 corresponds to the pattern ", pat_2, ". \nThe angle 270 corresponds to the pattern ", pat_3, "."))
## The angle 0 corresponds to the pattern 213. 
## The angle 90 corresponds to the pattern 132. 
## The angle 270 corresponds to the pattern 312.

Given an interval \(I=[a,b] \subset \mathbb{R}\) and a positive integer \(\mathcal{N}\), the \(\mathcal{N}\)-partition of \(I\) is defined as \[\begin{eqnarray} \label{eq:part} \pi_{\mathcal{N}}(I) = [a, a+\mathcal{L}] \cup \bigcup_{i=0}^{\mathcal{N}-1} (a+i \mathcal{L}, a+(i+1) \mathcal{L}], \end{eqnarray}\] where \(\mathcal{L}=(b-a)/\mathcal{N}\).

If \(D\) is the set of the euclidean norms of the GPT points, we need the \(\mathcal{N}\)-partition of the interval \(I=[\min(D), \max(D)]\) with \[\begin{eqnarray} \label{eq:optN} \mathcal{N}=\max\{\mathcal{M} \in \mathbb{Z}^+ / \forall A \text{ interval in } \pi_{\mathcal{M}}(I): A \cap D \ne \emptyset\}. \end{eqnarray}\]

# INPUT

# set: set to be partitioned
# n_part: number of parts

# OUTPUT

# Position in the partition.

partition <- function(set, n_part){
  
  m <- min(set) # minimum of the set
  M <- max(set) # maximum of the set
  part_length <- (M-m)/n_part # part length
  
  separation <- seq(m, M, part_length) # partition limits
  
  q <- length(separation) # number of parts
  
  part <- vector() # initialization
  
  # Partition assignment
  for (i in 1:length(set)){
    
    if (set[i] == m){
      part[i] <- 1
    } else if (set[i] == M){
      part[i] <- n_part
    } else {
      sorted <- sort(c(set[i], separation))
      if (length(sorted) == q + 1){
        part[i] <- min(which(sorted == set[i])) - 1
      } else {
        part[i] <- which(sep == set[i]) - 1
      }
    }
  }
  
  return(part)
}
# Example

A <- c(1, 2, 1.5, 3, 2.5, 4, 6) # set
sorted_set <- sort(A) # ordered set

two_part <- partition(A, 2) # 2-partition 
five_part <- partition(A, 5) # 5-partition in which the fourth part is empty

cat(paste0("A sorted is ", as.data.frame(sorted_set),  ". \nThe 2-partition of A is ", as.data.frame(two_part), ". \nThe 5-partition of A is ", as.data.frame(five_part), "."))
## A sorted is c(1, 1.5, 2, 2.5, 3, 4, 6). 
## The 2-partition of A is c(1, 1, 1, 1, 1, 2, 2). 
## The 5-partition of A is c(1, 1, 1, 2, 2, 3, 5).
# INPUT

# set: set to be partitioned

# OUTPUT

# Minimum number of partition length such that none part is empty.

partition_length <- function(set){
  
  nset <- length(set) # number of elements in the set
  opt_length <- 1 # initialization

  repeat{
    pp <- partition(set, opt_length)
    if (length(unique(pp)) == opt_length){
      opt_length <- opt_length + 1
    } else {
      break
    } 
  }

  return(opt_length-1)
  
}
# Example

A <- c(1, 2, 1.5, 3, 2.5, 4, 6) # set
optimum_length <- partition_length(A) 
A_partition <- partition(A, optimum_length)

cat(paste0("The optimal number of partitions of A is ", optimum_length,  ". \nThe corresponding partition of A sorted is ", as.data.frame(A_partition), "."))
## The optimal number of partitions of A is 4. 
## The corresponding partition of A sorted is c(1, 1, 1, 2, 2, 3, 4).

The normalized entropy of the pattern \(\rho\) is computed by \[\begin{eqnarray} \label{eq:ent_pat} H_{\rho} = \frac{1}{\ln(\mathcal{N})}\sum_{i=0}^{\mathcal{N}-1} p_i \ln(p_i), \end{eqnarray}\] where \(p_i\) is the probability that a GPT point in the region defined by the pattern \(\rho\) belongs to the \(i\)th element in the partition defined above.

Thus, a vector of entropy for a signal can be obtained by \[\begin{eqnarray} \label{eq:ent_vect} \overrightarrow{H} = (H_{213},H_{123},H_{132},H_{231},H_{321},H_{312}). \end{eqnarray}\]

# INPUT

# signal: vector with the values of the signal

# OUTPUT

# Vector of entropy of the signal.

# Required package
library(plyr)

voe <- function(signal){
  
  # GPT points
  gpt <- GPT(signal)
  
  # Patterns assignment
  gpt_pat <- vector()
  for (i in 1:dim(gpt)[1]){
    polar_coord <- cartesian2polar(gpt[i, 1], gpt[i, 2])
    gpt_pat[i] <- pattern(polar_coord[2])
  }
  
  # Data arrangement
  data_points <- as.data.frame(cbind(gpt, gpt_pat))
  colnames(data_points) <- c("Delta_M", "Delta_S", "Pattern")

  # Pattern splittings
  by_pattern <- split(data_points, data_points$Pattern)
  m <- length(by_pattern)

  # Vector of entropy computation
  vec_entropy <- rep(NA, 6)

  for (i in 1:m){
    
    # Data preparation
    data_points_splittings <- as.data.frame(by_pattern[i])
    colnames(data_points_splittings) <- colnames(data_points)
  
    if (data_points_splittings$Pattern[1] == 213){
      vec_coord <- 1
    } else if (data_points_splittings$Pattern[1] == 123){
      vec_coord <- 2
    } else if (data_points_splittings$Pattern[1] == 132){
      vec_coord <- 3
    } else if (data_points_splittings$Pattern[1] == 231){
      vec_coord <- 4
    } else if (data_points_splittings$Pattern[1] == 321){
      vec_coord <- 5
    } else if (data_points_splittings$Pattern[1] == 312){
      vec_coord <- 6
    }
  
    # Partition of the set of norms
    number_points <- dim(data_points_splittings)[1]
  
    norms <- vector()
    for (j in 1:number_points){
      norms[j] <- vec_norm(data_points_splittings[j, 1:2])
    }
    norms <- unlist(norms)
    
    number_part <- partition_length(norms)
    parts <- partition(norms, number_part)
  
    data_points_splittings <- cbind(data_points_splittings, parts)
  
    # Probability function
    frequencies <- count(data_points_splittings, vars = "parts")["freq"]
    prob <- frequencies / sum(frequencies)
  
    # Entropy
    H <- -sum(prob * log(prob))
    if (H > 0){
      vec_entropy[vec_coord] <- H / log(dim(prob)[1])
    } else {
      vec_entropy[vec_coord] <- 0
    }
  }

  return(vec_entropy)
}
# Example

x <- 1:10
linear <- x
quadratic <- (x-5)^2 
y <- c(1,3,2,5,6,8,2)

voe_1 <- voe(linear)
voe_2 <- voe(quadratic)
voe_3 <- voe(y)

cat(paste0("The vector of entropy of the linear signal is ", as.data.frame(voe_1),  ". \nThe vector of entropy of the quadratic signal is ", as.data.frame(voe_2), ". \nThe vector of entropy of the  signal y is ", as.data.frame(voe_3),"."))
## The vector of entropy of the linear signal is c(NA, 0, NA, NA, NA, NA). 
## The vector of entropy of the quadratic signal is c(0, 1, NA, NA, 1, NA). 
## The vector of entropy of the  signal y is c(0, 1, 0, 0, NA, NA).