-
Notifications
You must be signed in to change notification settings - Fork 76
/
sits_machine_learning.R
495 lines (484 loc) · 21.1 KB
/
sits_machine_learning.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
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
#' @title Train random forest models
#' @name sits_rfor
#'
#' @author Alexandre Ywata de Carvalho, \email{alexandre.ywata@@ipea.gov.br}
#' @author Rolf Simoes, \email{rolf.simoes@@inpe.br}
#' @author Gilberto Camara, \email{gilberto.camara@@inpe.br}
#'
#' @description Use Random Forest algorithm to classify samples.
#' This function is a front-end to the \code{randomForest} package.
#' Please refer to the documentation in that package for more details.
#'
#' @param samples Time series with the training samples
#' (tibble of class "sits").
#' @param num_trees Number of trees to grow. This should not be set to too
#' small a number, to ensure that every input
#' row gets predicted at least a few times (default: 100)
#' (integer, min = 50, max = 150).
#' @param mtry Number of variables randomly sampled as candidates at
#' each split (default: NULL - use default value of
#' \code{randomForest::randomForest()} function, i.e.
#' \code{floor(sqrt(features))}).
#' @param ... Other parameters to be passed
#' to `randomForest::randomForest` function.
#' @return Model fitted to input data
#' (to be passed to \code{\link[sits]{sits_classify}}).
#'
#' @examples
#' if (sits_run_examples()) {
#' # Example of training a model for time series classification
#' # Retrieve the samples for Mato Grosso
#' # train a random forest model
#' rf_model <- sits_train(samples_modis_ndvi,
#' ml_method = sits_rfor
#' )
#' # classify the point
#' point_ndvi <- sits_select(point_mt_6bands, bands = "NDVI")
#' # classify the point
#' point_class <- sits_classify(
#' data = point_ndvi, ml_model = rf_model
#' )
#' plot(point_class)
#' }
#' @export
#'
sits_rfor <- function(samples = NULL, num_trees = 100, mtry = NULL, ...) {
.check_set_caller("sits_rfor")
# Function that trains a random forest model
train_fun <- function(samples) {
# Verifies if 'randomForest' package is installed
.check_require_packages("randomForest")
# Checks 'num_trees'
.check_int_parameter(num_trees, min = 20)
# Get labels (used later to ensure column order in result matrix)
labels <- .samples_labels(samples)
# Get predictors features
train_samples <- .predictors(samples)
# Post condition: is predictor data valid?
.check_predictors(pred = train_samples, samples = samples)
# determine number of random forest
n_features <- ncol(train_samples) - 2
# Apply the 'mtry' default value of 'randomForest' package
if (.has(mtry)) {
# Checks 'mtry'
.check_int_parameter(mtry, min = 1, max = n_features)
} else {
# set the default values of `mtry`
mtry <- floor(sqrt(n_features))
}
# Train a random forest model
model <- randomForest::randomForest(
x = .pred_features(train_samples),
y = as.factor(.pred_references(train_samples)),
samples = NULL, ntree = num_trees, mtry = mtry,
nodesize = 1, localImp = TRUE, norm.votes = FALSE, ...,
na.action = stats::na.fail
)
# Function that predicts results
predict_fun <- function(values) {
# Verifies if randomForest package is installed
.check_require_packages("randomForest")
# Used to check values (below)
input_pixels <- nrow(values)
# Do classification
values <- stats::predict(
object = model, newdata = values, type = "prob"
)
# Are the results consistent with the data input?
.check_processed_values(values, input_pixels)
# Reorder matrix columns if needed
if (any(labels != colnames(values))) {
values <- values[, labels]
}
return(values)
}
# Set model class
predict_fun <- .set_class(
predict_fun, "rfor_model", "sits_model", class(predict_fun)
)
return(predict_fun)
}
# If samples is informed, train a model and return a predict function
# Otherwise give back a train function to train model further
result <- .factory_function(samples, train_fun)
return(result)
}
#' @title Train support vector machine models
#' @name sits_svm
#' @author Alexandre Ywata de Carvalho, \email{alexandre.ywata@@ipea.gov.br}
#' @author Rolf Simoes, \email{rolf.simoes@@inpe.br}
#' @author Gilberto Camara, \email{gilberto.camara@@inpe.br}
#'
#' @description This function receives a tibble with a set of attributes X
#' for each observation Y. These attributes are the values of the time series
#' for each band.
#' The SVM algorithm is used for multiclass-classification.
#' For this purpose, it uses the "one-against-one" approach,
#' in which k(k-1)/2 binary classifiers are trained;
#' the appropriate class is found by a voting scheme.
#' This function is a front-end to the "svm" method in the "e1071" package.
#' Please refer to the documentation in that package for more details.
#'
#' @param samples Time series with the training samples.
#' @param formula Symbolic description of the model to be fit.
#' (default: sits_formula_linear).
#' @param scale Logical vector indicating the variables to be scaled.
#' @param cachesize Cache memory in MB (default = 1000).
#' @param kernel Kernel used in training and predicting.
#' options: "linear", "polynomial", "radial", "sigmoid"
#' (default: "radial").
#' @param degree Exponential of polynomial type kernel (default: 3).
#' @param coef0 Parameter needed for kernels of type polynomial
#' and sigmoid (default: 0).
#' @param cost Cost of constraints violation (default: 10).
#' @param tolerance Tolerance of termination criterion (default: 0.001).
#' @param epsilon Epsilon in the insensitive-loss function
#' (default: 0.1).
#' @param cross Number of cross validation folds applied
#' to assess the quality of the model (default: 10).
#' @param ... Other parameters to be passed to e1071::svm function.
#' @return Model fitted to input data
#' (to be passed to \code{\link[sits]{sits_classify}})
#' @note
#' Please refer to the sits documentation available in
#' <https://e-sensing.github.io/sitsbook/> for detailed examples.
#' @examples
#' if (sits_run_examples()) {
#' # Example of training a model for time series classification
#' # Retrieve the samples for Mato Grosso
#' # train an SVM model
#' ml_model <- sits_train(samples_modis_ndvi, ml_method = sits_svm)
#' # classify the point
#' point_ndvi <- sits_select(point_mt_6bands, bands = "NDVI")
#' # classify the point
#' point_class <- sits_classify(
#' data = point_ndvi, ml_model = ml_model
#' )
#' plot(point_class)
#' }
#' @export
#'
sits_svm <- function(samples = NULL, formula = sits_formula_linear(),
scale = FALSE, cachesize = 1000, kernel = "radial",
degree = 3, coef0 = 0, cost = 10, tolerance = 0.001,
epsilon = 0.1, cross = 10, ...) {
.check_set_caller("sits_svm")
# Function that trains a support vector machine model
train_fun <- function(samples) {
# Verifies if e1071 package is installed
.check_require_packages("e1071")
# Get labels (used later to ensure column order in result matrix)
labels <- .samples_labels(samples)
# Get normalized training samples
ml_stats <- .samples_stats(samples)
# Get predictors features
train_samples <- .predictors(samples)
# Normalize predictors
train_samples <- .pred_normalize(pred = train_samples, stats = ml_stats)
# Post condition: is predictor data valid?
.check_predictors(pred = train_samples, samples = samples)
# Update formula parameter
if (inherits(formula, "function")) {
formula <- formula(train_samples)
}
# Train an svm model
model <- e1071::svm(
formula = formula, data = train_samples, scale = scale,
kernel = kernel, degree = degree, cost = cost, coef0 = coef0,
cachesize = cachesize, tolerance = tolerance, epsilon = epsilon,
cross = cross, probability = TRUE, ..., na.action = stats::na.fail
)
# Function that predicts labels of input values
predict_fun <- function(values) {
# Verifies if e1071 package is installed
.check_require_packages("e1071")
# Used to check values (below)
input_pixels <- nrow(values)
# Performs data normalization
values <- .pred_normalize(pred = values, stats = ml_stats)
# Do classification
values <- stats::predict(
object = model, newdata = values, probability = TRUE
)
# Get the predicted probabilities
values <- attr(values, "probabilities")
# Are the results consistent with the data input?
.check_processed_values(values, input_pixels)
# Reorder matrix columns if needed
if (any(labels != colnames(values))) {
values <- values[, labels]
}
return(values)
}
# Set model class
predict_fun <- .set_class(
predict_fun, "svm_model", "sits_model", class(predict_fun)
)
return(predict_fun)
}
# If samples is informed, train a model and return a predict function
# Otherwise give back a train function to train model further
result <- .factory_function(samples, train_fun)
return(result)
}
#' @title Train extreme gradient boosting models
#' @name sits_xgboost
#' @author Rolf Simoes, \email{rolf.simoes@@inpe.br}
#' @author Gilberto Camara, \email{gilberto.camara@@inpe.br}
#'
#' @description This function uses the extreme gradient boosting algorithm.
#' Boosting iteratively adds basis functions in a greedy fashion
#' so that each new basis function further reduces the selected loss function.
#' This function is a front-end to the methods in the "xgboost" package.
#' Please refer to the documentation in that package for more details.
#'
#' @references Tianqi Chen, Carlos Guestrin,
#' "XGBoost : Reliable Large-scale Tree Boosting System",
#' SIG KDD 2016.
#'
#' @param samples Time series with the training samples.
#' @param learning_rate Learning rate: scale the contribution
#' of each tree by a factor of 0 < lr < 1
#' when it is added to the current approximation.
#' Used to prevent overfitting. Default: 0.15
#' @param min_split_loss Minimum loss reduction to make a further
#' partition of a leaf. Default: 1.
#' @param max_depth Maximum depth of a tree.
#' Increasing this value makes the model more complex
#' and more likely to overfit. Default: 5.
#' @param min_child_weight If the leaf node has a minimum sum of instance
#' weights lower than min_child_weight,
#' tree splitting stops. The larger min_child_weight is,
#' the more conservative the algorithm is. Default: 1.
#' @param max_delta_step Maximum delta step we allow each leaf output to be.
#' If the value is set to 0, there is no constraint.
#' If it is set to a positive value, it can help making
#' the update step more conservative. Default: 1.
#' @param subsample Percentage of samples supplied to a tree.
#' Default: 0.8.
#' @param nfold Number of the subsamples for the cross-validation.
#' @param nrounds Number of rounds to iterate the cross-validation
#' (default: 100)
#' @param nthread Number of threads (default = 6)
#' @param early_stopping_rounds Training with a validation set will stop
#' if the performance doesn't improve for k rounds.
#' @param verbose Print information on statistics during the process
#' @return Model fitted to input data
#' (to be passed to \code{\link[sits]{sits_classify}})
#'
#' @note
#' Please refer to the sits documentation available in
#' <https://e-sensing.github.io/sitsbook/> for detailed examples.
#' @examples
#' if (sits_run_examples()) {
#' # Example of training a model for time series classification
#' # Retrieve the samples for Mato Grosso
#' # train a xgboost model
#' ml_model <- sits_train(samples_modis_ndvi, ml_method = sits_xgboost)
#' # classify the point
#' point_ndvi <- sits_select(point_mt_6bands, bands = "NDVI")
#' # classify the point
#' point_class <- sits_classify(
#' data = point_ndvi, ml_model = ml_model
#' )
#' plot(point_class)
#' }
#' @export
#'
sits_xgboost <- function(samples = NULL, learning_rate = 0.15,
min_split_loss = 1, max_depth = 5,
min_child_weight = 1, max_delta_step = 1,
subsample = 0.8, nfold = 5, nrounds = 100,
nthread = 6,
early_stopping_rounds = 20, verbose = FALSE) {
.check_set_caller("sits_xgboost")
# Function that trains a xgb model
train_fun <- function(samples) {
# verifies if xgboost package is installed
.check_require_packages("xgboost")
# Get labels (used later to ensure column order in result matrix)
labels <- .samples_labels(samples)
# Get predictors features
train_samples <- .predictors(samples)
# Post condition: is predictor data valid?
.check_predictors(pred = train_samples, samples = samples)
# Transform labels to integer code before train
code_labels <- seq_along(labels)
names(code_labels) <- labels
# Reference labels for each sample expressed as numerical values
references <-
unname(code_labels[.pred_references(train_samples)]) - 1
# Define the parameters of the model
params <- list(
booster = "gbtree", objective = "multi:softprob",
eval_metric = "mlogloss", eta = learning_rate,
gamma = min_split_loss, max_depth = max_depth,
min_child_weight = min_child_weight,
max_delta_step = max_delta_step, subsample = subsample,
nthread = nthread
)
if (verbose)
verbose <- 1
else
verbose <- 0
# transform predictors in a xgb.DMatrix
xgb_matrix <- xgboost::xgb.DMatrix(
data = as.matrix(.pred_features(train_samples)),
label = references)
# train the model
model <- xgboost::xgb.train(xgb_matrix,
num_class = length(labels), params = params,
nrounds = nrounds, verbose = verbose
)
# Get best ntreelimit
ntreelimit <- model[["best_ntreelimit"]]
# Function that predicts labels of input values
predict_fun <- function(values) {
# Verifies if xgboost package is installed
.check_require_packages("xgboost")
# Used to check values (below)
input_pixels <- nrow(values)
# Do classification
values <- stats::predict(
object = model, as.matrix(values), ntreelimit = ntreelimit,
reshape = TRUE
)
# Are the results consistent with the data input?
.check_processed_values(values, input_pixels)
# Update the columns names to labels
colnames(values) <- labels
return(values)
}
# Set model class
predict_fun <- .set_class(
predict_fun, "xgb_model", "sits_model", class(predict_fun)
)
return(predict_fun)
}
# If samples is informed, train a model and return a predict function
# Otherwise give back a train function to train model further
result <- .factory_function(samples, train_fun)
return(result)
}
#' @title Define a loglinear formula for classification models
#' @name sits_formula_logref
#'
#' @author Alexandre Ywata de Carvalho, \email{alexandre.ywata@@ipea.gov.br}
#' @author Rolf Simoes, \email{rolf.simoes@@inpe.br}
#'
#' @description A function to be used as a symbolic description
#' of some fitting models such as svm and random forest.
#' This function tells the models to do a log transformation of the inputs.
#' The `predictors_index` parameter informs
#' the positions of `tb` fields corresponding to formula independent variables.
#' If no value is given, the default is NULL,
#' a value indicating that all fields will be used as predictors.
#'
#' @param predictors_index Index of the valid columns
#' to compose formula (default: -2:0).
#' @return A function that computes a valid formula using a log function.
#' @examples
#' if (sits_run_examples()) {
#' # Example of training a model for time series classification
#' # Retrieve the samples for Mato Grosso
#' # train an SVM model
#' ml_model <- sits_train(samples_modis_ndvi,
#' ml_method = sits_svm(formula = sits_formula_logref())
#' )
#' # classify the point
#' point_ndvi <- sits_select(point_mt_6bands, bands = "NDVI")
#' # classify the point
#' point_class <- sits_classify(
#' data = point_ndvi, ml_model = ml_model
#' )
#' plot(point_class)
#' }
#' @export
#'
sits_formula_logref <- function(predictors_index = -2:0) {
# set caller to show in errors
.check_set_caller("sits_formula_logref")
# store configuration information about model formula
sits_env[["model_formula"]] <- "log"
# this function returns a formula like
# 'factor(reference~log(f1)+log(f2)+...+log(fn)' where f1, f2, ..., fn are
# the predictor fields given by the predictor index.
result_fun <- function(tb) {
.check_that(nrow(tb) > 0)
n_rows_tb <- nrow(tb)
# if no predictors_index are given, assume all tb's fields are used
if (!.has(predictors_index))
predictors_index <- 1:n_rows_tb
# get predictors names
categories <- names(tb)[c(predictors_index)]
# compute formula result
result_for <- stats::as.formula(paste0(
"factor(label)~",
paste0(paste0("log(`", categories, "`)"),
collapse = "+"
)
))
return(result_for)
}
return(result_fun)
}
#' @title Define a linear formula for classification models
#' @name sits_formula_linear
#' @author Gilberto Camara, \email{gilberto.camara@@inpe.br}
#' @author Alexandre Ywata de Carvalho, \email{alexandre.ywata@@ipea.gov.br}
#' @author Rolf Simoes, \email{rolf.simoes@@inpe.br}
#'
#' @description Provides a symbolic description of a fitting model.
#' Tells the model to do a linear transformation of the input values.
#' The `predictors_index` parameter informs the positions of fields
#' corresponding to formula independent variables.
#' If no value is given, that all fields will be used as predictors.
#'
#' @param predictors_index Index of the valid columns
#' whose names are used to compose formula (default: -2:0).
#' @return A function that computes a valid formula using a linear function.
#' @examples
#' if (sits_run_examples()) {
#' # Example of training a model for time series classification
#' # Retrieve the samples for Mato Grosso
#' # train an SVM model
#' ml_model <- sits_train(samples_modis_ndvi,
#' ml_method = sits_svm(formula = sits_formula_logref())
#' )
#' # classify the point
#' point_ndvi <- sits_select(point_mt_6bands, bands = "NDVI")
#' # classify the point
#' point_class <- sits_classify(
#' data = point_ndvi, ml_model = ml_model
#' )
#' plot(point_class)
#' }
#' @export
#'
sits_formula_linear <- function(predictors_index = -2:0) {
# set caller to show in errors
.check_set_caller("sits_formula_linear")
# store configuration information about model formula
sits_env[["model_formula"]] <- "linear"
# this function returns a formula
# 'factor(reference~log(f1)+log(f2)+...+log(fn)' where f1, f2, ..., fn are
# the predictor fields.
result_fun <- function(tb) {
.check_that(nrow(tb) > 0)
n_rows_tb <- nrow(tb)
# if no predictors_index are given, assume that all fields are used
if (!.has(predictors_index))
predictors_index <- 1:n_rows_tb
# get predictors names
categories <- names(tb)[c(predictors_index)]
# compute formula result
result_for <- stats::as.formula(paste0(
"factor(label)~",
paste0(paste0(categories,
collapse = "+"
))
))
return(result_for)
}
return(result_fun)
}