diff --git a/NAMESPACE b/NAMESPACE index de50687..43afa78 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -22,6 +22,7 @@ S3method(print,ggwrm) S3method(print,gtdrm) S3method(print,gwavgm) S3method(print,gwcorrm) +S3method(print,gwdam) S3method(print,gwlcrm) S3method(print,gwpcam) S3method(print,gwrm) @@ -39,6 +40,7 @@ export(gtdr_config) export(gwaverage) export(gwcorr_config) export(gwcorrelation) +export(gwda) export(gwpca) export(gwr_basic) export(gwr_generalized) diff --git a/R/RcppExports.R b/R/RcppExports.R index 61ea322..36d3bfd 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -13,6 +13,10 @@ gw_correlation_cal <- function(x1, x2, coords, bw, adaptive, kernel, longlat, p, .Call(`_GWmodel3_gw_correlation_cal`, x1, x2, coords, bw, adaptive, kernel, longlat, p, theta, initial_type, optim_bw_criterion, parallel_type, parallel_arg, variable_names, verbose) } +gwda_cal <- function(x, y, coords, bw, adaptive, kernel, longlat, p, theta, method, parallel_type, parallel_arg) { + .Call(`_GWmodel3_gwda_cal`, x, y, coords, bw, adaptive, kernel, longlat, p, theta, method, parallel_type, parallel_arg) +} + gwpca_cal <- function(x, coords, bw, adaptive, kernel, longlat, p, theta, keep_components) { .Call(`_GWmodel3_gwpca_cal`, x, coords, bw, adaptive, kernel, longlat, p, theta, keep_components) } diff --git a/R/gwda.R b/R/gwda.R new file mode 100644 index 0000000..3835453 --- /dev/null +++ b/R/gwda.R @@ -0,0 +1,180 @@ +#' Geographically Weighted Discriminant Analysis +#' +#' @param formula Regresison model. +#' @param data A `sf` objects. +#' @param bw Bandwidth value +#' @param adaptive Whether the bandwidth value is adaptive or not. +#' @param quantile Whether to calculate local quantiles. +#' @param kernel Kernel function used. +#' @param longlat Whether the coordinates +#' @param p Power of the Minkowski distance, +#' default to 2, i.e., Euclidean distance. +#' @param theta Angle in radian to roate the coordinate system, default to 0. +#' @param method Discriminant method for GWDA, can be set as +#' - `wqda`, weighted quadratic discriminant analysis +#' - `wlda`, weighted linear discriminant analysis +#' If not provided, this parameter will be set to `wqda`. +#' @param parallel_method Parallel method. +#' @param parallel_arg Parallel method argument. +#' +#' @return A `gwdam` object. +#' +#' @export +gwda <- function( + formula, + data, + bw, + adaptive = FALSE, + kernel = c("gaussian", "exp", "bisquare", "tricube", "boxcar"), + longlat = FALSE, + p = 2.0, + theta = 0.0, + method = c("wqda", "wlda"), + parallel_method = c("no", "omp"), + parallel_arg = c(0)) +{ + ### Check args + kernel <- match.arg(kernel) + parallel_method <- match.arg(parallel_method) + attr(data, "na.action") <- getOption("na.action") + + ### Extract coords + data <- do.call(na.action(data), args = list(data)) + coords <- as.matrix(sf::st_coordinates(sf::st_centroid(data))) + if (is.null(coords) || nrow(coords) != nrow(data)) { + stop("Missing coordinates.") + } + + ### Extract variables + mc <- match.call(expand.dots = FALSE) + mf <- model.frame(formula = formula(update(formula, ~ . - 1)), data = sf::st_drop_geometry(data)) + mt <- attr(mf, "terms") + y <- model.extract(mf, "response") + y <- as.character(y) + x <- model.matrix(mt, mf) + dep_var <- as.character(attr(terms(mf), "variables")[[2]]) + indep_vars <- colnames(x) + + ### Check whether bandwidth is valid. + if (missing(bw)) { + stop("Bandwidth is missing.") + } else if (is.numeric(bw) || is.integer(bw)) { + if (bw == 0) stop("Bandwidth must be non-zero.") + } else { + stop("Bandwidth must be a number.") + } + + method <- match.arg(method) + wqda_method <- ifelse(method == "wqda", TRUE, FALSE) + + c_results <- gwda_cal( + x, y, coords, bw, + adaptive, enum(kernel), + longlat, p, theta, + wqda_method, + enum_list(parallel_method, parallel_types), parallel_arg + ) + + sdf_data <- data.frame() + res <- c_results$res + group <- as.data.frame(c_results$group) + probs <- c_results$probs + pmax <- c_results$pmax + entropy <- c_results$entropy + + colnames(res) <- indep_vars + colnames(group) <- "groups" + colnames(probs) <- paste0(indep_vars, "_probs") + colnames(pmax) <- "pmax" + colnames(entropy) <- "entropy" + sdf_data <- as.data.frame(cbind(res, group, probs, pmax, entropy)) + + sdf_data$geometry <- sf::st_geometry(data) + sdf <- sf::st_sf(sdf_data) + + ### Return result + gwdam <- list( + SDF = sdf, + args = list( + x = x, + y = y, + correctRate = c_results$correctRate, + coords = coords, + bw = bw, + adaptive = adaptive, + kernel = kernel, + longlat = longlat, + p = p, + theta = theta, + method = method, + parallel_method = parallel_method, + parallel_arg = parallel_arg + ), + call = mc, + indep_vars = indep_vars, + dep_var = dep_var + ) + class(gwdam) <- "gwdam" + gwdam +} + +#' Print description of a `gwdam` object +#' +#' @param x An `gwdam` object returned by [gwr_basic()]. +#' @param decimal_fmt The format string passing to [base::sprintf()]. +#' @inheritDotParams print_table_md +#' +#' @method print gwdam +#' @rdname print +#' @export +print.gwdam <- function(x, ..., decimal_fmt) { + if (!inherits(x, "gwdam")) { + stop("It's not a 'gwdam' object.") + } + + cat(" ***********************************************************************\n") + cat(" * Package GWmodel3 *\n") + cat(" ***********************************************************************\n") + cat(" * Results of Geographically Discriminant Analysis *\n") + cat(" ***********************************************************************\n") + cat("\n *********************Model Calibration Information*********************\n") + + cat(" Formula:", deparse(x$call$formula), fill = T) + cat(" Data:", deparse(x$call$data), fill = T) + cat(" Method:", x$args$method, fill = T) + cat(" Number of summary points:", nrow(x$args$x), fill = T) + cat(" Kernel function:", x$args$kernel, fill = T) + cat(" Bandwidth:", x$args$bw, + ifelse(x$args$adaptive, "(Nearest Neighbours)", "(Meters)"), + fill = T + ) + distance_type <- "Euclidean" + if (x$args$longlat) { + distance_type <- "Geodetic" + } else if (x$args$p == 2) { + distance_type <- "Euclidean" + } else if (x$args$p == 1) { + distance_type <- "Manhattan" + } else if (is.infinite(x$args$p)) { + distance_type <- "Chebyshev" + } else { + distance_type <- "Generalized Minkowski" + } + distance_rotated <- (x$args$theta != 0 && x$args$p != 2 && !x$args$longlat) + cat(" Distance:", distance_type, ifelse(distance_rotated, " (rotated)", ""), fill = T) + res <- st_drop_geometry(x$SDF) + group_counts <- as.data.frame(table(res$groups)) + cat(" Discrimination Result:", fill = T) + colnames(group_counts) <- c("Group", "Count") + output <- capture.output(print(group_counts, row.names = FALSE)) + output <- paste0(" ", output) + cat(output, sep = "\n") + + cat("\n ***********************Local Summary Statistics************************", fill = T) + numeric_res <- res[sapply(res, is.numeric)] + df <- as.data.frame(t(sapply(numeric_res, summary))) + output <- capture.output(print(df)) + output <- paste0(" ", output) + cat(output, sep = "\n") + cat(" ***********************************************************************\n") +} \ No newline at end of file diff --git a/man/gwda.Rd b/man/gwda.Rd new file mode 100644 index 0000000..4eb683e --- /dev/null +++ b/man/gwda.Rd @@ -0,0 +1,57 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/gwda.R +\name{gwda} +\alias{gwda} +\title{Geographically Weighted Discriminant Analysis} +\usage{ +gwda( + formula, + data, + bw, + adaptive = FALSE, + kernel = c("gaussian", "exp", "bisquare", "tricube", "boxcar"), + longlat = FALSE, + p = 2, + theta = 0, + method = c("wqda", "wlda"), + parallel_method = c("no", "omp"), + parallel_arg = c(0) +) +} +\arguments{ +\item{formula}{Regresison model.} + +\item{data}{A \code{sf} objects.} + +\item{bw}{Bandwidth value} + +\item{adaptive}{Whether the bandwidth value is adaptive or not.} + +\item{kernel}{Kernel function used.} + +\item{longlat}{Whether the coordinates} + +\item{p}{Power of the Minkowski distance, +default to 2, i.e., Euclidean distance.} + +\item{theta}{Angle in radian to roate the coordinate system, default to 0.} + +\item{method}{Discriminant method for GWDA, can be set as +\itemize{ +\item \code{wqda}, weighted quadratic discriminant analysis +\item \code{wlda}, weighted linear discriminant analysis +If not provided, this parameter will be set to \code{wqda}. +}} + +\item{parallel_method}{Parallel method.} + +\item{parallel_arg}{Parallel method argument.} + +\item{quantile}{Whether to calculate local quantiles.} +} +\value{ +A \code{gwdam} object. +} +\description{ +Geographically Weighted Discriminant Analysis +} diff --git a/man/print.Rd b/man/print.Rd index 95a0c3c..f0d97e6 100644 --- a/man/print.Rd +++ b/man/print.Rd @@ -1,11 +1,12 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/gtdr.R, R/gwaverage.R, R/gwcorrelation.R, -% R/gwpca.R, R/gwr_basic.R, R/gwr_generalized.R, R/gwr_localcollinearity.R, -% R/gwr_multiscale.R +% R/gwda.R, R/gwpca.R, R/gwr_basic.R, R/gwr_generalized.R, +% R/gwr_localcollinearity.R, R/gwr_multiscale.R \name{print.gtdrm} \alias{print.gtdrm} \alias{print.gwavgm} \alias{print.gwcorrm} +\alias{print.gwdam} \alias{print.gwpcam} \alias{print.gwrm} \alias{print.ggwrm} @@ -19,6 +20,8 @@ \method{print}{gwcorrm}(x, ..., decimal_fmt = "\%.3f") +\method{print}{gwdam}(x, ..., decimal_fmt) + \method{print}{gwpcam}(x, decimal_fmt = "\%.3f", ...) \method{print}{gwrm}(x, decimal_fmt = "\%.3f", ...) @@ -35,7 +38,7 @@ \item{decimal_fmt}{The format string passing to \code{\link[base:sprintf]{base::sprintf()}}.} \item{...}{ - Arguments passed on to \code{\link[=print_table_md]{print_table_md}}, \code{\link[=print_table_md]{print_table_md}}, \code{\link[=print_table_md]{print_table_md}}, \code{\link[=print_table_md]{print_table_md}}, \code{\link[=print_table_md]{print_table_md}}, \code{\link[=print_table_md]{print_table_md}}, \code{\link[=print_table_md]{print_table_md}}, \code{\link[=print_table_md]{print_table_md}} + Arguments passed on to \code{\link[=print_table_md]{print_table_md}}, \code{\link[=print_table_md]{print_table_md}}, \code{\link[=print_table_md]{print_table_md}}, \code{\link[=print_table_md]{print_table_md}}, \code{\link[=print_table_md]{print_table_md}}, \code{\link[=print_table_md]{print_table_md}}, \code{\link[=print_table_md]{print_table_md}}, \code{\link[=print_table_md]{print_table_md}}, \code{\link[=print_table_md]{print_table_md}} \describe{ \item{\code{col.sep}}{Column seperator. Default to \code{""}.} \item{\code{header.sep}}{Header seperator. Default to \code{"-"}.} @@ -54,6 +57,8 @@ Print description of a \code{gwavgm} object Print description of a \code{gwcorrm} object +Print description of a \code{gwdam} object + Print description of a \code{gwpcam} object Print description of a \code{gwrm} object diff --git a/src/Makevars.in b/src/Makevars.in index 9201bef..ad401ef 100644 --- a/src/Makevars.in +++ b/src/Makevars.in @@ -50,6 +50,7 @@ OBJECTS_LIBGWMODEL = \ libgwmodel/src/gwmodelpp/GWRGeneralized.o \ libgwmodel/src/gwmodelpp/GWRMultiscale.o \ libgwmodel/src/gwmodelpp/GWAverage.o \ + libgwmodel/src/gwmodelpp/GWDA.o \ libgwmodel/src/gwmodelpp/GWCorrelation.o \ libgwmodel/src/gwmodelpp/GWRLocalCollinearity.o \ libgwmodel/src/gwmodelpp/GWPCA.o \ @@ -73,6 +74,7 @@ OBJECTS_GWMODEL = \ gwr_localcollinearity.o \ gtdr.o \ gwaverage.o \ + gwda.o \ gwcorrelation.o \ gwpca.o \ RcppExports.o diff --git a/src/Makevars.win b/src/Makevars.win index f657e3b..e782fd1 100644 --- a/src/Makevars.win +++ b/src/Makevars.win @@ -52,6 +52,7 @@ OBJECTS_LIBGWMODEL = \ libgwmodel/src/gwmodelpp/GWRGeneralized.o \ libgwmodel/src/gwmodelpp/GWRMultiscale.o \ libgwmodel/src/gwmodelpp/GWAverage.o \ + libgwmodel/src/gwmodelpp/GWDA.o \ libgwmodel/src/gwmodelpp/GWCorrelation.o \ libgwmodel/src/gwmodelpp/GWRLocalCollinearity.o \ libgwmodel/src/gwmodelpp/GWPCA.o \ @@ -75,6 +76,7 @@ OBJECTS_GWMODEL = \ gwr_localcollinearity.o \ gtdr.o \ gwaverage.o \ + gwda.o \ gwcorrelation.o \ gwpca.o \ RcppExports.o diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 91f4dc0..1419e0a 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -86,6 +86,28 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// gwda_cal +List gwda_cal(const arma::mat& x, std::vector& y, const arma::mat& coords, double bw, bool adaptive, size_t kernel, bool longlat, double p, double theta, bool method, size_t parallel_type, const IntegerVector& parallel_arg); +RcppExport SEXP _GWmodel3_gwda_cal(SEXP xSEXP, SEXP ySEXP, SEXP coordsSEXP, SEXP bwSEXP, SEXP adaptiveSEXP, SEXP kernelSEXP, SEXP longlatSEXP, SEXP pSEXP, SEXP thetaSEXP, SEXP methodSEXP, SEXP parallel_typeSEXP, SEXP parallel_argSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const arma::mat& >::type x(xSEXP); + Rcpp::traits::input_parameter< std::vector& >::type y(ySEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type coords(coordsSEXP); + Rcpp::traits::input_parameter< double >::type bw(bwSEXP); + Rcpp::traits::input_parameter< bool >::type adaptive(adaptiveSEXP); + Rcpp::traits::input_parameter< size_t >::type kernel(kernelSEXP); + Rcpp::traits::input_parameter< bool >::type longlat(longlatSEXP); + Rcpp::traits::input_parameter< double >::type p(pSEXP); + Rcpp::traits::input_parameter< double >::type theta(thetaSEXP); + Rcpp::traits::input_parameter< bool >::type method(methodSEXP); + Rcpp::traits::input_parameter< size_t >::type parallel_type(parallel_typeSEXP); + Rcpp::traits::input_parameter< const IntegerVector& >::type parallel_arg(parallel_argSEXP); + rcpp_result_gen = Rcpp::wrap(gwda_cal(x, y, coords, bw, adaptive, kernel, longlat, p, theta, method, parallel_type, parallel_arg)); + return rcpp_result_gen; +END_RCPP +} // gwpca_cal List gwpca_cal(const arma::mat& x, const arma::mat& coords, double bw, bool adaptive, size_t kernel, bool longlat, double p, double theta, size_t keep_components); RcppExport SEXP _GWmodel3_gwpca_cal(SEXP xSEXP, SEXP coordsSEXP, SEXP bwSEXP, SEXP adaptiveSEXP, SEXP kernelSEXP, SEXP longlatSEXP, SEXP pSEXP, SEXP thetaSEXP, SEXP keep_componentsSEXP) { @@ -279,6 +301,7 @@ static const R_CallMethodDef CallEntries[] = { {"_GWmodel3_gtdr_fit", (DL_FUNC) &_GWmodel3_gtdr_fit, 19}, {"_GWmodel3_gwaverage_fit", (DL_FUNC) &_GWmodel3_gwaverage_fit, 11}, {"_GWmodel3_gw_correlation_cal", (DL_FUNC) &_GWmodel3_gw_correlation_cal, 15}, + {"_GWmodel3_gwda_cal", (DL_FUNC) &_GWmodel3_gwda_cal, 12}, {"_GWmodel3_gwpca_cal", (DL_FUNC) &_GWmodel3_gwpca_cal, 9}, {"_GWmodel3_gwr_basic_fit", (DL_FUNC) &_GWmodel3_gwr_basic_fit, 22}, {"_GWmodel3_gwr_basic_predict", (DL_FUNC) &_GWmodel3_gwr_basic_predict, 14}, diff --git a/src/gwda.cpp b/src/gwda.cpp new file mode 100644 index 0000000..7265d64 --- /dev/null +++ b/src/gwda.cpp @@ -0,0 +1,90 @@ +// [[Rcpp::depends(RcppArmadillo)]] +#include +#include "utils.h" +#include "gwmodelpp/GWDA.h" + +using namespace std; +using namespace Rcpp; +using namespace arma; +using namespace gwm; + +// [[Rcpp::export]] +List gwda_cal( + const arma::mat &x, + std::vector &y, + const arma::mat &coords, + double bw, + bool adaptive, + size_t kernel, + bool longlat, + double p, + double theta, + bool method, + size_t parallel_type, + const IntegerVector ¶llel_arg) +{ + + // Make Spatial Weight + BandwidthWeight bandwidth(bw, adaptive, BandwidthWeight::KernelFunctionType((size_t)kernel)); + + Distance *distance = nullptr; + if (longlat) + { + distance = new CRSDistance(true); + } + else + { + if (p == 2.0 && theta == 0.0) + { + distance = new CRSDistance(false); + } + else + { + distance = new MinkwoskiDistance(p, theta); + } + } + + SpatialWeight spatial(&bandwidth, distance); + + GWDA algorithm; + algorithm.setCoords(coords); + algorithm.setVariables(x); + algorithm.setGroup(y); + algorithm.setSpatialWeight(spatial); + + vector vpar_args = as>(IntegerVector(parallel_arg)); + switch (ParallelType(size_t(parallel_type))) + { +#ifdef _OPENMP + case ParallelType::OpenMP: + algorithm.setParallelType(ParallelType::OpenMP); + algorithm.setOmpThreadNum(vpar_args[0]); + break; +#endif + default: + algorithm.setParallelType(ParallelType::SerialOnly); + break; + } + + algorithm.setIsWqda(method); + + try + { + algorithm.run(); + } + catch(const std::exception& e) + { + stop(e.what()); + } + + List results; + results = List::create( + Named("correctRate") = algorithm.correctRate(), + Named("group") = algorithm.group(), + Named("res") = algorithm.res(), + Named("probs") = algorithm.probs(), + Named("pmax") = algorithm.pmax(), + Named("entropy") = algorithm.entropy()); + + return results; +} \ No newline at end of file diff --git a/tests/testthat/test-gwda.R b/tests/testthat/test-gwda.R new file mode 100644 index 0000000..774c698 --- /dev/null +++ b/tests/testthat/test-gwda.R @@ -0,0 +1,19 @@ +data(LondonHP) + +test_that("GWDA: wqda", { + m1 <<- expect_no_error({ + gwda(TYPEDETCH~FLOORSZ+UNEMPLOY, LondonHP, 64, TRUE) + }) +}) + +test_that("GWDA: wlda", { + m2 <<- expect_no_error({ + gwda(TYPEDETCH~FLOORSZ+UNEMPLOY, LondonHP, 64, TRUE, method = "wlda") + }) +}) + +test_that("GWDA: omp", { + expect_no_error({ + gwda(TYPEDETCH~FLOORSZ+UNEMPLOY, LondonHP, 64, TRUE, parallel_method = "omp", parallel_arg = 2) + }) +}) \ No newline at end of file