From e13905b970705486d02cb9b88f349b993c342978 Mon Sep 17 00:00:00 2001 From: Yi Zhang Date: Wed, 9 Oct 2024 14:11:00 +0800 Subject: [PATCH] tmp update --- example/kde_ex.cpp | 43 +++++------ lib/algorithm/kde.cpp | 172 ++++++++++++------------------------------ lib/algorithm/kde.h | 68 +++++++++++++---- 3 files changed, 121 insertions(+), 162 deletions(-) diff --git a/example/kde_ex.cpp b/example/kde_ex.cpp index 5c29c36..4b9b6a2 100644 --- a/example/kde_ex.cpp +++ b/example/kde_ex.cpp @@ -31,7 +31,8 @@ using namespace gctl; int main(int argc, char const *argv[]) try -{ +{ +/* array x(201); x.sequence(-1.0, 0.01); kde k(0.02, x); @@ -50,7 +51,7 @@ int main(int argc, char const *argv[]) try array dm(201); array am(100000, 0.0); -/* + for (size_t i = 0; i < am.size(); i++) { k.get_gradient_at(i, a, dm); @@ -62,41 +63,41 @@ int main(int argc, char const *argv[]) try std::cout << a[i] << " " << am[i] << "\n"; } -*/ for (size_t i = 0; i < x.size(); i++) { std::cout << x[i] << " " << d[i] - g[i] << "\n"; } +*/ - /* - array a(100000), b(100000); - a.random(0, 0.2, RdNormal, 0); - b.random(1, 0.2, RdNormal, 0); - kde2d k(0.02, a, b); - - array x, y; - linespace(-1.0, 1.0, 201, x); - linespace(0.0, 2.0, 201, y); + array x(201), y(301); + x.sequence(-1.0, 0.01); + y.sequence(0.0, 0.01); + kde2d k(0.1, x, y); - gaussian_para2d g1(0, 1, 0.2, 0.2, 0); + array a(10000), b(10000); + a.random_float(0, 0.2, RdNormal, 0); + b.random_float(1.5, 0.3, RdNormal, 0); + + gaussian_para2d g1(0, 1.5, 0.2, 0.3, 0); + + array d(201*301); + //k.get_distribution(a, b, d); + + a[0] = 0; + b[0] = 1.5; + k.get_gradient_y_at(0, 0, a, b, d); double t, sum = 0; for (size_t i = 0; i < y.size(); i++) { for (size_t j = 0; j < x.size(); j++) { - t = k.get_density_at(x[j], y[i]); - sum += t*0.01*0.01; - std::cout << x[j] << " " << y[i] << " " - << gaussian_dist2d(x[j], y[i], g1) << " " - << t << "\n"; + //<< gaussian_dist2d(x[j], y[i], g1) << " " + << d[201*i + j] << "\n"; } } - - std::clog << "sum = " << sum << "\n"; - */ return 0; } catch(std::exception &e) diff --git a/lib/algorithm/kde.cpp b/lib/algorithm/kde.cpp index 1e89f01..d2ba75e 100644 --- a/lib/algorithm/kde.cpp +++ b/lib/algorithm/kde.cpp @@ -38,8 +38,8 @@ gctl::kde::kde(double h, const array &x) void gctl::kde::init(double h, const array &x) { - if (h <= 0) throw std::runtime_error("[GCTL Kernel Density Estimation] Invalid averaging width."); - if (x.size() < 2) throw std::runtime_error("[GCTL Kernel Density Estimation] Invalid sample size."); + if (h <= 0) throw std::runtime_error("[gctl::kde2d] Invalid averaging width."); + if (x.size() < 2) throw std::runtime_error("[gctl::kde2d] Invalid sample size."); h_ = h; x_ = x; @@ -193,167 +193,89 @@ gctl::kde2d::kde2d(double h, const array &x, const array &y) init(h, x, y); } -gctl::kde2d::kde2d(double h, const std::vector &x, const std::vector &y) -{ - init(h, x, y); -} - void gctl::kde2d::init(double h, const array &x, const array &y) { - if (h <= 0) throw std::runtime_error("[GCTL Kernel Density Estimation] Invalid averaging width."); - if (x.size() < 2 || y.size() < 2 || x.size() != y.size()) throw std::runtime_error("[GCTL Kernel Density Estimation] Invalid sample size."); + if (h <= 0) throw std::runtime_error("[gctl::kde2d] Invalid averaging width."); + if (x.size() < 2 || y.size() < 2) throw std::runtime_error("[gctl::kde2d] Invalid sample size."); h_ = h; + xs_ = x.size(); + ys_ = y.size(); x_ = x; y_ = y; return; } -void gctl::kde2d::init(double h, const std::vector &x, const std::vector &y) +void gctl::kde2d::get_distribution(const array &mx, + const array &my, array &dxy, kde_kernel_e k_type) { - if (h <= 0) throw std::runtime_error("[GCTL Kernel Density Estimation] Invalid averaging width."); - if (x.size() < 2 || y.size() < 2 || x.size() != y.size()) throw std::runtime_error("[GCTL Kernel Density Estimation] Invalid sample size."); + if (mx.size() != my.size()) throw std::runtime_error("[gctl::kde2d] Invalid evaluating size.");; + int ms = mx.size(); - h_ = h; - x_.import_vector(x); - y_.import_vector(y); - return; -} + dxy.resize(xs_*ys_); -double gctl::kde2d::get_density_at(double x, double y, kde_kernel_e k_type) -{ - double out = 0; - if (k_type == KDE_Gaussian) - { - for (size_t i = 0; i < x_.size(); i++) - { - out += gaussian_kernel((x - x_[i])/h_, (y - y_[i])/h_); - } - } - else throw std::runtime_error("[GCTL Kernel Density Estimation] Invalid kernel type."); - return out/(h_*h_*x_.size()); -} - -double gctl::kde2d::get_kernel_density_at(size_t k_id, double x, double y, kde_kernel_e k_type) -{ - if (k_id >= x_.size()) throw std::runtime_error("[gctl::kde2d::get_kernel_density_at(...)] Invalid kernel index."); - double out; - if (k_type == KDE_Gaussian) out = gaussian_kernel((x - x_[k_id])/h_, (y - y_[k_id])/h_); - else throw std::runtime_error("[gctl::kde2d::get_kernel_density_at(...)] Invalid kernel type."); - - return out/(h_*h_); -} - -void gctl::kde2d::get_gradient_at(double x, double y, double &gx, double &gy, kde_kernel_e k_type) -{ - double out_x = 0.0, out_y = 0.0; if (k_type == KDE_Gaussian) { - for (size_t i = 0; i < x_.size(); i++) + for (size_t i = 0; i < ys_; i++) { - out_x += ((x - x_[i])/h_)*gaussian_kernel((x - x_[i])/h_, (y - y_[i])/h_); - out_y += ((y - y_[i])/h_)*gaussian_kernel((x - x_[i])/h_, (y - y_[i])/h_); + for (size_t j = 0; j < xs_; j++) + { + out = 0.0; + for (size_t k = 0; k < ms; k++) + { + out += gaussian_kernel((x_[j] - mx[k])/h_, (y_[i] - my[k])/h_); + } + + dxy[i*xs_ + j] = out/(h_*h_*ms); + } } } - else throw std::runtime_error("[GCTL Kernel Density Estimation] Invalid kernel type."); - - gx = -1.0*out_x/(h_*h_*h_*x_.size()); - gy = -1.0*out_y/(h_*h_*h_*x_.size()); + else throw std::runtime_error("[gctl::kde2d] Invalid kernel type."); return; } -void gctl::kde2d::get_kernel_gradient_at(size_t k_id, double x, double y, double &gx, double &gy, kde_kernel_e k_type) +void gctl::kde2d::get_gradient_x_at(size_t mx_id, size_t my_id, + const array &mx, const array &my, + array &dmx, kde_kernel_e k_type) { - if (k_id >= x_.size()) throw std::runtime_error("[gctl::kde2d::get_kernel_gradient_at(...)] Invalid kernel index."); + if (mx.size() != my.size()) throw std::runtime_error("[gctl::kde2d] Invalid evaluating size.");; + int ms = mx.size(); + + dmx.resize(xs_*ys_); - double out_x, out_y; if (k_type == KDE_Gaussian) { - out_x = ((x - x_[k_id])/h_)*gaussian_kernel((x - x_[k_id])/h_, (y - y_[k_id])/h_); - out_y = ((y - y_[k_id])/h_)*gaussian_kernel((x - x_[k_id])/h_, (y - y_[k_id])/h_); - } - else throw std::runtime_error("[gctl::kde2d::get_kernel_gradient_at(...)] Invalid kernel type."); - - gx = -1.0*out_x/(h_*h_*h_); - gy = -1.0*out_y/(h_*h_*h_); - return; -} - -void gctl::kde2d::get_distribution(const array x, const array y, array &d, - array &dx, array &dy, kde_kernel_e k_type, kde_norm_e n_type, double norm) -{ - if (x.size() != y.size()) throw std::runtime_error("[GCTL Kernel Density Estimation] Invalid distribution size."); - if (norm < 0.0) throw std::runtime_error("[GCTL Kernel Density Estimation] Invalid normalization value."); - - size_t xnum = x.size(); - d.resize(xnum); - dx.resize(xnum); - dy.resize(xnum); - - double s = 0.0; - if (n_type == KDE_MAX2ONE) - { - for (size_t i = 0; i < xnum; i++) + for (size_t i = 0; i < ys_; i++) { - d[i] = get_density_at(x[i], y[i], k_type); - get_gradient_at(x[i], y[i], dx[i], dy[i], k_type); - - s = std::max(s, d[i]); + for (size_t j = 0; j < xs_; j++) + { + dmx[i*xs_ + j] = ((x_[j] - mx[mx_id])/h_)*gaussian_kernel((x_[j] - mx[mx_id])/h_, (y_[i] - my[my_id])/h_)/(h_*h_*h_*ms); + } } } - else if (n_type == KDE_SUM2ONE) - { - for (size_t i = 0; i < xnum; i++) - { - d[i] = get_density_at(x[i], y[i], k_type); - get_gradient_at(x[i], y[i], dx[i], dy[i], k_type); - - s += d[i]; - } - } - else - { - for (size_t i = 0; i < xnum; i++) - { - d[i] = get_density_at(x[i], y[i], k_type); - get_gradient_at(x[i], y[i], dx[i], dy[i], k_type); - } - - s = norm; - } - - for (size_t i = 0; i < xnum; i++) - { - d[i] /= s; - dx[i]/= s; - dy[i]/= s; - } return; } -void gctl::kde2d::save(double xmin, double xmax, double ymin, double ymax, int xnum, int ynum, std::string file) +void gctl::kde2d::get_gradient_y_at(size_t mx_id, size_t my_id, + const array &mx, const array &my, + array &dmy, kde_kernel_e k_type) { - std::string suffix_str = file.substr(file.find_last_of('.') + 1); - if (suffix_str != "nc") - { - throw std::runtime_error("[gctl::kde2d::save(...)] Invalid file extension type."); - } + if (mx.size() != my.size()) throw std::runtime_error("[gctl::kde2d] Invalid evaluating size.");; + int ms = mx.size(); - array dist(xnum*ynum, 0.0); - double dx = (xmax - xmin)/(xnum - 1); - double dy = (ymax - ymin)/(ynum - 1); + dmy.resize(xs_*ys_); - for (size_t i = 0; i < ynum; i++) + if (k_type == KDE_Gaussian) { - for (size_t j = 0; j < xnum; j++) + for (size_t i = 0; i < ys_; i++) { - dist[j + i*xnum] = get_density_at(xmin + dx*j, ymin + dy*i); + for (size_t j = 0; j < xs_; j++) + { + dmy[i*xs_ + j] = ((y_[i] - my[my_id])/h_)*gaussian_kernel((x_[j] - mx[mx_id])/h_, (y_[i] - my[my_id])/h_)/(h_*h_*h_*ms); + } } } - - if (suffix_str == "nc") save_netcdf_grid(file, dist, xnum, ynum, xmin, dx, ymin, dy, "x", "y", "probability density"); - return; } diff --git a/lib/algorithm/kde.h b/lib/algorithm/kde.h index 90b0fdf..8ac8918 100644 --- a/lib/algorithm/kde.h +++ b/lib/algorithm/kde.h @@ -107,31 +107,67 @@ namespace gctl { public: kde2d(); - kde2d(double h, const array &x, const array &y); - kde2d(double h, const std::vector &x, const std::vector &y); virtual ~kde2d(); + /** + * @brief Construct a new kde2d object + * + * @param h + * @param x + * @param y + */ + kde2d(double h, const array &x, const array &y); + + /** + * @brief + * + * @param h + * @param x + * @param y + */ void init(double h, const array &x, const array &y); - void init(double h, const std::vector &x, const std::vector &y); - double get_density_at(double x, double y, kde_kernel_e k_type = KDE_Gaussian); /** - * @brief Get the probability density of a single kernel. Note the value is not normalized by the kernel number. + * @brief Get the distribution object * - * @param k_id kernel index - * @param x inquiring location x - * @param y inquiring location y - * @param k_type kernel type - * @return kernel value + * @param mx + * @param my + * @param dxy + * @param k_type */ - double get_kernel_density_at(size_t k_id, double x, double y, kde_kernel_e k_type = KDE_Gaussian); - void get_gradient_at(double x, double y, double &gx, double &gy, kde_kernel_e k_type = KDE_Gaussian); - void get_kernel_gradient_at(size_t k_id, double x, double y, double &gx, double &gy, kde_kernel_e k_type = KDE_Gaussian); - void get_distribution(const array x, const array y, array &d, array &dx, - array &dy, kde_kernel_e k_type = KDE_Gaussian, kde_norm_e n_type = KDE_MAX2ONE, double norm = 1.0); - void save(double xmin, double xmax, double ymin, double ymax, int xnum, int ynum, std::string file); + void get_distribution(const array &mx, const array &my, + array &dxy, kde_kernel_e k_type = KDE_Gaussian); + + /** + * @brief Get the gradient at object + * + * @param mx_id + * @param my_id + * @param mx + * @param my + * @param dmx + * @param k_type + */ + void get_gradient_x_at(size_t mx_id, size_t my_id, + const array &mx, const array &my, + array &dmx, kde_kernel_e k_type = KDE_Gaussian); + + /** + * @brief Get the gradient y at object + * + * @param mx_id + * @param my_id + * @param mx + * @param my + * @param dmy + * @param k_type + */ + void get_gradient_y_at(size_t mx_id, size_t my_id, + const array &mx, const array &my, + array &dmy, kde_kernel_e k_type = KDE_Gaussian); private: + size_t xs_, ys_; double h_; array x_, y_;