Mercurial > repos > marie-tremblay-metatoul > nmr_annotation

comparison nmr_preprocessing/Int_Funct.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 beginTreatment <- function(name, Signal_data = NULL, Signal_info = NULL,
2 force.real = FALSE) {
3
4 cat("Begin", name, "\n")
5
6
7 # Formatting the Signal_data and Signal_info -----------------------
8
9 vec <- is.vector(Signal_data)
10 if (vec) {
11 Signal_data <- vec2mat(Signal_data)
12 }
13 if (is.vector(Signal_info)) {
14 Signal_info <- vec2mat(Signal_info)
15 }
16 if (!is.null(Signal_data)) {
17 if (!is.matrix(Signal_data)) {
18 stop("Signal_data is not a matrix.")
19 }
20 if (!is.complex(Signal_data) && !is.numeric(Signal_data)) {
21 stop("Signal_data contains non-numerical values.")
22 }
23 }
24 if (!is.null(Signal_info) && !is.matrix(Signal_info)) {
25 stop("Signal_info is not a matrix.")
26 }
27
28
29 Original_data <- Signal_data
30
31 # Extract the real part of the spectrum ---------------------------
32
33 if (force.real) {
34 if (is.complex(Signal_data)) {
35 Signal_data <- Re(Signal_data)
36 } else {
37 # The signal is numeric Im(Signal_data) is zero anyway so let's avoid
38 # using complex(real=...,imaginary=0) which would give a complex signal
39 # in endTreatment()
40 force.real <- FALSE
41 }
42 }
43
44
45 # Return the formatted data and metadata entries --------------------
46
47 return(list(start = proc.time(), vec = vec, force.real = force.real,
48 Original_data = Original_data, Signal_data = Signal_data, Signal_info = Signal_info))
49 }
50
51 binarySearch <- function(a, target, lower = TRUE) {
52 # search the index i in a such that a[i] == target
53 # if it doesn't exists and lower, it searches the closer a[i] such that a[i] < target
54 # if !lower, it seraches the closer a[i] such that a[i] > target
55 # a should be monotone but can be increasing or decreasing
56
57 # if a is increasing INVARIANT: a[amin] < target < a[amax]
58 N <- length(a)
59 if ((a[N] - target) * (a[N] - a[1]) <= 0) {
60 return(N)
61 }
62 if ((a[1] - target) * (a[N] - a[1]) >= 0) {
63 return(1)
64 }
65 amin <- 1
66 amax <- N
67 while (amin + 1 < amax) {
68 amid <- floor((amin + amax)/2)
69 if ((a[amid] - target) * (a[amax] - a[amid]) < 0) {
70 amin <- amid
71 } else if ((a[amid] - target) * (a[amax] - a[amid]) > 0) {
72 amax <- amid
73 } else {
74 # a[amid] == a[amax] or a[amid] == target In both cases, a[amid] ==
75 # target
76 return(amid)
77 }
78 }
79 if (xor(lower, a[amin] > a[amax])) {
80 # (lower && a[amin] < a[amax]) || (!lower && a[min] > a[max])
81 # If increasing and we want the lower, we take amin
82 # If decreasing and we want the bigger, we take amin too
83 return(amin)
84 } else {
85 return(amax)
86 }
87 }
88
89 checkArg <- function(arg, checks, can.be.null=FALSE) {
90 check.list <- list(bool=c(is.logical, "a boolean"),
91 int =c(function(x){x%%1==0}, "an integer"),
92 num =c(is.numeric, "a numeric"),
93 str =c(is.character, "a string"),
94 pos =c(function(x){x>0}, "positive"),
95 pos0=c(function(x){x>=0}, "positive or zero"),
96 l1 =c(function(x){length(x)==1}, "of length 1")
97 )
98 if (is.null(arg)) {
99 if (!can.be.null) {
100 stop(deparse(substitute(arg)), " is null.")
101 }
102 } else {
103 if (is.matrix(arg)) {
104 stop(deparse(substitute(arg)), " is not scalar.")
105 }
106 for (c in checks) {
107 if (!check.list[[c]][[1]](arg)) {
108 stop(deparse(substitute(arg)), " is not ", check.list[[c]][[2]], ".")
109 }
110 }
111 }
112 }
113
114 endTreatment <- function(name, begin_info, Signal_data) {
115
116 # begin_info: object outputted from beginTreatment
117
118
119 # Formatting the entries and printing process time -----------------------
120 end_time <- proc.time() # record it as soon as possible
121 start_time <- begin_info[["start"]]
122 delta_time <- end_time - start_time
123 delta <- delta_time[]
124 cat("End", name, "\n")
125 cat("It lasted", round(delta["user.self"], 3), "s user time,", round(delta["sys.self"],3),
126 "s system time and", round(delta["elapsed"], 3), "s elapsed time.\n")
127
128
129 if (begin_info[["force.real"]]) {
130 # The imaginary part is left untouched
131 i <- complex(real = 0, imaginary = 1)
132 Signal_data <- Signal_data + i * Im(begin_info[["Original_data"]])
133 }
134
135 if (begin_info[["vec"]]) {
136 Signal_data <- Signal_data[1, ]
137 }
138
139 # Return the formatted data and metadata entries --------------------
140 return(Signal_data)
141 }
142
143 Interpol <- function(t, y) {
144 # y: sample
145 # t : warping function
146
147 m <- length(y)
148 # t <= m-1
149 # because if t > m-1, y[ti+1] will be NA when we compute g
150 valid <- 1 <= t & t <= m-1 # FIXME it was '<' in Bubble v2
151 s <- (1:m)[valid]
152 ti <- floor(t[s])
153 tr <- t[s] - ti
154 g <- y[ti + 1] - y[ti]
155 f <- y[ti] + tr * g
156 list(f=f, s=s, g=g)
157 }
158
159 vec2mat <- function(vec) {
160
161 return(matrix(vec, nrow = 1, dimnames = list(c(1), names(vec))))
162
163 }
164
165 getArg <- function(arg, info, argname, can.be.absent=FALSE) {
166 if (is.null(arg)) {
167 start <- paste("impossible to get argument", argname, "it was not given directly and");
168 if (!is.matrix(info)) {
169 stop(paste(start, "the info matrix was not given"))
170 }
171 if (!(argname %in% colnames(info))) {
172 if (can.be.absent) {
173 return(NULL)
174 } else {
175 stop(paste(start, "is not in the info matrix"))
176 }
177 }
178 if (nrow(info) < 1) {
179 stop(paste(start, "the info matrix has no row"))
180 }
181 arg <- info[1,argname]
182 if (is.na(arg)) {
183 stop(paste(start, "it is NA in the info matrix"))
184 }
185 }
186 return(arg)
187 }
188
189 binarySearch <- function(a, target, lower = TRUE) {
190 # search the index i in a such that a[i] == target
191 # if it doesn't exists and lower, it searches the closer a[i] such that a[i] < target
192 # if !lower, it seraches the closer a[i] such that a[i] > target
193 # a should be monotone but can be increasing or decreasing
194
195 # if a is increasing INVARIANT: a[amin] < target < a[amax]
196 N <- length(a)
197 if ((a[N] - target) * (a[N] - a[1]) <= 0) {
198 return(N)
199 }
200 if ((a[1] - target) * (a[N] - a[1]) >= 0) {
201 return(1)
202 }
203 amin <- 1
204 amax <- N
205 while (amin + 1 < amax) {
206 amid <- floor((amin + amax)/2)
207 if ((a[amid] - target) * (a[amax] - a[amid]) < 0) {
208 amin <- amid
209 } else if ((a[amid] - target) * (a[amax] - a[amid]) > 0) {
210 amax <- amid
211 } else {
212 # a[amid] == a[amax] or a[amid] == target In both cases, a[amid] ==
213 # target
214 return(amid)
215 }
216 }
217 if (xor(lower, a[amin] > a[amax])) {
218 # (lower && a[amin] < a[amax]) || (!lower && a[min] > a[max])
219 # If increasing and we want the lower, we take amin
220 # If decreasing and we want the bigger, we take amin too
221 return(amin)
222 } else {
223 return(amax)
224 }
225 }
226
227 # returns the discrete borders of the interval for a numeric vector a
228
229 indexInterval <- function (a, from, to, inclusive=TRUE) {
230 # If inclusive and from <= to, we need to take the lower
231 # If not inclusive and from > to, we need to take the lower too
232 lowerFrom <- xor(inclusive, from > to)
233 fromIndex <- binarySearch(a, from, lowerFrom)
234 toIndex <- binarySearch(a, to, !lowerFrom)
235 return(fromIndex:toIndex)
236 }