Mercurial > repos > marie-tremblay-metatoul > nmr_annotation

comparison nmr_preprocessing/ptw.R @ 2:7304ec2c9ab7 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author marie-tremblay-metatoul
date Mon, 30 Jul 2018 10:33:03 -0400
parents
children
comparison
equal deleted inserted replaced
1:b55559a2854f 2:7304ec2c9ab7
1 ptw <- function (ref, samp, selected.traces, init.coef = c(0, 1, 0),
2 try = FALSE, warp.type = c("individual", "global"), optim.crit = c("WCC",
3 "RMS"), mode = c("forward", "backward"), smooth.param = ifelse(try,
4 0, 1e+05), trwdth = 20, trwdth.res = trwdth, verbose = FALSE,
5 ...)
6 {
7 optim.crit <- match.arg(optim.crit)
8 warp.type <- match.arg(warp.type)
9 mode <- match.arg(mode)
10 if (is.vector(ref))
11 ref <- matrix(ref, nrow = 1)
12 if (is.vector(samp))
13 samp <- matrix(samp, nrow = 1)
14 if (nrow(ref) > 1 && nrow(ref) != nrow(samp))
15 stop("The number of references does not equal the number of samples")
16 if (length(dim(ref)) > 2)
17 stop("Reference cannot be an array")
18 if (length(dim(samp)) > 2)
19 stop("Sample cannot be an array")
20 if (nrow(samp) == 1)
21 warp.type <- "individual"
22 r <- nrow(samp)
23 if (!missing(selected.traces)) {
24 samp <- samp[selected.traces, , drop = FALSE]
25 if (nrow(ref) > 1)
26 ref <- ref[selected.traces, , drop = FALSE]
27 }
28 if (is.vector(init.coef))
29 init.coef <- matrix(init.coef, nrow = 1)
30 if (warp.type == "global") {
31 if (nrow(init.coef) != 1)
32 stop("Only one warping function is allowed with global alignment.")
33 }
34 else {
35 if (nrow(init.coef) != nrow(samp))
36 if (nrow(init.coef) == 1) {
37 init.coef <- matrix(init.coef, byrow = TRUE,
38 nrow = nrow(samp), ncol = length(init.coef))
39 }
40 else {
41 stop("The number of warping functions does not match the number of samples")
42 }
43 }
44 if (warp.type == "individual") {
45 w <- matrix(0, nrow(samp), ncol(ref))
46 a <- matrix(0, nrow(samp), ncol(init.coef))
47 v <- rep(0, nrow(samp))
48 warped.sample <- matrix(NA, nrow = nrow(samp), ncol = ncol(samp))
49 for (i in 1:nrow(samp)) {
50 if (verbose & nrow(samp) > 1)
51 cat(ifelse(nrow(ref) == 1, paste("Warping sample",
52 i, "with the reference \n"), paste("Warping sample",
53 i, "with reference \n", i)))
54 if (nrow(ref) == 1) {
55 rfrnc <- ref
56 }
57 else {
58 rfrnc <- ref[i, , drop = FALSE]
59 }
60 quad.res <- pmwarp(rfrnc, samp[i, , drop = FALSE],
61 optim.crit, init.coef[i, ], try = try, mode = mode,
62 smooth.param = smooth.param, trwdth = trwdth,
63 trwdth.res = trwdth.res, ...)
64 w[i, ] <- quad.res$w
65 a[i, ] <- quad.res$a
66 v[i] <- quad.res$v
67 warped.sample[i, ] <- c(warp.sample(samp[i, , drop = FALSE],
68 w[i, ], mode = mode))
69 }
70 }
71 else {
72 if (nrow(ref) == 1)
73 ref <- matrix(ref, nrow = nrow(samp), ncol = ncol(ref),
74 byrow = TRUE)
75 if (verbose) {
76 if (nrow(ref) == 1) {
77 cat("Simultaneous warping of samples with reference... \n")
78 }
79 else {
80 cat("Simultaneous warping of samples with references... \n")
81 }
82 }
83 quad.res <- pmwarp(ref, samp, optim.crit, c(init.coef),
84 try = try, mode = mode, smooth.param = smooth.param,
85 trwdth = trwdth, trwdth.res = trwdth.res, ...)
86 w <- t(as.matrix(quad.res$w))
87 a <- t(as.matrix(quad.res$a))
88 v <- quad.res$v
89 warped.sample <- t(warp.sample(samp, w, mode))
90 }
91 if (verbose)
92 cat("\nFinished.\n")
93 result <- list(reference = ref, sample = samp, warped.sample = warped.sample,
94 warp.coef = a, warp.fun = w, crit.value = v, optim.crit = optim.crit,
95 mode = mode, warp.type = warp.type)
96 class(result) <- "ptw"
97 result
98 }