-
Notifications
You must be signed in to change notification settings - Fork 9
/
Copy pathcreate.power.basis.R
214 lines (207 loc) · 8.49 KB
/
create.power.basis.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
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
create.power.basis <- function(rangeval = c(0,1), nbasis = NULL,
exponents = NULL, dropind=NULL,
quadvals=NULL, values=NULL,
basisvalues=NULL, names='power',
axes=NULL)
{
# This function creates an power functional data basis
# Arguments
# RANGEVAL ... An array of length 2 containing the lower and upper
# boundaries for the rangeval of argument values
# NBASIS ... The number of basis functions. If this conflicts with
# the length of exponent, the latter is used.
# EXPONENT ... The rate parameters defining x^exponent[i]
# DROPIND ... A vector of integers specifying the basis functions to
# be dropped, if any.
# QUADVALS .. A NQUAD by 2 matrix. The firs t column contains quadrature
# points to be used in a fixed point quadrature. The second
# contains quadrature weights. For example, for (Simpson"s
# rule for (NQUAD = 7, the points are equally spaced and the
# weights are delta.*[1, 4, 2, 4, 2, 4, 1]/3. DELTA is the
# spacing between quadrature points. The default is
# matrix("numeric",0,0).
# VALUES ... A list, with entries containing the values of
# the basis function derivatives starting with 0 and
# going up to the highest derivative needed. The values
# correspond to quadrature points in QUADVALS and it is
# up to the user to decide whether or not to multiply
# the derivative values by the square roots of the
# quadrature weights so as to make numerical integration
# a simple matrix multiplication.
# Values are checked against QUADVALS to ensure the correct
# number of rows, and against NBASIS to ensure the correct
# number of columns.
# The default value of is VALUES is vector("list",0).
# VALUES contains values of basis functions and derivatives at
# quadrature points weighted by square root of quadrature weights.
# These values are only generated as required, and only if slot
# QUADVALS is not matrix("numeric",0,0).
# BASISVALUES ... A vector of lists, allocated by code such as
# vector("list",1).
# This field is designed to avoid evaluation of a
# basis system repeatedly at a set of argument values.
# Each list within the vector corresponds to a specific set
# of argument values, and must have at least two components,
# which may be tagged as you wish.
# The first component in an element of the list vector contains the
# argument values.
# The second component in an element of the list vector
# contains a matrix of values of the basis functions evaluated
# at the arguments in the first component.
# The third and subsequent components, if present, contain
# matrices of values their derivatives up to a maximum
# derivative order.
# Whenever function getbasismatrix is called, it checks
# the first list in each row to see, first, if the number of
# argument values corresponds to the size of the first dimension,
# and if this test succeeds, checks that all of the argument
# values match. This takes time, of course, but is much
# faster than re-evaluation of the basis system. Even this
# time can be avoided by direct retrieval of the desired
# array.
# For example, you might set up a vector of argument values
# called "evalargs" along with a matrix of basis function
# values for these argument values called "basismat".
# You might want too use tags like "args" and "values",
# respectively for these. You would then assign them
# to BASISVALUES with code such as
# basisobj$basisvalues <- vector("list",1)
# basisobj$basisvalues[[1]] <-
# list(args=evalargs, values=basismat)
# Returns
# BASISOBJ ... a functional data basis object of type "power"
# Last modified 9 November 2008 by Spencer Graves
# Previously modified 6 January 2008 by Jim Ramsay
# Default basis for missing arguments
##
## 1. check RANGEVAL
##
if(!is.numeric(rangeval))
stop('rangaval must be numeric; class(rangeval) = ',
class(rangeval) )
if(length(rangeval)<1)
stop('rangeval must be a numeric vector of length 2; ',
'length(rangeval) = 0.')
if (length(rangeval) == 1) {
if( rangeval <= 0)
stop("rangeval a single value that is not positive: ",
rangeval)
rangeval <- c(0,rangeval)
}
if(length(rangeval)>2)
stop('rangeval must be a vector of length 2; ',
'length(rangeval) = ', length(rangeval))
if(diff(rangeval)<=0)
stop('rangeval must cover a positive range; diff(rangeval) = ',
diff(rangeval) )
##
## 2. check 0 < nbasis == length(exponents)
##
{
if(is.null(nbasis)){
if(is.null(exponents)){
nbasis <- 2
exponents <- 0:1
}
else {
if(is.numeric(exponents)){
nbasis <- length(exponents)
if(length(unique(exponents)) != nbasis)
stop('duplicates found in exponents; not allowed.')
}
else
stop('exponents must be numeric; class(exponents) = ',
class(exponents) )
}
}
else {
if(is.numeric(nbasis)){
if(length(nbasis)!=1)
stop('nbasis must be a scalar; length(nbasis) = ',
length(nbasis) )
if((nbasis %%1) != 0)
stop('nbasis just be an integer; nbasis%%1 = ',
nbasis%%1)
{
if(is.null(exponents))
exponents <- 0:(nbasis-1)
else {
if(is.numeric(exponents)){
if(length(exponents) != nbasis)
stop('length(exponents) must = nbasis; ',
'length(exponents) = ', length(exponents),
' != nbasis = ', nbasis)
if(length(unique(exponents)) != nbasis)
stop('duplicates found in exponents; not allowed.')
}
else
stop('exponents must be numeric; class(exponents) = ',
class(exponents) )
}
}
}
else stop('nbasis must be numeric; class(nbasis) = ',
class(nbasis) )
}
}
# check whether exponents are negative,
# and if so whether the range includes nonpostive values
if (any(exponents < 0) && rangeval[1] <= 0)
stop("An exponent is negative and range contains",
" 0 or negative values.")
##
## 3. check DROPIND
##
if (length(dropind) == 0) dropind <- NULL
#
if (length(dropind) > 0){
if(!is.numeric(dropind))
stop('dropind must be numeric; is ', class(dropind))
doops <- which((dropind%%1)>0)
if(length(doops)>0)
stop('dropind must be integer; element ', doops[1],
" = ", dropind[doops[1]], '; fractional part = ',
dropind[doops[1]] %%1)
#
doops0 <- which(dropind<=0)
if(length(doops0)>0)
stop('dropind must be positive integers; element ',
doops0[1], ' = ', dropind[doops0[1]], ' is not.')
doops2 <- which(dropind>nbasis)
if(length(doops2)>0)
stop("dropind must not exceed nbasis = ", nbasis,
'; dropind[', doops2[1], '] = ', dropind[doops2[1]])
#
dropind <- sort(dropind)
if(length(dropind) > 1) {
if(min(diff(dropind)) == 0)
stop("Multiple index values in DROPIND.")
}
}
##
## 4. set up the basis object
##
type <- "power"
params <- sort(as.vector(exponents))
basisobj <- basisfd(type=type, rangeval=rangeval, nbasis=nbasis,
params=params, dropind=dropind, quadvals=quadvals,
values=values, basisvalues=basisvalues)
##
## 5. names
##
{
if(length(names) == nbasis)
basisobj$names <- names
else {
if(length(names)>1)
stop('length(names) = ', length(names), '; must be either ',
'1 or nbasis = ', nbasis)
basisobj$names <- paste(names, 0:(nbasis-1), sep="")
}
}
##
## 6. Done
##
if(!is.null(axes))basisobj$axes <- axes
basisobj
}