From 390b39f38cac9b4da00c58273d6c4f078ae2916a Mon Sep 17 00:00:00 2001 From: HPDell Date: Thu, 5 Oct 2023 15:08:30 +0100 Subject: [PATCH 1/4] edit: add a benchmark --- test/CMakeLists.txt | 1 + test/testGWRMultiscale.cpp | 45 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 46 insertions(+) diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index f5fe3bc2..f808657c 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -29,6 +29,7 @@ set(SAMPLE_DATA_DIR ${CMAKE_CURRENT_SOURCE_DIR}/data) add_definitions(-DSAMPLE_DATA_DIR="${SAMPLE_DATA_DIR}") set(ADDON_SOURCES + data/londonhp.cpp data/londonhp100.cpp TerminateCheckTelegram.cpp FileTelegram.cpp diff --git a/test/testGWRMultiscale.cpp b/test/testGWRMultiscale.cpp index 5d70a997..3e95438b 100644 --- a/test/testGWRMultiscale.cpp +++ b/test/testGWRMultiscale.cpp @@ -10,6 +10,7 @@ #include "gwmodelpp/spatialweight/CRSDistance.h" #include "gwmodelpp/spatialweight/BandwidthWeight.h" #include "gwmodelpp/spatialweight/SpatialWeight.h" +#include "londonhp.h" #include "londonhp100.h" #include "TerminateCheckTelegram.h" #include "FileTelegram.h" @@ -336,3 +337,47 @@ TEST_CASE("Multiscale GWR: cancel") } } + +TEST_CASE("Large Data Repeate") +{ + mat londonhp_coord, londonhp_data; + vector londonhp_fields; + if (!read_londonhp(londonhp_coord, londonhp_data, londonhp_fields)) + { + FAIL("Cannot load londonhp100 data."); + } + + vec y = londonhp_data.col(0); + mat x = join_rows(ones(londonhp_data.n_rows), londonhp_data.cols(uvec({1, 3}))); + uword nVar = 3; + + BENCHMARK("specified bandwidth") + { + vector spatials; + vector preditorCentered; + vector bandwidthInitialize; + vector bandwidthSelectionApproach; + for (size_t i = 0; i < nVar; i++) + { + CRSDistance distance; + BandwidthWeight bandwidth(64, true, BandwidthWeight::Bisquare); + spatials.push_back(SpatialWeight(&bandwidth, &distance)); + preditorCentered.push_back(i != 0); + bandwidthInitialize.push_back(GWRMultiscale::BandwidthInitilizeType::Specified); + bandwidthSelectionApproach.push_back(GWRMultiscale::BandwidthSelectionCriterionType::CV); + } + + GWRMultiscale algorithm; + algorithm.setCoords(londonhp_coord); + algorithm.setDependentVariable(y); + algorithm.setIndependentVariables(x); + algorithm.setSpatialWeights(spatials); + algorithm.setHasHatMatrix(true); + algorithm.setPreditorCentered(preditorCentered); + algorithm.setBandwidthInitilize(bandwidthInitialize); + algorithm.setBandwidthSelectionApproach(bandwidthSelectionApproach); + algorithm.setBandwidthSelectThreshold(vector(3, 1e-5)); + algorithm.fit(); + return algorithm.betas(); + }; +} From de8dbca65da40fe4b9ed10de044aeb58da34bc68 Mon Sep 17 00:00:00 2001 From: HPDell Date: Thu, 5 Oct 2023 15:14:11 +0100 Subject: [PATCH 2/4] edit: test set bandwidth as initial value --- test/testGWRMultiscale.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/testGWRMultiscale.cpp b/test/testGWRMultiscale.cpp index 3e95438b..02b71ede 100644 --- a/test/testGWRMultiscale.cpp +++ b/test/testGWRMultiscale.cpp @@ -363,7 +363,7 @@ TEST_CASE("Large Data Repeate") BandwidthWeight bandwidth(64, true, BandwidthWeight::Bisquare); spatials.push_back(SpatialWeight(&bandwidth, &distance)); preditorCentered.push_back(i != 0); - bandwidthInitialize.push_back(GWRMultiscale::BandwidthInitilizeType::Specified); + bandwidthInitialize.push_back(GWRMultiscale::BandwidthInitilizeType::Initial); bandwidthSelectionApproach.push_back(GWRMultiscale::BandwidthSelectionCriterionType::CV); } From e0ebbb9cad0405552113312480f8ba4bcb899e78 Mon Sep 17 00:00:00 2001 From: HPDell Date: Thu, 5 Oct 2023 15:20:58 +0100 Subject: [PATCH 3/4] edit: test add benchmark label --- test/testGWRMultiscale.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/testGWRMultiscale.cpp b/test/testGWRMultiscale.cpp index 02b71ede..5ceda765 100644 --- a/test/testGWRMultiscale.cpp +++ b/test/testGWRMultiscale.cpp @@ -338,7 +338,7 @@ TEST_CASE("Multiscale GWR: cancel") } -TEST_CASE("Large Data Repeate") +TEST_CASE("Large Data Repeate", "[benchmark]") { mat londonhp_coord, londonhp_data; vector londonhp_fields; From 8c1e0838c1d02a720292199ed6e7c1681c9a07bd Mon Sep 17 00:00:00 2001 From: HPDell Date: Thu, 5 Oct 2023 15:21:33 +0100 Subject: [PATCH 4/4] edit: optimise matrix inverse --- src/gwmodelpp/GWRMultiscale.cpp | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/src/gwmodelpp/GWRMultiscale.cpp b/src/gwmodelpp/GWRMultiscale.cpp index ed16facd..b19d0e42 100644 --- a/src/gwmodelpp/GWRMultiscale.cpp +++ b/src/gwmodelpp/GWRMultiscale.cpp @@ -426,17 +426,17 @@ vec GWRMultiscale::fitVarSerial(const vec &x, const vec &y, const uword var, mat if (mHasHatMatrix) { mat ci, si; - S = mat(mHasHatMatrix ? nDp : 1, nDp, fill::zeros); + S = mat(nDp, nDp, fill::zeros); for (uword i = 0; i < nDp ; i++) { GWM_LOG_STOP_BREAK(mStatus); vec w = mSpatialWeights[var].weightVector(i); mat xtw = trans(x % w); - mat xtwx = xtw * x; - mat xtwy = xtw * y; + double xtwx = as_scalar(xtw * x); + double xtwy = as_scalar(xtw * y); try { - mat xtwx_inv = inv_sympd(xtwx); + double xtwx_inv = 1.0 / xtwx; betas.col(i) = xtwx_inv * xtwy; ci = xtwx_inv * xtw; si = x(i) * ci; @@ -570,12 +570,12 @@ double GWRMultiscale::bandwidthSizeCriterionVarCVSerial(BandwidthWeight *bandwid vec w = bandwidthWeight->weight(d); w(i) = 0.0; mat xtw = trans(mXi % w); - mat xtwx = xtw * mXi; - mat xtwy = xtw * mYi; + double xtwx = as_scalar(xtw * mXi); + double xtwy = as_scalar(xtw * mYi); try { - mat xtwx_inv = inv_sympd(xtwx); - vec beta = xtwx_inv * xtwy; + double xtwx_inv = 1.0 / xtwx; + double beta = xtwx_inv * xtwy; double res = mYi(i) - det(mXi(i) * beta); cv += res * res; } @@ -607,11 +607,11 @@ double GWRMultiscale::bandwidthSizeCriterionVarAICSerial(BandwidthWeight *bandwi vec d = mSpatialWeights[var].distance()->distance(i); vec w = bandwidthWeight->weight(d); mat xtw = trans(mXi % w); - mat xtwx = xtw * mXi; - mat xtwy = xtw * mYi; + double xtwx = as_scalar(xtw * mXi); + double xtwy = as_scalar(xtw * mYi); try { - mat xtwx_inv = inv_sympd(xtwx); + double xtwx_inv = 1.0 / xtwx; betas.col(i) = xtwx_inv * xtwy; mat ci = xtwx_inv * xtw; mat si = mXi(i) * ci;