Skip to content

Commit

Permalink
Moves dist calculation to fortran code
Browse files Browse the repository at this point in the history
  • Loading branch information
@vwmaus
vwmaus committed Oct 1, 2021
1 parent c517a49 commit 085dee7
Show file tree
Hide file tree
Showing 4 changed files with 197 additions and 13 deletions.
66 changes: 53 additions & 13 deletions R/twdtw_reduce_time.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -82,21 +82,28 @@ twdtwReduceTime = function(x,
py <- py[,names(px),drop=FALSE]

# Compute local cost matrix
cm <- proxy::dist(py, px, method = dist.method)
#cm <- proxy::dist(py, px, method = dist.method)

if(!is.null(weight.fun)){
# Get day of the year for pattern and time series
doyy <- as.numeric(format(ty, "%j"))
doyx <- as.numeric(format(tx, "%j"))

# Compute time-weght matrix
w <- .g(proxy::dist(doyy, doyx, method = dist.method))

# Apply time-weight to local cost matrix
cm <- weight.fun(cm, w)
}
# Get day of the year for pattern and time series
doyy <- as.numeric(format(ty, "%j"))
doyx <- as.numeric(format(tx, "%j"))

# if(!is.null(weight.fun)){
# # Get day of the year for pattern and time series
# doyy <- as.numeric(format(ty, "%j"))
# doyx <- as.numeric(format(tx, "%j"))
#
# # Compute time-weght matrix
# w <- .g(proxy::dist(doyy, doyx, method = dist.method))
#
# # Apply time-weight to local cost matrix
# cm <- weight.fun(cm, w)
# }
# Compute accumulated DTW cost matrix
internals <- .computecost(cm = cm, step.matrix = step.matrix)
#internals <- dtwSat:::.computecost(cm = cm, step.matrix = step.matrix)
xm = na.omit(cbind(doyx, as.matrix(px)))
ym = na.omit(cbind(doyy, as.matrix(py)))
internals = .computecost_fast(xm, ym, step.matrix)

# Find all low cost candidates
a <- internals$startingMatrix[internals$N,1:internals$M]
Expand Down Expand Up @@ -147,3 +154,36 @@ twdtwReduceTime = function(x,
return(out)
}

# @useDynLib dtwSat computecost
.computecost_fast = function(xm, ym, step.matrix){

# cm = rbind(0, cm)
n = nrow(ym)
m = nrow(xm)
d = ncol(ym)

if(is.loaded("computecostfast", PACKAGE = "dtwSat", type = "Fortran")){
out = .Fortran(computecostfast,
XM = matrix(as.double(xm), m, d),
YM = matrix(as.double(ym), n, d),
CM = matrix(as.double(0), n+1, m),
DM = matrix(as.integer(0), n+1, m),
VM = matrix(as.integer(0), n+1, m),
SM = matrix(as.integer(step.matrix), nrow(step.matrix), ncol(step.matrix)),
N = as.integer(n),
M = as.integer(m),
D = as.integer(d),
NS = as.integer(nrow(step.matrix)))
} else {
stop("Fortran computecostfast lib is not loaded")
}
# sqrt(sum((ym[1,-1] - xm[1,-1])*(ym[1,-1] - xm[1,-1])))
res = list()
res$costMatrix = out$CM[-1,]
res$directionMatrix = out$DM[-1,]
res$startingMatrix = out$VM[-1,]
res$stepPattern = step.matrix
res$N = n
res$M = m
res
}
96 changes: 96 additions & 0 deletions src/computecost_fast.f
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,96 @@
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C C
C (c) Victor Maus <[email protected]> C
C Institute for Geoinformatics (IFGI) C
C University of Muenster (WWU), Germany C
C C
C Earth System Science Center (CCST) C
C National Institute for Space Research (INPE), Brazil C
C C
C C
C Efficient computation of DTW cost matrix - 2015-10-16 C
C C
C This function was adpted from the C function 'computeCM' C
C implemented in the R package 'dtw' by Toni Giorgino. C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C XM - matrix with the time series (N,D)
C YM - matrix with the temporal profile (M,D)
C CM - Output cumulative cost matrix
C DM - Direction matrix
C VM - Starting points matrix
C SM - Matrix of step patterns
C N - Number of rows in CM, DM, and VM - time series
C M - Number of columns CM, DM, and VM - temporal profile
C D - Number of spectral dimensions including time in XM and YM
C NS - Number of rows in SM
SUBROUTINE computecostfast(XM, YM, CM, DM, VM, SM, N, M, D, NS)
C I/O Variables
INTEGER N, M, D, NS, SM(NS,4), DM(N+1,M), VM(N+1,M)
DOUBLE PRECISION XM(M,D), YM(N,D), CM(N+1,M)
C Internals
DOUBLE PRECISION W, CP(NS), VMIN
INTEGER I, J, IL(NS), JL(NS), K, PK, KMIN, ZERO, ONE
PARAMETER(ZERO=0,ONE=1)
REAL NAN, INF
NAN = ZERO
NAN = NAN / NAN
INF = HUGE(ZERO)
VM(1,1) = 1
C Initialize the firt row and col of the matrices
DO 21 I = 2, N+1
CM(I,1) = CM(I-1,1) + distance(YM, XM, N, M, D, I-1, 1)
C WRITE (*,*) 'The distance ',I,',1 is ', CM(I,1)
DM(I,1) = 3
VM(I,1) = 1
21 CONTINUE
DO 31 J = 2, M
CM(2,J) = CM(2,J-1) + distance(YM, XM, N, M, D, 1, J)
C WRITE (*,*) 'The distance 2,',J,' is ', CM(2,J)
DM(1,J) = 2
VM(1,J) = J
31 CONTINUE
C Compute cumulative cost matrix
DO 32 J = 2, M
DO 22 I = 2, N+1
C Calculate local distance
CM(I,J) = distance(YM, XM, N, M, D, I-1, J)
C Initialize list of step cost
DO 10 K = 1, NS
CP(K) = NAN
10 CONTINUE
DO 11 K = 1, NS
PK = SM(K,1)
IL(K) = I - SM(K,2)
JL(K) = J - SM(K,3)
IF ((IL(K).GT.ZERO).AND.(JL(K).GT.ZERO)) THEN
W = SM(K,4)
IF (W.EQ.-ONE) THEN
CP(PK) = CM(IL(K),JL(K))
ELSE
CP(PK) = CP(PK) + CM(IL(K),JL(K))*W
ENDIF
ENDIF
11 CONTINUE
KMIN = -ONE
VMIN = INF
ILMIN = -ONE
JLMIN = -ONE
DO 12 K = 1, NS
PK = SM(K,1)
IF (CP(PK).EQ.CP(PK).AND.CP(PK).LT.VMIN) THEN
KMIN = PK
VMIN = CP(PK)
ILMIN = IL(K)
JLMIN = JL(K)
ENDIF
12 CONTINUE
IF (KMIN.GT.-ONE) THEN
CM(I,J) = VMIN
DM(I,J) = KMIN
VM(I,J) = VM(ILMIN, JLMIN)
ENDIF
22 CONTINUE
32 CONTINUE
99 CONTINUE
END
46 changes: 46 additions & 0 deletions src/distance.f
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,46 @@
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C C
C (c) Victor Maus <[email protected]> C
C Institute for Geoinformatics (IFGI) C
C University of Muenster (WWU), Germany C
C C
C Earth System Science Center (CCST) C
C National Institute for Space Research (INPE), Brazil C
C C
C C
C Computation allapsed time - 2016-01-22 C
C C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C XM - matrix with the time series (N,D)
C YM - matrix with the temporal profile (M,D)
C N - Number of rows in CM, DM, and VM - time series
C M - Number of columns CM, DM, and VM - temporal profile
C D - Number of spectral dimensions including time in XM and YM
C I - Single point in the time series to calculate the local distance
C J - Single point in the temporal profile to calculate the local distance
REAL FUNCTION distance(YM, XM, N, M, D, I, J)
INTEGER N, M, D, I, J, K
DOUBLE PRECISION XM(M,D), YM(N,D), PC, TD, BD, CD, HPC
PARAMETER(ZERO=0,ONE=1)
REAL NAN, INF
NAN = ZERO
NAN = NAN / NAN
PC = 366
HPC = PC/2
distance = NAN
C Compute ellapsed time difference
TD = SQRT((YM(I,1) - XM(J,1)) * (YM(I,1) - XM(J,1)))
IF (TD.GT.HPC) THEN
TD = ABS(PC - TD)
ENDIF
CD = ZERO
DO 30 K = 2, D
BD = YM(I,K) - XM(J,K)
CD = CD + (BD * BD)
30 CONTINUE
C WRITE (*,*) 'The value of X J', J, 'is', XM(J,D)
distance = SQRT(CD) + 1.0 / (1.0 + EXP(-0.1 * (TD - 50.0)))
RETURN
END

2 changes: 2 additions & 0 deletions src/init.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,12 +10,14 @@
/* .Fortran calls */
extern void F77_NAME(bestmatches)(void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *);
extern void F77_NAME(computecost)(void *, void *, void *, void *, void *, void *, void *);
extern void F77_NAME(computecostfast)(void *, void *, void *, void *, void *, void *, void *, void *, void *, void *);
extern void F77_NAME(g)(void *, void *, void *, void *);
extern void F77_NAME(tracepath)(void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *);

static const R_FortranMethodDef FortranEntries[] = {
{"bestmatches", (DL_FUNC) &F77_NAME(bestmatches), 11},
{"computecost", (DL_FUNC) &F77_NAME(computecost), 7},
{"computecostfast", (DL_FUNC) &F77_NAME(computecostfast), 10},
{"g", (DL_FUNC) &F77_NAME(g), 4},
{"tracepath", (DL_FUNC) &F77_NAME(tracepath), 11},
{NULL, NULL, 0}
Expand Down

0 comments on commit 085dee7

Please sign in to comment.