딥러닝이 발달함에 따라 비전과 음성인식 분야에서 좋은 결과를 보이고 있으나, 시계열 분야에서는 RNN과 1D-CNN으로는 만족스러운 결과를 보이지 못하고 있다. 오늘은 시계열 자료를 이미지로 인코딩하는 방법 중 하나인 Gramian Angular Field(GAF)에 대해 알아보고자 한다. GAF는 시계열 데이터를 극좌표계 변환을 통해 이미지로 표현한 것을 의미한다. 또한 가시화된 결과가 우상단에서 좌하단으로 시계열자료를 가시화하므로 시간적 종속성을 유지한다는 특징이 있다.

GAF 분류

GAF는 크게 GASF와 GADF로 구분된다고 한다. pyts패키지의 GADF 함수의 경우 $sin(-\phi_i+\phi_j)$로 계산되므로 역산의 편의성을 위해 아래와 같이 기술하였다.

$$GASF=\begin{pmatrix} cos(\phi_1+\phi_1) & cos(\phi_1+\phi_2) & \cdots & cos(\phi_1+\phi_n) \\ cos(\phi_2+\phi_1) & cos(\phi_2+\phi_2)& \cdots & cos(\phi_2+\phi_n) \\ \vdots & \vdots & \ddots &\vdots \\ cos(\phi_m + \phi_1) & cos(\phi_m+\phi_2) & \cdots & cos(\phi_m+\phi_n) \end{pmatrix} $$

$$GADF=\begin{pmatrix} sin(\phi_1+\phi_1) & sin(\phi_1+\phi_2) & \cdots & sin(\phi_1+\phi_n) \\ cos(\phi_2+\phi_1) & sin(\phi_2+\phi_2)& \cdots & sin(\phi_2+\phi_n) \\ \vdots & \vdots & \ddots &\vdots \\ sin(\phi_m + \phi_1) & sin(\phi_m+\phi_2) & \cdots & sin(\phi_m+\phi_n) \end{pmatrix} $$

GAF를 파이썬으로 구현하면 아래와 같다.

import numpy as np

def cos_sum(a, b):
    """To work with tabulate."""
    return(np.cos(a+b))

def sin_sum(a, b):
    """To work with tabulate."""
    return(np.sin(a+b))
    
x_min=np.min(x)
x_max=np.max(x)
scale=(2*x - x_max- x_min)/(x_max-x_min)
x = np.linspace(0, 1, len(scale))

phi=np.arccos(scale)

gasf=np.vectorize(cos_sum)(
  *np.meshgrid(phi, phi, sparse=True))
gadf=np.vectorize(sin_sum)(
  *np.meshgrid(phi, phi, sparse=True))

python의 구현

python에는 pyts패키지라는 유용한 패키지가 있다. pyts패키지 구현은 아래와 같다.

from pyts.image import GramianAngularField
gasf = GramianAngularField(method='summation')
gadf = GramianAngularField(method='difference')

plt.imshow(gasf.fit_transform(X)[0,:,:])
plt.imshow(gadf.fit_transform(X)[0,:,:])

import numpy as np
import math
import matplotlib.pyplot as plt
time_points = np.linspace(0, 4 * np.pi, 1000)
x = np.sin(time_points)
X = np.array([x])

def cos_sum(a, b):
    """To work with tabulate."""
    return(math.cos(a+b))

def sin_sum(a, b):
    """To work with tabulate."""
    return(math.sin(b-a))
x_min=np.min(x)
x_max=np.max(x)
scale=(2*x - x_max- x_min)/(x_max-x_min)

phi=np.arccos(scale)
r = np.linspace(0, 1, len(scale))

#gasf
plt.imshow(np.vectorize(cos_sum)(
  *np.meshgrid(phi, phi, sparse=True)))

#gadf
plt.imshow(np.vectorize(sin_sum)(
  *np.meshgrid(phi, phi, sparse=True)))

추가 내용은 아래 링크를 참고하기 바란다.

GAF의 역변환

수식을 통한 역변환

GAF는 아쉽게도 역변환이 안된다고 한다. GADF의 주 대각선은 $cos(2\Phi_i)$, GASF의 주 대각선은 $sin(2 \Phi_i)$이 되므로 arctan를 통해 각도를 계산할 수 있다. 다만 pyts패키지를 통해서는 주 대각선이 $sin(0)$으로 0이 되므로, 계산이 불가하였다.

import numpy as np
import math
import matplotlib.pyplot as plt
time_points = np.linspace(0, 4 * np.pi, 1000)
x = np.sin(time_points)
X = np.array([x])

def cos_sum(a, b):
    """To work with tabulate."""
    return(np.cos(a+b))

def sin_sum(a, b):
    """To work with tabulate."""
    return(np.sin(+a+b))
x_min=np.min(x)
x_max=np.max(x)
scale=(2*x - x_max- x_min)/(x_max-x_min)

phi=np.arccos(scale)
r = np.linspace(0, 1, len(scale))

diag_c=np.diag(np.vectorize(cos_sum)(
  *np.meshgrid(phi, phi, sparse=True)))
diag_s=np.diag(np.vectorize(sin_sum)(
  *np.meshgrid(phi, phi, sparse=True)))


angs=np.arctan2(diag_s,diag_c)/2

ang_=list()
for ang in angs:
    if ang>=0:
        ang_.append(ang)
    else:
        ang_.append(ang+np.pi)
plt.plot(phi)
plt.plot(np.arctan2(diag_s,diag_c)/2)
plt.plot(ang_)

패키지를 통한 역변환

역변환 연산이 복잡하여 딥러닝 모델을 활용하여 해결해 보았다. 적어도 기온 자료에 대해서는 적용이 가능해 보였다.

library(tensorflow)
library(reticulate)
library(keras)
library(data.table)
reticulate::py_config()
pyts=import('pyts')
tf$test$is_gpu_available()
# 디비 연결
con <- RMariaDB::dbConnect(RMariaDB::MariaDB(),
    port = 3307,
    host = "ducj3.iptime.org",
    user = "ducj", password = Sys.getenv('db_key'),
	dbname='mysql'
)
# DB 데이터 로드
RMariaDB::dbGetQuery(con,'select * from asos_h limit 1')
load(file='/home/ducj/shiny/shinyapps/asos_flexdashboard/station.rda')
df=RMariaDB::dbGetQuery(con,
  'select tm,ta from asos_h where stnId=100')
df[,'tm']=as.POSIXct(df[,'tm'])
RMariaDB::dbDisconnect(con)

times=df$tm
time=times[1]

# -1~1 사이로 자료 변환
df$scale_ta=(
  2*df$ta-max(df$ta)-min(df$ta))/
    (max(df$ta)-min(df$ta))
gasf=pyts$image$GramianAngularField(
  method='summation',sample_range =NULL)
gadf=pyts$image$GramianAngularField(
  method='difference',sample_range =NULL)

# 학습자료 생성 
## train_x
## (71시간 전부터 현재까지 자료를 GASF, GADF로 변환)
## train_y
## (71시간 전부터 현재까지 -1~1 사이로 스케일 변환한 기온 자료)
gasf_tas=list(); gadf_tas=list();series=list()
for(time in as.character(times[71:length(times)])){
  message(as.character(time))
  time=as.POSIXct(time)
  scale_ta=df[
    which(df$tm%in%(time-(71:0)*3600)),'scale_ta']
  gasf_ta=gasf$fit_transform(list(scale_ta))
  gadf_ta=gadf$fit_transform(list(scale_ta))
  gasf_ta=array(gasf_ta,c(1,72,72,1))
  gadf_ta=array(gadf_ta,c(1,72,72,1))
  scale_ta=array(scale_ta,c(1,72))
  gasf_tas[[as.character(time)]]=gasf_ta
  gadf_tas[[as.character(time)]]=gadf_ta
  series[[as.character(time)]]=scale_ta
  
  rm(list=c('gasf_ta','gadf_ta','scale_ta'))
}
gc()

series=abind::abind(series,along=1)
gasf_ta_df=abind::abind(gasf_tas,along=1)
gadf_ta_df=abind::abind(gadf_tas,along=1)
dfs=abind::abind(gasf_ta_df,gadf_ta_df,along=4)
rm(list=c('gasf','gadf','gasf_ta_df',
  'gadf_ta_df','gadf_tas','gasf_tas'))
gc()

time_idx=times[71:length(times)]

# 학습 모델 생성
height=72
width=72
channels=2

# 3일(72시간) 자료
set_random_seed(42)
model_input=layer_input(shape=c(height,width,channels))
model_output=model_input%>%layer_conv_2d(
  filters = 32, kernel_size = c(3,3),
  activation='relu')%>%
  layer_max_pooling_2d(pool_size=c(2,2))%>%
  layer_conv_2d(filters=64, kernel_size=c(3,3),
  activation='relu')%>%
  layer_max_pooling_2d(pool_size=c(2,2))%>%
  layer_conv_2d(filters=64,kernel_size = c(3,3),
  activation='relu')%>%
  layer_flatten()%>%
  layer_dense(units=512)%>%
  layer_dense(units=72)
model=keras_model(model_input,model_output)

model%>%compile(
  optimizer=optimizer_rmsprop(learning_rate =.001),
  loss=loss_mean_absolute_error()
)

library(lubridate)
min(which(year(times[71:length(times)])==2022)) # 8691
dim(dfs);dim(series)
history=model%>%fit(
  dfs[1:8691,,,1:2,drop=F],
  series[1:8691,],
  epochs=500,
  batch_size=256,
  validation_split=.2
)
pred=model%>%predict(dfs[8692:10682,,,1:2,drop=F])

R에서의 구현

R을 통해 GASF와 GADF를 구현하는 것은 아래와 같다.

library(forecast)
minmax=function(x){
  (2*x-max(x,na.rm=T)-min(x,na.rm=T))/
    (max(x,na.rm=T)-min(x,na.rm=T))
}
phi=acos(
  as.vector(minmax(AirPassengers)))

cos_sum=function(a,b){
  cos(a+b)
}
sin_sum=function(a,b){
  sin(b-a)
}

raster::plot(
  raster::raster(
    outer(phi,phi,FUN=Vectorize(cos_sum))))
raster::plot(
  raster::raster(
    outer(phi,phi,FUN=Vectorize(sin_sum))))

library(reticulate)
pyts=import('pyts')
gasf=pyts$image$GramianAngularField(
  method='summation',sample_range =NULL)
gadf=pyts$image$GramianAngularField(
  method='difference',sample_range =NULL)

gasf_ta=gasf$fit_transform(
  list(as.vector(minmax(AirPassengers))))
gadf_ta=gadf$fit_transform(
  list(as.vector(minmax(AirPassengers))))

raster::plot(
  raster::raster(gasf_ta[1,,]))
raster::plot(
  raster::raster(gadf_ta[1,,]))

# 역변환
library(raster)
gasf=outer(phi,phi,FUN=Vectorize(cos_sum))
plot(raster(gasf))
gadf=outer(phi,phi,FUN=Vectorize(sin_sum))
plot(phi,type='l')
plot(atan2(diag(gadf),diag(gasf)),type='l')
ang=atan2(diag(gadf),diag(gasf))/2
plot(ifelse(ang>=0,ang,ang+pi),type='l')
plot(cos(ifelse(ang>=0,ang,ang+pi)),type='l')