-
Notifications
You must be signed in to change notification settings - Fork 2
/
squareDisp.R
151 lines (115 loc) · 5.29 KB
/
squareDisp.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
## squareDisp-methods
##-----------------------------------------------------------------------------
## squareDisp
##' @name squareDisp
##' @aliases squareDisp squareDispCpp
##' @title squareDisp
##' @rdname squareDisp-methods
##' @docType methods
##' @description
##' calculate square displacement of a track/trajectory as a function of
##' time/step. data.frame has two column, x and y
##' also calculate dx, dy bivariate
##' squareDispCpp is the cpp version of squareDisp
##' @usage
##' squareDisp(track,dt=1,resolution=0.107)
##' squareDispCpp(track,dt=1,resolution=0.107)
##' @param track track dataframe with x and y coordinates.
##' @param dt time step size(in frames).
##' @param resolution resolution value, default is 0.107.
##'
##' @return
##' \itemize{
##' \item{list of square displacements(dx^2 + dy^2) for varying dt
##' values from 1 to dt}
##' }
##'
##' @examples
##' folder1=system.file('extdata','SWR1',package='sojourner')
##' folder2=system.file('extdata','HTZ1',package='sojourner')
##' trackll=compareFolder(folders=c(folder1,folder2), input=3)
##' #use default filter with min=7
##' filtered.trackll=filterTrack(trackll)
##' track.dt=squareDisp(filtered.trackll[[1]][[1]],dt=6)
##' @export squareDisp
## such a basic function that has been used in msd and cdf, need to
## export add help page later TODO
## returns a list of length dt. i.e. when dt=2, it returns a list of 2
## with two item, first corresponding to dt=1, second corresponding to
## dt=2. colnames are x,y, square.disp,dx,dy.
# it should be explicitly say that the dt is not at dt but from 1 to dt
# use to.dt or dt=1:1 or simply only calculate at dt, then use
# lapply(1:7, squareDisp) this way of confusing is because the intial
# msd caculation, that dt actually means 1:dt
# add this TODO, change it to specific dt or change its nomenclature.
# this will affect msd and displacementCDF, change accordingly.
# it is actally at dt
squareDisp = function(track, dt = 1, resolution = 0.107) {
# validity check, stop if dt greater than track length
if (dt > (dim(track)[1] - 1)) {
stop("\ntrack length:\t", dim(track)[1], "\ndt:\t\t", dt,
"\nTime interval (dt) greater than track length-1\n")
}
# store dt-wise sub tracks into a list
track.dt = list()
for (i in seq_len(dt)) {
# dt, time lags; the number of time interval used for measurement.
# default every one step. i, indexing through dt
# define start divide track into dt-wise-subtracks, store them in the
# form of list # at each specified dt, there are N-dt number of
# subtracks
track.dt[[i]] = track[seq(from = i, to = dim(track)[1], by = dt),
]
# define lag compute square.disp in dt-wise-subtracks, stored them in
# original data.frame. note this lag is a dplyr::lag, not base::lag,
# which only works for time series
x.disp = (track.dt[[i]]$x - dplyr::lag(track.dt[[i]]$x, n = 1)) *
resolution
y.disp = (track.dt[[i]]$y - dplyr::lag(track.dt[[i]]$y, n = 1)) *
resolution
# below is an alternative definition view displacement as displacement
# to initial position, how far it has moved, rather than to previous
# position
# x.disp=(track.dt[[i]]$x-track.dt[[i]]$x[1])*resolution
# y.disp=(track.dt[[i]]$y-track.dt[[i]]$y[1])*resolution
square.disp = x.disp^2 + y.disp^2
index = rownames(track.dt[[i]])
# track.dt[[i]]=dplyr::mutate(track.dt[[i]],index,square.disp,
# dx=x.disp,dy=y.disp) for performance purpose
track.dt[[i]] = cbind(track.dt[[i]], index, square.disp, dx = x.disp,
dy = y.disp)
}
return(track.dt)
}
# compile track wise computations to bytecode for performance
# compiler::enableJIT(3)
# squareDisp=compiler::cmpfun(squareDisp,options=list(optimize=3))
# compiler::enableJIT(0)
# no difference when enableJIT
# #----------------------------------------------------------------------------
# # TODO: calculate displacement variance variance for each trajectory
# move small steps, or varies hugely distribution of displacement,
# tells how centralized or spread the trajectory is, is a parameter to
# measure trajectory calculate square displacement for all tracks
# R requirements: Rcpp, RcppArmadillo
##-----------------------------------------------------------------------------
## squareDispCpp
## calculate square displacement of a track/trajectory as a function of
## time/step data.frame has two column, x and y also calculate dx, dy
## bivariate
##' @export squareDispCpp
##' @importFrom Rcpp sourceCpp
squareDispCpp = function(track, dt = 1, resolution = 0.107) {
# convert given track to matrix type
track.out = data.matrix(track)
# Compile source C++ file
file = system.file("cpp", "squareDispRcpp.cpp", package = "sojourner")
cpp.invisibles = Rcpp::sourceCpp(file)
squareDispRcpp = cpp.invisibles[[1]]
# run squareDispRcpp.cpp
track.dt = squareDispRcpp(track.out, dt, resolution)
return(track.dt)
}
# TODO: Need to rewrite make 'file' an option rather than default
# calling in the function, this can leads to unexpected errors when the
# package change name or file changes location.