Skip to content

Commit 5984edd

Browse files
JamesRamsay5JamesRamsay5
committed
update
1 parent e529b80 commit 5984edd

File tree

1 file changed

+192
-171
lines changed

1 file changed

+192
-171
lines changed

‎R/landmarkreg.R

+192-171
Original file line numberDiff line numberDiff line change
@@ -1,196 +1,217 @@
1-
landmarkreg <- function(unregfd, ximarks, x0marks=xmeanmarks,
2-
WfdPar=NULL, ylambda=1e-10) {
3-
# Arguments:
4-
# UNREGFD ... functional data object for curves to be registered
5-
# XIMARKS ... N by NL array of times of interior landmarks for
6-
# each observed curve
7-
# XOMARKS ... vector of length NL of times of interior landmarks for
8-
# target curve
9-
# WFDPAR ... a functional parameter object defining a warping function
10-
# YLAMBDA ... smoothing parameter to be used in computing the registered
11-
# functions. For high dimensional bases, local wiggles may be
12-
# found in the registered functions or its derivatives that are
13-
# not seen in the unregistered functions. In this event, this
14-
# parameter should be increased.
15-
# Returns:
16-
# FDREG ... a functional data object for the registered curves
17-
# WARPFD ... a functional data object for the warping functions
18-
# WFD ... a functional data object for the W functions defining the
19-
# warping functions
20-
21-
# Warning: As of March 2022, landmark registration cannot be done using
22-
# function smooth.basis instead of function smooth.morph. The
23-
# warping function must be strictly monotonic, and we have found that using
24-
# smooth.basis too often violates this contraint. Function
25-
# smooth.morph ensures monotonicity.
26-
27-
# Last modified 29 March 2022 by Jim Ramsay
28-
29-
# check unregfd
30-
31-
if (!(inherits(unregfd, "fd"))) stop(
32-
"Argument unregfd not a functional data object.")
33-
34-
basisobj <- unregfd$basis
35-
nbasis <- basisobj$nbasis
36-
rangeval <- basisobj$rangeval
37-
38-
# ---------------------------------------------------------------------
39-
# check ximarks and x0marks
40-
# ---------------------------------------------------------------------
41-
42-
# check ximarks being matrix with ncurve rows and nmarks columns
43-
44-
if (is.numeric(ximarks)) {
45-
nximarks <- length(ximarks)
46-
# if ximarks is a vector, coerce it to a single row matrix
47-
if (is.vector(ximarks)) ximarks <- matrix(ximarks,1,nximarks)
48-
# if ximarks is a data.frame, coerce it to a matrix
49-
if (is.data.frame(ximarks)) ximarks <- as.matrix(ximarks)
1+
landmarkreg <- function(fdobj, ximarks, x0marks=xmeanmarks,
2+
WfdPar=NULL, monwrd=FALSE, ylambda=1e-10)
3+
{
4+
# Arguments:
5+
# FDOBJ ... functional data object for curves to be registered
6+
# XIMARKS ... N by NL array of times of interior landmarks for
7+
# each observed curve
8+
# XOMARKS ... vector of length NL of times of interior landmarks for
9+
# target curve
10+
# WFDPAR ... a functional parameter object defining a warping function
11+
# MONWRD ... If TRUE, warping functions are estimated by monotone smoothing,
12+
# otherwise by regular smoothing. The latter is faster, but
13+
# not guaranteed to produce a strictly monotone warping
14+
# function. If MONWRD is 0 and an error message results
15+
# indicating nonmonotonicity, rerun with MONWRD = 1.
16+
# Default: TRUE
17+
# YLAMBDA ... smoothing parameter to be used in computing the registered
18+
# functions. For high dimensional bases, local wiggles may be
19+
# found in the registered functions or its derivatives that are
20+
# not seen in the unregistered functions. In this event, this
21+
# parameter should be increased.
22+
# Returns:
23+
# FDREG ... a functional data object for the registered curves
24+
# WARPFD ... a functional data object for the warping functions
25+
# WFD ... a functional data object for the W functions defining the
26+
# warping functions
27+
28+
# Last modified 16 November 2021 by Jim Ramsay
29+
30+
# check FDOBJ
31+
32+
if (!(inherits(fdobj, "fd"))) stop(
33+
"Argument fdobj not a functional data object.")
34+
35+
# extract information from curve functional data object and its basis
36+
37+
coef <- fdobj$coefs
38+
coefd <- dim(coef)
39+
ndim <- length(coefd)
40+
ncurve <- coefd[2]
41+
if (ndim > 2) {
42+
nvar <- coefd[3]
5043
} else {
51-
stop("Argument ximarks is not numeric.")
44+
nvar <- 1
5245
}
53-
54-
# compute row-wise mean of ximarks to serve as x0marks if
55-
# needed
56-
46+
47+
basisobj <- fdobj$basis
48+
type <- basisobj$type
49+
nbasis <- basisobj$nbasis
50+
rangeval <- basisobj$rangeval
51+
# fdobj <- fd(matrix(0,nbasis,ncurve),basisobj)
52+
fdParobj <- fdPar(basisobj, 2, ylambda)
53+
54+
# check landmarks
55+
56+
if (is.vector(ximarks) | is.data.frame(ximarks) ) ximarks = as.matrix(ximarks)
57+
ximarksd <- dim(ximarks)
58+
if (ximarksd[1] != ncurve) stop(
59+
"Number of rows of XIMARKS is incorrect.")
60+
if (any(ximarks <= rangeval[1]) || any(ximarks >= rangeval[2])) stop(
61+
"Some landmark values are not within the range.")
62+
nlandm <- dim(ximarks)[2]
5763
xmeanmarks <- apply(ximarks,2,mean)
58-
59-
# check x0marks and coerce it to be a one-row matrix
60-
61-
if (is.numeric(x0marks)) {
62-
nx0marks <- length(x0marks)
63-
if (is.vector(x0marks)) x0marks <- matrix(x0marks,1,nx0marks)
64-
} else {
65-
stop("Argument x0marks is not numeric.")
66-
}
67-
68-
# check that ximarks and x0marks have same number of columns
69-
70-
if (ncol(ximarks) != length(x0marks))
71-
stop("The number of columns in ximarks is not equal to length of x0marks.")
72-
73-
# check that ximarks are within range of unregfd
74-
75-
if (any(ximarks <= rangeval[1]) || any(ximarks >= rangeval[2]))
76-
stop("Argument ximarks has values outside of range of unregfd.")
77-
78-
# check that x0marks are within range of unregfd
79-
80-
if (any(x0marks <= rangeval[1]) || any(x0marks >= rangeval[2]))
81-
stop("Argument x0marks has values outside of range of unregfd.")
82-
83-
# determine the number of curves to be registered
84-
85-
ncurve <- dim(ximarks)[1]
86-
87-
# ---------------------------------------------------------------------
88-
# check WFDPAR
89-
# ---------------------------------------------------------------------
90-
64+
65+
if (length(x0marks) != nlandm) stop(
66+
"Number of target landmarks not equal to number of curve landmarks.")
67+
9168
# set up default WfdPar
92-
69+
9370
if (is.null(WfdPar)) {
94-
Wnbasis <- length(x0marks) + 2
95-
Wbasis <- create.bspline.basis(rangeval, Wnbasis)
96-
Wfd <- fd(matrix(0,Wnbasis,ncurve), Wbasis)
97-
WfdPar <- fdPar(Wfd, 2, 1e-10)
71+
basisobj <- fdobj$basis
72+
rangex <- basisobj$rangeval
73+
wnbasis <- length(x0marks) + 2
74+
wbasis <- create.bspline.basis(rangex, wnbasis)
75+
wfd <- fd(matrix(0,wnbasis,ncurve))
76+
WfdParobj <- fdPar(wfd)
9877
}
99-
78+
79+
# check WFDPAR
80+
10081
WfdPar <- fdParcheck(WfdPar, ncurve)
101-
82+
10283
# set up WFD0 and WBASIS
103-
104-
Wfd0 <- WfdPar$fd
105-
WLfd <- WfdPar$Lfd
106-
Wbasis <- Wfd0$basis
107-
Wrange <- Wbasis$rangeval
108-
Wlambda <- WfdPar$lambda
109-
110-
# ---------------------------------------------------------------------
111-
# set up analysis
112-
# ---------------------------------------------------------------------
113-
114-
nfine <- min(c(101,10*nbasis))
115-
xfine <- seq(rangeval[1],rangeval[2],length=nfine)
116-
yfine <- eval.fd(xfine, unregfd)
117-
yregmat <- yfine
118-
hfunmat <- matrix(0,nfine,ncurve)
119-
hinvmat <- matrix(0,nfine,ncurve)
120-
Wlambda <- max(Wlambda,1e-10)
121-
84+
85+
Wfd0 <- WfdPar$fd
86+
wLfd <- WfdPar$Lfd
87+
wbasis <- Wfd0$basis
88+
89+
# set up WLAMBDA
90+
91+
wlambda <- WfdPar$lambda
92+
93+
# check landmark target values
94+
95+
wrange <- wbasis$rangeval
96+
if (any(rangeval != wrange)) stop(
97+
"Ranges for FD and WFDPAR do not match.")
98+
99+
# set up analysis
100+
101+
n <- min(c(101,10*nbasis))
102+
x <- seq(rangeval[1],rangeval[2],length=n)
103+
104+
y <- eval.fd(x, fdobj)
105+
yregmat <- y
106+
hfunmat <- matrix(0,n,ncurve)
107+
wlambda <- max(wlambda,1e-10)
108+
122109
xval <- c(rangeval[1],x0marks,rangeval[2])
123-
Wnbasis <- Wbasis$nbasis
124-
Wcoef <- matrix(0,Wnbasis,ncurve)
110+
nwbasis <- wbasis$nbasis
111+
Wcoef <- matrix(0,nwbasis,ncurve)
125112
nval <- length(xval)
126-
113+
wval <- rep(1,nval)
114+
127115
# --------------------------------------------------------------------
128116
# Iterate through curves to register
129117
# --------------------------------------------------------------------
130-
131-
if (ncurve > 1) cat("Progress: Each dot is a curve\n")
132-
118+
119+
cat("Progress: Each dot is a curve\n")
120+
133121
for (icurve in 1:ncurve) {
134-
if (ncurve > 1) cat(".")
122+
cat(".")
135123
# set up landmark times for this curve
136124
yval <- c(rangeval[1],ximarks[icurve,],rangeval[2])
137125
# smooth relation between this curve"s values and target"s values
138-
# use monotone smoother
139-
Wfd <- smooth.morph(xval, yval, Wrange, WfdPar)$Wfdobj
140-
hfun <- monfn(xfine, Wfd)
141-
b <- (rangeval[2]-rangeval[1])/(hfun[nfine]-hfun[1])
142-
a <- rangeval[1] - b*hfun[1]
143-
hfun <- a + b*hfun
144-
hfun[c(1,nfine)] <- rangeval
145-
Wcoefi <- Wfd$coef
146-
Wcoef[,icurve] <- Wcoefi
147-
hfunmat[,icurve] <- hfun
148-
126+
if (monwrd) {
127+
# use monotone smoother
128+
Wfds <- smooth.morph(xval, yval, wrange, WfdPar)
129+
Wfd <- Wfds$Wfdobj
130+
h <- monfn(x, Wfd)
131+
b <- (rangeval[2]-rangeval[1])/(h[n]-h[1])
132+
a <- rangeval[1] - b*h[1]
133+
h <- a + b*h
134+
h[c(1,n)] <- rangeval
135+
wcoefi <- Wfd$coef
136+
Wcoef[,icurve] <- wcoefi
137+
} else {
138+
# use unconstrained smoother
139+
warpfd <- smooth.basis(xval, yval, WfdPar, wval)$fd
140+
# set up warping function by evaluating at sampling values
141+
h <- as.vector(eval.fd(x, warpfd))
142+
b <- (rangeval[2]-rangeval[1])/(h[n]-h[1])
143+
a <- rangeval[1] - b*h[1]
144+
h <- a + b*h
145+
h[c(1,n)] <- rangeval
146+
# check for monotonicity
147+
deltah <- diff(h)
148+
if (any(deltah <= 0)) stop(
149+
paste("Non-increasing warping function estimated for curve",icurve,
150+
" Try setting MONWRD to TRUE."))
151+
wcoefi <- warpfd$coef
152+
Wcoef[,icurve] <- wcoefi
153+
}
154+
hfunmat[,icurve] <- h
155+
149156
# compute h-inverse in order to register curves
150-
151-
Wcoefi <- Wfd$coefs
152-
Wfdinv <- fd(-Wcoefi,Wbasis)
153-
WfdParinv <- fdPar(Wfdinv, WLfd, Wlambda)
154-
Wfdinv <- smooth.morph(hfun, xfine, Wrange, WfdParinv)$Wfdobj
155-
hinv <- monfn(xfine, Wfdinv)
156-
b <- (rangeval[2]-rangeval[1])/(hinv[nfine]-hinv[1])
157-
a <- rangeval[1] - b*hinv[1]
158-
hinv <- a + b*hinv
159-
hinv[c(1,nfine)] <- rangeval
160-
157+
158+
if (monwrd) {
159+
wcoef <- Wfd$coefs
160+
Wfdinv <- fd(-wcoef,wbasis)
161+
WfdParinv <- fdPar(Wfdinv, wLfd, wlambda)
162+
Wfdinv <- smooth.morph(h, x, wrange, WfdParinv)$Wfdobj
163+
hinv <- monfn(x, Wfdinv)
164+
b <- (rangeval[2]-rangeval[1])/(hinv[n]-hinv[1])
165+
a <- rangeval[1] - b*hinv[1]
166+
hinv <- a + b*hinv
167+
hinv[c(1,n)] <- rangeval
168+
} else {
169+
hinvfd <- smooth.basis(h, x, WfdPar)$fd
170+
hinv <- as.vector(eval.fd(x, hinvfd))
171+
b <- (rangeval[2]-rangeval[1])/(hinv[n]-hinv[1])
172+
a <- rangeval[1] - b*hinv[1]
173+
hinv <- a + b*hinv
174+
hinv[c(1,n)] <- rangeval
175+
deltahinv <- diff(hinv)
176+
if (any(deltahinv <= 0)) stop(
177+
paste("Non-increasing warping function estimated for curve",icurve))
178+
}
179+
161180
# compute registered curves
162-
163-
yregfd <- smooth.basis(hinv, yfine[,icurve], basisobj)$fd
164-
yregmat[,icurve] <- eval.fd(xfine, yregfd, 0)
181+
182+
if (length(dim(coef)) == 2) {
183+
# single variable case
184+
yregfd <- smooth.basis(hinv, y[,icurve], fdParobj)$fd
185+
yregmat[,icurve] <- eval.fd(x, yregfd, 0)
186+
}
187+
if (length(dim(coef)) == 3) {
188+
# multiple variable case
189+
for (ivar in 1:nvar) {
190+
# evaluate curve as a function of h at sampling points
191+
yregfd <- smooth.basis(hinv, y[,icurve,ivar], fdParobj)$fd
192+
yregmat[,icurve,ivar] <- eval.fd(x, yregfd)
193+
}
194+
}
165195
}
166-
167-
if (ncurve > 1) cat("\n")
168-
196+
197+
cat("\n")
198+
169199
# create functional data objects for the registered curves
170-
171-
regfdPar <- fdPar(basisobj, 2, ylambda)
172-
regfd <- smooth.basis(xfine, yregmat, regfdPar)$fd
173-
regnames <- unregfd$fdnames
200+
201+
fdParobj <- fdPar(basisobj, 2, ylambda)
202+
regfdobj <- smooth.basis(x, yregmat, fdParobj)$fd
203+
regnames <- fdobj$fdnames
174204
names(regnames)[3] <- paste("Registered",names(regnames)[3])
175-
regfd$fdnames <- regnames
176-
205+
regfdobj$fdnames <- regnames
206+
177207
# create functional data objects for the warping functions
178-
179-
warpfd <- smooth.basis(xfine, hfunmat, basisobj)$fd
180-
warpfdnames <- unregfd$fdnames
181-
names(warpfdnames)[3] <- paste("Warped",names(regnames)[1])
182-
warpfd$fdnames <- warpfdnames
183-
184-
# create functional data objects for the inverse warping functions
185-
186-
warpfdinv <- smooth.basis(xfine, hinvmat, basisobj)$fd
187-
warpfdnames <- unregfd$fdnames
208+
209+
warpfdobj <- smooth.basis(x, hfunmat, fdParobj)$fd
210+
warpfdnames <- fdobj$fdnames
188211
names(warpfdnames)[3] <- paste("Warped",names(regnames)[1])
189-
warpfdinv$fdnames <- warpfdnames
190-
191-
# The core function defining the strictly monotone warping
192-
193-
Wfd <- fd(Wcoef, Wbasis)
194-
195-
return( list(regfd=regfd, warpfd=warpfd, warpfdinv=warpfdinv, Wfd = Wfd) )
212+
warpfdobj$fdnames <- warpfdnames
213+
214+
Wfd <- fd(Wcoef, wbasis)
215+
216+
return( list("regfd" = regfdobj, "warpfd" = warpfdobj, "Wfd" = Wfd) )
196217
}

0 commit comments

Comments
 (0)