tmp update

example/kde_ex.cpp

@@ -32,6 +32,7 @@ using namespace gctl;
 
 int main(int argc, char const *argv[]) try
 {
+/*
     array<double> 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<double> dm(201);
     array<double> 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<double> x(201), y(301);
-    array<double> a(100000), b(100000);
+    x.sequence(-1.0, 0.01);
-    a.random(0, 0.2, RdNormal, 0);
+    y.sequence(0.0, 0.01);
-    b.random(1, 0.2, RdNormal, 0);
+    kde2d k(0.1, x, y);
-    kde2d k(0.02, a, b);
 
-    array<double> x, y;
+    array<double> a(10000), b(10000);
-    linespace(-1.0, 1.0, 201, x);
+    a.random_float(0, 0.2, RdNormal, 0);
-    linespace(0.0, 2.0, 201, y);
+    b.random_float(1.5, 0.3, RdNormal, 0);
 
-    gaussian_para2d g1(0, 1, 0.2, 0.2, 0);
+    gaussian_para2d g1(0, 1.5, 0.2, 0.3, 0);
+
+    array<double> 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) << " "
+                //<< gaussian_dist2d(x[j], y[i], g1) << " "
-                << t << "\n";
+                << d[201*i + j] << "\n";
         }
     }
-
-    std::clog << "sum = " << sum << "\n";
-    */
     return 0;
 }
 catch(std::exception &e)

lib/algorithm/kde.cpp

@@ -38,8 +38,8 @@ gctl::kde::kde(double h, const array<double> &x)
 
 void gctl::kde::init(double h, const array<double> &x)
 {
-    if (h <= 0) throw std::runtime_error("[GCTL Kernel Density Estimation] Invalid averaging width.");
+    if (h <= 0) throw std::runtime_error("[gctl::kde2d] Invalid averaging width.");
-    if (x.size() < 2) throw std::runtime_error("[GCTL Kernel Density Estimation] Invalid sample size.");
+    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<double> &x, const array<double> &y)
     init(h, x, y);
 }
 
-gctl::kde2d::kde2d(double h, const std::vector<double> &x, const std::vector<double> &y)
-{
-    init(h, x, y);
-}
-
 void gctl::kde2d::init(double h, const array<double> &x, const array<double> &y)
 {
-    if (h <= 0) throw std::runtime_error("[GCTL Kernel Density Estimation] Invalid averaging width.");
+    if (h <= 0) throw std::runtime_error("[gctl::kde2d] 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 (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<double> &x, const std::vector<double> &y)
+void gctl::kde2d::get_distribution(const array<double> &mx, 
+    const array<double> &my, array<double> &dxy, kde_kernel_e k_type)
 {
-    if (h <= 0) throw std::runtime_error("[GCTL Kernel Density Estimation] Invalid averaging width.");
+    if (mx.size() != my.size()) throw std::runtime_error("[gctl::kde2d] Invalid evaluating size.");;
-    if (x.size() < 2 || y.size() < 2 || x.size() != y.size()) throw std::runtime_error("[GCTL Kernel Density Estimation] Invalid sample size.");
+    int ms = mx.size();
 
-    h_ = h;
+    dxy.resize(xs_*ys_);
-    x_.import_vector(x);
-    y_.import_vector(y);
-    return;
-}
-
-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_);
+            for (size_t j = 0; j < xs_; j++)
-            out_y += ((y - y_[i])/h_)*gaussian_kernel((x - x_[i])/h_, (y - y_[i])/h_);
+            {
+                out = 0.0;
+                for (size_t k = 0; k < ms; k++)
+                {
+                    out += gaussian_kernel((x_[j] - mx[k])/h_, (y_[i] - my[k])/h_);
                 }
-    }
-    else throw std::runtime_error("[GCTL Kernel Density Estimation] Invalid kernel type.");
                 
-    gx = -1.0*out_x/(h_*h_*h_*x_.size());
+                dxy[i*xs_ + j] = out/(h_*h_*ms);
-    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<double> &mx, const array<double> &my, 
+    array<double> &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_);
+        for (size_t i = 0; i < ys_; i++)
-        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<double> x, const array<double> y, array<double> &d, 
-    array<double> &dx, array<double> &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 j = 0; j < xs_; j++)
             {
-            d[i] = get_density_at(x[i], y[i], k_type);
+                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);
-            get_gradient_at(x[i], y[i], dx[i], dy[i], k_type);
-
-            s = std::max(s, d[i]);
             }
         }
-    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<double> &mx, const array<double> &my, 
+    array<double> &dmy, kde_kernel_e k_type)
 {
-    std::string suffix_str = file.substr(file.find_last_of('.') + 1);
+    if (mx.size() != my.size()) throw std::runtime_error("[gctl::kde2d] Invalid evaluating size.");;
-    if (suffix_str != "nc")
+    int ms = mx.size();
-    {
-        throw std::runtime_error("[gctl::kde2d::save(...)] Invalid file extension type.");
-    }
 
-    array<double> dist(xnum*ynum, 0.0);
+    dmy.resize(xs_*ys_);
-    double dx = (xmax - xmin)/(xnum - 1);
-    double dy = (ymax - ymin)/(ynum - 1);
 
-    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;
 }
 

lib/algorithm/kde.h

@@ -107,31 +107,67 @@ namespace gctl
     {
     public:
         kde2d();
-        kde2d(double h, const array<double> &x, const array<double> &y);
-        kde2d(double h, const std::vector<double> &x, const std::vector<double> &y);
         virtual ~kde2d();
 
-        void init(double h, const array<double> &x, const array<double> &y);
+        /**
-        void init(double h, const std::vector<double> &x, const std::vector<double> &y);
+         * @brief Construct a new kde2d object
-        double get_density_at(double x, double y, kde_kernel_e k_type = KDE_Gaussian);
+         * 
+         * @param h 
+         * @param x 
+         * @param y 
+         */
+        kde2d(double h, const array<double> &x, const array<double> &y);
         
         /**
-         * @brief Get the probability density of a single kernel. Note the value is not normalized by the kernel number.
+         * @brief 
          * 
-         * @param k_id kernel index
+         * @param h 
-         * @param x inquiring location x
+         * @param x 
-         * @param y inquiring location y
+         * @param y 
-         * @param k_type kernel type
-         * @return kernel value
          */
-        double get_kernel_density_at(size_t k_id, double x, double y, kde_kernel_e k_type = KDE_Gaussian);
+        void init(double h, const array<double> &x, const array<double> &y);
-        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<double> x, const array<double> y, array<double> &d, array<double> &dx, 
+         * @brief Get the distribution object
-            array<double> &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);
+         * @param mx 
+         * @param my 
+         * @param dxy 
+         * @param k_type 
+         */
+        void get_distribution(const array<double> &mx, const array<double> &my, 
+            array<double> &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<double> &mx, const array<double> &my, 
+            array<double> &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<double> &mx, const array<double> &my, 
+            array<double> &dmy, kde_kernel_e k_type = KDE_Gaussian);
     
     private:
+        size_t xs_, ys_;
         double h_;
         array<double> x_, y_;