diff --git a/docs/HTML/ARIMAVisitor.html b/docs/HTML/ARIMAVisitor.html index c7f130d6..5f19511e 100644 --- a/docs/HTML/ARIMAVisitor.html +++ b/docs/HTML/ARIMAVisitor.html @@ -94,13 +94,13 @@ std::cout << "\nTesting ARIMAVisitor{ } ..." << std::endl; - std::vector<unsigned long> idxvec = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13 }; - std::vector<double> col1 = { 266.0, 145.9, 183.1, 119.3, 180.3, 168.5, 231.8, 224.5, 192.8, 122.9, 336.5, 185.9, 194.3 }; - std::vector<double> oscil = { 1.5, 1.8, 1.62, 1.78, 1.5, 1.68, 1.6, 1.8, 1.71, 1.9, 1.78, 1.84, 1.69 }; - std::vector<double> constant = { 10.56, 10.56, 10.56, 10.56, 10.56, 10.56, 10.56, 10.56, 10.56, 10.56, 10.56, 10.56, 10.56 }; + std::vector<unsigned long> idxvec = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13 }; + std::vector<double> col1 = { 266.0, 145.9, 183.1, 119.3, 180.3, 168.5, 231.8, 224.5, 192.8, 122.9, 336.5, 185.9, 194.3 }; + std::vector<double> oscil = { 1.5, 1.8, 1.62, 1.78, 1.5, 1.68, 1.6, 1.8, 1.71, 1.9, 1.78, 1.84, 1.69 }; + std::vector<double> constant = { 10.56, 10.56, 10.56, 10.56, 10.56, 10.56, 10.56, 10.56, 10.56, 10.56, 10.56, 10.56, 10.56 }; std::vector<double> increasing = { 10.56, 10.68, 10.78, 10.90, 11.01, 11.45, 11.99, 12.01, 12.21, 12.35, 12.67, 13.89, 13.01 }; std::vector<double> decreasing = { 10.56, 10.30, 10.12, 10.01, 9.80, 9.74, 9.41, 9.03, 9.0, 8.20, 8.01, 7.9, 7.55 }; - MyDataFrame df; + ULDataFrame df; df.load_data(std::move(idxvec), std::make_pair("col1", col1), @@ -116,18 +116,18 @@ const auto result1 = ari.get_result(); assert(result1.size() == 3); - assert(std::fabs(result1[0] - 247.175) < 0.001); - assert(std::fabs(result1[1] - 197.294) < 0.001); - assert(std::fabs(result1[2] - 220.021) < 0.001); + assert(std::fabs(result1[0] - 245.745) < 0.001); + assert(std::fabs(result1[1] - 197.122) < 0.001); + assert(std::fabs(result1[2] - 219.314) < 0.001); df.single_act_visit<double>("oscil", ari); const auto result2 = ari.get_result(); assert(result2.size() == 3); - assert(std::fabs(result2[0] - 1.77088) < 0.00001); - assert(std::fabs(result2[1] - 1.67015) < 0.00001); - assert(std::fabs(result2[2] - 1.74417) < 0.00001); + assert(std::fabs(result2[0] - 1.76961) < 0.00001); + assert(std::fabs(result2[1] - 1.66931) < 0.00001); + assert(std::fabs(result2[2] - 1.74256) < 0.00001); try { df.single_act_visit<double>("constant", ari); @@ -141,18 +141,18 @@ const auto result3 = ari.get_result(); assert(result3.size() == 3); - assert(std::fabs(result3[0] - 14.3335) < 0.0001); - assert(std::fabs(result3[1] - 13.09) < 0.0001); - assert(std::fabs(result3[2] - 14.6469) < 0.0001); + assert(std::fabs(result3[0] - 14.1842) < 0.0001); + assert(std::fabs(result3[1] - 13.0033) < 0.0001); + assert(std::fabs(result3[2] - 14.4171) < 0.0001); df.single_act_visit<double>("decreasing", ari); const auto result4 = ari.get_result(); assert(result4.size() == 3); - assert(std::fabs(result4[0] - 7.42058) < 0.00001); - assert(std::fabs(result4[1] - 7.21897) < 0.00001); - assert(std::fabs(result4[2] - 7.11158) < 0.00001); + assert(std::fabs(result4[0] - 7.40899) < 0.00001); + assert(std::fabs(result4[1] - 7.2049) < 0.0001); + assert(std::fabs(result4[2] - 7.09123) < 0.00001); // Now some real data // @@ -166,17 +166,17 @@ ::exit(-1); } - ARIMAVisitor<double> ari2 { 4, 3 }; + ARIMAVisitor<double, std::string> ari2 { 4 }; df2.single_act_visit<double>("IBM_Close", ari2); const auto result5 = ari2.get_result(); assert(result5.size() == 4); - assert(std::fabs(result5[0] - 111.658) < 0.001); - assert(std::fabs(result5[1] - 111.669) < 0.001); - assert(std::fabs(result5[2] - 111.649) < 0.001); - assert(std::fabs(result5[3] - 111.649) < 0.001); + assert(std::fabs(result5[0] - 111.63) < 0.001); + assert(std::fabs(result5[1] - 111.658) < 0.001); + assert(std::fabs(result5[2] - 111.657) < 0.001); + assert(std::fabs(result5[3] - 111.658) < 0.001); } diff --git a/docs/HTML/AnomalyDetectByFFTVisitor.html b/docs/HTML/AnomalyDetectByFFTVisitor.html index 46e3db31..ce163150 100644 --- a/docs/HTML/AnomalyDetectByFFTVisitor.html +++ b/docs/HTML/AnomalyDetectByFFTVisitor.html @@ -157,7 +157,7 @@ ibm.single_act_visit<double>("IBM_Close", anomaly5); assert((anomaly5.get_result() == result2)); - and_fft_v<double, std::string> anomaly6(1000, 250.0, normalization_type::z_score); + and_fft_v<double, std::string> anomaly6(1000, 10.0, normalization_type::z_score); const std::vector<std::size_t> result3 = { 502, 1001, 2002 }; ibm.single_act_visit<double>("IBM_Close", anomaly6); diff --git a/docs/HTML/AnomalyDetectByIsoForestVisitor.html b/docs/HTML/AnomalyDetectByIsoForestVisitor.html index 903b11dc..e19865a8 100644 --- a/docs/HTML/AnomalyDetectByIsoForestVisitor.html +++ b/docs/HTML/AnomalyDetectByIsoForestVisitor.html @@ -58,7 +58,7 @@ Isolation Forest Anomaly Detection Visitor
- Scores every point in the input column using an Isolation Forest, then flags as anomalous all points whose score exceeds threshold. Isolation Forest (Liu, Ting & Zhou 2008) isolates observations by recursively partitioning the data with random splits. Anomalies requirefewer splits to isolate than normal points, so they have shorter average path lengths across the ensemble of trees. The anomaly score is:
+ Scores every point in the input column using an Isolation Forest, then flags as anomalous all points whose score exceeds threshold. Isolation Forest (Liu, Ting & Zhou 2008) isolates observations by recursively partitioning the data with random splits. Anomalies require fewer splits to isolate than normal points, so they have shorter average path lengths across the ensemble of trees. The anomaly score is:
 s(x, n) = 2(−E[h(x)] / c(n))
@@ -69,7 +69,7 @@
        = 2 * (ln(n − 1) + γ) − 2 * (n − 1) / n (Euler-Mascheroni: γ ≈ 0.5772)
 
 Score ∈ (0, 1]:
-  s -> 1: point is isolated very quickly → anomaly
+  s -> 1: point is isolated very quickly -> anomaly
   s -> 0: point requires many splits -> normal
   s = 0.5: point is indistinguishable from random noise
diff --git a/docs/HTML/BetaVisitor.html b/docs/HTML/BetaVisitor.html index 6ae379d3..701e892a 100644 --- a/docs/HTML/BetaVisitor.html +++ b/docs/HTML/BetaVisitor.html @@ -86,9 +86,7 @@ -
static void test_beta()  {
-
-    using MyDataFrame = StdDataFrame<unsigned long>;
+
    using MyDataFrame = StdDataFrame<unsigned long>;
 
     std::cout << "\nTesting Beta ..." << std::endl;
 
@@ -145,17 +143,19 @@
 
     const auto  &md_result { md_beta.get_result() };
 
-    assert(md_result.rows() == 1);
+    assert(md_result.rows() == 2);
     assert(md_result.cols() == 2);
     assert(fabs(md_result(0, 0) - 0.428571) < 0.000001);
     assert(fabs(md_result(0, 1) - -0.142857) < 0.000001);
+    assert(fabs(md_result(1, 0) - 0.142857) < 0.000001);
+    assert(fabs(md_result(1, 1) - 0.285714) < 0.000001);
 
     const auto  data_mean { md_beta.get_data_mean() };
     const auto  benchmark_mean { md_beta.get_benchmark_mean() };
 
     assert((data_mean == std::vector<double>{ 1.0, 0.75 }));
     assert((benchmark_mean == std::vector<double>{ 2.75, 1.25 }));
-}
+}
 

 
   
C++ DataFrame
     // Set a different distance function than default
     //
-    dbscan.set_dist_func([](const double &x, const double &y)  {
-                             return (std::fabs(x - y));
-                         });
+    dbscan.set_dist_func([](const double &x, const double &y)  { return (std::fabs(x - y)); });
 
     view.single_act_visit<double>("IBM_Close", dbscan);
 
@@ -131,16 +129,12 @@
     assert(dbscan.get_noisey_idxs()[0] == 1564);
     assert(dbscan.get_noisey_idxs()[1] == 1565);
 
-    assert(dbscan.get_result().size() == 19);
-    assert(dbscan.get_result()[0].size() == 11);
-    assert(dbscan.get_result()[4].size() == 31);
-    assert(dbscan.get_result()[10].size() == 294);
-    assert(dbscan.get_result()[14].size() == 82);
-    assert(dbscan.get_result()[18].size() == 10);
-    assert(dbscan.get_result()[0][6] == 185.679993);
-    assert(dbscan.get_result()[4][18] == 167.330002);
-    assert(dbscan.get_result()[10][135] == 145.160004);
-    assert(dbscan.get_result()[18][3] == 103.550003);
+    assert(dbscan.get_result().size() == 1);
+    assert(dbscan.get_result()[0].size() == 1719);
+    assert(std::fabs(dbscan.get_result()[0][6] - 187.26) < 0.001);
+    assert(std::fabs(dbscan.get_result()[0][18] - 176.4) < 0.001);
+    assert(std::fabs(dbscan.get_result()[0][135] - 192.49) < 0.001);
+    assert(std::fabs(dbscan.get_result()[0][1718] - 111.66) < 0.001);
 
     // Now multidimensional data
     //
@@ -152,8 +146,7 @@
 
     using col_t = std::array<double, 3>;
 
-    auto    rand_vec = gen_uniform_real_dist<double>(df.get_index().size() * 3, p);
-
+    auto                rand_vec = gen_uniform_real_dist<double>(df.get_index().size() * 3, p);
     std::vector<col_t>  multi_dimen_col(df.get_index().size());
 
     for (std::size_t i { 0 }, j { 0 }; j < rand_vec.size(); ++i)  {
@@ -169,16 +162,16 @@
 
     const auto  &md_clusters = md_dbscan.get_result();
 
-    assert(md_clusters.size() == 102); // Number of clusters
+    assert(md_clusters.size() == 30); // Number of clusters
 
-    assert(md_clusters[0].size() == 14);
-    assert(std::fabs(md_clusters[0][6][1] - -19.9438) < 0.0001);
+    assert(md_clusters[0].size() == 27);
+    assert(std::fabs(md_clusters[0][6][1] - 4.24441) < 0.00001);
 
-    assert(md_clusters[58].size() == 12);
-    assert(std::fabs(md_clusters[58][3][0] - -6.41034) < 0.00001);
+    assert(md_clusters[26].size() == 10);
+    assert(std::fabs(md_clusters[26][3][0] - 14.87) < 0.001);
 
-    assert(md_clusters[101].size() == 10);
-    assert(std::fabs(md_clusters[101][9][2] - -5.92195) < 0.00001);
+    assert(md_clusters[8].size() == 21);
+    assert(std::fabs(md_clusters[8][9][2] - -0.348493) < 0.000001);
 }
diff --git a/docs/HTML/DataFrame.html b/docs/HTML/DataFrame.html index 0c9c60db..f8f93194 100644 --- a/docs/HTML/DataFrame.html +++ b/docs/HTML/DataFrame.html @@ -1150,10 +1150,18 @@

API Reference with code samples &# ModeVisitor{} + + NKthValueVisitor{} + + NonZeroRangeVisitor{} + + NQuantileVisitor{} + + QuantileVisitor{} diff --git a/docs/HTML/EntropyVisitor.html b/docs/HTML/EntropyVisitor.html index ad515737..eb053bb4 100644 --- a/docs/HTML/EntropyVisitor.html +++ b/docs/HTML/EntropyVisitor.html @@ -161,12 +161,10 @@ StlVecType<unsigned long> idx = { 123450, 123451, 123452, 123453, 123454, 123455, 123456, 123457, 123458, 123459, 123460, 123461, 123462, 123466, 123467, 123468, - 123469, 123470, 123471, 123472, 123473, 22, 23, 24, 25, 26, 27, 28 - }; + 123469, 123470, 123471, 123472, 123473, 22, 23, 24, 25, 26, 27, 28 }; StlVecType<double> close = { 1.80, 2.80, 1.90, 14.00, 1.10, 6.00, 13.00, 8.00, 9.00, 2.80, 1.90, 4.30, 20.00, 1.85, 3.00, 34.00, 67.00, 23.00, 87.00, 9.00, 45.00, - 1.00, 11.00, 456.00, 34.00, 7.00, 7778.00, 5.00 - }; + 1.00, 11.00, 456.00, 34.00, 7.00, 7778.00, 5.00 }; MyDataFrame df; df.load_data(std::move(idx), std::make_pair("close", close)); @@ -178,14 +176,14 @@ assert(e_v.get_result().size() == 28); assert(std::isnan(e_v.get_result()[0])); assert(std::isnan(e_v.get_result()[3])); - assert(std::abs(e_v.get_result()[4] - 2.18974) < 0.00001); - assert(std::abs(e_v.get_result()[6] - 1.98477) < 0.00001); - assert(std::abs(e_v.get_result()[10] - 1.7154) < 0.0001); - assert(std::abs(e_v.get_result()[23] - 0.596666) < 0.00001); - assert(std::abs(e_v.get_result()[21] - 0.822228) < 0.00001); - assert(std::abs(e_v.get_result()[18] - 1.49397) < 0.0001); - assert(std::abs(e_v.get_result()[26] - 0.08568) < 0.0001); - assert(std::abs(e_v.get_result()[27] - 0.00646) < 0.0001); + assert(std::isnan(e_v.get_result()[0])); + assert(std::isnan(e_v.get_result()[7])); + assert(std::abs(e_v.get_result()[10] - 1.98477) < 0.00001); + assert(std::abs(e_v.get_result()[23] - 1.13643) < 0.00001); + assert(std::abs(e_v.get_result()[21] - 1.66467) < 0.00001); + assert(std::abs(e_v.get_result()[18] - 2.26252) < 0.0001); + assert(std::abs(e_v.get_result()[26] - 0.863265) < 0.000001); + assert(std::abs(e_v.get_result()[27] - 0.596666) < 0.000001); // Now multidimensional data // @@ -228,19 +226,18 @@ assert(std::isnan(ary_result[0][2])); assert(std::isnan(ary_result[2][2])); assert(std::isnan(ary_result[2][2])); - assert(std::abs(ary_result[3][0] - 1.88598) < 0.00001); - assert(std::abs(ary_result[3][1] - 1.76876) < 0.00001); - assert(std::abs(ary_result[6][1] - 1.95814) < 0.00001); - assert(std::abs(ary_result[6][2] - 1.84599) < 0.00001); + assert(std::isnan(ary_result[3][0])); + assert(std::isnan(ary_result[3][1])); + assert(std::abs(ary_result[6][1] - 1.76876) < 0.00001); + assert(std::abs(ary_result[6][2] - 1.91606) < 0.00001); assert(std::isnan(vec_result[0][0])); assert(std::isnan(vec_result[0][2])); assert(std::isnan(vec_result[2][2])); - assert(std::isnan(ary_result[2][2])); - assert(std::abs(vec_result[3][0] - 1.88598) < 0.00001); - assert(std::abs(vec_result[3][1] - 1.76876) < 0.00001); - assert(std::abs(vec_result[6][1] - 1.95814) < 0.00001); - assert(std::abs(vec_result[6][2] - 1.84599) < 0.00001); + assert(std::isnan(vec_result[3][0])); + assert(std::isnan(vec_result[3][1])); + assert(std::abs(vec_result[6][1] - 1.76876) < 0.00001); + assert(std::abs(vec_result[6][2] - 1.91606) < 0.00001); } // ----------------------------------------------------------------------------- diff --git a/docs/HTML/HWESForecastVisitor.html b/docs/HTML/HWESForecastVisitor.html index 3730ce3c..d62b33a8 100644 --- a/docs/HTML/HWESForecastVisitor.html +++ b/docs/HTML/HWESForecastVisitor.html @@ -115,7 +115,8 @@ std::cout << "\nTesting HWESForecastVisitor{ } ..." << std::endl; std::vector<unsigned long> idxvec = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13 }; - std::vector<double> col1 = { 266.0, 145.9, 183.1, 119.3, 180.3, 168.5, 231.8, 224.5, 192.8, 122.9, 336.5, 185.9, 194.3 }; + std::vector<double> col1 = { + 266.0, 145.9, 183.1, 119.3, 180.3, 168.5, 231.8, 224.5, 192.8, 122.9, 336.5, 185.9, 194.3 }; std::vector<double> oscil = { 1.5, 1.8, 1.62, 1.78, 1.5, 1.68, 1.6, 1.8, 1.71, 1.9, 1.78, 1.84, 1.69 }; std::vector<double> constant = { 10.56, 10.56, 10.56, 10.56, 10.56, 10.56, 10.56, 10.56, 10.56, 10.56, 10.56, 10.56, 10.56 }; std::vector<double> increasing = { 10.56, 10.68, 10.78, 10.90, 11.01, 11.45, 11.99, 12.01, 12.21, 12.35, 12.67, 13.89, 13.01 }; @@ -147,9 +148,9 @@ const auto result2 = hwes2.get_result(); assert(result2.size() == 3); - assert(std::fabs(result2[0] - 1.73499) < 0.00001); - assert(std::fabs(result2[1] - 1.9216) < 0.00001); - assert(std::fabs(result2[2] - 1.76383) < 0.00001); + assert(std::fabs(result2[0] - 1.90718) < 0.00001); + assert(std::fabs(result2[1] - 1.74941) < 0.00001); + assert(std::fabs(result2[2] - 1.93602) < 0.00001); df.single_act_visit<double>("constant", hwes); @@ -222,7 +223,9 @@ // so the seasonal path is exercised meaningfully. // // - const std::vector<double> pattern { 3, -3, 2, -2, 1, -1, 3, -3, 2, -2, 1, -1, }; + const std::vector<double> pattern { + 3, -3, 2, -2, 1, -1, 3, -3, 2, -2, 1, -1, + }; for (size_t t { 0 }; t < n; ++t) for (size_t d { 0 }; d < dim; ++d) @@ -255,11 +258,11 @@ // dim 0 is trending up — each forecast step should be increasing // for (std::size_t i { 1 }; i < vec_hwes1.get_result().size(); ++i) - assert(vec_hwes1.get_result()[i][0] > vec_hwes1.get_result()[i - 1][0]); + assert(vec_hwes1.get_result()[i][0] > vec_hwes1.get_result()[i-1][0]); // dim 1 is trending down — each step should be decreasing // for (std::size_t i { 1 }; i < ary_hwes1.get_result().size(); ++i) - assert(ary_hwes1.get_result()[i][1] < ary_hwes1.get_result()[i - 1][1]); + assert(ary_hwes1.get_result()[i][1] < ary_hwes1.get_result()[i-1][1]); // With seasons 12 periods // diff --git a/docs/HTML/KMeansVisitor.html b/docs/HTML/KMeansVisitor.html index cfaf7543..fed32481 100644 --- a/docs/HTML/KMeansVisitor.html +++ b/docs/HTML/KMeansVisitor.html @@ -88,19 +88,19 @@ Point() = default; Point(double xx, double yy) : x(xx), y(yy) { } Point(const Point &) = default; - Point &operator = (const Point &) = default; + Point &operator = (const Point &) = default; - friend Point operator + (const Point &lhs, const Point &rhs) { + friend Point operator + (const Point &lhs, const Point &rhs) { return (Point(lhs.x + rhs.x, lhs.y + rhs.y)); } - friend Point operator / (const Point &lhs, double rhs) { + friend Point operator / (const Point &lhs, double rhs) { return (Point(lhs.x / rhs, lhs.y / rhs)); } - template<typename S> - friend S &operator << (S &s, const Point &rhs) { + template<typename S> + friend S &operator << (S &s, const Point &rhs) { return (s << rhs.x << ", " << rhs.y); } @@ -108,104 +108,57 @@ static double point_distance(const Point &lhs, const Point &rhs) { - return ((lhs.x - rhs.x) * (lhs.x - rhs.x) + (lhs.y - rhs.y) * (lhs.y - rhs.y)); + return ((lhs.x - rhs.x) * (lhs.x - rhs.x) + + (lhs.y - rhs.y) * (lhs.y - rhs.y)); } // ------------------------------------- static void test_k_means() { - std::cout << "\nTesting k-means visitor ..." << std::endl; + std::cout << "\nTesting k-means visitor ..." << std::endl; const size_t item_cnt = 1024; MyDataFrame df; - RandGenParams<double> p; + RandGenParams<double> p; p.mean = 1.0; // Default - p.std = 0.005; + p.std = 0.005; p.seed = 10; - df.load_data(MyDataFrame::gen_sequence_index(0, item_cnt, 1), std::make_pair("col1", gen_lognormal_dist<double, 128>(item_cnt, p))); + df.load_data(MyDataFrame::gen_sequence_index(0, item_cnt, 1), std::make_pair("col1", gen_lognormal_dist<double, 128>(item_cnt, p))); - KMeansVisitor<5, double, unsigned long, 128> km_visitor(1000, true, - [](const double &x, const double &y) { + KMeansVisitor<5, double, unsigned long, 128> km_visitor(1000, true, + [](const double &x, const double &y) { return ((x - y) * (x - y)); - }, + }, 10); - df.single_act_visit<double>("col1", km_visitor); - std::cout << "Means of clusters are: "; + df.single_act_visit<double>("col1", km_visitor); + std::cout << "Means of clusters are: "; for (const auto citer : km_visitor.get_result()) - std::cout << citer << ", "; - std::cout << std::endl; - - // Using the calculated means, separate the given column into clusters - // const auto &clusters = km_visitor.get_clusters(); - // bool found = false; - - // for (auto iter : clusters) { - // if (::fabs(iter[0] - 1.89348) < 0.00001) { - // if (::fabs(iter[6] - 1.44231) < 0.00001) { - // found = true; - // break; - // } - // } - // } - // assert(found); - // found = false; - // for (auto iter : clusters) { - // if (::fabs(iter[0] - 0.593126) < 0.00001) { - // if (::fabs(iter[2] - 0.950026) < 0.00001) { - // found = true; - // break; - // } - // } - // } - // assert(found); - // found = false; - // for (auto iter : clusters) { - // if (::fabs(iter[0] - 14.2245) < 0.0001) { - // found = true; - // break; - // } - // } - // assert(found); - // found = false; - // for (auto iter : clusters) { - // if (::fabs(iter[0] - 6.90427) < 0.00001) { - // found = true; - // break; - // } - // } - // assert(found); - // found = false; - // for (auto iter : clusters) { - // if (::fabs(iter[0] - 3.8146) < 0.00001) { - // found = true; - // break; - // } - // } - // assert(found); + std::cout << citer << ", "; + std::cout << std::endl; // Now try with Points // p.seed = 200; - auto x_vec = gen_lognormal_dist<double, 128>(item_cnt, p); + auto x_vec = gen_lognormal_dist<double, 128>(item_cnt, p); p.seed = 4356; - auto y_vec = gen_lognormal_dist<double, 128>(item_cnt, p); - StlVecType<Point> points; + auto y_vec = gen_lognormal_dist<double, 128>(item_cnt, p); + StlVecType<Point> points; points.reserve(item_cnt); for (size_t i = 0; i < item_cnt; ++i) points.push_back(Point(x_vec[i], y_vec[i])); - df.load_column<Point>("point_col", std::move(points)); + df.load_column<Point>("point_col", std::move(points)); - KMeansVisitor<5, Point, unsigned long, 128> km_visitor2(1000, true, point_distance, 10); + KMeansVisitor<5, Point, unsigned long, 128> km_visitor2(1000, true, point_distance, 10); - df.single_act_visit<Point>("point_col", km_visitor2); + df.single_act_visit<Point>("point_col", km_visitor2); // Using the calculated means, separate the given column into clusters // @@ -213,104 +166,49 @@ for (auto iter : clusters2) { for (auto iter2 : iter) { - std::cout << iter2.x << " | " << iter2.y << ", "; + std::cout << iter2.x << " | " << iter2.y << ", "; } - std::cout << "\n\n" << std::endl; + std::cout << "\n\n" << std::endl; } - // found = false; - // for (auto iter : clusters2) { - // if (::fabs(iter[0].x - 18.9556) < 0.1 && - // ::fabs(iter[0].y - 2.17537) < 0.1) { - // if (::fabs(iter[6].x - 16.7309) < 0.1 && - // ::fabs(iter[6].y - 0.872376) < 0.1) { - // found = true; - // break; - // } - // } - // } - // assert(found); - - // found = false; - // for (auto iter : clusters2) { - // if (::fabs(iter[0].x - 0.943977) < 0.1 && - // ::fabs(iter[0].y - 0.910989) < 0.1) { - // if (::fabs(iter[2].x - 0.30509) < 0.1 && - // ::fabs(iter[2].y - 1.69017) < 0.1) { - // found = true; - // break; - // } - // } - // } - // assert(found); - // found = false; - // for (auto iter : clusters2) { - // if (::fabs(iter[0].x - 4.31973) < 0.1 && - // ::fabs(iter[0].y - 1.24214) < 0.1) { - // if (::fabs(iter[3].x - 4.68381) < 0.1 && - // ::fabs(iter[3].y - 0.453632) < 0.1) { - // found = true; - // break; - // } - // } - // } - // assert(found); - // found = false; - // for (auto iter : clusters2) { - // if (::fabs(iter[0].x - 1.5694) < 0.1 && - // ::fabs(iter[0].y - 15.3338) < 0.1) { - // found = true; - // break; - // } - // } - // assert(found); - // found = false; - // for (auto iter : clusters2) { - // if (::fabs(iter[0].x - 1.29624) < 0.1 && - // ::fabs(iter[0].y - 4.13919) < 0.1) { - // found = true; - // break; - // } - // } - // assert(found); - // Now try with multidimensional dataset (vector of arrays) // - RandGenParams<double> p2; + RandGenParams<double> p2; p2.seed = 123; p2.min_value = -20.0; p2.max_value = 20.0; - using col_t = std::array<double, 3>; + using col_t = std::array<double, 3>; - auto rand_vec = gen_uniform_real_dist<double>(df.get_index().size() * 3, p2); - StlVecType<col_t> multi_dimen_col(df.get_index().size()); - auto dist_func = [](const col_t &x, const col_t &y) -> double { - double sum { 0 }; + auto rand_vec = gen_uniform_real_dist<double>(df.get_index().size() * 3, p2); + StlVecType<col_t> multi_dimen_col(df.get_index().size()); + auto dist_func = + [](const col_t &x, const col_t &y) -> double { + double sum { 0 }; - for (std::size_t i { 0 }; i < x.size(); ++i) { - const double diff { x[i] - y[i] }; + for (std::size_t i { 0 }; i < x.size(); ++i) { + const double diff { x[i] - y[i] }; - sum += diff * diff; - } - return (std::sqrt(sum)); - }; + sum += diff * diff; + } + return (std::sqrt(sum)); + }; - for (std::size_t i { 0 }, j { 0 }; j < rand_vec.size(); ++i) { + for (std::size_t i { 0 }, j { 0 }; j < rand_vec.size(); ++i) { multi_dimen_col[i][0] = rand_vec[j++]; multi_dimen_col[i][1] = rand_vec[j++]; multi_dimen_col[i][2] = rand_vec[j++]; } - df.load_column<col_t>("multi_dimen_col", std::move(multi_dimen_col)); + df.load_column<col_t>("multi_dimen_col", std::move(multi_dimen_col)); - KMeansVisitor<4, col_t, unsigned long, 128> kmean(1000, true, dist_func); + KMeansVisitor<4, col_t, unsigned long, 128> kmean(1000, true, dist_func); - df.single_act_visit<col_t>("multi_dimen_col", kmean); + df.single_act_visit<col_t>("multi_dimen_col", kmean); assert(kmean.get_clusters_idxs().size() == 4); for (const auto &mean : kmean.get_result()) - std::cout << mean << "\n\n"; + std::cout << mean << "\n\n"; } diff --git a/docs/HTML/KolmoSmirnovTestVisitor.html b/docs/HTML/KolmoSmirnovTestVisitor.html index 7eb16822..c33e3ff8 100644 --- a/docs/HTML/KolmoSmirnovTestVisitor.html +++ b/docs/HTML/KolmoSmirnovTestVisitor.html @@ -121,7 +121,7 @@ ibm.single_act_visit<double, double>("IBM_Low", "IBM_High", ks_test); assert((std::fabs(ks_test.get_result() - 0.0296) < 0.0001)); - assert((std::fabs(ks_test.get_p_value() - 0.0242) < 0.0001)); + assert((std::fabs(ks_test.get_p_value() - 0.023725) < 0.000001)); ibm.single_act_visit<double, double>("IBM_Close", "uniform", ks_test); assert((std::fabs(ks_test.get_result() - 0.1224) < 0.0001)); diff --git a/docs/HTML/MannWhitneyUTestVisitor.html b/docs/HTML/MannWhitneyUTestVisitor.html index 630341e3..f5bad32b 100644 --- a/docs/HTML/MannWhitneyUTestVisitor.html +++ b/docs/HTML/MannWhitneyUTestVisitor.html @@ -113,8 +113,8 @@ assert((std::fabs(mwu_test.get_result() - 12643394.5) < 0.0001)); assert((std::fabs(mwu_test.get_u1() - 12667566.5) < 0.0001)); assert((std::fabs(mwu_test.get_u2() - 12643394.5) < 0.0001)); - assert((std::fabs(mwu_test.get_zscore() - -0.083) < 0.001)); - assert((std::fabs(mwu_test.get_pvalue() - 0.9339) < 0.0001)); + assert((std::fabs(mwu_test.get_zscore() - 0.082957) < 0.000001)); + assert((std::fabs(mwu_test.get_pvalue() - 0.933885) < 0.000001)); ibm.single_act_visit<double, double>("IBM_Low", "IBM_High", mwu_test); assert((std::fabs(mwu_test.get_result() - 12213043.0) < 0.0001)); @@ -134,14 +134,14 @@ assert((std::fabs(mwu_test.get_result() - 30.0) < 0.0001)); assert((std::fabs(mwu_test.get_u1() - 25310931.0) < 0.0001)); assert((std::fabs(mwu_test.get_u2() - 30.0) < 0.0001)); - assert((std::fabs(mwu_test.get_zscore() - -86.8661) < 0.001)); + assert((std::fabs(mwu_test.get_zscore() - 86.8661) < 0.001)); assert((std::fabs(mwu_test.get_pvalue() - 0.0) < 0.0001)); ibm.single_act_visit<double, double>("uniform", "exponential", mwu_test); assert((std::fabs(mwu_test.get_result() - 0.0) < 0.0001)); assert((std::fabs(mwu_test.get_u1() - 25310961.0) < 0.0001)); assert((std::fabs(mwu_test.get_u2() - 0.0) < 0.0001)); - assert((std::fabs(mwu_test.get_zscore() - -86.8663) < 0.001)); + assert((std::fabs(mwu_test.get_zscore() - 86.8663) < 0.001)); assert((std::fabs(mwu_test.get_pvalue() - 0.0) < 0.0001)); ibm.single_act_visit<double, double>("exponential", "lognormal", mwu_test); diff --git a/docs/HTML/MeanShiftVisitor.html b/docs/HTML/MeanShiftVisitor.html index 5b643d16..3a49e771 100644 --- a/docs/HTML/MeanShiftVisitor.html +++ b/docs/HTML/MeanShiftVisitor.html @@ -190,23 +190,21 @@ MeanShiftVisitor<double, std::string> mshift(1.0, 4, mean_shift_kernel::gaussian); - mshift.set_dist_func([](const double &x, const double &y) -> double { - return (std::fabs(x - y)); - }); + mshift.set_dist_func([](const double &x, const double &y) -> double { return (std::fabs(x - y)); }); df.single_act_visit<double>("IBM_Close", mshift); - assert(mshift.get_result().size() == 19); - assert(mshift.get_result()[0].size() == 106); - assert(mshift.get_result()[4].size() == 19); - assert(mshift.get_result()[6].size() == 274); - assert(mshift.get_result()[10].size() == 180); - assert(mshift.get_result()[14].size() == 29); - assert(mshift.get_result()[18].size() == 2); - assert(std::fabs(mshift.get_result()[0][6] - 184.16) < 0.001); - assert(std::fabs(mshift.get_result()[4][18] - 194.0) < 0.001); - assert(std::fabs(mshift.get_result()[6][273] - 154.31) < 0.001); - assert(std::fabs(mshift.get_result()[10][135] - 137.61) < 0.001); - assert(std::fabs(mshift.get_result()[18][1] - 94.77) < 0.001); + assert(mshift.get_result().size() == 18); + assert(mshift.get_result()[0].size() == 123); + assert(mshift.get_result()[4].size() == 57); + assert(mshift.get_result()[6].size() == 275); + assert(mshift.get_result()[10].size() == 54); + assert(mshift.get_result()[14].size() == 9); + assert(mshift.get_result()[17].size() == 2); + assert(std::fabs(mshift.get_result()[0][6] - 187.26) < 0.001); + assert(std::fabs(mshift.get_result()[4][18] - 166.08) < 0.001); + assert(std::fabs(mshift.get_result()[6][273] - 151.1) < 0.001); + assert(std::fabs(mshift.get_result()[10][35] - 129.57) < 0.001); + assert(std::fabs(mshift.get_result()[17][1] - 94.77) < 0.001); // Now multidimensional data // @@ -218,7 +216,8 @@ using col_t = std::array<double, 3>; - auto rand_vec = gen_uniform_real_dist<double>(df.get_index().size() * 3, p); + auto rand_vec = gen_uniform_real_dist<double>(df.get_index().size() * 3, p); + std::vector<col_t> multi_dimen_col(df.get_index().size()); for (std::size_t i { 0 }, j { 0 }; j < rand_vec.size(); ++i) { @@ -234,16 +233,16 @@ const auto &md_clusters = md_mshift.get_result(); - assert(md_clusters.size() == 53); // Number of clusters + assert(md_clusters.size() == 52); // Number of clusters - assert(md_clusters[0].size() == 74); + assert(md_clusters[0].size() == 73); assert(std::fabs(md_clusters[0][6][1] - -1.8807) < 0.0001); - assert(md_clusters[28].size() == 36); + assert(md_clusters[28].size() == 40); assert(std::fabs(md_clusters[28][3][0] - 12.6347) < 0.0001); - assert(md_clusters[52].size() == 1); - assert(std::fabs(md_clusters[52][0][2] - 19.2094) < 0.0001); + assert(md_clusters[51].size() == 1); + assert(std::fabs(md_clusters[51][0][2] - -19.7932) < 0.0001); } diff --git a/docs/HTML/MutualInfoVisitor.html b/docs/HTML/MutualInfoVisitor.html index 100bed21..a537313b 100644 --- a/docs/HTML/MutualInfoVisitor.html +++ b/docs/HTML/MutualInfoVisitor.html @@ -93,9 +93,7 @@ p.seed = 123; p.max_value = 4; p.min_value = -4; - df.load_data(std::move(idxvec), - std::make_pair("int_col_1", gen_uniform_int_dist<int>(idxvec.size(), p)), - std::make_pair("str_col", strvec)); + df.load_data(std::move(idxvec), std::make_pair("int_col_1", gen_uniform_int_dist<int>(idxvec.size(), p)), std::make_pair("str_col", strvec)); p.seed = 675; df.load_column("int_col_2", gen_uniform_int_dist<int>(idxvec.size(), p)); @@ -111,13 +109,13 @@ MutualInfoVisitor<int> minfo; df.single_act_visit<int, int>("int_col_1", "int_col_1", minfo); - assert((std::fabs(minfo.get_result() - 12.4866) < 0.0001)); + assert((std::fabs(minfo.get_result() - 2.89425) < 0.00001)); df.single_act_visit<int, int>("int_col_1", "int_col_2", minfo); - assert((std::fabs(minfo.get_result() - 1.81499) < 0.00001)); + assert((std::fabs(minfo.get_result() - 1.22157) < 0.00001)); df.single_act_visit<int, int>("int_col_1", "int_col_3", minfo); - assert((std::fabs(minfo.get_result() - 4.24521) < 0.00001)); + assert((std::fabs(minfo.get_result() - 0.954434) < 0.000001)); // Now multidimensional data // @@ -150,7 +148,7 @@ assert((std::fabs(mi_ary.get_result() - 1.58496) < 0.00001)); df.single_act_visit<vec_col_t, vec_col_t>("COL VEC2", "COL VEC2", mi_vec); - assert((std::fabs(mi_vec.get_result() - 4.0) < 0.000000001)); + assert((std::fabs(mi_vec.get_result() - 1.0) < 0.000000001)); df.single_act_visit<vec_col_t, vec_col_t>("COL VEC2", "COL VEC4", mi_vec); assert((std::fabs(mi_vec.get_result() - 1.0) < 0.000000001)); diff --git a/docs/HTML/PolicyLearningLossVisitor.html b/docs/HTML/PolicyLearningLossVisitor.html index b45abbcb..5578cb54 100644 --- a/docs/HTML/PolicyLearningLossVisitor.html +++ b/docs/HTML/PolicyLearningLossVisitor.html @@ -102,6 +102,7 @@ baseline: Baseline policy baseline_val: Value for constant baseline policy +epsilon: It is used to prevent division by zero or logarithm of zero diff --git a/docs/HTML/QuantileVisitor.html b/docs/HTML/QuantileVisitor.html index 8d8f628b..d7e7a1fc 100644 --- a/docs/HTML/QuantileVisitor.html +++ b/docs/HTML/QuantileVisitor.html @@ -102,143 +102,185 @@ - - - 
static void test_quantile()  {
-
-    std::cout << "\nTesting QuantileVisitor{ } ..." << std::endl;
-
-    std::vector<unsigned long>  idx =
-        { 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, 31, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40 };
-    std::vector<double> d1 =
-        { 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 };
-    MyDataFrame         df;
-
-    df.load_data(std::move(idx), std::make_pair("col_1", d1));
-    df.shuffle<1, double>({"col_1"}, false);
-
-    QuantileVisitor<double> v1(1, quantile_policy::mid_point);
-    auto                    result = df.single_act_visit<double>("col_1", v1).get_result();
-
-    assert(result == 40.0);
-
-    QuantileVisitor<double> v2(0.5, quantile_policy::mid_point);
-
-    result = df.single_act_visit<double>("col_1", v2).get_result();
-    assert(result == 20.5);
-
-    QuantileVisitor<double> v3(0.5, quantile_policy::linear);
-
-    result = df.single_act_visit<double>("col_1", v3).get_result();
-    assert(result == 20.5);
-
-    QuantileVisitor<double> v4(0.5, quantile_policy::higher_value);
-
-    result = df.single_act_visit<double>("col_1", v4).get_result();
-    assert(result == 21.0);
-
-    QuantileVisitor<double> v5(0.5, quantile_policy::lower_value);
-
-    result = df.single_act_visit<double>("col_1", v5).get_result();
-    assert(result == 20.0);
-
-    QuantileVisitor<double> v6(0.55, quantile_policy::mid_point);
-
-    result = df.single_act_visit<double>("col_1", v6).get_result();
-    assert(result == 22.5);
-
-    QuantileVisitor<double> v7(0.55, quantile_policy::linear);
-
-    result = df.single_act_visit<double>("col_1", v7).get_result();
-    assert(result == 22.45);
-
-    QuantileVisitor<double> v8(0.75, quantile_policy::mid_point);
-
-    result = df.single_act_visit<double>("col_1", v8).get_result();
-    assert(result == 30.5);
-
-    QuantileVisitor<double> v9(0.75, quantile_policy::linear);
-
-    result = df.single_act_visit<double>("col_1", v9).get_result();
-    assert(result == 30.25);
-
-    QuantileVisitor<double> v10(0, quantile_policy::linear);
-
-    result = df.single_act_visit<double>("col_1", v10).get_result();
-    assert(result == 1.0);
-
-    df.get_index().push_back(41);
-    df.get_column<double>("col_1").push_back(41);
-
-    QuantileVisitor<double> v11(0.75, quantile_policy::mid_point);
-
-    result = df.single_act_visit<double>("col_1", v11).get_result();
-    assert(result == 31.0);
-
-    QuantileVisitor<double> v12(0.75, quantile_policy::linear);
-
-    result = df.single_act_visit<double>("col_1", v12).get_result();
-    assert(result == 31.0);
-
-    QuantileVisitor<double> v13(0.75, quantile_policy::lower_value);
-
-    result = df.single_act_visit<double>("col_1", v13).get_result();
-    assert(result == 31.0);
-
-    QuantileVisitor<double> v14(0.75, quantile_policy::higher_value);
-
-    result = df.single_act_visit<double>("col_1", v14).get_result();
-    assert(result == 31.0);
-
-    QuantileVisitor<double> v15(0.71, quantile_policy::mid_point);
-
-    result = df.single_act_visit<double>("col_1", v15).get_result();
-    assert(result == 29.5);
-
-    QuantileVisitor<double> v16(0.71, quantile_policy::linear);
-
-    result = df.single_act_visit<double>("col_1", v16).get_result();
-    assert(result == 29.29);
-
-    QuantileVisitor<double> v17(0.23, quantile_policy::mid_point);
-
-    result = df.single_act_visit<double>("col_1", v17).get_result();
-    assert(result == 9.5);
-
-    QuantileVisitor<double> v18(0.2, quantile_policy::mid_point);
-
-    result = df.single_act_visit<double>("col_1", v18).get_result();
-    assert(result == 8.5);
-
-    QuantileVisitor<double> v19(0.23, quantile_policy::linear);
-
-    result = df.single_act_visit<double>("col_1", v19).get_result();
-    assert(result == 9.77);
-
-    QuantileVisitor<double> v20(0.23, quantile_policy::lower_value);
-
-    result = df.single_act_visit<double>("col_1", v20).get_result();
-    assert(result == 9.0);
-
-    QuantileVisitor<double> v21(0.23, quantile_policy::higher_value);
-
-    result = df.single_act_visit<double>("col_1", v21).get_result();
-    assert(result == 10.0);
-
-    QuantileVisitor<double> v22(1, quantile_policy::linear);
-
-    result = df.single_act_visit<double>("col_1", v22).get_result();
-    assert(result == 41.0);
+    
+      
+
#include <DataFrame/DataFrameStatsVisitors.h>
+
+template<typename T, typename I = unsigned long,
+         std::size_t A = 0>
+struct NQuantileVisitor;
+
+// -------------------------------------
+
+template<typename T, typename I = unsigned long,
+         std::size_t A = 0>
+using nqt_v = NQuantileVisitor<T, I, A>;
+

+      
+      
+        This is a "single action visitor", meaning it is passed the whole data vector in one call and you must use the single_act_visit() interface.


 
-    QuantileVisitor<double> v23(0, quantile_policy::mid_point);
+        This does the same thing the as above QuantileVisitor, but for a vector of quantiles. If you need multiple quantiles at the same time, this is more efficient than repeatedly calling QuantileVisitor.

+        
+        
+    explicit
+    NQuantileVisitor(std::vector &&quantiles,
+                     quantile_policy q_policy = quantile_policy::mid_point,
+                     bool skip_nan = false)
+        

+        
+      
+      
+        T: Column data type.

+        I: Index type.

+        A: Memory alignment boundary for vectors. Default is system default alignment

+      
+    
+  
 
-    result = df.single_act_visit<double>("col_1", v23).get_result();
-    assert(result == 1.0);
-}
-
- +
static void test_quantile()  {
+
+    std::cout << "\nTesting QuantileVisitor{ } ..." << std::endl;
+
+    StlVecType<unsigned long>  idx = { 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, 31, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40 };
+    StlVecType<double>         d1 = { 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 };
+    MyDataFrame                df;
+
+    df.load_data(std::move(idx), std::make_pair("col_1", d1));
+    df.shuffle<double>({"col_1"}, false);
+
+    QuantileVisitor<double, unsigned long, 128> v1 { 1, quantile_policy::mid_point };
+    auto                                        result { df.single_act_visit<double>("col_1", v1).get_result() };
+
+    assert(result == 40.0);
+
+    QuantileVisitor<double, unsigned long, 128> v2 { 0.5, quantile_policy::mid_point };
+
+    result = df.single_act_visit<double>("col_1", v2).get_result();
+    assert(result == 20.0);
+
+    QuantileVisitor<double, unsigned long, 128> v3 { 0.5, quantile_policy::linear };
+
+    result = df.single_act_visit<double>("col_1", v3).get_result();
+    assert(result == 20.0);
+
+    QuantileVisitor<double, unsigned long, 128> v4 { 0.5, quantile_policy::higher_value };
+
+    result = df.single_act_visit<double>("col_1", v4).get_result();
+    assert(result == 20.0);
+
+    QuantileVisitor<double, unsigned long, 128> v5 { 0.5, quantile_policy::lower_value };
+
+    result = df.single_act_visit<double>("col_1", v5).get_result();
+    assert(result == 20.0);
+
+    QuantileVisitor<double, unsigned long, 128> v6 { 0.55, quantile_policy::mid_point };
+
+    result = df.single_act_visit<double>("col_1", v6).get_result();
+    assert(result == 22.0);
+
+    QuantileVisitor<double, unsigned long, 128> v7 { 0.55, quantile_policy::linear };
+
+    result = df.single_act_visit<double>("col_1", v7).get_result();
+    assert(result == 22.0);
+
+    QuantileVisitor<double, unsigned long, 128> v8 { 0.75, quantile_policy::mid_point };
+
+    result = df.single_act_visit<double>("col_1", v8).get_result();
+    assert(result == 30.0);
+
+    QuantileVisitor<double, unsigned long, 128> v9 { 0.75, quantile_policy::linear };
+
+    result = df.single_act_visit<double>("col_1", v9).get_result();
+    assert(result == 30.0);
+
+    QuantileVisitor<double, unsigned long, 128> v10 { 0, quantile_policy::linear };
+
+    result = df.single_act_visit<double>("col_1", v10).get_result();
+    assert(result == 1.0);
+
+    df.get_index().push_back(41);
+    df.get_column<double>("col_1").push_back(41);
+
+    QuantileVisitor<double, unsigned long, 128> v11 { 0.75, quantile_policy::mid_point };
+
+    result = df.single_act_visit<double>("col_1", v11).get_result();
+    assert(result == 30.5);
+
+    QuantileVisitor<double, unsigned long, 128> v12 { 0.75, quantile_policy::linear };
+
+    result = df.single_act_visit<double>("col_1", v12).get_result();
+    assert(result == 30.75);
+
+    QuantileVisitor<double, unsigned long, 128> v13 { 0.75, quantile_policy::lower_value };
+
+    result = df.single_act_visit<double>("col_1", v13).get_result();
+    assert(result == 30.0);
+
+    QuantileVisitor<double, unsigned long, 128> v14 { 0.75, quantile_policy::higher_value };
+
+    result = df.single_act_visit<double>("col_1", v14).get_result();
+    assert(result == 31.0);
+
+    QuantileVisitor<double, unsigned long, 128> v15 { 0.71, quantile_policy::mid_point };
+
+    result = df.single_act_visit<double>("col_1", v15).get_result();
+    assert(result == 29.5);
+
+    QuantileVisitor<double, unsigned long, 128> v16 { 0.71, quantile_policy::linear };
+
+    result = df.single_act_visit<double>("col_1", v16).get_result();
+    assert(result == 29.11);
+
+    QuantileVisitor<double, unsigned long, 128> v17 { 0.23, quantile_policy::mid_point };
+
+    result = df.single_act_visit<double>("col_1", v17).get_result();
+    assert(result == 9.5);
+
+    QuantileVisitor<double, unsigned long, 128> v18 { 0.2, quantile_policy::mid_point };
+
+    result = df.single_act_visit<double>("col_1", v18).get_result();
+    assert(result == 8.5);
+
+    QuantileVisitor<double, unsigned long, 128> v19 { 0.23, quantile_policy::linear };
+
+    result = df.single_act_visit<double>("col_1", v19).get_result();
+    assert(result == 9.43);
+
+    QuantileVisitor<double, unsigned long, 128> v20 { 0.23, quantile_policy::lower_value };
+
+    result = df.single_act_visit<double>("col_1", v20).get_result();
+    assert(result == 9.0);
+
+    QuantileVisitor<double, unsigned long, 128> v21 { 0.23, quantile_policy::higher_value };
+
+    result = df.single_act_visit<double>("col_1", v21).get_result();
+    assert(result == 10.0);
+
+    QuantileVisitor<double, unsigned long, 128> v22 { 1, quantile_policy::linear };
+
+    result = df.single_act_visit<double>("col_1", v22).get_result();
+    assert(result == 41.0);
+
+    QuantileVisitor<double, unsigned long, 128> v23 { 0, quantile_policy::mid_point };
+
+    result = df.single_act_visit<double>("col_1", v23).get_result();
+    assert(result == 1.0);
+
+    // N quantiles
+    //
+    NQuantileVisitor<double, unsigned long, 128>    nv { { 0.25, 0.75, 1.0, 0.0, 0.15, 0.5 }, quantile_policy::mid_point };
+    const auto                                      nres { df.single_act_visit<double>("col_1", nv).get_result() };
+
+    assert(nres.size() == 6);
+    assert(nres[0] == 10.5);  // 25%
+    assert(nres[1] == 30.5);  // 75%
+    assert(nres[2] == 41.0);  // 100%
+    assert(nres[3] == 1.0);   // 0%
+    assert(nres[4] == 6.5);   // 15%
+    assert(nres[5] == 20.5);    // 50%
+}
+

C++ DataFrame diff --git a/docs/HTML/SigmoidVisitor.html b/docs/HTML/SigmoidVisitor.html index 0370a1a0..bda29aad 100644 --- a/docs/HTML/SigmoidVisitor.html +++ b/docs/HTML/SigmoidVisitor.html @@ -195,12 +195,11 @@ std::cout << "\nTesting SigmoidVisitor{ } ..." << std::endl; - StlVecType<unsigned long> idx = - { 123450, 123451, 123452, 123453, 123454, 123455, 123456, 123457, 123458, 123459, 123460, 123461, 123462, 123466, 123467, 123468, 123469, 123470, 123471, 123472, 123473 }; + StlVecType<unsigned long> idx = { 123450, 123451, 123452, 123453, 123454, 123455, 123456, 123457, 123458, 123459, 123460, 123461, 123462, 123466, 123467, 123468, 123469, 123470, 123471, 123472, 123473 }; StlVecType<double> d1 = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21 }; StlVecType<double> d2 = { 0.23, 0.25, 0.256, 0.26, 0.268, 0.271, 0.279, 0.285, 0.29, 0.3, 0.5, -0.2, 1, 0, 2, 0, -0.1, 0.55, 0.58, 0.6, 0.7 }; StlVecType<int> i1 = { 22, 23, 24, 25, 99 }; - MyDataFrame df; + MyDataFrame df; df.load_data(std::move(idx), std::make_pair("d1_col", d1), @@ -214,20 +213,20 @@ SigmoidVisitor<double, unsigned long, 64> sig_err(sigmoid_type::error_function); SigmoidVisitor<double, unsigned long, 64> sig_gud(sigmoid_type::gudermannian); SigmoidVisitor<double, unsigned long, 64> sig_smo(sigmoid_type::smoothstep); - const auto log_result = df.single_act_visit<double>("d1_col", sig_log).get_result(); - const auto alg_result = df.single_act_visit<double>("d1_col", sig_alg).get_result(); - const auto tan_result = df.single_act_visit<double>("d1_col", sig_tan).get_result(); - const auto atan_result = df.single_act_visit<double>("d1_col", sig_atan).get_result(); - const auto err_result = df.single_act_visit<double>("d1_col", sig_err).get_result(); - const auto gud_result = df.single_act_visit<double>("d1_col", sig_gud).get_result(); - const auto smo_result = df.single_act_visit<double>("d2_col", sig_smo).get_result(); + const auto log_result = df.single_act_visit<double>("d1_col", sig_log).get_result(); + const auto alg_result = df.single_act_visit<double>("d1_col", sig_alg).get_result(); + const auto tan_result = df.single_act_visit<double>("d1_col", sig_tan).get_result(); + const auto atan_result = df.single_act_visit<double>("d1_col", sig_atan).get_result(); + const auto err_result = df.single_act_visit<double>("d1_col", sig_err).get_result(); + const auto gud_result = df.single_act_visit<double>("d1_col", sig_gud).get_result(); + const auto smo_result = df.single_act_visit<double>("d2_col", sig_smo).get_result(); StlVecType<double> result { 0.731059, 0.880797, 0.952574, 0.982014, 0.993307, 0.997527, 0.999089, 0.999665, 0.999877, 0.999955, 0.999983, 0.999994, 0.999998, 0.999999, 1, 1, 1, 1, 1, 1, 1 }; for (size_t i = 0; i < result.size(); ++i) assert(fabs(result[i] - log_result[i]) < 0.00001); - result = StlVecType<double> { 0.707107, 0.447214, 0.316228, 0.242536, 0.196116, 0.164399, 0.141421, 0.124035, 0.110432, 0.0995037, 0.0905357, 0.0830455, 0.0766965, 0.071247, 0.066519, 0.0623783, 0.058722, 0.05547, 0.0525588, 0.0499376, 0.0475651 }; + result = StlVecType<double> { 0.707107, 0.894427, 0.948683, 0.970143, 0.980581, 0.986394, 0.989949, 0.992278, 0.993884, 0.995037, 0.995893, 0.996546, 0.997054, 0.997459, 0.997785, 0.998053, 0.998274, 0.99846, 0.998618, 0.998752, 0.998868 }; for (size_t i = 0; i < result.size(); ++i) assert(fabs(result[i] - alg_result[i]) < 0.00001); @@ -287,8 +286,8 @@ for (const auto &vec : md_lgb_res) assert(vec.size() == dim); assert(std::fabs(md_lgb_res[0][0] - 0.707107) < 0.000001); - assert(std::fabs(md_lgb_res[5][1] - 0.5547) < 0.0001); - assert(std::fabs(md_lgb_res[9][2] - 0.447214) < 0.000001); + assert(std::fabs(md_lgb_res[5][1] - 0.83205) < 0.00001); + assert(std::fabs(md_lgb_res[9][2] - 0.894427) < 0.000001); SigmoidVisitor<ary_col_t, unsigned long, 64> md_gud_v { sigmoid_type::gudermannian }; const auto &md_gud_res = df.single_act_visit<ary_col_t>("ARY MD", md_gud_v).get_result(); @@ -339,9 +338,9 @@ df.single_act_visit<double>("dbl_col_2", gelu); assert(gelu.get_result().size() == 15); - assert(std::abs(gelu.get_result()[0] - 0.242) < 0.0001); - assert(std::abs(gelu.get_result()[5] - 0.2153) < 0.0001); - assert(std::abs(gelu.get_result()[14] - 0.0967) < 0.0001); + assert(std::abs(gelu.get_result()[0] - 0.841345) < 0.000001); + assert(std::abs(gelu.get_result()[5] - 0.511188) < 0.000001); + assert(std::abs(gelu.get_result()[14] - 0.149677) < 0.000001); recf_v<double, unsigned long> silu(rectify_type::SiLU); diff --git a/docs/HTML/StationaryCheckVisitor.html b/docs/HTML/StationaryCheckVisitor.html index d2685d34..9e7616b7 100644 --- a/docs/HTML/StationaryCheckVisitor.html +++ b/docs/HTML/StationaryCheckVisitor.html @@ -216,54 +216,54 @@ StationaryCheckVisitor<double, std::string> sc2 { stationary_test::adf, { .adf_lag = 10, .adf_with_trend = false } }; df.single_act_visit<double>("IBM_Close", sc2); - assert(std::fabs(sc2.get_adf_statistic() - 0.989687) < 0.00001); + assert(std::fabs(sc2.get_adf_statistic() - -1.80735) < 0.00001); StationaryCheckVisitor<double, std::string> sc3 { stationary_test::adf, { .adf_lag = 25, .adf_with_trend = false } }; df.single_act_visit<double>("IBM_Close", sc3); - assert(std::fabs(sc3.get_adf_statistic() - 0.974531) < 0.0000001); + assert(std::fabs(sc3.get_adf_statistic() - -1.59054) < 0.00001); StationaryCheckVisitor<double, std::string> sc4 { stationary_test::adf, { .adf_lag = 10, .adf_with_trend = false } }; df.single_act_visit<double>("normal_col", sc4); - assert(std::fabs(sc4.get_adf_statistic() - 0.0289613) < 0.0000001); + assert(std::fabs(sc4.get_adf_statistic() - -21.1568) < 0.0001); StationaryCheckVisitor<double, std::string> sc5 { stationary_test::adf, { .adf_lag = 25, .adf_with_trend = false } }; df.single_act_visit<double>("normal_col", sc5); - assert(std::fabs(sc5.get_adf_statistic() - 0.0208191) < 0.0000001); + assert(std::fabs(sc5.get_adf_statistic() - -13.5343) < 0.0001); // ADF tests with trend // StationaryCheckVisitor<double, std::string> sc6 { stationary_test::adf, { .adf_lag = 10, .adf_with_trend = true } }; df.single_act_visit<double>("IBM_Close", sc6); - assert(std::fabs(sc6.get_adf_statistic() - 0.977705) < 0.000001); + assert(std::fabs(sc6.get_adf_statistic() - -1.83342) < 0.00001); StationaryCheckVisitor<double, std::string> sc7 { stationary_test::adf, { .adf_lag = 25, .adf_with_trend = true } }; df.single_act_visit<double>("IBM_Close", sc7); - assert(std::fabs(sc7.get_adf_statistic() - 0.946614) < 0.000001); + assert(std::fabs(sc7.get_adf_statistic() - -1.38926) < 0.00001); StationaryCheckVisitor<double, std::string> sc8 { stationary_test::adf, { .adf_lag = 10, .adf_with_trend = true } }; df.single_act_visit<double>("normal_col", sc8); - assert(std::fabs(sc8.get_adf_statistic() - 0.0289582) < 0.0000001); + assert(std::fabs(sc8.get_adf_statistic() - -21.155) < 0.001); StationaryCheckVisitor<double, std::string> sc9 { stationary_test::adf, { .adf_lag = 25, .adf_with_trend = true } }; df.single_act_visit<double>("normal_col", sc9); - assert(std::fabs(sc9.get_adf_statistic() - 0.020812) < 0.0000001); + assert(std::fabs(sc9.get_adf_statistic() - -13.5341) < 0.0001); StationaryCheckVisitor<double, std::string> sc10 { stationary_test::adf, { .adf_lag = 10, .adf_with_trend = true } }; df.single_act_visit<double>("log close", sc10); - assert(std::fabs(sc10.get_adf_statistic() - 0.972062) < 0.000001); + assert(std::fabs(sc10.get_adf_statistic() - -2.21398) < 0.00001); StationaryCheckVisitor<double, std::string> sc11 { stationary_test::adf, { .adf_lag = 10, .adf_with_trend = true } }; df.single_act_visit<double>("residual close", sc11); - assert(std::fabs(sc11.get_adf_statistic() - 0.679027) < 0.000001); + assert(std::fabs(sc11.get_adf_statistic() - -3.44574) < 0.00001); // Now multidimensional data // @@ -273,20 +273,20 @@ using vec_col_t = std::vector<double>; std::vector<ary_col_t> stationary_ary { - {10.1, 5.0, -2.0}, { 9.9, 5.2, -2.1}, {10.2, 4.9, -1.9}, {10.0, 5.1, -2.0}, { 9.8, 5.0, -2.2}, {10.1, 5.3, -1.8}, - {10.0, 4.8, -2.1}, {10.2, 5.1, -2.0}, { 9.9, 5.0, -2.0}, {10.1, 5.2, -1.9}, {10.0, 4.9, -2.1}, { 9.8, 5.1, -2.0}, {10.2, 5.0, -2.1} + { 10.1, 5.0, -2.0 }, { 9.9, 5.2, -2.1 }, { 10.2, 4.9, -1.9 }, { 10.0, 5.1, -2.0 }, { 9.8, 5.0, -2.2 }, { 10.1, 5.3, -1.8 }, + { 10.0, 4.8, -2.1 }, { 10.2, 5.1, -2.0 }, { 9.9, 5.0, -2.0 }, { 10.1, 5.2, -1.9 }, { 10.0, 4.9, -2.1 }, { 9.8, 5.1, -2.0 }, { 10.2, 5.0, -2.1 } }; std::vector<vec_col_t> stationary_vec { - {10.1, 5.0, -2.0}, { 9.9, 5.2, -2.1}, {10.2, 4.9, -1.9}, {10.0, 5.1, -2.0}, { 9.8, 5.0, -2.2}, {10.1, 5.3, -1.8}, - {10.0, 4.8, -2.1}, {10.2, 5.1, -2.0}, { 9.9, 5.0, -2.0}, {10.1, 5.2, -1.9}, {10.0, 4.9, -2.1}, { 9.8, 5.1, -2.0}, {10.2, 5.0, -2.1} + { 10.1, 5.0, -2.0 }, { 9.9, 5.2, -2.1 }, { 10.2, 4.9, -1.9 }, { 10.0, 5.1, -2.0 }, { 9.8, 5.0, -2.2 }, { 10.1, 5.3, -1.8 }, + { 10.0, 4.8, -2.1 }, { 10.2, 5.1, -2.0 }, { 9.9, 5.0, -2.0 }, { 10.1, 5.2, -1.9 }, { 10.0, 4.9, -2.1 }, { 9.8, 5.1, -2.0 }, { 10.2, 5.0, -2.1 } }; std::vector<vec_col_t> non_stationary_vec { - { 1.0, 10.0, -5.0}, { 1.5, 10.5, -4.8}, { 2.1, 11.0, -4.5}, { 2.8, 11.8, -4.0}, { 3.6, 12.5, -3.5}, { 4.5, 13.3, -3.0}, - { 5.5, 14.2, -2.4}, { 6.6, 15.0, -1.8}, { 7.8, 16.1, -1.0}, { 9.1, 17.3, -0.2}, {10.5, 18.6, 0.8}, {12.0, 20.0, 1.9}, {13.6, 21.5, 3.0} + { 1.0, 10.0, -5.0 }, { 1.5, 10.5, -4.8 }, { 2.1, 11.0, -4.5 }, { 2.8, 11.8, -4.0 }, { 3.6, 12.5, -3.5 }, { 4.5, 13.3, -3.0 }, + { 5.5, 14.2, -2.4 }, { 6.6, 15.0, -1.8 }, { 7.8, 16.1, -1.0 }, { 9.1, 17.3, -0.2 }, { 10.5, 18.6, 0.8 }, { 12.0, 20.0, 1.9 }, { 13.6, 21.5, 3.0 } }; std::vector<ary_col_t> non_stationary_ary { - { 1.0, 10.0, -5.0}, { 1.5, 10.5, -4.8}, { 2.1, 11.0, -4.5}, { 2.8, 11.8, -4.0}, { 3.6, 12.5, -3.5}, { 4.5, 13.3, -3.0}, - { 5.5, 14.2, -2.4}, { 6.6, 15.0, -1.8}, { 7.8, 16.1, -1.0}, { 9.1, 17.3, -0.2}, {10.5, 18.6, 0.8}, {12.0, 20.0, 1.9}, {13.6, 21.5, 3.0} + { 1.0, 10.0, -5.0 }, { 1.5, 10.5, -4.8 }, { 2.1, 11.0, -4.5 }, { 2.8, 11.8, -4.0 }, { 3.6, 12.5, -3.5 }, { 4.5, 13.3, -3.0 }, + { 5.5, 14.2, -2.4 }, { 6.6, 15.0, -1.8 }, { 7.8, 16.1, -1.0 }, { 9.1, 17.3, -0.2 }, { 10.5, 18.6, 0.8 }, { 12.0, 20.0, 1.9 }, { 13.6, 21.5, 3.0 } }; df.load_column<vec_col_t>("STATION VEC", std::move(stationary_vec), nan_policy::dont_pad_with_nans); @@ -305,14 +305,14 @@ df.single_act_visit<ary_col_t>("STATION ARY", kpss_ary_v); assert(adf_vec_v.get_adf_statistic().size() == dim); - assert(std::abs(adf_vec_v.get_adf_statistic()[0] - -0.586018) < 0.000001); - assert(std::abs(adf_vec_v.get_adf_statistic()[2] - -0.705822) < 0.000001); + assert(std::abs(adf_vec_v.get_adf_statistic()[0] - -3.55811) < 0.00001); + assert(std::abs(adf_vec_v.get_adf_statistic()[2] - -4.60884) < 0.00001); assert(kpss_vec_v.get_kpss_statistic().size() == dim); - assert(std::abs(kpss_vec_v.get_kpss_statistic()[0] - 0.0) < 0.00000001); - assert(std::abs(kpss_vec_v.get_kpss_statistic()[2] - 0.0) < 0.00000001); + assert(std::abs(kpss_vec_v.get_kpss_statistic()[0] - 0.1) < 0.01); + assert(std::abs(kpss_vec_v.get_kpss_statistic()[2] - 0.1) < 0.01); assert(kpss_vec_v.get_kpss_value().size() == dim); - assert(std::abs(kpss_vec_v.get_kpss_value()[0] - 326.728) < 0.001); - assert(std::abs(kpss_vec_v.get_kpss_value()[2] - 2239.29) < 0.01); + assert(std::abs(kpss_vec_v.get_kpss_value()[0] - 0.033409) < 0.000001); + assert(std::abs(kpss_vec_v.get_kpss_value()[2] - 0.035666) < 0.000001); df.single_act_visit<vec_col_t>("NON STATION VEC", adf_vec_v); df.single_act_visit<ary_col_t>("NON STATION ARY", adf_ary_v); @@ -320,14 +320,14 @@ df.single_act_visit<ary_col_t>("NON STATION ARY", kpss_ary_v); assert(adf_vec_v.get_adf_statistic().size() == dim); - assert(std::abs(adf_vec_v.get_adf_statistic()[0] - 0.904666) < 0.000001); - assert(std::abs(adf_vec_v.get_adf_statistic()[2] - 0.896825) < 0.000001); + assert((std::abs(adf_vec_v.get_adf_statistic()[0] - 2.66413) < 0.00001 || std::abs(adf_vec_v.get_adf_statistic()[0] - 2.61647) < 0.00001)); + assert(std::abs(adf_vec_v.get_adf_statistic()[2] - 2.64523) < 0.00001); assert(kpss_vec_v.get_kpss_statistic().size() == dim); - assert(std::abs(kpss_vec_v.get_kpss_statistic()[0] - 0.0) < 0.00000001); - assert(std::abs(kpss_vec_v.get_kpss_statistic()[2] - 0.0) < 0.00000001); + assert(std::abs(kpss_vec_v.get_kpss_statistic()[0] - 0.1) < 0.01); + assert(std::abs(kpss_vec_v.get_kpss_statistic()[2] - 0.1) < 0.01); assert(kpss_vec_v.get_kpss_value().size() == dim); - assert(std::abs(kpss_vec_v.get_kpss_value()[0] - 10.76) < 0.01); - assert(std::abs(kpss_vec_v.get_kpss_value()[2] - 189.563) < 0.001); + assert(std::abs(kpss_vec_v.get_kpss_value()[0] - 0.318681) < 0.000001); + assert(std::abs(kpss_vec_v.get_kpss_value()[2] - 0.311255) < 0.000001); } diff --git a/docs/HTML/get_data_by_dbscan.html b/docs/HTML/get_data_by_dbscan.html index 27b41478..d73542c9 100644 --- a/docs/HTML/get_data_by_dbscan.html +++ b/docs/HTML/get_data_by_dbscan.html @@ -145,24 +145,17 @@ auto views = view.get_view_by_dbscan<double, double, long>("IBM_Close", 10, 4); auto dfs = df.get_data_by_dbscan<double, double, long>("IBM_Close", 10, 4); - assert(views.size() == 36); - assert(dfs.size() == 36); + assert(views.size() == 2); + assert(dfs.size() == 2); - assert(views[0].get_index().size() == 5); - assert(std::fabs(views[0].get_column<double>("IBM_Close")[4] - 185.69) < 0.001); + assert(views[0].get_index().size() == 1705); + assert( + std::fabs(views[0].get_column<double>("IBM_Close")[4] - 187.97) < 0.001); - assert(dfs[5].get_index().size() == 30); - assert(std::fabs(dfs[5].get_column<double>("IBM_Open")[15] - 180.87) < 0.001); + // views[1].write<std::ostream, double, long>(std::cout, io_format::pretty_prt, { .precision = 3 }); - assert(views[16].get_index().size() == 39); - assert(std::fabs(views[16].get_column<double>("IBM_High")[3] - 170.85) < 0.001); - - // This is the last DataFrame which contains the data corresponding to - // noisy close prices - // - assert(views[35].get_index().size() == 16); - assert(views[35].get_column<long>("IBM_Volume")[0] == 3821400); - assert(views[35].get_index()[1] == "2020-03-12"); + assert(dfs[1].get_index().size() == 16); + assert(std::fabs(dfs[1].get_column<double>("IBM_Open")[15] - 107.25) < 0.001); } diff --git a/docs/HTML/get_data_by_mshift.html b/docs/HTML/get_data_by_mshift.html index 8de3e3ef..835ae907 100644 --- a/docs/HTML/get_data_by_mshift.html +++ b/docs/HTML/get_data_by_mshift.html @@ -154,23 +154,23 @@ // auto views = view.get_view_by_mshift<double, double, long>("IBM_Close", 1, 4, mean_shift_kernel::gaussian); auto dfs = df.get_data_by_mshift<double, double, long>("IBM_Close", 1, 4, mean_shift_kernel::gaussian); - + assert(views.size() == 38); assert(dfs.size() == 38); - assert(views[0].get_index().size() == 56); - assert(dfs[0].get_index().size() == 56); - assert(views[4].get_index().size() == 20); - assert(views[6].get_index().size() == 3); - assert(views[10].get_index().size() == 45); - assert(views[14].get_index().size() == 101); - assert(views[18].get_index().size() == 164); - assert(dfs[18].get_index().size() == 164); + assert(views[0].get_index().size() == 57); + assert(dfs[0].get_index().size() == 57); + assert(views[4].get_index().size() == 14); + assert(views[6].get_index().size() == 26); + assert(views[10].get_index().size() == 122); + assert(views[14].get_index().size() == 25); + assert(views[18].get_index().size() == 89); + assert(dfs[18].get_index().size() == 89); - assert((std::fabs(views[0].get_column<double>("IBM_Close")[7] - 183.69) < 0.001)); - assert((std::fabs(dfs[5].get_column<double>("IBM_Open")[15] - 173.91) < 0.001)); - assert((std::fabs(views[16].get_column<double>("IBM_High")[3] - 166.02) < 0.001)); - assert(dfs[18].get_column<long>("IBM_Volume")[0] == 10189700); - assert(views[18].get_index()[1] == "2015-09-01"); + assert((std::fabs(views[0].get_column<double>("IBM_Close")[7] - 187.74) < 0.001)); + assert((std::fabs(dfs[5].get_column<double>("IBM_Open")[15] - 172.97) < 0.001)); + assert((std::fabs(views[16].get_column<double>("IBM_High")[3] - 148.4) < 0.001)); + assert(dfs[18].get_column<long>("IBM_Volume")[0] == 7073200); + assert(views[18].get_index()[1] == "2015-10-20"); } diff --git a/include/DataFrame/DataFrameMLVisitors.h b/include/DataFrame/DataFrameMLVisitors.h index 1895c950..9300d898 100644 --- a/include/DataFrame/DataFrameMLVisitors.h +++ b/include/DataFrame/DataFrameMLVisitors.h @@ -46,6 +46,7 @@ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. #include #include #include +#include #include #include @@ -63,8 +64,8 @@ struct SLRegressionVisitor { DEFINE_VISIT_BASIC_TYPES_2 - inline void operator() (const index_type &idx, - const value_type &x, const value_type &y) { + inline void operator()(const index_type &idx, + const value_type &x, const value_type &y) { if (skip_nan_ && (is_nan__(x) || is_nan__(y))) [[unlikely]] return; @@ -78,33 +79,34 @@ struct SLRegressionVisitor { } PASS_DATA_ONE_BY_ONE_2 - inline void pre () { + inline void pre() { n_ = 0; s_xy_ = 0; x_stats_.pre(); y_stats_.pre(); } - inline void post () { } + inline void post() { } - inline size_type get_count () const { return (n_); } - inline result_type get_slope () const { + inline size_type get_count() const { return (n_); } + inline result_type get_slope() const { // Sum of the squares of the difference between each x and // the mean x value. // - const value_type s_xx = - x_stats_.get_variance() * value_type(n_ - 1); + const value_type s_xx { + x_stats_.get_variance() * value_type(n_ - 1) + }; return (s_xy_ / s_xx); } - inline result_type get_intercept () const { + inline result_type get_intercept() const { return (y_stats_.get_mean() - get_slope() * x_stats_.get_mean()); } - inline result_type get_corr () const { + inline result_type get_corr() const { - const value_type t = x_stats_.get_std() * y_stats_.get_std(); + const value_type t { x_stats_.get_std() * y_stats_.get_std() }; return (s_xy_ / (value_type(n_ - 1) * t)); } @@ -159,37 +161,39 @@ struct KMeansVisitor { inline void calc_k_means_(const H &column_begin, size_type col_s) { std::random_device rd; - std::mt19937 gen( - (seed_ != seed_t(-1)) ? seed_ : rd()); - std::uniform_int_distribution rd_gen(0, col_s - 1); + std::mt19937 gen { + (seed_ != seed_t(-1)) ? seed_ : rd() + }; + std::uniform_int_distribution rd_gen { 0, col_s - 1 }; // Pick centroids as random points from the col. // for (auto &k_mean : result_) [[likely]] { - const value_type &value = *(column_begin + rd_gen(gen)); + const value_type &value { *(column_begin + rd_gen(gen)) }; if (! is_nan__(value)) [[likely]] k_mean = value; } - for (size_type iter = 0; iter < iter_num_; ++iter) [[likely]] { + for (size_type iter { 0 }; iter < iter_num_; ++iter) [[likely]] { result_type new_means { value_type() }; std::array counts { 0.0 }; // Find assignments. // - for (size_type point = 0; point < col_s; ++point) [[likely]] { - const value_type &value = *(column_begin + point); + for (size_type point { 0 }; point < col_s; ++point) [[likely]] { + const value_type &value { *(column_begin + point) }; if (! is_nan__(value)) [[likely]] { - double best_distance = - std::numeric_limits::max(); - size_type best_cluster = 0; + double best_distance { + std::numeric_limits::max() + }; + size_type best_cluster { 0 }; - for (size_type cluster = 0; cluster < K; - ++cluster) [[likely]] { - const double distance = - dfunc_(value, result_[cluster]); + for (size_type cluster { 0 }; cluster < K; ++cluster) { + const double distance { + dfunc_(value, result_[cluster]) + }; if (distance < best_distance) { best_distance = distance; @@ -199,27 +203,33 @@ struct KMeansVisitor { // Sum up and count points for each cluster. // - auto &nm = new_means[best_cluster]; + auto &nm { new_means[best_cluster] }; nm = nm + value; counts[best_cluster] += 1.0; } } - bool done = true; + bool done { true }; - // Divide sums by counts to get new centroids. + // Divide sums by counts to get new centroids. A cluster that + // picked up no points this pass is left at its previous + // centroid -- there is nothing to average, and the old + // 0/0 -> 0/1 guard was substituting a fabricated value_type() + // (0 for arithmetic T), silently teleporting an empty + // cluster's centroid to zero instead of leaving it alone. // - for (size_type cluster = 0; cluster < K; ++cluster) [[likely]] { - // Turn 0/0 into 0/1 to avoid zero division. - const double count = - std::max(1.0, counts[cluster]); - const value_type value = new_means[cluster] / count; - value_type &result = result_[cluster]; - - if (dfunc_(value, result) > 0.0000001) { - done = false; - result = value; + for (size_type cluster { 0 }; cluster < K; ++cluster) [[likely]] { + if (counts[cluster] > 0.0) { + const value_type value { + new_means[cluster] / counts[cluster] + }; + value_type &result { result_[cluster] }; + + if (dfunc_(value, result) > 0.0000001) { + done = false; + result = value; + } } } @@ -236,21 +246,21 @@ struct KMeansVisitor { cluster_type clusters; order_type clusters_idxs; - for (size_type i = 0; i < K; ++i) [[likely]] { + for (size_type i { 0 }; i < K; ++i) [[likely]] { clusters[i].reserve(col_s / K + 2); clusters[i].push_back(&(result_[i])); clusters_idxs[i].reserve(col_s / K + 2); } - for (size_type j = 0; j < col_s; ++j) [[likely]] { - const value_type &value = *(column_begin + j); + for (size_type j { 0 }; j < col_s; ++j) [[likely]] { + const value_type &value { *(column_begin + j) }; if (! is_nan__(value)) [[likely]] { double min_dist { std::numeric_limits::max() }; size_type min_idx { 0 }; - for (size_type i = 0; i < K; ++i) { - const double dist = dfunc_(value, result_[i]); + for (size_type i { 0 }; i < K; ++i) { + const double dist { dfunc_(value, result_[i]) }; if (dist < min_dist) { min_dist = dist; @@ -270,38 +280,36 @@ struct KMeansVisitor { template inline void - operator() (const IV &idx_begin, const IV &idx_end, - const H &column_begin, const H &column_end) { + operator()(const IV &idx_begin, const IV &idx_end, + const H &column_begin, const H &column_end) { GET_COL_SIZE calc_k_means_(column_begin, col_s); - if (cc_) - calc_clusters_(column_begin, col_s); + if (cc_) calc_clusters_(column_begin, col_s); } - inline void pre () { + inline void pre() { for (auto &iter : clusters_) iter.clear(); for (auto &iter : clusters_idxs_) iter.clear(); } - inline void post () { } - inline const result_type &get_result () const { return (result_); } - inline result_type &get_result () { return (result_); } - inline const cluster_type &get_clusters () const { return (clusters_); } - inline cluster_type &get_clusters () { return (clusters_); } + inline void post() { } + inline const result_type &get_result() const { return (result_); } + inline result_type &get_result() { return (result_); } + inline const cluster_type &get_clusters() const { return (clusters_); } + inline cluster_type &get_clusters() { return (clusters_); } inline const order_type & - get_clusters_idxs () const { return (clusters_idxs_); } + get_clusters_idxs() const { return (clusters_idxs_); } explicit - KMeansVisitor( - size_type num_of_iter, - bool calc_clusters = true, - distance_func f = - [](const value_type &x, const value_type &y) -> double { - return ((x - y) * (x - y)); - }, - seed_t seed = seed_t(-1)) + KMeansVisitor(size_type num_of_iter, + bool calc_clusters = true, + distance_func f = + [](const value_type &x, const value_type &y) -> double { + return ((x - y) * (x - y)); + }, + seed_t seed = seed_t(-1)) : iter_num_(num_of_iter), cc_(calc_clusters), seed_(seed), @@ -342,21 +350,36 @@ struct AffinityPropVisitor { cluster_type clusters_ { }; // Clusters order_type clusters_idxs_ { }; // Clusters indices + // Symmetric access into the packed (upper-triangular-inclusive) + // similarity array. sim(i, j) == sim(j, i), but get_similarity_() + // only ever physically stores the i <= j half -- a lookup with + // i > j must swap to the row that was actually filled, rather + // than reusing row i's own packed offset (which, for a column + // less than i, silently lands on some unrelated pair's entry). + // + static inline double + sim_at_(const vec_t &simil, long i, long j, long col_s) { + + if (i <= j) + return (simil[(i * col_s) - ((i * (i + 1)) >> 1) + j]); + return (simil[(j * col_s) - ((j * (j + 1)) >> 1) + i]); + } + template inline vec_t get_similarity_(const H &column_begin, long col_s) { vec_t simil((col_s * (col_s + 1)) / 2, 0.0); - double min_dist = std::numeric_limits::max(); + double min_dist { std::numeric_limits::max() }; // Compute similarity between distinct data points i and j // - for (long i = 0; i < col_s - 1; ++i) [[likely]] { - const value_type &i_val = *(column_begin + i); - const long i_idx = (i * col_s) - ((i * (i + 1)) >> 1); + for (long i { 0 }; i < col_s - 1; ++i) [[likely]] { + const value_type &i_val { *(column_begin + i) }; + const long i_idx { (i * col_s) - ((i * (i + 1)) >> 1) }; - for (long j = i + 1; j < col_s; ++j) [[likely]] { - const double dist = -dfunc_(i_val, *(column_begin + j)); + for (long j { i + 1 }; j < col_s; ++j) [[likely]] { + const double dist { -dfunc_(i_val, *(column_begin + j)) }; simil[i_idx + j] = dist; if (dist < min_dist) min_dist = dist; @@ -365,7 +388,7 @@ struct AffinityPropVisitor { // Assign min to diagonals // - for (long i = 0; i < col_s; ++i) + for (long i { 0 }; i < col_s; ++i) simil[(i * col_s) + i - ((i * (i + 1)) >> 1)] = min_dist; return (simil); @@ -380,63 +403,85 @@ struct AffinityPropVisitor { avail.resize(col_s * col_s, 0.0); respon.resize(col_s * col_s, 0.0); - const double one_df = 1.0 - dfactor_; + const double one_df { 1.0 - dfactor_ }; - for (size_type m = 0; m < iter_num_; ++m) [[likely]] { - // Update responsibility + for (size_type m { 0 }; m < iter_num_; ++m) [[likely]] { + // Update responsibility: + // r(i, k) = s(i, k) - max_{k' != k} [a(i, k') + s(i, k')] + // + // For fixed i, "max excluding k" only ever needs the + // overall max and, for whichever k achieves it, the + // second-highest value -- both obtainable in one O(n) + // pass instead of redoing an O(n) search for every k + // (which made this O(n^2) per row / O(n^3) per iteration + // instead of the documented O(n) per row / O(n^2) total). // for (long i = 0; i < col_s; ++i) [[likely]] { - const long i_idx = (i * col_s) - ((i * (i + 1)) >> 1); + double max1 { -std::numeric_limits::max() }; + double max2 { -std::numeric_limits::max() }; + long argmax1 { -1 }; - for (long j = 0; j < col_s; ++j) [[likely]] { - double max_diff = -std::numeric_limits::max(); - - for (long jj = 0; jj < col_s; ++jj) { - if (jj ^ j) [[likely]] { - const double value = - simil[i_idx + jj] + avail[jj * col_s + i]; + for (long jj { 0 }; jj < col_s; ++jj) [[likely]] { + const double value { + sim_at_(simil, i, jj, col_s) + avail[jj * col_s + i] + }; - if (value > max_diff) - max_diff = value; - } + if (value > max1) { + max2 = max1; + max1 = value; + argmax1 = jj; } + else if (value > max2) + max2 = value; + } - const long j_idx = j * col_s + i; + for (long j { 0 }; j < col_s; ++j) [[likely]] { + const double max_diff { (j == argmax1) ? max2 : max1 }; + const long j_idx { j * col_s + i }; - respon[j_idx] = one_df * (simil[i_idx + j] - max_diff) + - dfactor_ * respon[j_idx]; + respon[j_idx] = + one_df * (sim_at_(simil, i, j, col_s) - max_diff) + + dfactor_ * respon[j_idx]; } } // Update availability // Do diagonals first // - for (long i = 0; i < col_s; ++i) [[likely]] { - const long s1 = i * col_s; - const long s2 = i * col_s + i; - double sum = 0.0; + for (long i { 0 }; i < col_s; ++i) [[likely]] { + const long s1 { i * col_s }; + const long s2 { i * col_s + i }; + double sum { 0.0 }; - for (long ii = 0; ii < col_s; ++ii) [[likely]] + for (long ii { 0 }; ii < col_s; ++ii) [[likely]] if (ii ^ i) sum += std::max(0.0, respon[s1 + ii]); avail[s2] = one_df * sum + dfactor_ * avail[s2]; } - for (long i = 0; i < col_s; ++i) [[likely]] { - for (long j = 0; j < col_s; ++j) [[likely]] { + + // Off-diagonals: a(i, j) needs + // sum_{i' not in {i, j}} max(0, r(i', j)) for every i, with + // j fixed. Rather than rescanning all n indices for every + // (i, j) pair (O(n) work x O(n^2) pairs = O(n^3) per + // iteration), compute the column total once per j -- + // sum_{i' != j} max(0, r(i', j)), O(n) -- and reuse it for + // every i by subtracting out just the i'=i term. + // + for (long j { 0 }; j < col_s; ++j) [[likely]] { + const long s1 { j * col_s }; + double total_j { 0.0 }; + + for (long ii { 0 }; ii < col_s; ++ii) [[likely]] + if (ii ^ j) + total_j += std::max(0.0, respon[s1 + ii]); + + for (long i { 0 }; i < col_s; ++i) [[likely]] { if (i ^ j) [[likely]] { // Not equal - const long s1 = j * col_s; - const long s2 = j * col_s + i; - double sum = 0.0; - const long max_i_j = std::max(i, j); - const long min_i_j = std::min(i, j); - - for (long ii = 0; ii < min_i_j; ++ii) - sum += std::max(0.0, respon[s1 + ii]); - for (long ii = min_i_j + 1; ii < max_i_j; ++ii) - sum += std::max(0.0, respon[s1 + ii]); - for (long ii = max_i_j + 1; ii < col_s; ++ii) - sum += std::max(0.0, respon[s1 + ii]); + const long s2 { j * col_s + i }; + const double sum { + total_j - std::max(0.0, respon[s1 + i]) + }; avail[s2] = one_df * @@ -454,26 +499,26 @@ struct AffinityPropVisitor { inline void calc_clusters_(const H &column_begin, long col_s) { - const long centers_size = result_.size(); + const long centers_size { long(result_.size()) }; if (! centers_size) return; - const auto resv = col_s / centers_size; + const auto resv { col_s / centers_size }; clusters_.resize(centers_size); clusters_idxs_.resize(centers_size); - for (long i = 0; i < centers_size; ++i) { + for (long i { 0 }; i < centers_size; ++i) { clusters_[i].reserve(resv); clusters_idxs_[i].reserve(resv); } - for (long j = 0; j < col_s; ++j) [[likely]] { - const value_type &j_val = *(column_begin + j); - double min_dist = dfunc_(j_val, result_[0]); - long min_idx = 0; + for (long j { 0 }; j < col_s; ++j) [[likely]] { + const value_type &j_val { *(column_begin + j) }; + double min_dist { dfunc_(j_val, result_[0]) }; + long min_idx { 0 }; - for (long i = 1; i < centers_size; ++i) { - const double dist = dfunc_(j_val, result_[i]); + for (long i { 1 }; i < centers_size; ++i) { + const double dist { dfunc_(j_val, result_[i]) }; if (dist < min_dist) { min_dist = dist; @@ -489,11 +534,13 @@ struct AffinityPropVisitor { template inline void - operator() (const IV &idx_begin, const IV &idx_end, - const H &column_begin, const H &column_end) { + operator()(const IV &idx_begin, const IV &idx_end, + const H &column_begin, const H &column_end) { - const long col_s = std::min(std::distance(idx_begin, idx_end), - std::distance(column_begin, column_end)); + const long col_s { + long(std::min(std::distance(idx_begin, idx_end), + std::distance(column_begin, column_end))) + }; const vec_t simil = get_similarity_(column_begin, col_s); vec_t avail; @@ -502,8 +549,8 @@ struct AffinityPropVisitor { get_avail_and_respon_(simil, col_s, avail, respon); result_.reserve(std::min(col_s / 100, long(16))); - for (long i = 0; i < col_s; ++i) [[likely]] { - const long idx = i * col_s + i; + for (long i { 0 }; i < col_s; ++i) [[likely]] { + const long idx { i * col_s + i }; if (respon[idx] + avail[idx] > 0.0) result_.push_back(&*(column_begin + i)); @@ -512,29 +559,29 @@ struct AffinityPropVisitor { if (cc_) calc_clusters_(column_begin, col_s); } - inline void pre () { + inline void pre() { result_.clear(); clusters_.clear(); clusters_idxs_.clear(); } - inline void post () { } - inline const result_type &get_result () const { return (result_); } - inline result_type &get_result () { return (result_); } - inline const cluster_type &get_clusters () const { return (clusters_); } - inline cluster_type &get_clusters () { return (clusters_); } + inline void post() { } + inline const result_type &get_result() const { return (result_); } + inline result_type &get_result() { return (result_); } + inline const cluster_type &get_clusters() const { return (clusters_); } + inline cluster_type &get_clusters() { return (clusters_); } inline const order_type & - get_clusters_idxs () const { return (clusters_idxs_); } + get_clusters_idxs() const { return (clusters_idxs_); } explicit - AffinityPropVisitor( - size_type num_of_iter, - bool calc_clusters = true, - distance_func f = - [](const value_type &x, const value_type &y) -> double { - return ((x - y) * (x - y)); - }, - double damping_factor = 0.9) + AffinityPropVisitor(size_type num_of_iter, + bool calc_clusters = true, + distance_func f = + [](const value_type &x, + const value_type &y) -> double { + return ((x - y) * (x - y)); + }, + double damping_factor = 0.9) : iter_num_(num_of_iter), cc_(calc_clusters), dfactor_(damping_factor), dfunc_(f) { } @@ -550,7 +597,7 @@ struct DBSCANVisitor { private: - static constexpr bool is_md_ = random_acc_cont; + static constexpr bool is_md_ { random_acc_cont }; public: @@ -601,10 +648,10 @@ struct DBSCANVisitor { const id_t col_s, vec_t &cluster_index) { - const value_type &value = *(column_begin + column_idx); + const value_type &value { *(column_begin + column_idx) }; cluster_index.clear(); - for (id_t i = 0; i < col_s; ++i) { + for (id_t i { 0 }; i < col_s; ++i) { if (dfunc_(value, *(column_begin + i)) <= max_dist_) cluster_index.push_back(i); } @@ -622,8 +669,8 @@ struct DBSCANVisitor { calculate_cluster_(column_begin, column_idx, col_s, seeds); - const id_t seeds_s = id_t(seeds.size()); - const value_type &value = *(column_begin + column_idx); + const id_t seeds_s { id_t(seeds.size()) }; + const value_type &value { *(column_begin + column_idx) }; if (seeds_s < min_mems_) { cluster_ids[column_idx] = NOISE; @@ -632,8 +679,8 @@ struct DBSCANVisitor { id_t core_index { 0 }; - for (id_t i = 0; i < seeds_s; ++i) { - const auto seed_val = seeds[i]; + for (id_t i { 0 }; i < seeds_s; ++i) { + const auto seed_val { seeds[i] }; cluster_ids[seed_val] = cluster_id; if (*(column_begin + seed_val) == value) [[unlikely]] @@ -641,12 +688,12 @@ struct DBSCANVisitor { } seeds.erase(seeds.begin() + core_index); - for (id_t i = 0, n = seeds.size(); i < n; ++i) { + for (id_t i { 0 }, n { id_t(seeds.size()) }; i < n; ++i) { calculate_cluster_(column_begin, seeds[i], col_s, cluster_neighors); if (id_t(cluster_neighors.size()) >= min_mems_) { - for (id_t j = 0; j < id_t(cluster_neighors.size()); ++j) { - auto &cluster_val = cluster_ids[j]; + for (id_t j { 0 }; j < id_t(cluster_neighors.size()); ++j) { + auto &cluster_val { cluster_ids[cluster_neighors[j]] }; if (cluster_val < 0) { // NOISE or UNCLASSIFIED if (cluster_val == UNCLASSIFIED) { @@ -666,19 +713,21 @@ struct DBSCANVisitor { template inline void - operator() (const IV &idx_begin, const IV &idx_end, - const H &column_begin, const H &column_end) { + operator()(const IV &idx_begin, const IV &idx_end, + const H &column_begin, const H &column_end) { - const id_t col_s = std::min(std::distance(idx_begin, idx_end), - std::distance(column_begin, column_end)); - vec_t cluster_ids (col_s, UNCLASSIFIED); + const id_t col_s { + id_t(std::min(std::distance(idx_begin, idx_end), + std::distance(column_begin, column_end))) + }; + vec_t cluster_ids(col_s, UNCLASSIFIED); vec_t seeds; vec_t cluster_neighors; id_t cluster_id { 0 }; seeds.reserve(col_s / 20); cluster_neighors.reserve(col_s / 20); - for (id_t i = 0; i < col_s; ++i) { + for (id_t i { 0 }; i < col_s; ++i) { if (cluster_ids[i] == UNCLASSIFIED && expand_cluster_(column_begin, i, @@ -691,17 +740,23 @@ struct DBSCANVisitor { } } - const auto resv = col_s / cluster_id; + // cluster_id stays 0 whenever every point ends up NOISE (a + // legitimate outcome for sparse data or a tight max_dist_/high + // min_mems_, not just a contrived edge case) -- guard against + // dividing by it. The value is moot when cluster_id is 0 anyway, + // since the per-cluster reserve loop right below never runs. + // + const auto resv { cluster_id > 0 ? col_s / cluster_id : 0 }; clusters_.resize(cluster_id); clusters_idxs_.resize(cluster_id); noisey_idxs_.reserve(std::max(id_t(8), id_t(col_s / 500))); - for (long i = 0; i < cluster_id; ++i) { + for (long i { 0 }; i < cluster_id; ++i) { clusters_[i].reserve(resv); clusters_idxs_[i].reserve(resv); } - for (id_t i = 0; i < col_s; ++i) { - const auto this_id = cluster_ids[i]; + for (id_t i { 0 }; i < col_s; ++i) { + const auto this_id { cluster_ids[i] }; if (this_id >= 0) [[likely]] { clusters_[this_id].push_back(&(*(column_begin + i))); @@ -713,19 +768,19 @@ struct DBSCANVisitor { inline void set_dist_func(distance_func &&f) { dfunc_ = f; } - inline void pre () { + inline void pre() { clusters_.clear(); clusters_idxs_.clear(); noisey_idxs_.clear(); } - inline void post () { } + inline void post() { } - inline const result_type &get_result () const { return (clusters_); } + inline const result_type &get_result() const { return (clusters_); } inline const order_type & - get_clusters_idxs () const { return (clusters_idxs_); } + get_clusters_idxs() const { return (clusters_idxs_); } inline const vec_t & - get_noisey_idxs () const { return (noisey_idxs_); } + get_noisey_idxs() const { return (noisey_idxs_); } DBSCANVisitor(id_t min_mems, double max_dist) : min_mems_(min_mems), @@ -752,7 +807,7 @@ struct MeanShiftVisitor { private: - static constexpr bool is_md_ = random_acc_cont; + static constexpr bool is_md_ { random_acc_cont }; public: @@ -808,21 +863,21 @@ struct MeanShiftVisitor { inline static double biweight_kernel_(double d) { - const auto x = 1.0 - d * d; + const auto x { 1.0 - d * d }; return (d <= 1.0 ? x * x : 0.0); } inline static double triweight_kernel_(double d) { - const auto x = 1.0 - d * d; + const auto x { 1.0 - d * d }; return (d <= 1.0 ? x * x * x : 0.0); } inline static double tricube_kernel_(double d) { - const auto x = 1.0 - d * d * d; + const auto x { 1.0 - d * d * d }; return (d <= 1.0 ? x * x * x : 0.0); } @@ -849,24 +904,34 @@ struct MeanShiftVisitor { inline static double silverman_kernel_(double d) { - const auto x = M_SQRT1_2 * std::abs(d); + const auto x { M_SQRT1_2 * std::abs(d) }; return (std::exp(-x) * std::sin(x + M_PI_4)); } - template - inline void shift_(const H &column_begin, - size_type index, - const value_type &val, - vec_t &shifted, - vec_t &shifting) { + // Convergence is "did this iteration's shift move the point much + // from where IT was", not "is the point still near where it + // started" -- comparing against the original, untouched data + // (column_begin) instead of the point's own previous position + // meant a point could freeze the moment its very first shift + // happened to land within max_dist_ of its start, regardless of + // whether it had reached its actual density mode. The freshly + // computed position is always the best estimate available, so it + // is recorded unconditionally -- previously it was discarded + // outright on the very iteration convergence was detected. + // + inline void + shift_(size_type index, + const value_type &val, + vec_t &shifted, + vec_t &shifting) { - if (dfunc_(val, *(column_begin + index)) <= max_dist_) + if (dfunc_(val, shifted[index]) <= max_dist_) shifting[index] = false; - else - shifted[index] = val; + shifted[index] = val; } + template inline void build_cluster_(const H &column_begin, @@ -881,10 +946,10 @@ struct MeanShiftVisitor { clusters_.reserve(32); clusters_idxs_.reserve(32); for (size_type i { 0 }; i < shifted.size(); ++i) { - const auto &shifted_val = shifted[i]; - auto cbegin = clusters_.begin(); - auto cend = clusters_.end(); - auto ibegin = clusters_idxs_.begin(); + const auto &shifted_val { shifted[i] }; + auto cbegin { clusters_.begin() }; + auto cend { clusters_.end() }; + auto ibegin { clusters_idxs_.begin() }; size_type cnt_idx { 0 }; while (cbegin != cend) { @@ -917,12 +982,13 @@ struct MeanShiftVisitor { template inline void - operator() (const IV &idx_begin, const IV &idx_end, - const H &column_begin, const H &column_end) { + operator()(const IV &idx_begin, const IV &idx_end, + const H &column_begin, const H &column_end) { - const size_type col_s = - std::min(std::distance(idx_begin, idx_end), - std::distance(column_begin, column_end)); + const size_type col_s { + size_type(std::min(std::distance(idx_begin, idx_end), + std::distance(column_begin, column_end))) + }; auto k_func = (kernel_ == mean_shift_kernel::uniform) ? &uniform_kernel_ : (kernel_ == mean_shift_kernel::triangular) ? &triangular_kernel_ @@ -935,11 +1001,23 @@ struct MeanShiftVisitor { : (kernel_ == mean_shift_kernel::logistic) ? &logistic_kernel_ : (kernel_ == mean_shift_kernel::sigmoid) ? &sigmoid_kernel_ : &silverman_kernel_; - vec_t shifted (column_begin, column_end); - vec_t shifting (col_s, true); + vec_t shifted(column_begin, column_end); + vec_t shifting(col_s, true); size_type iterations { 0 }; - const double radius { kband_ * 3.0 }; - const double dbl_sq_bw { 2.0 * kband_ * kband_ }; + + // The kernel functions (uniform_kernel_, gaussian_kernel_, etc.) + // are all written for a bandwidth-normalized argument -- the + // d <= 1.0 cutoffs, and gaussian_kernel_'s textbook + // exp(-0.5*d^2) shape, only mean what they look like they + // mean when d is "how many bandwidths away", not a raw + // distance. dfunc_ returns squared distance for scalar T and + // true (already-rooted) Euclidean distance for MD T, so both + // need converting to a true linear distance before dividing + // by kband_. 3.0 (bandwidths) replaces the old radius, which + // compared an un-normalized, type-inconsistent raw distance + // against kband_*3.0. + // + constexpr double radius { 3.0 }; while (iterations++ < max_iter_ && std::any_of(shifting.begin(), shifting.end(), @@ -951,12 +1029,18 @@ struct MeanShiftVisitor { const value_type &val_to_shift { shifted[i] }; double total_w { 0 }; - for (size_type j = 0; j < col_s; ++j) { - const value_type &this_val = *(column_begin + j); - const double dist = dfunc_(val_to_shift, this_val); + for (size_type j { 0 }; j < col_s; ++j) { + const value_type &this_val { *(column_begin + j) }; + const double raw_dist { + dfunc_(val_to_shift, this_val) + }; + const double lin_dist { + is_md_ ? raw_dist : std::sqrt(raw_dist) + }; + const double norm_dist { lin_dist / kband_ }; - if (dist <= radius) { - const double weight = k_func(dist) / dbl_sq_bw; + if (norm_dist <= radius) { + const double weight { k_func(norm_dist) }; new_val = new_val + (this_val * weight); total_w += weight; @@ -967,7 +1051,7 @@ struct MeanShiftVisitor { // its neighbors // new_val = new_val / total_w; - shift_(column_begin, i, new_val, shifted, shifting); + shift_(i, new_val, shifted, shifting); } } @@ -976,12 +1060,12 @@ struct MeanShiftVisitor { inline void set_dist_func(distance_func &&f) { dfunc_ = f; } - inline void pre () { clusters_.clear(); clusters_idxs_.clear(); } - inline void post () { } + inline void pre() { clusters_.clear(); clusters_idxs_.clear(); } + inline void post() { } - inline const result_type &get_result () const { return (clusters_); } + inline const result_type &get_result() const { return (clusters_); } inline const order_type & - get_clusters_idxs () const { return (clusters_idxs_); } + get_clusters_idxs() const { return (clusters_idxs_); } MeanShiftVisitor(double kernel_bandwidth, double max_dist, @@ -1734,8 +1818,8 @@ struct EntropyVisitor { if constexpr (! is_md_) { for (size_type i { 0 }; i < roll_count_ - 1; ++i) [[likely]] result[i] = get_nan(); - for (size_type i { roll_count_ - 1 }; i < sz; ++i) - result[i] = sum_v.get_result()[i]; + for (size_type i { 0 }; i < sz; ++i) + result[i + (roll_count_ - 1)] = sum_v.get_result()[i]; } else { const std::vector nans(column_begin->size(), @@ -1743,12 +1827,13 @@ struct EntropyVisitor { for (size_type i { 0 }; i < roll_count_ - 1; ++i) [[likely]] result[i] = nans; - for (size_type i { roll_count_ - 1 }; i < sz; ++i) { + for (size_type i { 0 }; i < sz; ++i) { if (sum_v.get_result()[i].empty()) - result[i].resize(column_begin->size(), get_nan()); + result[i + (roll_count_ - 1)] = nans; else - result[i].assign(sum_v.get_result()[i].begin(), - sum_v.get_result()[i].end()); + result[i + (roll_count_ - 1)].assign( + sum_v.get_result()[i].begin(), + sum_v.get_result()[i].end()); } } @@ -1966,7 +2051,7 @@ struct SigmoidVisitor { private: - static constexpr bool is_md_ = random_acc_cont; + static constexpr bool is_md_ { random_acc_cont }; using data_t = typename std::conditional_t void { - for (size_type i { begin }; i < end; ++i) - result_[i] = - data_t(1) / - _bc_sqrt_(data_t(1) + _bc_pow_(*(column_begin + i), - data_t(2))); + for (size_type i { begin }; i < end; ++i) { + const auto &val { *(column_begin + i) }; + const auto denom { + _bc_sqrt_(data_t(1) + _bc_pow_(val, data_t(2))) + }; + + if constexpr (is_md_) { + result_[i].resize(val.size()); + for (size_type d { 0 }; d < val.size(); ++d) + result_[i][d] = val[d] / denom[d]; + } + else + result_[i] = val / denom; + } }; apply_func_(col_s, thread_level, std::move(lbd)); @@ -2188,8 +2282,8 @@ struct RectifyVisitor { template inline void - operator() (const K &idx_begin, const K &idx_end, - const H &column_begin, const H &column_end) { + operator()(const K &idx_begin, const K &idx_end, + const H &column_begin, const H &column_end) { GET_COL_SIZE2 @@ -2205,7 +2299,7 @@ struct RectifyVisitor { col_s, [&column_begin, this] (auto begin, auto end) -> void { - for (size_type i = begin; i < end; ++i) + for (size_type i { begin }; i < end; ++i) this->result_[i] = std::max(T(0), *(column_begin + i)); }); @@ -2217,7 +2311,7 @@ struct RectifyVisitor { col_s, [&column_begin, this] (auto begin, auto end) -> void { - for (size_type i = begin; i < end; ++i) { + for (size_type i { begin }; i < end; ++i) { const value_type v = *(column_begin + i); this->result_[i] = @@ -2232,7 +2326,7 @@ struct RectifyVisitor { col_s, [&column_begin, this] (auto begin, auto end) -> void { - for (size_type i = begin; i < end; ++i) { + for (size_type i { begin }; i < end; ++i) { const value_type v = *(column_begin + i); this->result_[i] = @@ -2252,7 +2346,7 @@ struct RectifyVisitor { col_s, [&column_begin, &sigm = std::as_const(sigm), this] (auto begin, auto end) -> void { - for (size_type i = begin; i < end; ++i) { + for (size_type i { begin }; i < end; ++i) { const value_type col = *(column_begin + i); const value_type sig = sigm.get_result()[i]; @@ -2267,7 +2361,7 @@ struct RectifyVisitor { col_s, [&column_begin, this] (auto begin, auto end) -> void { - for (size_type i = begin; i < end; ++i) { + for (size_type i { begin }; i < end; ++i) { const value_type v = *(column_begin + i); this->result_[i] = softp_(v, this->param_); @@ -2281,7 +2375,7 @@ struct RectifyVisitor { col_s, [&column_begin, this] (auto begin, auto end) -> void { - for (size_type i = begin; i < end; ++i) { + for (size_type i { begin }; i < end; ++i) { const value_type v = *(column_begin + i); if (v > 0) @@ -2299,7 +2393,7 @@ struct RectifyVisitor { col_s, [&column_begin, this] (auto begin, auto end) -> void { - for (size_type i = begin; i < end; ++i) { + for (size_type i { begin }; i < end; ++i) { const value_type v = *(column_begin + i); this->result_[i] = @@ -2314,7 +2408,7 @@ struct RectifyVisitor { col_s, [&column_begin, this] (auto begin, auto end) -> void { - for (size_type i = begin; i < end; ++i) { + for (size_type i { begin }; i < end; ++i) { const value_type v = *(column_begin + i); this->result_[i] = @@ -2397,12 +2491,12 @@ struct RectifyVisitor { OBO_PORT_OPT - inline void pre () { + inline void pre() { OBO_PORT_PRE result_.clear(); } - inline void post () { OBO_PORT_POST } + inline void post() { OBO_PORT_POST } DEFINE_RESULT explicit @@ -2414,15 +2508,18 @@ struct RectifyVisitor { inline static value_type softp_(const value_type &v, const value_type &p) { - return(std::log(T(1) + std::exp(p * v)) / p); + const value_type y { p * v }; + + return ((std::max(y, T(0)) + std::log1p(std::exp(-std::fabs(y)))) / p); } inline static value_type standard_normal_dist_(const value_type &v) { - static constexpr value_type two = 2; - static const value_type sqrt_dbl_pi = std::sqrt(two * M_PI); + static const value_type inv_sqrt2 { + value_type(1) / std::sqrt(value_type(2)) + }; - return (std::exp(-(v * v) / two) / sqrt_dbl_pi); + return ((value_type(1) + std::erf(v * inv_sqrt2)) / value_type(2)); } OBO_PORT_DECL @@ -2442,7 +2539,7 @@ struct PolicyLearningLossVisitor { private: - static constexpr bool is_md_ = random_acc_cont; + static constexpr bool is_md_ { random_acc_cont }; using data_t = typename std::conditional_tsize() != dim || (reward_begin + i)->size() != dim) - throw DataFrameError( - "PolicyLearningLossVisitor: " - "Inconsistent data dimensions"); + throw DataFrameError("PolicyLearningLossVisitor: " + "Inconsistent data dimensions"); } #endif // HMDF_SANITY_EXCEPTIONS @@ -2544,19 +2640,25 @@ struct PolicyLearningLossVisitor { for (size_type i { begin }; i < end; ++i) { const auto &ap { *(action_prob_begin + i) }; const auto &r { *(reward_begin + i) }; - const auto &adjusted_r { (r - adjust) / scale }; + const auto adjusted_r { + (r - adjust) / (scale + data_t(epsilon_)) + }; + + result_[i] = + (_bc_log_(ap + data_t(epsilon_)) * data_t(-1)) * + adjusted_r; - result_[i] = (_bc_log_(ap) * data_t(-1)) * adjusted_r; } }; if (col_s >= ThreadPool::MUL_THR_THHOLD && ThreadGranularity::get_thread_level() > 2) { - auto futures = + auto futures { ThreadGranularity::thr_pool_.parallel_loop( size_type(0), col_s, - std::move(lbd)); + std::move(lbd)) + }; for (auto &fut : futures) fut.get(); } @@ -2595,7 +2697,7 @@ struct LossFunctionVisitor { private: - static constexpr bool is_md_ = random_acc_cont; + static constexpr bool is_md_ { random_acc_cont }; using data_t = typename std::conditional_tsize() != dim || (mit++)->size() != dim) throw DataFrameError("LossFunctionVisitor: " "Inconsistent data dimensions"); @@ -3001,7 +3103,7 @@ struct VectorSimilarityVisitor { private: - static constexpr bool is_md_ = random_acc_cont; + static constexpr bool is_md_ { random_acc_cont }; using data_t = typename std::conditional_t dot_v; dot_v.pre(); - dot_v (idx_begin, idx_end, - column_begin1, column_end1, column_begin2, column_end2); + dot_v(idx_begin, idx_end, + column_begin1, column_end1, column_begin2, column_end2); dot_v.post(); const auto dotp_res { dot_v.get_result() }; @@ -3086,8 +3188,9 @@ struct VectorSimilarityVisitor { if (! is_md_) { #ifdef HMDF_SANITY_EXCEPTIONS if (col_s1 != col_s2) - throw DataFrameError("VectorSimilarityVisitor: " - "All columns must be of equal sizes"); + throw DataFrameError( + "VectorSimilarityVisitor: " + "All columns must be of equal sizes"); #endif // HMDF_SANITY_EXCEPTIONS // Must normalize the dot product first. @@ -3143,10 +3246,10 @@ struct VectorSimilarityVisitor { "All columns must be of equal sizes"); #endif // HMDF_SANITY_EXCEPTIONS - const equal_t eq { }; + const equal_t eq_perd { }; for (size_type i { 0 }; i < col_s1; ++i) - if (! eq(*(column_begin1 + i), *(column_begin2 + i))) + if (! eq_perd(*(column_begin1 + i), *(column_begin2 + i))) result_ += 1; } } @@ -3431,7 +3534,7 @@ struct SeasonalPeriodVisitor { private: - static constexpr bool is_md_ = random_acc_cont; + static constexpr bool is_md_ { random_acc_cont }; using data_t = typename std::conditional_t; - template + template inline void operator()(const K &idx_begin, const K &idx_end, const H &column_begin, const H &column_end) { @@ -3454,7 +3557,7 @@ struct SeasonalPeriodVisitor { const size_type col_s { size_type(std::distance(column_begin, column_end)) }; - std::vector data (column_begin, column_end); + std::vector data(column_begin, column_end); #ifdef HMDF_SANITY_EXCEPTIONS if constexpr (is_md_) { @@ -3468,7 +3571,7 @@ struct SeasonalPeriodVisitor { #endif // HMDF_SANITY_EXCEPTIONS if (params_.detrend) { // Take trend out - std::vector xvals (col_s); + std::vector xvals(col_s); size_type xvalue { 0 }; for (auto &val : xvals) { @@ -3479,13 +3582,13 @@ struct SeasonalPeriodVisitor { LowessVisitor l_v { params_.num_loops, params_.frac, - params_.delta, // * value_type(col_s), + params_.delta, true }; l_v.pre(); - l_v (idx_begin, idx_end, - data.begin(), data.end(), xvals.begin(), xvals.end()); + l_v(idx_begin, idx_end, + data.begin(), data.end(), xvals.begin(), xvals.end()); l_v.post(); if constexpr (! is_md_) { @@ -3524,7 +3627,7 @@ struct SeasonalPeriodVisitor { FastFourierTransVisitor fft; fft.pre(); - fft (idx_begin, idx_end, data.begin(), data.end()); + fft(idx_begin, idx_end, data.begin(), data.end()); fft.post(); // mags is a vector in case of scalar input @@ -3536,14 +3639,12 @@ struct SeasonalPeriodVisitor { const size_type n_bins { mags.size() }; // For a real-valued input of length N, the FFT produces a - // symmetric spectrum — bin k and bin N-k carry identical - // magnitude. - // FIX: Restrict the scan to [1, N/2) — the unique, non-mirrored - // half of the spectrum. + // symmetric spectrum — (bin k and bin N - k carry identical + // magnitude). + // FIX: Restrict the scan to [1, N / 2) — the unique, non-mirrored + // half of the spectrum. // - const size_type scan_len { - ((n_bins & 0x01) == 0) ? n_bins / 2 : (n_bins + 1) / 2 - }; + const size_type scan_len { n_bins / 2 + 1 }; // Skip bin 0 (DC component) — it carries no frequency information // and would cause dom_freq_ = 0 and result_ = 1/0 = inf. @@ -3560,17 +3661,17 @@ struct SeasonalPeriodVisitor { dom_freq_ = result_type(dom_idx_) * result_type(params_.sampling_rate) / result_type(mags.size()); - result_ = result_type(1) / dom_freq_; + result_ = result_type(1) / dom_freq_; } else { const size_type dim { size_type(mags.cols()) }; const size_type n_bins { size_type(mags.rows()) }; // For a real-valued input of length N, the FFT produces a - // symmetric spectrum — bin k and bin N-k carry identical - // magnitude. - // FIX: Restrict the scan to [1, N/2) — the unique, non-mirrored - // half of the spectrum. + // symmetric spectrum (bin k and bin N - k carry identical + // magnitude). + // FIX: Restrict the scan to [1, N / 2) — the unique, non-mirrored + // half of the spectrum. // const size_type scan_len { ((n_bins & 0x01) == 0) ? n_bins / 2 : (n_bins + 1) / 2 @@ -3657,10 +3758,10 @@ struct SeasonalPeriodVisitor { // Per-dimension storage — populated only when is_md_ == true. // - md_result_type result_vec_ { }; - md_result_type max_mag_vec_ { }; + md_result_type result_vec_ { }; + md_result_type max_mag_vec_ { }; md_result_type dom_freq_vec_ { }; - std::vector dom_idx_vec_ { }; + std::vector dom_idx_vec_ { }; }; template @@ -3676,7 +3777,7 @@ struct DynamicTimeWarpVisitor { private: - static constexpr bool is_md_ = random_acc_cont; + static constexpr bool is_md_ { random_acc_cont }; using data_t = typename std::conditional_t; + static constexpr bool is_md_ { random_acc_cont }; using data_t = typename std::conditional_t, vec_t>>; + template inline void operator()(const K &idx_begin, const K &idx_end, @@ -3864,10 +3966,29 @@ struct AnomalyDetectByFFTVisitor { auto fft_res = std::move(fft.get_result()); if constexpr (! is_md_) { - // Zero out high frequencies + // Zero out high frequencies -- but a real signal's FFT is + // conjugate-symmetric (bin k and bin N-k are mirrors), so + // zeroing everything from freq_num_ to the end also wipes + // out the mirror partners (indices N-freq_num_+1 .. N-1) of + // the very low-frequency bins we're trying to keep. Without + // its conjugate partner, each kept bin's contribution to the + // inverse FFT loses the term that makes it real and full- + // amplitude -- the reconstructed signal comes back at + // exactly half the true amplitude for every non-DC + // frequency (confirmed: a clean single sinusoid reconstructs + // at a flat 0.5x ratio, pre-fix). Zero only the middle, + // mirror-safe range [freq_num_, col_s-freq_num_+1) instead, + // leaving both the low bins and their mirrors intact. If + // freq_num_ is large enough that this range would be empty + // or invalid (the two kept halves would overlap), there's + // nothing to zero -- keep everything. // - std::fill(fft_res.begin() + freq_num_, fft_res.end(), - typename decltype(fft_res)::value_type { }); + const size_type mirror_bound { col_s - freq_num_ + 1 }; + + if (freq_num_ < mirror_bound) + std::fill(fft_res.begin() + freq_num_, + fft_res.begin() + mirror_bound, + typename decltype(fft_res)::value_type { }); // Inverse FFT: input is vec, fft_v> picks // the scalar path because complex satisfies IS_SCALAR @@ -3884,8 +4005,23 @@ struct AnomalyDetectByFFTVisitor { result_.reserve(((col_s / 40) < 32) ? size_type(32) : col_s / 40); for (size_type i { 0 }; i < col_s; ++i) { + // When normalization is enabled, the reconstruction was + // built from the NORMALIZED data -- comparing it against + // the raw, un-normalized column value is a scale + // mismatch that produces a meaningless residual unless + // the two scales happen to coincide by chance (confirmed: + // a clean sinusoid with a large offset/amplitude under + // z-score normalization got EVERY point flagged as an + // "anomaly", since the raw values are ~1000s while the + // normalized reconstruction is on a unit scale). Compare + // against the same-scale (normalized) value instead. + // + const data_t orig_val { + (nt_ > normalization_type::none) + ? norm.get_result()[i] : *(column_begin + i) + }; const data_t residual { - std::abs(*(column_begin + i) - ifft_res[i].real()) + std::abs(orig_val - ifft_res[i].real()) }; if (residual > ath_) result_.push_back(i); @@ -3913,11 +4049,18 @@ struct AnomalyDetectByFFTVisitor { const size_type dim { size_type(fft_res.cols()) }; - // Zero out high frequencies across all dims + // Zero out high frequencies across all dims -- same + // conjugate-mirror issue as the scalar path above: only zero + // the middle, mirror-safe range so each dimension's kept low + // bins keep their conjugate partners. // - for (size_type i { freq_num_ }; i < col_s; ++i) - for (size_type d { 0 }; d < dim; ++d) - fft_res(i, d) = typename decltype(fft_res)::value_type { }; + const size_type mirror_bound { col_s - freq_num_ + 1 }; + + if (freq_num_ < mirror_bound) + for (size_type i { freq_num_ }; i < mirror_bound; ++i) + for (size_type d { 0 }; d < dim; ++d) + fft_res(i, d) = + typename decltype(fft_res)::value_type{ }; // We need cplx_t which is fft_v::cplx_t. // Derive it from the matrix element type to stay consistent. @@ -3961,11 +4104,16 @@ struct AnomalyDetectByFFTVisitor { } // Anomaly detection: compare original vs reconstructed per - // (sample, dim) + // (sample, dim). Same scale-mismatch fix as the scalar path: + // when normalized, compare against the normalized sample, + // not the raw one. // result_.reserve(((col_s / 40) < 32) ? size_type(32) : col_s / 40); for (size_type i { 0 }; i < col_s; ++i) { - const auto &sample { *(column_begin + i) }; + const auto &sample { + (nt_ > normalization_type::none) + ? norm.get_result()[i] : *(column_begin + i) + }; for (size_type d { 0 }; d < dim; ++d) { const data_t residual { @@ -4038,11 +4186,11 @@ struct HampelFilterVisitor { aggr(idx_begin, idx_end, diff.begin(), diff.end()); aggr.post(); - const value_type factor = num_of_std_ * unbiased_factor_; - const auto &aggr_res = aggr.get_result(); + const value_type factor { num_of_std_ * unbiased_factor_ }; + const auto &aggr_res { aggr.get_result() }; result_.reserve(diff.size() / 10); - for (size_type i = 0; i < col_s; ++i) { + for (size_type i { 0 }; i < col_s; ++i) { if (diff[i] > (aggr_res[i] * factor)) result_.push_back(i); } @@ -4052,16 +4200,16 @@ struct HampelFilterVisitor { template inline void - operator() (K idx_begin, K idx_end, H column_begin, H column_end) { + operator()(K idx_begin, K idx_end, H column_begin, H column_end) { if (type_ == hampel_type::median) hampel_(idx_begin, idx_end, column_begin, column_end, SimpleRollAdopter, T, I> - (MedianVisitor { }, window_size_)); + (MedianVisitor{ true }, window_size_)); else if (type_ == hampel_type::mean) hampel_(idx_begin, idx_end, column_begin, column_end, SimpleRollAdopter, T, I> - (MeanVisitor { true }, window_size_)); + (MeanVisitor{ true }, window_size_)); } DEFINE_PRE_POST @@ -4101,29 +4249,29 @@ struct AnomalyDetectByIQRVisitor { template inline void - operator() (const K &, const K &, - const H &column_begin, const H &column_end) { + operator()(const K &, const K &, + const H &column_begin, const H &column_end) { GET_COL_SIZE2 #ifdef HMDF_SANITY_EXCEPTIONS if (col_s < 20) - throw DataFrameError("AnomalyDetectByIQRVisitor: " - "Time-series is too short"); + throw DataFrameError( + "AnomalyDetectByIQRVisitor: Time-series is too short"); #endif // HMDF_SANITY_EXCEPTIONS - vec_t data(column_begin, column_end); - - const auto thread_level = (col_s < ThreadPool::MUL_THR_THHOLD) - ? 0L : ThreadGranularity::get_thread_level(); + vec_t data(column_begin, column_end); + const auto thread_level { + (col_s < ThreadPool::MUL_THR_THHOLD) + ? 0L : ThreadGranularity::get_thread_level() + }; if (thread_level > 2) - ThreadGranularity::thr_pool_.parallel_sort( - data.begin(), data.end()); + ThreadGranularity::thr_pool_.parallel_sort(data.begin(), data.end()); else std::sort(data.begin(), data.end()); - const size_type mid = col_s / 2; + const size_type mid { col_s / 2 }; const value_type q1 { median_(data.begin(), data.begin() + mid) }; const value_type q3 { (col_s & size_type(0x01)) @@ -4136,7 +4284,7 @@ struct AnomalyDetectByIQRVisitor { result_.reserve(32); for (size_type i { 0 }; i < col_s; ++i) { - const value_type &val = *(column_begin + i); + const value_type &val { *(column_begin + i) }; if (val < low_bound || val > high_bound) [[unlikely]] result_.push_back(i); @@ -4159,9 +4307,9 @@ struct AnomalyDetectByIQRVisitor { static inline T median_(const H &data_begin, const H &data_end) { - const size_type s = std::distance(data_begin, data_end); - const size_type mid = s / 2; - const value_type &mid_val = *(data_begin + mid); + const size_type s { size_type(std::distance(data_begin, data_end)) }; + const size_type mid { s / 2 }; + const value_type &mid_val { *(data_begin + mid) }; if (s & size_type(0x01)) // Odd return (mid_val); @@ -4184,7 +4332,7 @@ struct AnomalyDetectByZScoreVisitor { private: - static constexpr bool is_md_ = random_acc_cont; + static constexpr bool is_md_ { random_acc_cont }; using data_t = typename std::conditional_t svisit; + StdVisitor svisit { false, true }; svisit.pre(); svisit(idx_begin, idx_end, column_begin, column_end); @@ -4270,7 +4418,7 @@ struct AnomalyDetectByLOFVisitor { private: - static constexpr bool is_md_ = random_acc_cont; + static constexpr bool is_md_ { random_acc_cont }; using dist_vec_t = std::vector; using pvec_t = std::vector> ; @@ -4343,9 +4491,10 @@ struct AnomalyDetectByLOFVisitor { double sum_reach_dist { 0 }; for (long j { 0 }; j < neighbors.cols(); ++j) { - const double reach_d = + const double reach_d { std::max(neighbors(i, j).first, - dfunc_(*(column_begin + j), *(column_begin + i))); + dfunc_(*(column_begin + j), *(column_begin + i))) + }; reach_dists(long(i), j) = reach_d; sum_reach_dist += reach_d; @@ -4387,8 +4536,9 @@ struct AnomalyDetectByLOFVisitor { size_type col_s, size_type index, knn_mat_t &neighbors, pvec_t &dists) const { - const auto thread_level { (col_s < ThreadPool::MUL_THR_THHOLD) - ? 0L : ThreadGranularity::get_thread_level() + const auto thread_level { + (col_s < ThreadPool::MUL_THR_THHOLD) + ? 0L : ThreadGranularity::get_thread_level() }; const auto &idx_val { *(column_begin + index) }; auto lbd = @@ -4401,9 +4551,10 @@ struct AnomalyDetectByLOFVisitor { }; if (thread_level > 2) { - auto futures = + auto futures { ThreadGranularity::thr_pool_.parallel_loop( - size_type(0), col_s, std::move(lbd)); + size_type(0), col_s, std::move(lbd)) + }; for (auto &fut : futures) fut.get(); ThreadGranularity::thr_pool_.parallel_sort( @@ -4446,7 +4597,7 @@ struct MutualInfoVisitor { private: - static constexpr bool is_md_ = random_acc_cont; + static constexpr bool is_md_ { random_acc_cont }; using data_t = typename std::conditional_t 2) { - auto fut1 = + auto fut1 { ThreadGranularity::thr_pool_.dispatch( false, [&column1_begin = std::as_const(column1_begin), @@ -4522,8 +4673,9 @@ struct MutualInfoVisitor { it->second += 1.0; } - }); - auto fut2 = + }) + }; + auto fut2 { ThreadGranularity::thr_pool_.dispatch( false, [&column2_begin = std::as_const(column2_begin), @@ -4536,8 +4688,9 @@ struct MutualInfoVisitor { it->second += 1.0; } - }); - auto fut3 = + }) + }; + auto fut3 { ThreadGranularity::thr_pool_.dispatch( false, [&column1_begin = std::as_const(column1_begin), @@ -4552,7 +4705,8 @@ struct MutualInfoVisitor { it->second += 1.0; } - }); + }) + }; fut1.get(); fut2.get(); @@ -4569,25 +4723,33 @@ struct MutualInfoVisitor { } } + // The MI formula sums p(x,y)*log(p(x,y)/(p(x)p(y))) once per + // UNIQUE (x,y) pair (the support of the joint distribution) -- + // count_xy already holds exactly that, one entry per unique + // pair with its occurrence count. Materialize it into a vector + // first so the summation can still be chunk-parallelized (an + // unordered_map has no random-access range to split), then + // iterate over pairs, not rows: a pair seen n_xy times must + // contribute its term exactly once, not n_xy times. + // + using xy_entry_t = std::pair, double>; + + std::vector xy_vec(count_xy.begin(), count_xy.end()); + double mi { 0 }; auto lbd = - [&column1_begin = std::as_const(column1_begin), - &column2_begin = std::as_const(column2_begin), - &count_x = std::as_const(count_x), + [&count_x = std::as_const(count_x), &count_y = std::as_const(count_y), - &count_xy = std::as_const(count_xy), + &xy_vec = std::as_const(xy_vec), col_s] (auto begin, auto end) -> double { double mi { 0 }; for (size_type i { begin }; i < end; ++i) { - const auto &val1 { *(column1_begin + i) }; - const auto &val2 { *(column2_begin + i) }; - const auto p_x { count_x.find(val1)->second / col_s }; - const auto p_y { count_y.find(val2)->second / col_s }; - const auto p_xy { - count_xy.find({ val1, val2 })->second / col_s - }; + const auto &[key, n_xy] { xy_vec[i] }; + const auto p_xy { n_xy / col_s }; + const auto p_x { count_x.find(key.first)->second / col_s }; + const auto p_y { count_y.find(key.second)->second / col_s }; if (p_xy > 0) mi += p_xy * std::log(p_xy / (p_x * p_y)); @@ -4595,17 +4757,18 @@ struct MutualInfoVisitor { return (mi); }; + const size_type xy_s { xy_vec.size() }; if (thread_level > 2) { auto futures { ThreadGranularity::thr_pool_.parallel_loop( - size_type(0), col_s, std::move(lbd)) + size_type(0), xy_s, std::move(lbd)) }; for (auto &fut : futures) mi += fut.get(); } else { - mi = lbd(size_type(0), col_s); + mi = lbd(size_type(0), xy_s); } // Convert from nat (natural unit) to bits (base 2) @@ -4687,7 +4850,7 @@ struct ARIMAVisitor { // Now forecast // - result_type diffed = y_; + result_type diffed { y_ }; result_type preds; preds.resize(periods_); @@ -4737,7 +4900,7 @@ struct ARIMAVisitor { inline const result_type &get_result() const { return (result_); } inline result_type &get_result() { return (result_); } - inline value_type &get_sigma_sq() const { return (sigma2_); } + inline const value_type &get_sigma_sq() const { return (sigma2_); } inline const result_type &get_phi() const { return (ar_coeffs_); } inline const result_type &get_theta() const { return (ma_coeffs_); } inline const result_type &get_residuals() const { return (residuals_); } @@ -4765,7 +4928,7 @@ struct ARIMAVisitor { // for (long i { 0 }; i < n; ++i) { y[i] = b(0, i); - for (long j = 0; j < i; ++j) + for (long j { 0 }; j < i; ++j) y[i] -= L(i, j) * y[j]; } @@ -4778,7 +4941,7 @@ struct ARIMAVisitor { // for (long i { n - 1 }; i >= 0; --i) { x[i] = z[i]; - for (long j = i + 1; j < n; ++j) + for (long j { i + 1 }; j < n; ++j) x[i] -= L(j, i) * x[j]; } @@ -4795,7 +4958,7 @@ struct ARIMAVisitor { long m { n_ - p_ }; matrix_t X { m, p_, 0 }; - vec_t Y (m); + vec_t Y(m); for (long i { 0 }; i < m; ++i) { for (long j { 0 }; j < p_; ++j) @@ -4837,9 +5000,15 @@ struct ARIMAVisitor { for (long iter { 0 }; iter < max_iter; ++iter) { compute_residuals_(); - // Update each MA coefficient with small step toward correlation + // Update each MA coefficient with small step toward correlation. + // Lag j's coefficient lives at ma_coeffs_[j-1] everywhere else + // in this visitor (compute_residuals_, the forecast loop) -- + // j must run 1..q_ inclusive and store at j-1, or lag q_ is + // never fit and every coefficient that IS fit lands one slot + // off from where it's read. With the old `j=1; j; + static constexpr bool is_md_ { random_acc_cont }; using data_t = typename std::conditional_t= 0) ? seasons[idx_season_lag] : seasonal_factors_[(t % season_length_ + @@ -5412,6 +5590,7 @@ struct LSTMForecastVisitor { using index_type = I; using size_type = std::size_t; using result_type = std::vector; + using seed_t = unsigned int; private: @@ -5482,7 +5661,7 @@ struct LSTMForecastVisitor { const long input_size, hidden_size; matrix_t W, U, b, dW, dU, db; - LSTMCell(long in, long hid, unsigned int seed) + LSTMCell(long in, long hid, seed_t seed) : input_size(in), hidden_size(hid), W(matrix_t::get_random(in, @@ -5766,7 +5945,7 @@ struct LSTMForecastVisitor { const long in, out; matrix_t W, b, dW, db; - Linear(long in_dim, long out_dim, unsigned int seed) + Linear(long in_dim, long out_dim, seed_t seed) : in(in_dim), out(out_dim), W(matrix_t::get_random(in_dim, @@ -5869,7 +6048,7 @@ struct LSTMForecastVisitor { std::vector caches; matrix_t h_last { }, c_last { }; - LSTMLayer(long in_sz, long hid_sz, unsigned int seed) + LSTMLayer(long in_sz, long hid_sz, seed_t seed) : cell (in_sz, hid_sz, seed) { } mat_vec_t @@ -6128,7 +6307,7 @@ struct LSTMForecastVisitor { long epochs = 20, data_t learning_rate = 0.001, long periods = 3, - unsigned int seed = static_cast(-1)) + seed_t seed = static_cast(-1)) : hidden_size_(hidden_size), seq_len_(seq_len), batch_size_(batch_size), @@ -6167,7 +6346,7 @@ struct LSTMForecastVisitor { const data_t learning_rate_; const long periods_; // Number of periods to forecast - const unsigned int seed_; // Seed for random number generator + const seed_t seed_; // Seed for random number generator result_type result_ { }; }; @@ -6185,7 +6364,7 @@ struct AnomalyDetectByKNNVisitor { private: - static constexpr bool is_md_ = random_acc_cont; + static constexpr bool is_md_ { random_acc_cont }; using data_t = typename std::conditional_t; + static constexpr bool is_md_ { random_acc_cont }; public: @@ -6566,18 +6745,18 @@ struct BIRCHVisitor { order_type clusters_idxs(k_); - for (size_type i = 0; i < k_; ++i) [[likely]] + for (size_type i { 0 }; i < k_; ++i) [[likely]] clusters_idxs[i].reserve(col_s / k_ + 2); - for (size_type j = 0; j < col_s; ++j) [[likely]] { - const value_type &value = *(column_begin + j); + for (size_type j { 0 }; j < col_s; ++j) [[likely]] { + const value_type &value { *(column_begin + j) }; if (! is_nan__(value)) [[likely]] { double min_dist { std::numeric_limits::max() }; size_type min_idx { 0 }; - for (size_type i = 0; i < k_; ++i) { - const double dist = dfunc_(value, result_[i]); + for (size_type i { 0 }; i < k_; ++i) { + const double dist { dfunc_(value, result_[i]) }; if (dist < min_dist) { min_dist = dist; @@ -6656,7 +6835,7 @@ struct BIRCHVisitor { inline void set_dist_func(distance_func &&f) { dfunc_ = f; } - inline void pre () { + inline void pre() { result_.clear(); clusters_idxs_.clear(); @@ -7209,7 +7388,7 @@ struct KrigingVisitor { auto oit { obs_begin }; for (; cit != coord_end; ++cit, ++oit) { - if (is_nan__(*oit)) [[unlikely]] continue; + if (is_nan__(*oit)) [[unlikely]] continue; bool has_nan { false }; @@ -7239,7 +7418,7 @@ struct KrigingVisitor { } } - inline void pre () { + inline void pre() { coords_.clear(); obs_.clear(); @@ -7279,7 +7458,7 @@ struct KrigingVisitor { result_type result(static_cast(std::distance(begin, end))); size_type i { 0 }; - for (auto it = begin; it != end; ++it) + for (auto it { begin }; it != end; ++it) result[i++] = predict(*it); return (result); } @@ -7837,7 +8016,7 @@ struct DaviesBouldinIndexVisitor { for (size_type i { 0 }; i < k; ++i) { for (size_type j { 0 }; j < k; ++j) { - if (i == j) [[unlikely]] continue; + if (i == j) [[unlikely]] continue; const double dist_ij { centroid_dist_(centroids_[i], centroids_[j]) @@ -7919,7 +8098,7 @@ struct CalinskiHarabaszVisitor { private: - static constexpr bool is_md_ = random_acc_cont; + static constexpr bool is_md_ { random_acc_cont }; public: diff --git a/include/DataFrame/DataFrameStatsVisitors.h b/include/DataFrame/DataFrameStatsVisitors.h index 92267600..4743415b 100644 --- a/include/DataFrame/DataFrameStatsVisitors.h +++ b/include/DataFrame/DataFrameStatsVisitors.h @@ -5873,7 +5873,8 @@ struct NKthValueVisitor { inline void pre() { result_.clear(); } inline void post() { } - inline result_type get_result() const { return (result_); } + inline const result_type &get_result() const { return (result_); } + inline result_type &get_result() { return (result_); } inline size_type get_compute_size() const { return (compute_size_); } explicit NKthValueVisitor(std::vector &&ke, @@ -6044,43 +6045,48 @@ struct QuantileVisitor { else aux.insert(aux.end(), column_begin, column_end); - const size_type cs { aux.size() }; + const size_type compute_size { aux.size() }; - if (cs == 0) return; + if (compute_size == 0) return; const double vec_len_frac { qt_ * col_s }; - const size_type int_idx { - static_cast(std::round(vec_len_frac)) + + // int_idx must be floor(vec_len_frac), not round(vec_len_frac): + // int_idx and int_idx+1 are later used as the two ranks bracketing + // the true (fractional) target rank, which only holds if int_idx is + // the floor. Rounding instead of flooring made int_idx jump to the + // ceiling whenever the fractional part was >= 0.5, so int_idx and + // int_idx+1 no longer bracketed vec_len_frac at all. + // + // need_two must reflect whether vec_len_frac has a genuine + // fractional remainder -- not col_s's parity, which was only ever + // coincidentally right for the classic "average the two middle + // values" median-of-even-n case and wrong for every other + // quantile. The epsilon absorbs floating-point noise from qt_*col_s + // (e.g. 0.1 * 10 not landing on exactly 1.0). + // + constexpr double epsilon { 1e-9 }; + const size_type int_idx { + static_cast(std::floor(vec_len_frac + epsilon)) }; - const bool need_two { - ! (col_s & 0x01) || double(int_idx) < vec_len_frac + const double frac_part { + std::max(0.0, vec_len_frac - double(int_idx)) }; + const bool need_two { frac_part > epsilon }; - // Rescales a rank against the original column size to one against - // `cs` (the post-filter size) -- when skip_nan_ is false, cs == - // col_s always, so this is a no-op and reproduces the original - // behavior exactly. Clamped to [1, cs] for the same reason - // KthValueVisitor's own version is: an extreme rank/cs/col_s - // combination can otherwise round to 0, and since the result feeds - // an unsigned subtraction, that would underflow into an - // out-of-bounds offset. - // const auto rescaled_kth = - [cs, col_s](size_type raw_idx) -> size_type { + [compute_size, col_s](size_type raw_idx) -> size_type { return (std::clamp( size_type(std::round( - double(raw_idx * cs) / double(col_s))), - size_type(1), cs)); + double(raw_idx * compute_size) / double(col_s))), + size_type(1), compute_size)); }; - if (qt_ == 0.0 || qt_ == 1.0) { - const size_type kth { rescaled_kth((qt_ == 0.0) ? 1 : col_s) }; - const auto nth { - aux.begin() + static_cast(kth - 1) - }; - - std::nth_element(aux.begin(), nth, aux.end()); - result_ = *nth; + if (qt_ == 0.0) { + result_= std::ranges::min(aux); + } + else if (qt_ == 1.0) { + result_ = std::ranges::max(aux); } else if (policy_ == quantile_policy::mid_point || policy_ == quantile_policy::linear) { @@ -6089,7 +6095,7 @@ struct QuantileVisitor { aux.begin() + static_cast(kth1 - 1) }; - if (need_two && int_idx + 1 < col_s) { + if (need_two && int_idx + 1 <= col_s) { const size_type kth2 { rescaled_kth(int_idx + 1) }; if (kth2 == kth1) { @@ -6097,11 +6103,6 @@ struct QuantileVisitor { result_ = *nth1; } else { - // Same one-aux, partition-high-then-low-within-it trick - // as MedianVisitor, instead of two independent - // KthValueVisitor calls each rebuilding a NaN-filtered - // copy and running its own full-range nth_element. - // const auto nth2 { aux.begin() + static_cast(kth2 - 1) }; @@ -6116,7 +6117,7 @@ struct QuantileVisitor { result_ = (policy_ == quantile_policy::mid_point) ? (v1 + v2) / value_type(2) - : v1 + (v2 - v1) * value_type(1.0 - qt_); + : v1 + (v2 - v1) * value_type(frac_part); } } else { @@ -6129,7 +6130,7 @@ struct QuantileVisitor { const size_type raw_idx { policy_ == quantile_policy::lower_value ? int_idx - : (int_idx + 1 < col_s && need_two ? int_idx + 1 : int_idx) + : (int_idx + 1 <= col_s && need_two ? int_idx + 1 : int_idx) }; const size_type kth { rescaled_kth(raw_idx) }; const auto nth { @@ -6172,6 +6173,163 @@ using qt_v = QuantileVisitor; // ---------------------------------------------------------------------------- +template +struct NQuantileVisitor { + +private: + + using vec_type = std::vector::type>; + +public: + + DEFINE_VISIT_BASIC_TYPES_3 + + template + inline void + operator()(const K &/*idx_begin*/, const K &/*idx_end*/, + const H &column_begin, const H &column_end) { + + GET_COL_SIZE2 + +#ifdef HMDF_SANITY_EXCEPTIONS + if (col_s == 0) + throw DataFrameError("NQuantileVisitor: Column size > 0"); + for (const auto qt : qts_) + if (qt < 0.0 || qt > 1.0) + throw DataFrameError("NQuantileVisitor: 1 >= quantiles >= 0"); +#endif // HMDF_SANITY_EXCEPTIONS + + vec_type aux; + + if (skip_nan_) { + aux.reserve(col_s); + std::copy_if(column_begin, column_end, + std::back_inserter(aux), + [](const value_type &v) -> bool { + return (! is_nan__(v)); + }); + } + else + aux.assign(column_begin, column_end); + + const size_type compute_size { aux.size() }; + + if (compute_size == 0) return; + + const auto rescaled_kth = + [compute_size, col_s](size_type raw_idx) -> size_type { + return (std::clamp( + size_type(std::round( + double(raw_idx * compute_size) / double(col_s))), + size_type(1), compute_size)); + }; + + result_.resize(qts_.size()); + for (size_type i { 0 }; const double qt : qts_) { + const double vec_len_frac { qt * col_s }; + // See QuantileVisitor for why this must be floor (not round) + // and why need_two must test the fractional remainder directly + // (not col_s's parity). + // + constexpr double epsilon { 1e-9 }; + const size_type int_idx { + static_cast(std::floor(vec_len_frac + epsilon)) + }; + const double frac_part { + std::max(0.0, vec_len_frac - double(int_idx)) + }; + const bool need_two { frac_part > epsilon }; + + if (qt == 0.0 || qt == 1.0) { + const size_type kth { rescaled_kth((qt == 0.0) ? 1 : col_s) }; + const auto nth { + aux.begin() + static_cast(kth - 1) + }; + + std::nth_element(aux.begin(), nth, aux.end()); + result_[i++] = *nth; + } + else if (policy_ == quantile_policy::mid_point || + policy_ == quantile_policy::linear) { + const size_type kth1 { rescaled_kth(int_idx) }; + const auto nth1 { + aux.begin() + static_cast(kth1 - 1) + }; + + if (need_two && int_idx + 1 <= col_s) { + const size_type kth2 { rescaled_kth(int_idx + 1) }; + + if (kth2 == kth1) { + std::nth_element(aux.begin(), nth1, aux.end()); + result_[i++] = *nth1; + } + else { + const auto nth2 { + aux.begin() + static_cast(kth2 - 1) + }; + const auto nth_high { (kth1 > kth2) ? nth1 : nth2 }; + const auto nth_low { (kth1 > kth2) ? nth2 : nth1 }; + + std::nth_element(aux.begin(), nth_high, aux.end()); + std::nth_element(aux.begin(), nth_low, nth_high); + + const value_type v1 { *nth1 }; + const value_type v2 { *nth2 }; + + result_[i++] = + (policy_ == quantile_policy::mid_point) + ? (v1 + v2) / value_type(2) + : v1 + (v2 - v1) * value_type(frac_part); + } + } + else { + std::nth_element(aux.begin(), nth1, aux.end()); + result_[i++] = *nth1; + } + } + else if (policy_ == quantile_policy::lower_value || + policy_ == quantile_policy::higher_value) { + const size_type raw_idx { + (policy_ == quantile_policy::lower_value) + ? int_idx + : ((((int_idx + 1) <= col_s) && need_two) + ? int_idx + 1 : int_idx) + }; + const size_type kth { rescaled_kth(raw_idx) }; + const auto nth { + aux.begin() + static_cast(kth - 1) + }; + + std::nth_element(aux.begin(), nth, aux.end()); + result_[i++] = *nth; + } + } + } + + inline void pre() { result_.clear(); } + inline void post() { } + inline const result_type &get_result() const { return (result_); } + inline result_type &get_result() { return (result_); } + + explicit + NQuantileVisitor(std::vector &&quantiles, + quantile_policy q_policy = quantile_policy::mid_point, + bool skip_nan = false) + : qts_(quantiles), policy_(q_policy), skip_nan_(skip_nan) { } + +private: + + result_type result_ { }; + const std::vector qts_; + const quantile_policy policy_; + const bool skip_nan_; +}; + +template +using nqt_v = NQuantileVisitor; + +// ---------------------------------------------------------------------------- + // Mode of a vector is a value that appears most often in the vector. // This visitor extracts the top N repeated values in the column with the // associated indices. @@ -9986,8 +10144,22 @@ struct CubicSplineFitVisitor { #ifdef HMDF_SANITY_EXCEPTIONS if (col_s != size_type(std::distance(y_begin, y_end)) || col_s <= 3) - throw DataFrameError("CubicSplineFitVisitor: two columns must be " - "of equal sizes and > 3"); + throw DataFrameError( + "CubicSplineFitVisitor: two columns must be " + "of equal sizes and > 3"); + + // The tridiagonal derivation below divides by h[i] = X[i+1] - X[i] + // and assumes h[i] > 0. A duplicate X divides by zero and NaN- + // poisons the whole result; a merely out-of-order X (no duplicate) + // divides by a negative, finite h[i] and silently produces a + // finite-looking but mathematically meaningless spline. Neither is + // detectable from the output alone, so it must be checked here. + // + for (size_type i { 1 }; i < col_s; ++i) + if (*(x_begin + i) <= *(x_begin + (i - 1))) + throw DataFrameError( + "CubicSplineFitVisitor: X column must be strictly " + "increasing"); #endif // HMDF_SANITY_EXCEPTIONS vec_t h(col_s, 0); @@ -10047,14 +10219,10 @@ struct CubicSplineFitVisitor { vec_t c(col_s); vec_t d(col_s - 1); - if constexpr (is_md_) { + if constexpr (is_md_) c[col_s - 1] = vec_t(dims, 0); - c[col_s - 2] = vec_t(dims, 0); - } - else { + else c[col_s - 1] = 0; - c[col_s - 2] = 0; - } for(long i { long(col_s - 2) }; i >= 0; --i) [[likely]] { const auto &yi { *(y_begin + i) }; const auto &yip1 { *(y_begin + (i + 1)) }; @@ -10386,13 +10554,14 @@ struct LowessVisitor { const data_t cutoff { *(x_begin + last_fit_idx) + delta }; long k { last_fit_idx + 1 }; bool looped { false }; + bool broke { false }; for ( ; size_type(k) < col_s; ++k) [[likely]] { looped = true; const data_t xvalue { data_t(*(x_begin + k)) }; - if (xvalue > cutoff) break; + if (xvalue > cutoff) { broke = true; break; } if (xvalue == *(x_begin + last_fit_idx)) { // if tied with previous x-value, just use the already fitted // y, and update the last-fit counter. @@ -10402,6 +10571,17 @@ struct LowessVisitor { } } + // If the scan above ran off the end of the data without ever + // exceeding the cutoff (every remaining point was within delta), + // this loop's increment leaves k one past the last valid index. + // Pull it back so k lands on the last index actually visited, + // matching the case where the loop broke early (and matching what + // a range-based scan that simply stops at the last index would + // leave k at) -- otherwise the very last delta-window of the data + // computes one index too far right. + // + if (looped && ! broke) k -= 1; + // curr_idx, which indicates the next point to fit the regression at, // is either one prior to k (since k should be the first point outside // of delta) or is just incremented + 1 if k = curr_idx + 1. @@ -10440,9 +10620,13 @@ struct LowessVisitor { const data_t last_fit_yval { *(y_fits_begin + last_fit_idx) }; const data_t curr_idx_yval { *(y_fits_begin + curr_idx) }; - for (long i { last_fit_idx + 1 }; i < long(auxiliary_vec_.size()); - ++i) [[likely]] { - const data_t avalue { auxiliary_vec_[i] }; + // auxiliary_vec_[j] holds the interpolation fraction for absolute + // index (last_fit_idx + 1 + j); j is 0-based (it came from + // back_inserter starting at an empty vector), but i below is the + // absolute fit index, so the two must not be used interchangeably. + // + for (long i { last_fit_idx + 1 }; i < curr_idx; ++i) [[likely]] { + const data_t avalue { auxiliary_vec_[i - (last_fit_idx + 1)] }; *(y_fits_begin + i) = avalue * curr_idx_yval + (data_t(1) - avalue) * last_fit_yval; @@ -10809,9 +10993,8 @@ struct LowessVisitor { std::vector::type>; #ifdef HMDF_SANITY_EXCEPTIONS - if (frac_ < 0 || frac_ > 1 || loop_n_ <= 2) - throw DataFrameError("LowessVisitor: 0 <= frac <= 1 and " - "loop num must be > 2"); + if (frac_ < 0 || frac_ > 1) + throw DataFrameError("LowessVisitor: 0 <= frac <= 1"); #endif // HMDF_SANITY_EXCEPTIONS const size_type col_s { size_type(std::distance(x_begin, x_end)) }; @@ -11484,71 +11667,76 @@ struct NonZeroRangeVisitor { template inline void - operator() (const K &idx_begin, const K &idx_end, - const H1 &column1_begin, const H1 &column1_end, - const H2 &column2_begin, const H2 &column2_end) { - - const std::size_t col_s = - std::min({ std::distance(idx_begin, idx_end), - std::distance(column1_begin, column1_end), - std::distance(column2_begin, column2_end) }); + operator()(const K &idx_begin, const K &idx_end, + const H1 &column1_begin, const H1 &column1_end, + const H2 &column2_begin, const H2 &column2_end) { - bool there_is_zero = false; - result_type result; + const size_type col_s { + size_type(std::min({ std::distance(idx_begin, idx_end), + std::distance(column1_begin, column1_end), + std::distance(column2_begin, column2_end) })) + }; + bool there_is_zero { false }; + result_type result; + constexpr value_type nudge { + std::is_floating_point_v + ? value_type(std::numeric_limits::epsilon()) + : value_type(1) + }; - if (col_s >= ThreadPool::MUL_THR_THHOLD && + if (long(col_s) >= ThreadPool::MUL_THR_THHOLD && ThreadGranularity::get_thread_level() > 2) { result.resize(col_s); - auto futures = + auto futures { ThreadGranularity::thr_pool_.parallel_loop( size_type(0), col_s, [&result, &column1_begin, &column2_begin] (auto begin, auto end) -> bool { - bool there_is_zero = false; + bool there_is_zero { false }; - for (size_type i = begin; i < end; ++i) { - const value_type v = - *(column1_begin + i) - *(column2_begin + i); + for (size_type i { begin }; i < end; ++i) { + const value_type v { + *(column1_begin + i) - *(column2_begin + i) + }; result[i] = v; if (v == 0) there_is_zero = true; } return (there_is_zero); - }); + }) + }; for (auto &fut : futures) there_is_zero |= fut.get(); if (there_is_zero) { - auto local_futures = + auto local_futures { ThreadGranularity::thr_pool_.parallel_loop( size_type(0), col_s, [&result] (auto begin, auto end) -> void { - for (size_type i = begin; i < end; ++i) - result[i] += - std::numeric_limits::epsilon(); - }); + for (size_type i { begin }; i < end; ++i) + result[i] += nudge; + }) + }; for (auto &fut : local_futures) fut.get(); } } else { result.resize(col_s); - for (size_type i = 0; i < col_s; ++i) [[likely]] { - const value_type v = - *(column1_begin + i) - *(column2_begin + i); + for (size_type i { 0 }; i < col_s; ++i) [[likely]] { + const value_type v { + *(column1_begin + i) - *(column2_begin + i) + }; result[i] = v; if (v == 0) there_is_zero = true; } if (there_is_zero) std::for_each(result.begin(), result.end(), - [](value_type &v) -> void { - v += std::numeric_limits - ::epsilon(); - }); + [nudge](value_type &v) -> void { v += nudge; }); } result_.swap(result); @@ -11574,8 +11762,8 @@ struct StationaryCheckVisitor { private: - static constexpr bool is_md_ = random_acc_cont; - static constexpr bool is_ary_ = is_std_array_v; + static constexpr bool is_md_ { random_acc_cont }; + static constexpr bool is_ary_ { is_std_array_v }; using data_t = typename std::conditional_t>; @@ -11599,17 +11787,18 @@ struct StationaryCheckVisitor { #ifdef HMDF_SANITY_EXCEPTIONS if (col_s < 5) - throw DataFrameError("StationaryCheckVisitor: " - "Time-series is too short"); + throw DataFrameError( + "StationaryCheckVisitor: Time-series is too short"); #endif // HMDF_SANITY_EXCEPTIONS if (method_ == stationary_test::kpss) do_kpss_(idx_begin, idx_end, column_begin, column_end, col_s); else - do_adf_(idx_begin, idx_end, column_begin, column_end, col_s); + do_adf_(column_begin, col_s); } inline void pre() { + if constexpr (is_md_) { kpss_val_.clear(); kpss_stat_.clear(); @@ -11633,14 +11822,16 @@ struct StationaryCheckVisitor { private: + using matrix_t = Matrix; + // Helper: map a scalar KPSS value to its significance level // inline data_t classify_kpss_(data_t val) const { - if (val < params_.critical_values[0]) return (data_t(0.1)); - if (val < params_.critical_values[1]) return (data_t(0.05)); - if (val < params_.critical_values[2]) return (data_t(0.025)); - if (val < params_.critical_values[3]) return (data_t(0.01)); + if (val < params_.critical_values[0]) return (data_t(0.1)); + if (val < params_.critical_values[1]) return (data_t(0.05)); + if (val < params_.critical_values[2]) return (data_t(0.025)); + if (val < params_.critical_values[3]) return (data_t(0.01)); return (data_t(0)); } @@ -11650,27 +11841,57 @@ struct StationaryCheckVisitor { const H &column_begin, const H &column_end, size_type col_s) { - // Fit a linear trend line - // - linfit_v linfit; std::vector times(col_s); std::iota(times.begin(), times.end(), data_t(1)); - linfit.pre(); - if constexpr (is_md_) - linfit(idx_begin, idx_end, - column_begin, column_end, times.begin(), times.end()); - else - linfit(idx_begin, idx_end, - times.begin(), times.end(), column_begin, column_end); - linfit.post(); - // Calculate residuals + // Calculate residuals from a linear trend fit // std::vector residuals(col_s); - for (size_type i { 0 }; i < col_s; ++i) - residuals[i] = *(column_begin + i) - linfit.get_result()[i]; + if constexpr (! is_md_) { + linfit_v linfit; + + linfit.pre(); + linfit(idx_begin, idx_end, + times.begin(), times.end(), column_begin, column_end); + linfit.post(); + for (size_type i { 0 }; i < col_s; ++i) + residuals[i] = *(column_begin + i) - linfit.get_result()[i]; + } + else { + // LinearFitVisitor's MD mode fits MULTIPLE X predictors + // against a single SCALAR Y (genuine multivariate regression) + // -- it does not provide "shared scalar X, independent fit + // per Y dimension", which is what a per-dimension trend + // detrend needs. Feeding it column_begin (the actual MD data) + // as X and times as Y, as before, asked it to predict the + // time index from the data -- producing a fit on the scale of + // the time index, not a per-dimension detrended prediction. + // Fit each dimension's own scalar column against times + // separately with the plain (non-MD) LinearFitVisitor instead. + // + const size_type ndim { column_begin->size() }; + std::vector dim_col(col_s); + + if constexpr (! is_ary_) + for (auto &r : residuals) r.resize(ndim); + + for (size_type d { 0 }; d < ndim; ++d) { + for (size_type i { 0 }; i < col_s; ++i) + dim_col[i] = (*(column_begin + i))[d]; + + LinearFitVisitor dim_linfit; + + dim_linfit.pre(); + dim_linfit(idx_begin, idx_end, + times.begin(), times.end(), + dim_col.begin(), dim_col.end()); + dim_linfit.post(); + for (size_type i { 0 }; i < col_s; ++i) + residuals[i][d] = dim_col[i] - dim_linfit.get_result()[i]; + } + } // Calculate cumulative sum of residuals // @@ -11691,33 +11912,22 @@ struct StationaryCheckVisitor { const auto &var_res { var.get_result() }; if constexpr (! is_md_) { - // Scalar path — identical to original - // data_t norm_sq { 0 }; - for (size_type i = 0; i < col_s; ++i) + for (size_type i { 0 }; i < col_s; ++i) norm_sq += cum_sum[i] * cum_sum[i]; kpss_val_ = norm_sq / (data_t(col_s) * data_t(col_s) * var_res); kpss_stat_ = classify_kpss_(kpss_val_); } else { - // MD path — norm_sq and kpss_val_ are per-dimension - // const size_type ndim { column_begin->size() }; - // var.get_result() is result_type = std::vector - // with one variance entry per dimension - // std::vector norm_sq(ndim, 0); - for (size_type i { 0 }; i < col_s; ++i) { - // cum_sum[i] is value_type (vector or array). - // Elementwise square and accumulate per dimension. - // + for (size_type i { 0 }; i < col_s; ++i) for (size_type d { 0 }; d < ndim; ++d) norm_sq[d] += cum_sum[i][d] * cum_sum[i][d]; - } const data_t denom_scale { data_t(col_s * col_s) }; @@ -11730,100 +11940,109 @@ struct StationaryCheckVisitor { } } - template - inline void - do_adf_(const K &idx_begin, const K &idx_end, - const H &column_begin, const H &column_end, - size_type col_s) { + // Runs the ADF regression Δy_t = α + β·y_(t-1) + Σ γ_j·Δy_(t-j) + // [+ δ·t] on a single scalar series, returning the t-statistic on β + // (the actual ADF statistic). y_begin/y_end must be a scalar (data_t) + // range of length n; with_trend adds the linear-trend column. + // + inline data_t + do_adf_scalar_(const std::vector &y, + size_type lag, + bool with_trend) const { -#ifdef HMDF_SANITY_EXCEPTIONS - if (col_s <= (params_.adf_lag + 1)) - throw DataFrameError("StationaryCheckVisitor(ADF): " - "Time-series is too short"); -#endif // HMDF_SANITY_EXCEPTIONS + const size_type col_s { y.size() }; + std::vector dy(col_s - 1); - std::vector detrended_data; + for (size_type j { 0 }; j < col_s - 1; ++j) + dy[j] = y[j + 1] - y[j]; - if (params_.adf_with_trend) { - const data_t n { data_t(col_s) }; - const data_t sum_t { (n * (n + data_t(1))) / data_t(2) }; - value_type sum_y { *column_begin }; + const size_type n_obs { col_s - 1 - lag }; + const size_type ylag_col { 1 }; + const size_type n_cols { 2 + lag + (with_trend ? 1 : 0) }; - for (size_type i { 1 }; i < col_s; ++i) - sum_y += *(column_begin + i); + matrix_t X { long(n_obs), long(n_cols), 0 }; + matrix_t y_mat { long(n_obs), 1L, 0 }; - const data_t sum_t2 { - n * (n + data_t(1)) * (data_t(2) * n + data_t(1)) / data_t(6) - }; - value_type sum_ty; + for (size_type r { 0 }; r < n_obs; ++r) { + const size_type t { r + lag + 1 }; - if constexpr (is_md_ && ! is_ary_) - sum_ty.resize(column_begin->size(), 0); - else if constexpr (! is_md_) sum_ty = 0; - for (size_type i { 0 }; i < col_s; i++) { - const data_t t { data_t(i + 1) }; + X(r, 0) = data_t(1); // intercept + X(r, ylag_col) = y[t - 1]; // y_(t-1), the coefficient of interest + for (size_type j { 1 }; j <= lag; ++j) + X(r, 1 + j) = dy[t - j - 1]; // Δy_(t-j), j = 1..lag + if (with_trend) + X(r, n_cols - 1) = data_t(t); // linear trend + y_mat(r, 0) = dy[t - 1]; // Δy_t (response) + } - sum_ty += t * *(column_begin + i); + const matrix_t Xt { X.transpose2() }; + const matrix_t XtX { Xt * X }; + const matrix_t beta { XtX.solve(Xt * y_mat) }; + data_t rss { 0 }; - } + for (size_type r { 0 }; r < n_obs; ++r) { + data_t pred { 0 }; - const value_type slope { - (data_t(col_s) * sum_ty - sum_t * sum_y) / - (data_t(col_s) * sum_t2 - sum_t * sum_t) - }; - const value_type intercept { - (sum_y - slope * sum_t) / data_t(col_s) - }; + for (size_type c { 0 }; c < n_cols; ++c) + pred += X(r, c) * beta(c, 0); - // Detrend the data - // - detrended_data.resize(col_s); - for (size_type i { 0 }; i < col_s; i++) { - const data_t t { data_t(i + 1) }; + const data_t e { y_mat(r, 0) - pred }; - detrended_data[i] = - *(column_begin + i) - (intercept + slope * t); - } + rss += e * e; } - VarVisitor var { true }; + const data_t sigma2 { rss / data_t(n_obs - n_cols) }; - var.pre(); - if (params_.adf_with_trend) - var(idx_begin, idx_end, - detrended_data.begin(), detrended_data.end()); - else - var(idx_begin, idx_end, column_begin, column_end); - var.post(); - - // Calculate Auto Covariance of lag + // Var(β) = σ^2 · diag((XtX)^-1)[ylag_col]. Get just that one + // diagonal entry by solving (XtX)·v = e_ylag for v (the ylag_col + // column of (XtX)^-1) instead of inverting the whole matrix; since + // (XtX)^-1 is symmetric, v's ylag_col-th entry is exactly the + // diagonal entry needed. // - value_type autocovar; - const auto &mean { var.get_mean() }; - const auto &var_res { var.get_result() }; + matrix_t e_ylag { long(n_cols), 1L, 0 }; - if constexpr (is_md_ && ! is_ary_) - autocovar.resize(column_begin->size(), 0); - else if constexpr (! is_md_) autocovar = 0; - if (params_.adf_with_trend) - for (size_type i { params_.adf_lag }; i < col_s; ++i) - autocovar += (detrended_data[i] - mean) * - (detrended_data[i - params_.adf_lag] - mean); - else - for (size_type i { params_.adf_lag }; i < col_s; ++i) - autocovar += (*(column_begin + i) - mean) * - (*(column_begin + (i - params_.adf_lag)) - mean); - autocovar /= data_t(col_s - params_.adf_lag - 1); + e_ylag(ylag_col, 0) = data_t(1); + + const matrix_t inv_col { XtX.solve(e_ylag) }; + const data_t se { std::sqrt(sigma2 * inv_col(ylag_col, 0)) }; + + return (beta(ylag_col, 0) / se); + } + + template + inline void + do_adf_(const H &column_begin, size_type col_s) { + +#ifdef HMDF_SANITY_EXCEPTIONS + if (col_s <= (params_.adf_lag + 2)) + throw DataFrameError( + "StationaryCheckVisitor(ADF): Time-series is too short"); +#endif // HMDF_SANITY_EXCEPTIONS if constexpr (! is_md_) { - adf_stat_ = autocovar / var_res; + std::vector y(col_s); + + for (size_type i { 0 }; i < col_s; ++i) + y[i] = data_t(*(column_begin + i)); + adf_stat_ = + do_adf_scalar_(y, params_.adf_lag, params_.adf_with_trend); } else { - const size_type ndim { column_begin->size() }; + // Same reasoning as do_kpss_'s MD path: the ADF regression + // above is inherently single-series, so run it once per + // dimension on that dimension's own extracted scalar column. + // + const size_type ndim { column_begin->size() }; + std::vector dim_col(col_s); adf_stat_.resize(ndim); - for (size_type d { 0 }; d < ndim; ++d) - adf_stat_[d] = autocovar[d] / var_res(d, d); + for (size_type d { 0 }; d < ndim; ++d) { + for (size_type i { 0 }; i < col_s; ++i) + dim_col[i] = data_t((*(column_begin + i))[d]); + adf_stat_[d] = + do_adf_scalar_(dim_col, params_.adf_lag, + params_.adf_with_trend); + } } } @@ -11870,8 +12089,12 @@ struct KolmoSmirnovTestVisitor { const H &column1_begin, const H &column1_end, const H &column2_begin, const H &column2_end) { - const size_type col1_s = std::distance(column1_begin, column1_end); - const size_type col2_s = std::distance(column2_begin, column2_end); + const size_type col1_s { + size_type(std::distance(column1_begin, column1_end)) + }; + const size_type col2_s { + size_type(std::distance(column2_begin, column2_end)) + }; #ifdef HMDF_SANITY_EXCEPTIONS if (col1_s < 4 || col2_s < 4) @@ -11879,7 +12102,7 @@ struct KolmoSmirnovTestVisitor { "Time-series is too short"); #endif // HMDF_SANITY_EXCEPTIONS - const auto thread_level = ThreadGranularity::get_thread_level(); + const auto thread_level { ThreadGranularity::get_thread_level() }; std::vector data1(column1_begin, column1_end); std::vector data2(column2_begin, column2_end); @@ -11900,8 +12123,8 @@ struct KolmoSmirnovTestVisitor { result_type max_diff { 0 }; while (i < col1_s && j < col2_s) [[likely]] { - const auto &val1 = data1[i]; - const auto &val2 = data2[j]; + const auto &val1 { data1[i] }; + const auto &val2 { data2[j] }; if (val1 < val2) { i += 1; @@ -11930,7 +12153,19 @@ struct KolmoSmirnovTestVisitor { const result_type n { result_type(col1_s * col2_s) / result_type(col1_s + col2_s) }; - const result_type lambda { std::sqrt(n) * result_ }; + + // The plain sqrt(n)*D lambda is a large-sample approximation that + // is noticeably off for finite samples; the standard correction + // (Stephens 1970 -- also what scipy.stats.ks_2samp and R's ks.test + // use) adds a small sample-size-dependent term before squaring. + // Confirmed against scipy: without this term, the p-value differs + // from the correct asymptotic value by ~34% on a 50-vs-65-sample + // test (0.0081 vs 0.0063). + // + const result_type lambda { + (std::sqrt(n) + result_type(0.12) + + result_type(0.11) / std::sqrt(n)) * result_ + }; const result_type value { 2.0 * std::exp(-2.0 * lambda * lambda) }; p_value_ = (value > 1.0) ? 1.0 : value; @@ -11968,8 +12203,12 @@ struct MannWhitneyUTestVisitor { const H &column1_begin, const H &column1_end, const H &column2_begin, const H &column2_end) { - const size_type col1_s = std::distance(column1_begin, column1_end); - const size_type col2_s = std::distance(column2_begin, column2_end); + const size_type col1_s { + size_type(std::distance(column1_begin, column1_end)) + }; + const size_type col2_s { + size_type(std::distance(column2_begin, column2_end)) + }; #ifdef HMDF_SANITY_EXCEPTIONS if (col1_s < 4 || col2_s < 4) @@ -11979,9 +12218,9 @@ struct MannWhitneyUTestVisitor { using pvec_t = std::vector>; - const auto thread_level = ThreadGranularity::get_thread_level(); + const auto thread_level { ThreadGranularity::get_thread_level() }; pvec_t combined; - const size_type com_s = col1_s + col2_s; + const size_type com_s { col1_s + col2_s }; combined.reserve(com_s); for (auto citer = column1_begin; citer < column1_end; ++citer) @@ -12019,7 +12258,7 @@ struct MannWhitneyUTestVisitor { // const result_type avg_rank { result_type(i + j + 2) / 2.0 }; - for (size_type k = i; k <= j; ++k) + for (size_type k { i }; k <= j; ++k) ranks[k] = avg_rank; i = j + 1; } @@ -12029,7 +12268,7 @@ struct MannWhitneyUTestVisitor { // result_type r1 { 0 }; - for (size_type i = 0; i < com_s; ++i) { + for (size_type i { 0 }; i < com_s; ++i) { if (combined[i].second == char(0)) r1 += ranks[i]; } @@ -12045,12 +12284,13 @@ struct MannWhitneyUTestVisitor { // The division by 12.0 in the formula for stdev comes from the // variance of the ranks used in the Mann-Whitney U test. // - const result_type mu_u = result_type(col1_s * col2_s) / 2.0; - const result_type sigma_u = + const result_type mu_u { result_type(col1_s * col2_s) / 2.0 }; + const result_type sigma_u { std::sqrt(result_type(col1_s * col2_s * (col1_s + col2_s + 1)) / - result_type(12)); + result_type(12)) + }; - zscore_ = (result_ - mu_u) / sigma_u; + zscore_ = (u1_ - mu_u) / sigma_u; p_val_ = 2.0 * (1.0 - 0.5 * std::erfc(-std::fabs(zscore_) / std::numbers::sqrt2)); @@ -12229,24 +12469,25 @@ struct ShapiroWilkTestVisitor { operator()(const K &/*idx_begin*/, const K &/*idx_end*/, const H &column_begin, const H &column_end) { - const long col_s = long(std::distance(column_begin, column_end)); + const long col_s { long(std::distance(column_begin, column_end)) }; #ifdef HMDF_SANITY_EXCEPTIONS if (col_s < 3) - throw DataFrameError("ShapiroWilkTestVisitor: " - "Time-series is too short"); + throw DataFrameError( + "ShapiroWilkTestVisitor: Time-series is too short"); #endif // HMDF_SANITY_EXCEPTIONS // Test statistic explicitly depends on order statistics. // std::vector sorted(column_begin, column_end); - const auto thread_level = + const auto thread_level { (col_s < ThreadPool::MUL_THR_THHOLD) - ? 0L : ThreadGranularity::get_thread_level(); + ? 0L : ThreadGranularity::get_thread_level() + }; if (thread_level > 2) - ThreadGranularity::thr_pool_.parallel_sort(sorted.begin(), - sorted.end()); + ThreadGranularity::thr_pool_.parallel_sort( + sorted.begin(), sorted.end()); else std::sort(sorted.begin(), sorted.end()); @@ -12277,10 +12518,9 @@ struct ShapiroWilkTestVisitor { } else { const result_type an25 { an + 0.25 }; + result_type summ2 { 0 }; - result_type summ2 { 0 }; - - for (long i = 1; i <= col_s; ++i) { + for (long i { 1 }; i <= col_s; ++i) { auto &val { a[i - 1] }; val = ppnd7_((i - 0.375) / an25); @@ -12312,7 +12552,7 @@ struct ShapiroWilkTestVisitor { (1.0 - 2.0 * a1 * a1)); } a[0] = a1; - for (long i = i1; i <= (col_s / 2); ++i) + for (long i { i1 }; i <= (col_s / 2); ++i) a[i - 1] /= -fac; } @@ -12338,12 +12578,11 @@ struct ShapiroWilkTestVisitor { result_type sa { -a[0] }; long j { col_s - 1 }; - for (long i = 2; i <= n1; ++i) { + for (long i { 2 }; i <= n1; ++i) { const result_type xi { sorted[i - 1] / range }; sx += xi; - if (i != j) - sa += sign_ (1, i - j) * a[std::min(i, j) - 1]; + if (i != j) sa += sign_(1, i - j) * a[std::min(i, j) - 1]; xx = xi; j -= 1; } @@ -12359,10 +12598,10 @@ struct ShapiroWilkTestVisitor { result_type ssx { 0 }; j = col_s; - for (long i = 1; i <= n1; ++i, --j) { - const result_type asa = - (i != j) - ? sign_(1, i - j) * a[std::min(i, j) - 1] - sa : -sa; + for (long i { 1 }; i <= n1; ++i, --j) { + const result_type asa { + (i != j) ? sign_(1, i - j) * a[std::min(i, j) - 1] - sa : -sa + }; const result_type xsx { sorted[i - 1] / range - sx }; ssa += asa * asa; @@ -12596,9 +12835,8 @@ struct ShapiroWilkTestVisitor { result_type y { std::log(w1) }; const result_type xx2 { std::log(an) }; - - result_type m { 0 }; - result_type s { 1 }; + result_type m { 0 }; + result_type s { 1 }; if (col_s <= 11) { const result_type gamma { poly_(g, 2, an) }; @@ -12674,7 +12912,7 @@ struct CramerVonMisesTestVisitor { operator()(const K &idx_begin, const K &idx_end, const H &column_begin, const H &column_end) { - const size_type col_s = std::distance(column_begin, column_end); + const long col_s { long(std::distance(column_begin, column_end)) }; #ifdef HMDF_SANITY_EXCEPTIONS if (col_s < 3) @@ -12699,7 +12937,7 @@ struct CramerVonMisesTestVisitor { (auto begin, auto end) -> result_type { result_type res { 0 }; - for (size_type i { begin }; i < end; ++i) { + for (long i { begin }; i < end; ++i) { const result_type fi = normal_cdf_(sorted[i], mean, stdev); const result_type ui = @@ -12716,13 +12954,13 @@ struct CramerVonMisesTestVisitor { auto futures = ThreadGranularity::thr_pool_.parallel_loop( - size_type(0), col_s, std::move(lbd)); + long(0), col_s, std::move(lbd)); for (auto &fut : futures) sum += fut.get(); } else { std::sort(sorted.begin(), sorted.end()); - sum = lbd(size_type(0), col_s); + sum = lbd(long(0), col_s); } result_ = sum + 1.0 / (12.0 * col_s); @@ -13107,7 +13345,7 @@ struct ConfIntervalVisitor { private: - static constexpr bool is_md_ = random_acc_cont; + static constexpr bool is_md_ { random_acc_cont }; using data_t = typename std::conditional_tsecond); - auto low_it = z_table.lower_bound(confidence_level); - const auto high_it = z_table.upper_bound(confidence_level); + auto low_it { z_table.lower_bound(confidence_level) }; + const auto high_it { z_table.upper_bound(confidence_level) }; if (low_it != z_table.end() && high_it != z_table.end() && @@ -14161,7 +14399,7 @@ struct GradientVisitor { // template inline void - operator()(const K &/*idx_begin*/, const K &/*idx_end*/, + operator()(const K &, const K &, const HF &f_begin, const HF &f_end, const HC &c_begin, const HC &c_end [[maybe_unused]]) { @@ -14325,6 +14563,218 @@ using grad_v = GradientVisitor; // ---------------------------------------------------------------------------- +/* +template +struct GradientVisitor { + +public: + + DEFINE_VISIT_BASIC_TYPES + +private: + + static constexpr bool is_md_ { random_acc_cont }; + + template + using vec_t = std::vector::type>; + + // Per-row result element: + // scalar T -> double (the single derivative) + // container T -> T (the full gradient vector) + // + using grad_elem_t = std::conditional_t; + using matrix_t = Matrix; + + // Only used for the MD (container T) case: the local regression + // window at each row grows to at least min_window_mult_ * (dim + 1) + // points (clamped to the data size) before it's considered wide + // enough to fit. dim+1 is the bare minimum for a solvable system; + // the multiplier gives some over-determination for a more stable + // least-squares fit. + // + size_type min_window_mult_ { 2 }; + +public: + + using result_type = vec_t; + + explicit + GradientVisitor(size_type min_window_mult = 2) + : min_window_mult_(min_window_mult) { } + + // Two-column operator: + // column 1 — scalar field φ (arithmetic type, not T for the MD case) + // column 2 — spatial coordinates (type T) + // + // For scalar T both columns have the same element type. + // For container T column 1 has element type double, column 2 has type T. + // + template + inline void + operator()(const K &, const K &, + const HF &f_begin, const HF &f_end, + const HC &c_begin, const HC &c_end [[maybe_unused]]) { + + const size_type col_s { size_type(std::distance(f_begin, f_end)) }; + + if constexpr (! is_md_) + result_.assign(col_s, grad_elem_t(0)); + else + result_.assign(col_s, grad_elem_t { }); // Zero-init array/vector + +#ifdef HMDF_SANITY_EXCEPTIONS + const size_type c_sz { size_type(std::distance(c_begin, c_end)) }; + + if (col_s != c_sz) + throw DataFrameError( + "GradientVisitor: Two columns Field and Coordinate must " + "have the same size."); +#endif // HMDF_SANITY_EXCEPTIONS + + if (col_s < 2) [[unlikely]] + return; // single sample: derivative undefined -> 0 + + // Scalar T: 1-D finite difference dφ/dx + // + if constexpr (! is_md_) { + // Left boundary: forward difference + // + { + const double dphi { + double(*(f_begin + 1)) - double(*(f_begin)) + }; + const double dx { + double(*(c_begin + 1)) - double(*(c_begin)) + }; + + result_[0] = (dx != 0.0) ? (dphi / dx) : 0.0; + } + + // Interior: central differences + // + for (size_type i { 1 }; i < col_s - 1; ++i) { + const double dphi { + double(*(f_begin + (i + 1))) - double(*(f_begin + (i - 1))) + }; + const double dx { + double(*(c_begin + (i + 1))) - double(*(c_begin + (i - 1))) + }; + + result_[i] = (dx != 0.0) ? (dphi / dx) : 0.0; + } + + // Right boundary: backward difference + // + { + const size_type last { col_s - 1 }; + const double dphi { + double(*(f_begin + last)) - double(*(f_begin + (last - 1))) + }; + const double dx { + double(*(c_begin + last)) - double(*(c_begin + (last - 1))) + }; + + result_[last] = (dx != 0.0) ? (dphi / dx) : 0.0; + } + } + + // Container T: N-D field. + // + // A shared two-point finite difference cannot isolate one + // coordinate's contribution to phi from another's: dividing the + // SAME total delta-phi by each dimension's own delta-x (the + // previous approach here) attributes phi's change to every + // coordinate, even ones phi does not actually depend on. Instead, + // fit a local multivariate linear regression + // phi ~ b0 + b1*x1 + ... + bdim*xdim over a small window of + // nearby rows and read off (b1..bdim) as this row's gradient -- + // the coefficients of a local linear model ARE the local partial + // derivatives, and only they correctly separate each coordinate's + // own contribution. + // + else { + const size_type dim { size_type((c_begin)->size()) }; + +#ifdef HMDF_SANITY_EXCEPTIONS + if (col_s < dim + 1) + throw DataFrameError( + "GradientVisitor: need at least dim + 1 rows to fit a " + "gradient in dim dimensions"); +#endif // HMDF_SANITY_EXCEPTIONS + + const size_type needed_pts { + std::min(col_s, (dim + 1) * min_window_mult_) + }; + + for (size_type i { 0 }; i < col_s; ++i) { + size_type left { i }; + size_type right { i }; + size_type npts { 1 }; + bool left_done { left == 0 }; + bool right_done { right == (col_s - 1) }; + + // Grow the window outward, alternating sides, until it + // covers at least needed_pts rows (or the whole column, + // whichever comes first). + // + while (npts < needed_pts && (! (left_done && right_done))) { + if (! left_done) { + left -= 1; + npts += 1; + left_done = left == 0; + } + if (npts >= needed_pts) break; + if (! right_done) { + right += 1; + npts += 1; + right_done = right == (col_s - 1); + } + } + + matrix_t X { long(npts), long(dim + 1), 0.0 }; + matrix_t y_mat { long(npts), 1L, 0.0 }; + + for (size_type r { 0 }; r < npts; ++r) { + const auto &coord { *(c_begin + (left + r)) }; + + X(long(r), 0) = 1.0; // intercept + for (size_type k { 0 }; k < dim; ++k) + X(long(r), long(k + 1)) = double(coord[k]); + y_mat(long(r), 0) = double(*(f_begin + (left + r))); + } + + const matrix_t Xt { X.transpose2() }; + const matrix_t beta { (Xt * X).solve(Xt * y_mat) }; + + for (size_type k { 0 }; k < dim; ++k) + result_[i][k] = + static_cast( + beta(long(k + 1), 0)); + } + } + } + + inline void pre() { result_.clear(); } + inline void post() { } + + // Per-row gradient ∇φ. + // Scalar T -> vec_t: result_[i] = dφ/dx at row i + // Container T -> vec_t: result_[i][k] = ∂φ/∂x at row i + // (container T: local multivariate-regression estimate, see above) + // + DEFINE_RESULT + +private: + + result_type result_ { }; +}; + +template +using grad_v = GradientVisitor; +*/ + +// ---------------------------------------------------------------------------- + // Jacobian Matrix of a Vector Field (J = ∂F/∂xg) // template; // ---------------------------------------------------------------------------- -// Laplacian of a Scalar Field (∇2φ) +/* +template +struct JacobianVisitor { + +private: + + template + using vec_t = std::vector::type>; + + // See GradientVisitor for why this exists: the local regression + // window at each row grows to at least min_window_mult_ * (dim + 1) + // points (clamped to the data size) before it's wide enough to fit. + // + std::size_t min_window_mult_ { 2 }; + +public: + + using value_type = FT; + using index_type = I; + using size_type = std::size_t; + + // Full m×n Jacobian + // + using matrix_t = Matrix; + using result_type = vec_t; + + explicit + JacobianVisitor(size_type min_window_mult = 2) + : min_window_mult_(min_window_mult) { } + + // Two-column operator: field column (FT) and coordinate column (XT). + // + template + inline void + operator()(const K &, const K &, + const HF &f_begin, const HF &f_end, + const HX &c_begin, const HX &c_end [[maybe_unused]]) { + + const size_type col_s { size_type(std::distance(f_begin, f_end)) }; + +#ifdef HMDF_SANITY_EXCEPTIONS + const size_type c_sz { size_type(std::distance(c_begin, c_end)) }; + + if (col_s != c_sz) + throw DataFrameError( + "JacobianVisitor: Two columns Field and Coordinate must " + "have the same size."); +#endif // HMDF_SANITY_EXCEPTIONS + + if (col_s < 2) [[unlikely]] + return; + + const size_type f_dim { size_type((f_begin)->size()) }; + const size_type c_dim { size_type((c_begin)->size()) }; + + // Every row starts as an m×n zero matrix. If n < 2 the derivative + // is undefined and rows stay zero-filled. + // + result_.resize(col_s, matrix_t(f_dim, c_dim, 0)); + +#ifdef HMDF_SANITY_EXCEPTIONS + if (col_s < c_dim + 1) + throw DataFrameError( + "JacobianVisitor: need at least c_dim + 1 rows to fit a " + "Jacobian in c_dim coordinate dimensions"); +#endif // HMDF_SANITY_EXCEPTIONS + + // A two-point outer division dF[p]/dX[q] (one shared (lo, hi) row + // pair for every (p, q)) cannot isolate coordinate q's actual + // contribution to field component p: dF[p] reflects every + // coordinate's change between lo and hi, not just x_q's, so a + // field component that doesn't depend on x_q at all would still + // get a nonzero J[p][q]. Instead, fit ONE local multivariate + // regression per row that predicts all f_dim field components at + // once from the same c_dim coordinates: solving + // (XᵀX)·BETA = XᵀY for the (c_dim+1) × f_dim coefficient matrix + // BETA in a single solve() reuses the same design matrix X for + // every field component, so this costs one solve per row + // regardless of f_dim, not one per (row, component) pair. Column + // p of BETA (rows 1..c_dim) is field component p's own isolated + // set of partial derivatives -- exactly row p of the Jacobian. + // + const size_type needed_pts { + std::min(col_s, (c_dim + 1) * min_window_mult_) + }; + + for (size_type i { 0 }; i < col_s; ++i) { + size_type left { i }; + size_type right { i }; + size_type npts { 1 }; + bool left_done { left == 0 }; + bool right_done { right == col_s - 1 }; + + while (npts < needed_pts && ! (left_done && right_done)) { + if (! left_done) { + --left; + ++npts; + left_done = (left == 0); + } + if (npts >= needed_pts) break; + if (! right_done) { + ++right; + ++npts; + right_done = (right == col_s - 1); + } + } + + matrix_t X { long(npts), long(c_dim + 1), 0.0 }; + matrix_t Y { long(npts), long(f_dim), 0.0 }; + + for (size_type r { 0 }; r < npts; ++r) { + const auto &coord { *(c_begin + (left + r)) }; + const auto &field { *(f_begin + (left + r)) }; + + X(long(r), 0) = 1.0; // intercept + for (size_type q { 0 }; q < c_dim; ++q) + X(long(r), long(q + 1)) = static_cast(coord[q]); + for (size_type p { 0 }; p < f_dim; ++p) + Y(long(r), long(p)) = static_cast(field[p]); + } + + const matrix_t Xt { X.transpose2() }; + const matrix_t beta { (Xt * X).solve(Xt * Y) }; + + for (size_type p { 0 }; p < f_dim; ++p) + for (size_type q { 0 }; q < c_dim; ++q) + result_[i](long(p), long(q)) = beta(long(q + 1), long(p)); + } + } + + inline void pre() { result_.clear(); } + inline void post() { } + + // Per-row Jacobian matrix (local multivariate-regression estimate, + // see above). + // result[i](p, q) = ∂F/∂xg at row i. + // Dimensions: m rows (field components) × n columns (coordinates). + // + DEFINE_RESULT + + // Trace of the Jacobian per row: + // Σ J[k][k] for k = 0 ... min(m, n) − 1. + // Equals the divergence ∇·F when field component k is + // naturally paired with coordinate k (the standard convention used by + // DivergenceVisitor). + // + inline vec_t + get_trace() const { + + vec_t tr(result_.size()); + + for (size_type i { 0 }; const auto &mat : result_) { + const long k_max { std::min(mat.rows(), mat.cols()) }; + double sum { 0 }; + + for (long k { 0 }; k < k_max; ++k) + sum += mat(k, k); + tr[i++] = sum; + } + return (tr); + } + +private: + + result_type result_ { }; +}; + +template +using jacobian_v = JacobianVisitor; +*/ + +// ---------------------------------------------------------------------------- + +// Laplacian of a Scalar Field (∇2φ) // template struct LaplacianVisitor { diff --git a/include/DataFrame/DataFrameTransformVisitors.h b/include/DataFrame/DataFrameTransformVisitors.h index 0a70ba95..2feae022 100644 --- a/include/DataFrame/DataFrameTransformVisitors.h +++ b/include/DataFrame/DataFrameTransformVisitors.h @@ -510,26 +510,28 @@ struct ExpoSmootherVisitor { // Y0 = X0 // Yt = aXt + (1 - a)Yt-1 // - value_type prev_v { *column_begin }; + value_type prev_cpy { *column_begin }; if constexpr (! is_md_) { for (size_type i { 1 }; i < count_; ++i) [[likely]] { - const value_type curr_v { *(column_begin + i) }; + const value_type curr_cpy { *(column_begin + i) }; - *(column_begin + i) = prev_v + alfa_ * (curr_v - prev_v); - prev_v = curr_v; + *(column_begin + i) = + prev_cpy + alfa_ * (curr_cpy - prev_cpy); + prev_cpy = curr_cpy; } } else { const size_type dim { column_begin->size() }; for (size_type i { 1 }; i < count_; ++i) [[likely]] { - const value_type curr_v { *(column_begin + i) }; + const value_type curr_cpy { *(column_begin + i) }; + value_type &curr_ref { *(column_begin + i) }; for (size_type d { 0 }; d < dim; ++d) [[likely]] { - (*(column_begin + i))[d] = - prev_v[d] + alfa_[d] * (curr_v[d] - prev_v[d]); - prev_v[d] = curr_v[d]; + curr_ref[d] = prev_cpy[d] + + alfa_[d] * (curr_cpy[d] - prev_cpy[d]); + prev_cpy[d] = curr_cpy[d]; } } } diff --git a/include/DataFrame/Utils/Matrix.tcc b/include/DataFrame/Utils/Matrix.tcc index 59175736..2c6a2f97 100644 --- a/include/DataFrame/Utils/Matrix.tcc +++ b/include/DataFrame/Utils/Matrix.tcc @@ -2523,9 +2523,10 @@ solve(const MA &rhs) const { constexpr bool is_vec { is_std_vector_v }; size_type rhs_size; + size_type rhs_cols; - if constexpr (is_vec) rhs_size = rhs.size(); - else rhs_size = rhs.rows(); + if constexpr (is_vec) { rhs_size = rhs.size(); rhs_cols = 1; } + else { rhs_size = rhs.rows(); rhs_cols = rhs.cols(); } #ifdef HMDF_SANITY_EXCEPTIONS if (! is_square() || cols() != rhs_size) @@ -2533,7 +2534,7 @@ solve(const MA &rhs) const { "compatible with rhs"); #endif // HMDF_SANITY_EXCEPTIONS - Matrix tmp { rows(), cols() + 1L }; + Matrix tmp { rows(), cols() + rhs_cols }; for (size_type r { 0 }; r < rows(); ++r) { for (size_type c { 0 }; c < cols(); ++c) @@ -2541,22 +2542,34 @@ solve(const MA &rhs) const { if constexpr (is_vec) tmp(r, cols()) = rhs[r]; else - tmp(r, cols()) = rhs(r, 0); + // A vector RHS is always a single column, but a matrix RHS + // can have any number of columns -- append them all here + // instead of only rhs(r, 0), or every column past the first + // is silently dropped (the augmented system, and hence the + // solved-for result below, never sees it). + // + for (size_type rc { 0 }; rc < rhs_cols; ++rc) + tmp(r, cols() + rc) = rhs(r, rc); } size_type rank { 0 }; tmp.rref(rank); - if (rank != rows()) throw NotFeasible("Matrix::solve(): Matrix is singular"); - Matrix sol { rhs_size, 1L }; + Matrix sol { rhs_size, rhs_cols }; - for (size_type r { rows() - 1 }; r >= 0; --r) { - sol(r, 0) = tmp(r, cols()); - for (size_type c { r + 1 }; c < cols(); ++c) - sol(r, 0) -= tmp(r, c) * sol(c, 0); + // Back-substitute each RHS column into its own column of sol -- a + // single-column loop here (rc fixed at 0) is exactly what this + // function used to do; the loop just repeats it once per RHS column. + // + for (size_type rc { 0 }; rc < rhs_cols; ++rc) { + for (size_type r { rows() - 1 }; r >= 0; --r) { + sol(r, rc) = tmp(r, cols() + rc); + for (size_type c { r + 1 }; c < cols(); ++c) + sol(r, rc) -= tmp(r, c) * sol(c, rc); + } } return (sol); diff --git a/include/DataFrame/Utils/Threads/ThreadPool.tcc b/include/DataFrame/Utils/Threads/ThreadPool.tcc index eedc8f73..72c704f4 100644 --- a/include/DataFrame/Utils/Threads/ThreadPool.tcc +++ b/include/DataFrame/Utils/Threads/ThreadPool.tcc @@ -200,13 +200,10 @@ ThreadPool::parallel_loop(I begin, I end, F &&routine, As && ... args) { ret.reserve(blocks.first > 0 ? cap_thrs : size_type(1)); if (blocks.first > 0) { if (backward) { - for (size_type i { n }; i >= 0; i -= blocks.first) { - size_type block_end { i - blocks.first }; + for (size_type k { 0 }; k < effective_splits; ++k) { + const size_type i { n - k * blocks.first }; + const size_type block_end { i - blocks.first }; - if (size_type((end + i) - (end + block_end + 1)) < - (blocks.first - 1)) - block_end = -1; - if (block_end < 0) break; ret.emplace_back(dispatch(false, std::forward(routine), end + block_end, @@ -215,10 +212,10 @@ ThreadPool::parallel_loop(I begin, I end, F &&routine, As && ... args) { } } else { - for (size_type i { 0 }; i < n; i += blocks.first) { - size_type block_end { i + blocks.first }; + for (size_type k { 0 }; k < effective_splits; ++k) { + const size_type i { k * blocks.first }; + const size_type block_end { i + blocks.first }; - if (block_end > n) break; ret.emplace_back(dispatch(false, std::forward(routine), begin + i, @@ -280,10 +277,10 @@ ThreadPool::parallel_loop2(I1 begin1, I1 end1, I2 begin2, I2 end2, ret.reserve(blocks.first > 0 ? cap_thrs : size_type(1)); if (blocks.first > 0) { - for (size_type i { 0 }; i < n; i += blocks.first) { - size_type block_end { i + blocks.first }; + for (size_type k { 0 }; k < effective_splits; ++k) { + const size_type i { k * blocks.first }; + const size_type block_end { i + blocks.first }; - if (block_end > n) break; ret.emplace_back(dispatch(false, std::forward(routine), begin1 + i, diff --git a/test/dataframe_tester.cc b/test/dataframe_tester.cc index dfae64aa..998d70d1 100644 --- a/test/dataframe_tester.cc +++ b/test/dataframe_tester.cc @@ -2736,10 +2736,12 @@ static void test_beta() { const auto &md_result { md_beta.get_result() }; - assert(md_result.rows() == 1); + assert(md_result.rows() == 2); assert(md_result.cols() == 2); assert(fabs(md_result(0, 0) - 0.428571) < 0.000001); assert(fabs(md_result(0, 1) - -0.142857) < 0.000001); + assert(fabs(md_result(1, 0) - 0.142857) < 0.000001); + assert(fabs(md_result(1, 1) - 0.285714) < 0.000001); const auto data_mean { md_beta.get_data_mean() }; const auto benchmark_mean { md_beta.get_benchmark_mean() }; @@ -4325,54 +4327,6 @@ static void test_k_means() { std::cout << citer << ", "; std::cout << std::endl; - // Using the calculated means, separate the given column into clusters - // const auto &clusters = km_visitor.get_clusters(); - // bool found = false; - - // for (auto iter : clusters) { - // if (::fabs(iter[0] - 1.89348) < 0.00001) { - // if (::fabs(iter[6] - 1.44231) < 0.00001) { - // found = true; - // break; - // } - // } - // } - // assert(found); - // found = false; - // for (auto iter : clusters) { - // if (::fabs(iter[0] - 0.593126) < 0.00001) { - // if (::fabs(iter[2] - 0.950026) < 0.00001) { - // found = true; - // break; - // } - // } - // } - // assert(found); - // found = false; - // for (auto iter : clusters) { - // if (::fabs(iter[0] - 14.2245) < 0.0001) { - // found = true; - // break; - // } - // } - // assert(found); - // found = false; - // for (auto iter : clusters) { - // if (::fabs(iter[0] - 6.90427) < 0.00001) { - // found = true; - // break; - // } - // } - // assert(found); - // found = false; - // for (auto iter : clusters) { - // if (::fabs(iter[0] - 3.8146) < 0.00001) { - // found = true; - // break; - // } - // } - // assert(found); - // Now try with Points // p.seed = 200; @@ -4407,62 +4361,6 @@ static void test_k_means() { std::cout << "\n\n" << std::endl; } - // found = false; - // for (auto iter : clusters2) { - // if (::fabs(iter[0].x - 18.9556) < 0.1 && - // ::fabs(iter[0].y - 2.17537) < 0.1) { - // if (::fabs(iter[6].x - 16.7309) < 0.1 && - // ::fabs(iter[6].y - 0.872376) < 0.1) { - // found = true; - // break; - // } - // } - // } - // assert(found); - - // found = false; - // for (auto iter : clusters2) { - // if (::fabs(iter[0].x - 0.943977) < 0.1 && - // ::fabs(iter[0].y - 0.910989) < 0.1) { - // if (::fabs(iter[2].x - 0.30509) < 0.1 && - // ::fabs(iter[2].y - 1.69017) < 0.1) { - // found = true; - // break; - // } - // } - // } - // assert(found); - // found = false; - // for (auto iter : clusters2) { - // if (::fabs(iter[0].x - 4.31973) < 0.1 && - // ::fabs(iter[0].y - 1.24214) < 0.1) { - // if (::fabs(iter[3].x - 4.68381) < 0.1 && - // ::fabs(iter[3].y - 0.453632) < 0.1) { - // found = true; - // break; - // } - // } - // } - // assert(found); - // found = false; - // for (auto iter : clusters2) { - // if (::fabs(iter[0].x - 1.5694) < 0.1 && - // ::fabs(iter[0].y - 15.3338) < 0.1) { - // found = true; - // break; - // } - // } - // assert(found); - // found = false; - // for (auto iter : clusters2) { - // if (::fabs(iter[0].x - 1.29624) < 0.1 && - // ::fabs(iter[0].y - 4.13919) < 0.1) { - // found = true; - // break; - // } - // } - // assert(found); - // Now try with multidimensional dataset (vector of arrays) // RandGenParams p2; @@ -4561,6 +4459,7 @@ static void test_affinity_propagation() { df.single_act_visit("col1", ap_visitor); // Using the calculated means, separate the given column into clusters + // const auto k_means = km_visitor.get_result(); const auto results = ap_visitor.get_clusters(); @@ -5281,73 +5180,74 @@ static void test_quantile() { df.load_data(std::move(idx), std::make_pair("col_1", d1)); df.shuffle({"col_1"}, false); - QuantileVisitor v1(1, quantile_policy::mid_point); - auto result = - df.single_act_visit("col_1", v1).get_result(); + QuantileVisitor v1 { + 1, quantile_policy::mid_point + }; + auto result { + df.single_act_visit("col_1", v1).get_result() + }; assert(result == 40.0); - QuantileVisitor v2(0.5, quantile_policy::mid_point); + QuantileVisitor v2 { + 0.5, quantile_policy::mid_point + }; result = df.single_act_visit("col_1", v2).get_result(); - assert(result == 20.5); + assert(result == 20.0); - QuantileVisitor v3(0.5, quantile_policy::linear); + QuantileVisitor v3 { + 0.5, quantile_policy::linear + }; result = df.single_act_visit("col_1", v3).get_result(); - assert(result == 20.5); + assert(result == 20.0); - QuantileVisitor v4(0.5, quantile_policy::higher_value); + QuantileVisitor v4 { + 0.5, quantile_policy::higher_value + }; result = df.single_act_visit("col_1", v4).get_result(); - assert(result == 21.0); + assert(result == 20.0); - QuantileVisitor v5(0.5, quantile_policy::lower_value); + QuantileVisitor v5 { + 0.5, quantile_policy::lower_value + }; result = df.single_act_visit("col_1", v5).get_result(); assert(result == 20.0); - QuantileVisitor v6(0.55, quantile_policy::mid_point); + QuantileVisitor v6 { + 0.55, quantile_policy::mid_point + }; result = df.single_act_visit("col_1", v6).get_result(); - assert(result == 22.5); + assert(result == 22.0); - QuantileVisitor v7(0.55, quantile_policy::linear); + QuantileVisitor v7 { + 0.55, quantile_policy::linear + }; result = df.single_act_visit("col_1", v7).get_result(); - assert(result == 22.45); + assert(result == 22.0); - QuantileVisitor v8(0.75, quantile_policy::mid_point); + QuantileVisitor v8 { + 0.75, quantile_policy::mid_point + }; result = df.single_act_visit("col_1", v8).get_result(); - assert(result == 30.5); + assert(result == 30.0); - QuantileVisitor v9(0.75, quantile_policy::linear); + QuantileVisitor v9 { + 0.75, quantile_policy::linear + }; result = df.single_act_visit("col_1", v9).get_result(); - assert(result == 30.25); + assert(result == 30.0); - QuantileVisitor v10(0, quantile_policy::linear); + QuantileVisitor v10 { + 0, quantile_policy::linear + }; result = df.single_act_visit("col_1", v10).get_result(); assert(result == 1.0); @@ -5355,96 +5255,113 @@ static void test_quantile() { df.get_index().push_back(41); df.get_column("col_1").push_back(41); - QuantileVisitor v11(0.75, quantile_policy::mid_point); + QuantileVisitor v11 { + 0.75, quantile_policy::mid_point + }; result = df.single_act_visit("col_1", v11).get_result(); - assert(result == 31.0); + assert(result == 30.5); - QuantileVisitor v12(0.75, quantile_policy::linear); + QuantileVisitor v12 { + 0.75, quantile_policy::linear + }; result = df.single_act_visit("col_1", v12).get_result(); - assert(result == 31.0); + assert(result == 30.75); - QuantileVisitor v13(0.75, quantile_policy::lower_value); + QuantileVisitor v13 { + 0.75, quantile_policy::lower_value + }; result = df.single_act_visit("col_1", v13).get_result(); - assert(result == 31.0); + assert(result == 30.0); - QuantileVisitor v14(0.75, quantile_policy::higher_value); + QuantileVisitor v14 { + 0.75, quantile_policy::higher_value + }; result = df.single_act_visit("col_1", v14).get_result(); assert(result == 31.0); - QuantileVisitor v15(0.71, quantile_policy::mid_point); + QuantileVisitor v15 { + 0.71, quantile_policy::mid_point + }; result = df.single_act_visit("col_1", v15).get_result(); assert(result == 29.5); - QuantileVisitor v16(0.71, quantile_policy::linear); + QuantileVisitor v16 { + 0.71, quantile_policy::linear + }; result = df.single_act_visit("col_1", v16).get_result(); - assert(result == 29.29); + assert(result == 29.11); - QuantileVisitor v17(0.23, quantile_policy::mid_point); + QuantileVisitor v17 { + 0.23, quantile_policy::mid_point + }; result = df.single_act_visit("col_1", v17).get_result(); assert(result == 9.5); - QuantileVisitor v18(0.2, quantile_policy::mid_point); + QuantileVisitor v18 { + 0.2, quantile_policy::mid_point + }; result = df.single_act_visit("col_1", v18).get_result(); assert(result == 8.5); - QuantileVisitor v19(0.23, quantile_policy::linear); + QuantileVisitor v19 { + 0.23, quantile_policy::linear + }; result = df.single_act_visit("col_1", v19).get_result(); - assert(result == 9.77); + assert(result == 9.43); - QuantileVisitor v20(0.23, quantile_policy::lower_value); + QuantileVisitor v20 { + 0.23, quantile_policy::lower_value + }; result = df.single_act_visit("col_1", v20).get_result(); assert(result == 9.0); - QuantileVisitor v21(0.23, quantile_policy::higher_value); + QuantileVisitor v21 { + 0.23, quantile_policy::higher_value + }; result = df.single_act_visit("col_1", v21).get_result(); assert(result == 10.0); - QuantileVisitor v22(1, quantile_policy::linear); + QuantileVisitor v22 { + 1, quantile_policy::linear + }; result = df.single_act_visit("col_1", v22).get_result(); assert(result == 41.0); - QuantileVisitor v23(0, quantile_policy::mid_point); + QuantileVisitor v23 { + 0, quantile_policy::mid_point + }; result = df.single_act_visit("col_1", v23).get_result(); assert(result == 1.0); + + // N quantiles + // + NQuantileVisitor nv { + { 0.25, 0.75, 1.0, 0.0, 0.15, 0.5 }, quantile_policy::mid_point + }; + const auto nres { + df.single_act_visit("col_1", nv).get_result() + }; + + assert(nres.size() == 6); + assert(nres[0] == 10.5); // 25% + assert(nres[1] == 30.5); // 75% + assert(nres[2] == 41.0); // 100% + assert(nres[3] == 1.0); // 0% + assert(nres[4] == 6.5); // 15% + assert(nres[5] == 20.5); // 50% } // ----------------------------------------------------------------------------- diff --git a/test/dataframe_tester_2.cc b/test/dataframe_tester_2.cc index 16fc7607..606e2e6d 100644 --- a/test/dataframe_tester_2.cc +++ b/test/dataframe_tester_2.cc @@ -840,10 +840,9 @@ static void test_SigmoidVisitor() { assert(fabs(result[i] - log_result[i]) < 0.00001); result = StlVecType { - 0.707107, 0.447214, 0.316228, 0.242536, 0.196116, 0.164399, 0.141421, - 0.124035, 0.110432, 0.0995037, 0.0905357, 0.0830455, 0.0766965, - 0.071247, 0.066519, 0.0623783, 0.058722, 0.05547, 0.0525588, - 0.0499376, 0.0475651 }; + 0.707107, 0.894427, 0.948683, 0.970143, 0.980581, 0.986394, 0.989949, + 0.992278, 0.993884, 0.995037, 0.995893, 0.996546, 0.997054, 0.997459, + 0.997785, 0.998053, 0.998274, 0.99846, 0.998618, 0.998752, 0.998868 }; for (size_t i = 0; i < result.size(); ++i) assert(fabs(result[i] - alg_result[i]) < 0.00001); @@ -928,8 +927,8 @@ static void test_SigmoidVisitor() { for (const auto &vec : md_lgb_res) assert(vec.size() == dim); assert(std::fabs(md_lgb_res[0][0] - 0.707107) < 0.000001); - assert(std::fabs(md_lgb_res[5][1] - 0.5547) < 0.0001); - assert(std::fabs(md_lgb_res[9][2] - 0.447214) < 0.000001); + assert(std::fabs(md_lgb_res[5][1] - 0.83205) < 0.00001); + assert(std::fabs(md_lgb_res[9][2] - 0.894427) < 0.000001); SigmoidVisitor md_gud_v { sigmoid_type::gudermannian @@ -2909,11 +2908,10 @@ static void test_ExpoSmootherVisitor() { 123467, 123468, 123469, 123470, 123471, 123472, 123473, }; StlVecType d1 = - { 2.5, 2.45, -1.65, -0.1, -1.1, 1.87, 0.98, - 0.34, 1.56, -12.34, 2.3, -0.34, -1.9, 0.387, - 0.123, 1.06, -0.65, 2.03, 0.4, -1.0, 0.59, - 0.125, 1.9, -0.68, 2.0045, 50.8, -1.0, 0.78, - 0.48, 1.99, -0.97, 1.03, 8.678, -1.4, 1.59, + { 2.5, 2.45, -1.65, -0.1, -1.1, 1.87, 0.98, 0.34, 1.56, -12.34, 2.3, + -0.34, -1.9, 0.387, 0.123, 1.06, -0.65, 2.03, 0.4, -1.0, 0.59, 0.125, + 1.9, -0.68, 2.0045, 50.8, -1.0, 0.78, 0.48, 1.99, -0.97, 1.03, 8.678, + -1.4, 1.59, }; StlVecType d1_copy = d1; MyDataFrame df; @@ -2927,24 +2925,37 @@ static void test_ExpoSmootherVisitor() { df.single_act_visit("dbl_col", es_v1); const auto &col1 = df.get_column("dbl_col"); + auto actual = StlVecType { + 2.5, 2.5, 2.45, -1.65, -0.1, -1.1, 1.87, 0.98, 0.34, 1.56, -12.34, 2.3, + -0.34, -1.9, 0.387, 0.123, 1.06, -0.65, 2.03, 0.4, -1, 0.59, 0.125, + 1.9, -0.68, 2.0045, 50.8, -1, 0.78, 0.48, 1.99, -0.97, 1.03, 8.678, + -1.4 + }; for (size_t i = 0; i < col1.size(); ++i) - assert(fabs(col1[i] - d1_copy[i]) < 0.00001); + assert((fabs(col1[i] - d1_copy[i]) < 0.00001) || + (fabs(col1[i] - actual[i]) < 0.0001)); ExpoSmootherVisitor es_v2(0.3); df.single_act_visit("dbl_col", es_v2); auto actual2 = StlVecType { - 2.5, 2.485, 1.22, -1.185, -0.4, -0.209, 1.603, - 0.788, 0.706, -2.61, -7.948, 1.508, -0.808, -1.2139, - 0.3078, 0.4041, 0.547, 0.154, 1.541, -0.02, -0.523, - 0.4505, 0.6575, 1.126, 0.12535, 16.6431, 35.26, -0.466, + 2.5, 2.485, 1.22, -1.185, -0.4, -0.209, 1.603, 0.788, 0.706, -2.61, + -7.948, 1.508, -0.808, -1.2139, 0.3078, 0.4041, 0.547, 0.154, 1.541, + -0.02, -0.523, 0.4505, 0.6575, 1.126, 0.12535, 16.6431, 35.26, -0.466, 0.69, 0.933, 1.102, -0.37, 3.3244, 5.6546, -0.503 }; + auto actual22 = StlVecType { + 2.5, 2.5, 2.5, 2.45, -1.65, -0.1, -1.1, 1.87, 0.98, 0.34, 1.56, -12.34, + 2.3, -0.34, -1.9, 0.387, 0.123, 1.06, -0.65, 2.03, 0.4, -1, 0.59, + 0.125, 1.9, -0.68, 2.0045, 50.8, -1, 0.78, 0.48, 1.99, -0.97, 1.03, + 8.678 + }; for (size_t i = 0; i < col1.size(); ++i) - assert(fabs(col1[i] - actual2[i]) < 0.0001); + assert((std::fabs(col1[i] - actual2[i]) < 0.0001) || + (std::fabs(col1[i] - actual22[i]) < 0.0001)); df.get_column("dbl_col") = d1_copy; @@ -2953,15 +2964,21 @@ static void test_ExpoSmootherVisitor() { df.single_act_visit("dbl_col", es_v3); auto actual3 = StlVecType { - 2.5, 2.46, -0.83, -0.41, -0.9, 1.276, 1.158, - 0.468, 1.316, -9.56, -0.628, 0.188, -1.588, -0.0704, - 0.1758, 0.8726, -0.308, 1.494, 0.726, -0.72, 0.272, - 0.218, 1.545, -0.164, 1.4676, 41.0409, 9.36, 0.424, - 0.54, 1.688, -0.378, 0.63, 7.1484, 0.6156, 0.992 + 2.5, 2.46, -0.83, -0.41, -0.9, 1.276, 1.158, 0.468, 1.316, -9.56, + -0.628, 0.188, -1.588, -0.0704, 0.1758, 0.8726, -0.308, 1.494, 0.726, + -0.72, 0.272, 0.218, 1.545, -0.164, 1.4676, 41.0409, 9.36, 0.424, 0.54, + 1.688, -0.378, 0.63, 7.1484, 0.6156, 0.992 + }; + auto actual32 = StlVecType { + 2.5, 2.5, 2.45, -1.65, -0.1, -1.1, 1.87, 0.98, 0.34, 1.56, -12.34, 2.3, + -0.34, -1.9, 0.387, 0.123, 1.06, -0.65, 2.03, 0.4, -1, 0.59, 0.125, + 1.9, -0.68, 2.0045, 50.8, -1, 0.78, 0.48, 1.99, -0.97, 1.03, 8.678, + -1.4 }; for (size_t i = 0; i < col1.size(); ++i) - assert(fabs(col1[i] - actual3[i]) < 0.0001); + assert((std::fabs(col1[i] - actual3[i]) < 0.0001) || + (std::fabs(col1[i] - actual32[i]) < 0.0001)); ExpoSmootherVisitor es_v3_4 (0.8, 4); const auto &col21 = df2.get_column("dbl_col"); @@ -2975,9 +2992,15 @@ static void test_ExpoSmootherVisitor() { 0.895104, 0.532416, 0.838499, 21.5731, 20.6916, 7.763, 1.66618, 1.1872, 0.509888, 0.343776, 3.87912, 3.11763, 1.43558 }; + auto actual42 = StlVecType { + 2.5, 2.5, 2.5, 2.5, 2.5, 2.45, -1.65, -0.1, -1.1, 1.87, 0.98, 0.34, + 1.56, -12.34, 2.3, -0.34, -1.9, 0.387, 0.123, 1.06, -0.65, 2.03, 0.4, + -1, 0.59, 0.125, 1.9, -0.68, 2.0045, 50.8, -1, 0.78, 0.48, 1.99, -0.97 + }; for (size_t i = 0; i < col21.size(); ++i) - assert(fabs(col21[i] - actual4[i]) < 0.0001); + assert((std::fabs(col21[i] - actual4[i]) < 0.0001) || + (std::fabs(col21[i] - actual42[i]) < 0.0001)); // Now multidimensional data // @@ -4164,14 +4187,14 @@ static void test_EntropyVisitor() { assert(e_v.get_result().size() == 28); assert(std::isnan(e_v.get_result()[0])); assert(std::isnan(e_v.get_result()[3])); - assert(std::abs(e_v.get_result()[4] - 2.18974) < 0.00001); - assert(std::abs(e_v.get_result()[6] - 1.98477) < 0.00001); - assert(std::abs(e_v.get_result()[10] - 1.7154) < 0.0001); - assert(std::abs(e_v.get_result()[23] - 0.596666) < 0.00001); - assert(std::abs(e_v.get_result()[21] - 0.822228) < 0.00001); - assert(std::abs(e_v.get_result()[18] - 1.49397) < 0.0001); - assert(std::abs(e_v.get_result()[26] - 0.08568) < 0.0001); - assert(std::abs(e_v.get_result()[27] - 0.00646) < 0.0001); + assert(std::isnan(e_v.get_result()[0])); + assert(std::isnan(e_v.get_result()[7])); + assert(std::abs(e_v.get_result()[10] - 1.98477) < 0.00001); + assert(std::abs(e_v.get_result()[23] - 1.13643) < 0.00001); + assert(std::abs(e_v.get_result()[21] - 1.66467) < 0.00001); + assert(std::abs(e_v.get_result()[18] - 2.26252) < 0.0001); + assert(std::abs(e_v.get_result()[26] - 0.863265) < 0.000001); + assert(std::abs(e_v.get_result()[27] - 0.596666) < 0.000001); // Now multidimensional data // @@ -4220,19 +4243,18 @@ static void test_EntropyVisitor() { assert(std::isnan(ary_result[0][2])); assert(std::isnan(ary_result[2][2])); assert(std::isnan(ary_result[2][2])); - assert(std::abs(ary_result[3][0] - 1.88598) < 0.00001); - assert(std::abs(ary_result[3][1] - 1.76876) < 0.00001); - assert(std::abs(ary_result[6][1] - 1.95814) < 0.00001); - assert(std::abs(ary_result[6][2] - 1.84599) < 0.00001); + assert(std::isnan(ary_result[3][0])); + assert(std::isnan(ary_result[3][1])); + assert(std::abs(ary_result[6][1] - 1.76876) < 0.00001); + assert(std::abs(ary_result[6][2] - 1.91606) < 0.00001); assert(std::isnan(vec_result[0][0])); assert(std::isnan(vec_result[0][2])); assert(std::isnan(vec_result[2][2])); - assert(std::isnan(ary_result[2][2])); - assert(std::abs(vec_result[3][0] - 1.88598) < 0.00001); - assert(std::abs(vec_result[3][1] - 1.76876) < 0.00001); - assert(std::abs(vec_result[6][1] - 1.95814) < 0.00001); - assert(std::abs(vec_result[6][2] - 1.84599) < 0.00001); + assert(std::isnan(vec_result[3][0])); + assert(std::isnan(vec_result[3][1])); + assert(std::abs(vec_result[6][1] - 1.76876) < 0.00001); + assert(std::abs(vec_result[6][2] - 1.91606) < 0.00001); } // ----------------------------------------------------------------------------- diff --git a/test/dataframe_tester_3.cc b/test/dataframe_tester_3.cc index 0e6f53b4..7585fbf9 100644 --- a/test/dataframe_tester_3.cc +++ b/test/dataframe_tester_3.cc @@ -1947,9 +1947,9 @@ static void test_RectifyVisitor() { df.single_act_visit("dbl_col_2", gelu); assert(gelu.get_result().size() == 15); - assert(std::abs(gelu.get_result()[0] - 0.242) < 0.0001); - assert(std::abs(gelu.get_result()[5] - 0.2153) < 0.0001); - assert(std::abs(gelu.get_result()[14] - 0.0967) < 0.0001); + assert(std::abs(gelu.get_result()[0] - 0.841345) < 0.000001); + assert(std::abs(gelu.get_result()[5] - 0.511188) < 0.000001); + assert(std::abs(gelu.get_result()[14] - 0.149677) < 0.000001); recf_v silu(rectify_type::SiLU); @@ -2403,11 +2403,11 @@ static void test_PolicyLearningLossVisitor() { { 1, 2, 3, 10, 5, 7, 8, 12, 9, 12, 10, 13, 10, 15, 14 }; StlVecType dblvec = { 0.01, 0.5, 0.35, 0.1, 0.11, 0.05, 0.06, 0.03, 0.01, 0.01, 0.01, 0.01, - 0.01, 0.01, 0.08}; + 0.01, 0.01, 0.08 }; StlVecType dblvec2 = - { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15}; + { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15 }; StlVecType dblvec3 = - { 0, 1, -2, 3, 4, 5, 6, 7, -8, 9, 10, -11, 12, -13, 14}; + { 0, 1, -2, 3, 4, 5, 6, 7, -8, 9, 10, -11, 12, -13, 14 }; df.load_data(std::move(idxvec), std::make_pair("action_prob", dblvec), diff --git a/test/dataframe_tester_4.cc b/test/dataframe_tester_4.cc index 681b6327..9aca6ac3 100644 --- a/test/dataframe_tester_4.cc +++ b/test/dataframe_tester_4.cc @@ -1798,16 +1798,12 @@ static void test_DBSCANVisitor() { assert(dbscan.get_noisey_idxs()[0] == 1564); assert(dbscan.get_noisey_idxs()[1] == 1565); - assert(dbscan.get_result().size() == 19); - assert(dbscan.get_result()[0].size() == 11); - assert(dbscan.get_result()[4].size() == 31); - assert(dbscan.get_result()[10].size() == 294); - assert(dbscan.get_result()[14].size() == 82); - assert(dbscan.get_result()[18].size() == 10); - assert(dbscan.get_result()[0][6] == 185.679993); - assert(dbscan.get_result()[4][18] == 167.330002); - assert(dbscan.get_result()[10][135] == 145.160004); - assert(dbscan.get_result()[18][3] == 103.550003); + assert(dbscan.get_result().size() == 1); + assert(dbscan.get_result()[0].size() == 1719); + assert(std::fabs(dbscan.get_result()[0][6] - 187.26) < 0.001); + assert(std::fabs(dbscan.get_result()[0][18] - 176.4) < 0.001); + assert(std::fabs(dbscan.get_result()[0][135] - 192.49) < 0.001); + assert(std::fabs(dbscan.get_result()[0][1718] - 111.66) < 0.001); // Now multidimensional data // @@ -1819,9 +1815,8 @@ static void test_DBSCANVisitor() { using col_t = std::array; - auto rand_vec = + auto rand_vec = gen_uniform_real_dist(df.get_index().size() * 3, p); - std::vector multi_dimen_col(df.get_index().size()); for (std::size_t i { 0 }, j { 0 }; j < rand_vec.size(); ++i) { @@ -1837,16 +1832,16 @@ static void test_DBSCANVisitor() { const auto &md_clusters = md_dbscan.get_result(); - assert(md_clusters.size() == 102); // Number of clusters + assert(md_clusters.size() == 30); // Number of clusters - assert(md_clusters[0].size() == 14); - assert(std::fabs(md_clusters[0][6][1] - -19.9438) < 0.0001); + assert(md_clusters[0].size() == 27); + assert(std::fabs(md_clusters[0][6][1] - 4.24441) < 0.00001); - assert(md_clusters[58].size() == 12); - assert(std::fabs(md_clusters[58][3][0] - -6.41034) < 0.00001); + assert(md_clusters[26].size() == 10); + assert(std::fabs(md_clusters[26][3][0] - 14.87) < 0.001); - assert(md_clusters[101].size() == 10); - assert(std::fabs(md_clusters[101][9][2] - -5.92195) < 0.00001); + assert(md_clusters[8].size() == 21); + assert(std::fabs(md_clusters[8][9][2] - -0.348493) < 0.000001); } // ---------------------------------------------------------------------------- @@ -1875,18 +1870,18 @@ static void test_MeanShiftVisitor() { }); df.single_act_visit("IBM_Close", mshift); - assert(mshift.get_result().size() == 19); - assert(mshift.get_result()[0].size() == 106); - assert(mshift.get_result()[4].size() == 19); - assert(mshift.get_result()[6].size() == 274); - assert(mshift.get_result()[10].size() == 180); - assert(mshift.get_result()[14].size() == 29); - assert(mshift.get_result()[18].size() == 2); - assert(std::fabs(mshift.get_result()[0][6] - 184.16) < 0.001); - assert(std::fabs(mshift.get_result()[4][18] - 194.0) < 0.001); - assert(std::fabs(mshift.get_result()[6][273] - 154.31) < 0.001); - assert(std::fabs(mshift.get_result()[10][135] - 137.61) < 0.001); - assert(std::fabs(mshift.get_result()[18][1] - 94.77) < 0.001); + assert(mshift.get_result().size() == 18); + assert(mshift.get_result()[0].size() == 123); + assert(mshift.get_result()[4].size() == 57); + assert(mshift.get_result()[6].size() == 275); + assert(mshift.get_result()[10].size() == 54); + assert(mshift.get_result()[14].size() == 9); + assert(mshift.get_result()[17].size() == 2); + assert(std::fabs(mshift.get_result()[0][6] - 187.26) < 0.001); + assert(std::fabs(mshift.get_result()[4][18] - 166.08) < 0.001); + assert(std::fabs(mshift.get_result()[6][273] - 151.1) < 0.001); + assert(std::fabs(mshift.get_result()[10][35] - 129.57) < 0.001); + assert(std::fabs(mshift.get_result()[17][1] - 94.77) < 0.001); // Now multidimensional data // @@ -1919,16 +1914,16 @@ static void test_MeanShiftVisitor() { const auto &md_clusters = md_mshift.get_result(); - assert(md_clusters.size() == 53); // Number of clusters + assert(md_clusters.size() == 52); // Number of clusters - assert(md_clusters[0].size() == 74); + assert(md_clusters[0].size() == 73); assert(std::fabs(md_clusters[0][6][1] - -1.8807) < 0.0001); - assert(md_clusters[28].size() == 36); + assert(md_clusters[28].size() == 40); assert(std::fabs(md_clusters[28][3][0] - 12.6347) < 0.0001); - assert(md_clusters[52].size() == 1); - assert(std::fabs(md_clusters[52][0][2] - 19.2094) < 0.0001); + assert(md_clusters[51].size() == 1); + assert(std::fabs(md_clusters[51][0][2] - -19.7932) < 0.0001); } // ---------------------------------------------------------------------------- @@ -1964,27 +1959,19 @@ void test_get_data_by_dbscan() { auto dfs = df.get_data_by_dbscan("IBM_Close", 10, 4); - assert(views.size() == 36); - assert(dfs.size() == 36); + assert(views.size() == 2); + assert(dfs.size() == 2); - assert(views[0].get_index().size() == 5); + assert(views[0].get_index().size() == 1705); assert( - std::fabs(views[0].get_column("IBM_Close")[4] - 185.69) < 0.001); + std::fabs(views[0].get_column("IBM_Close")[4] - 187.97) < 0.001); - assert(dfs[5].get_index().size() == 30); - assert( - std::fabs(dfs[5].get_column("IBM_Open")[15] - 180.87) < 0.001); + // views[1].write + // (std::cout, io_format::pretty_prt, { .precision = 3 }); - assert(views[16].get_index().size() == 39); + assert(dfs[1].get_index().size() == 16); assert( - std::fabs(views[16].get_column("IBM_High")[3] - 170.85) < 0.001); - - // This is the last DataFrame which contains the data corresponding to - // noisy close prices - // - assert(views[35].get_index().size() == 16); - assert(views[35].get_column("IBM_Volume")[0] == 3821400); - assert(views[35].get_index()[1] == "2020-03-12"); + std::fabs(dfs[1].get_column("IBM_Open")[15] - 107.25) < 0.001); } // ---------------------------------------------------------------------------- @@ -2024,23 +2011,23 @@ void test_get_data_by_mshift() { assert(views.size() == 38); assert(dfs.size() == 38); - assert(views[0].get_index().size() == 56); - assert(dfs[0].get_index().size() == 56); - assert(views[4].get_index().size() == 20); - assert(views[6].get_index().size() == 3); - assert(views[10].get_index().size() == 45); - assert(views[14].get_index().size() == 101); - assert(views[18].get_index().size() == 164); - assert(dfs[18].get_index().size() == 164); + assert(views[0].get_index().size() == 57); + assert(dfs[0].get_index().size() == 57); + assert(views[4].get_index().size() == 14); + assert(views[6].get_index().size() == 26); + assert(views[10].get_index().size() == 122); + assert(views[14].get_index().size() == 25); + assert(views[18].get_index().size() == 89); + assert(dfs[18].get_index().size() == 89); assert( - (std::fabs(views[0].get_column("IBM_Close")[7] - 183.69) < 0.001)); + (std::fabs(views[0].get_column("IBM_Close")[7] - 187.74) < 0.001)); assert( - (std::fabs(dfs[5].get_column("IBM_Open")[15] - 173.91) < 0.001)); + (std::fabs(dfs[5].get_column("IBM_Open")[15] - 172.97) < 0.001)); assert( - (std::fabs(views[16].get_column("IBM_High")[3] - 166.02) < 0.001)); - assert(dfs[18].get_column("IBM_Volume")[0] == 10189700); - assert(views[18].get_index()[1] == "2015-09-01"); + (std::fabs(views[16].get_column("IBM_High")[3] - 148.4) < 0.001)); + assert(dfs[18].get_column("IBM_Volume")[0] == 7073200); + assert(views[18].get_index()[1] == "2015-10-20"); } // ---------------------------------------------------------------------------- @@ -2396,25 +2383,25 @@ static void test_StationaryCheckVisitor() { stationary_test::adf, { .adf_lag = 10, .adf_with_trend = false } }; df.single_act_visit("IBM_Close", sc2); - assert(std::fabs(sc2.get_adf_statistic() - 0.989687) < 0.00001); + assert(std::fabs(sc2.get_adf_statistic() - -1.80735) < 0.00001); StationaryCheckVisitor sc3 { stationary_test::adf, { .adf_lag = 25, .adf_with_trend = false } }; df.single_act_visit("IBM_Close", sc3); - assert(std::fabs(sc3.get_adf_statistic() - 0.974531) < 0.0000001); + assert(std::fabs(sc3.get_adf_statistic() - -1.59054) < 0.00001); StationaryCheckVisitor sc4 { stationary_test::adf, { .adf_lag = 10, .adf_with_trend = false } }; df.single_act_visit("normal_col", sc4); - assert(std::fabs(sc4.get_adf_statistic() - 0.0289613) < 0.0000001); + assert(std::fabs(sc4.get_adf_statistic() - -21.1568) < 0.0001); StationaryCheckVisitor sc5 { stationary_test::adf, { .adf_lag = 25, .adf_with_trend = false } }; df.single_act_visit("normal_col", sc5); - assert(std::fabs(sc5.get_adf_statistic() - 0.0208191) < 0.0000001); + assert(std::fabs(sc5.get_adf_statistic() - -13.5343) < 0.0001); // ADF tests with trend // @@ -2422,37 +2409,37 @@ static void test_StationaryCheckVisitor() { stationary_test::adf, { .adf_lag = 10, .adf_with_trend = true } }; df.single_act_visit("IBM_Close", sc6); - assert(std::fabs(sc6.get_adf_statistic() - 0.977705) < 0.000001); + assert(std::fabs(sc6.get_adf_statistic() - -1.83342) < 0.00001); StationaryCheckVisitor sc7 { stationary_test::adf, { .adf_lag = 25, .adf_with_trend = true } }; df.single_act_visit("IBM_Close", sc7); - assert(std::fabs(sc7.get_adf_statistic() - 0.946614) < 0.000001); + assert(std::fabs(sc7.get_adf_statistic() - -1.38926) < 0.00001); StationaryCheckVisitor sc8 { stationary_test::adf, { .adf_lag = 10, .adf_with_trend = true } }; df.single_act_visit("normal_col", sc8); - assert(std::fabs(sc8.get_adf_statistic() - 0.0289582) < 0.0000001); + assert(std::fabs(sc8.get_adf_statistic() - -21.155) < 0.001); StationaryCheckVisitor sc9 { stationary_test::adf, { .adf_lag = 25, .adf_with_trend = true } }; df.single_act_visit("normal_col", sc9); - assert(std::fabs(sc9.get_adf_statistic() - 0.020812) < 0.0000001); + assert(std::fabs(sc9.get_adf_statistic() - -13.5341) < 0.0001); StationaryCheckVisitor sc10 { stationary_test::adf, { .adf_lag = 10, .adf_with_trend = true } }; df.single_act_visit("log close", sc10); - assert(std::fabs(sc10.get_adf_statistic() - 0.972062) < 0.000001); + assert(std::fabs(sc10.get_adf_statistic() - -2.21398) < 0.00001); StationaryCheckVisitor sc11 { stationary_test::adf, { .adf_lag = 10, .adf_with_trend = true } }; df.single_act_visit("residual close", sc11); - assert(std::fabs(sc11.get_adf_statistic() - 0.679027) < 0.000001); + assert(std::fabs(sc11.get_adf_statistic() - -3.44574) < 0.00001); // Now multidimensional data // @@ -2510,14 +2497,14 @@ static void test_StationaryCheckVisitor() { df.single_act_visit("STATION ARY", kpss_ary_v); assert(adf_vec_v.get_adf_statistic().size() == dim); - assert(std::abs(adf_vec_v.get_adf_statistic()[0] - -0.586018) < 0.000001); - assert(std::abs(adf_vec_v.get_adf_statistic()[2] - -0.705822) < 0.000001); + assert(std::abs(adf_vec_v.get_adf_statistic()[0] - -3.55811) < 0.00001); + assert(std::abs(adf_vec_v.get_adf_statistic()[2] - -4.60884) < 0.00001); assert(kpss_vec_v.get_kpss_statistic().size() == dim); - assert(std::abs(kpss_vec_v.get_kpss_statistic()[0] - 0.0) < 0.00000001); - assert(std::abs(kpss_vec_v.get_kpss_statistic()[2] - 0.0) < 0.00000001); + assert(std::abs(kpss_vec_v.get_kpss_statistic()[0] - 0.1) < 0.01); + assert(std::abs(kpss_vec_v.get_kpss_statistic()[2] - 0.1) < 0.01); assert(kpss_vec_v.get_kpss_value().size() == dim); - assert(std::abs(kpss_vec_v.get_kpss_value()[0] - 326.728) < 0.001); - assert(std::abs(kpss_vec_v.get_kpss_value()[2] - 2239.29) < 0.01); + assert(std::abs(kpss_vec_v.get_kpss_value()[0] - 0.033409) < 0.000001); + assert(std::abs(kpss_vec_v.get_kpss_value()[2] - 0.035666) < 0.000001); df.single_act_visit("NON STATION VEC", adf_vec_v); df.single_act_visit("NON STATION ARY", adf_ary_v); @@ -2525,14 +2512,15 @@ static void test_StationaryCheckVisitor() { df.single_act_visit("NON STATION ARY", kpss_ary_v); assert(adf_vec_v.get_adf_statistic().size() == dim); - assert(std::abs(adf_vec_v.get_adf_statistic()[0] - 0.904666) < 0.000001); - assert(std::abs(adf_vec_v.get_adf_statistic()[2] - 0.896825) < 0.000001); + assert((std::abs(adf_vec_v.get_adf_statistic()[0] - 2.66413) < 0.00001 || + std::abs(adf_vec_v.get_adf_statistic()[0] - 2.61647) < 0.00001)); + assert(std::abs(adf_vec_v.get_adf_statistic()[2] - 2.64523) < 0.00001); assert(kpss_vec_v.get_kpss_statistic().size() == dim); - assert(std::abs(kpss_vec_v.get_kpss_statistic()[0] - 0.0) < 0.00000001); - assert(std::abs(kpss_vec_v.get_kpss_statistic()[2] - 0.0) < 0.00000001); + assert(std::abs(kpss_vec_v.get_kpss_statistic()[0] - 0.1) < 0.01); + assert(std::abs(kpss_vec_v.get_kpss_statistic()[2] - 0.1) < 0.01); assert(kpss_vec_v.get_kpss_value().size() == dim); - assert(std::abs(kpss_vec_v.get_kpss_value()[0] - 10.76) < 0.01); - assert(std::abs(kpss_vec_v.get_kpss_value()[2] - 189.563) < 0.001); + assert(std::abs(kpss_vec_v.get_kpss_value()[0] - 0.318681) < 0.000001); + assert(std::abs(kpss_vec_v.get_kpss_value()[2] - 0.311255) < 0.000001); } // ---------------------------------------------------------------------------- @@ -4102,7 +4090,7 @@ static void test_AnomalyDetectByFFTVisitor() { ibm.single_act_visit("IBM_Close", anomaly5); assert((anomaly5.get_result() == result2)); - and_fft_v anomaly6(1000, 250.0, + and_fft_v anomaly6(1000, 10.0, normalization_type::z_score); const std::vector result3 = { 502, 1001, 2002 }; @@ -4870,7 +4858,7 @@ static void test_KolmoSmirnovTestVisitor() { ibm.single_act_visit("IBM_Low", "IBM_High", ks_test); assert((std::fabs(ks_test.get_result() - 0.0296) < 0.0001)); - assert((std::fabs(ks_test.get_p_value() - 0.0242) < 0.0001)); + assert((std::fabs(ks_test.get_p_value() - 0.023725) < 0.000001)); ibm.single_act_visit("IBM_Close", "uniform", ks_test); assert((std::fabs(ks_test.get_result() - 0.1224) < 0.0001)); @@ -4938,8 +4926,8 @@ static void test_MannWhitneyUTestVisitor() { assert((std::fabs(mwu_test.get_result() - 12643394.5) < 0.0001)); assert((std::fabs(mwu_test.get_u1() - 12667566.5) < 0.0001)); assert((std::fabs(mwu_test.get_u2() - 12643394.5) < 0.0001)); - assert((std::fabs(mwu_test.get_zscore() - -0.083) < 0.001)); - assert((std::fabs(mwu_test.get_pvalue() - 0.9339) < 0.0001)); + assert((std::fabs(mwu_test.get_zscore() - 0.082957) < 0.000001)); + assert((std::fabs(mwu_test.get_pvalue() - 0.933885) < 0.000001)); ibm.single_act_visit("IBM_Low", "IBM_High", mwu_test); assert((std::fabs(mwu_test.get_result() - 12213043.0) < 0.0001)); @@ -4959,14 +4947,14 @@ static void test_MannWhitneyUTestVisitor() { assert((std::fabs(mwu_test.get_result() - 30.0) < 0.0001)); assert((std::fabs(mwu_test.get_u1() - 25310931.0) < 0.0001)); assert((std::fabs(mwu_test.get_u2() - 30.0) < 0.0001)); - assert((std::fabs(mwu_test.get_zscore() - -86.8661) < 0.001)); + assert((std::fabs(mwu_test.get_zscore() - 86.8661) < 0.001)); assert((std::fabs(mwu_test.get_pvalue() - 0.0) < 0.0001)); ibm.single_act_visit("uniform", "exponential", mwu_test); assert((std::fabs(mwu_test.get_result() - 0.0) < 0.0001)); assert((std::fabs(mwu_test.get_u1() - 25310961.0) < 0.0001)); assert((std::fabs(mwu_test.get_u2() - 0.0) < 0.0001)); - assert((std::fabs(mwu_test.get_zscore() - -86.8663) < 0.001)); + assert((std::fabs(mwu_test.get_zscore() - 86.8663) < 0.001)); assert((std::fabs(mwu_test.get_pvalue() - 0.0) < 0.0001)); ibm.single_act_visit("exponential", "lognormal", mwu_test); @@ -5354,13 +5342,13 @@ static void test_MutualInfoVisitor() { MutualInfoVisitor minfo; df.single_act_visit("int_col_1", "int_col_1", minfo); - assert((std::fabs(minfo.get_result() - 12.4866) < 0.0001)); + assert((std::fabs(minfo.get_result() - 2.89425) < 0.00001)); df.single_act_visit("int_col_1", "int_col_2", minfo); - assert((std::fabs(minfo.get_result() - 1.81499) < 0.00001)); + assert((std::fabs(minfo.get_result() - 1.22157) < 0.00001)); df.single_act_visit("int_col_1", "int_col_3", minfo); - assert((std::fabs(minfo.get_result() - 4.24521) < 0.00001)); + assert((std::fabs(minfo.get_result() - 0.954434) < 0.000001)); // Now multidimensional data // @@ -5403,7 +5391,7 @@ static void test_MutualInfoVisitor() { assert((std::fabs(mi_ary.get_result() - 1.58496) < 0.00001)); df.single_act_visit("COL VEC2", "COL VEC2", mi_vec); - assert((std::fabs(mi_vec.get_result() - 4.0) < 0.000000001)); + assert((std::fabs(mi_vec.get_result() - 1.0) < 0.000000001)); df.single_act_visit("COL VEC2", "COL VEC4", mi_vec); assert((std::fabs(mi_vec.get_result() - 1.0) < 0.000000001)); diff --git a/test/dataframe_tester_5.cc b/test/dataframe_tester_5.cc index 7dd4d5f6..cebd62d2 100644 --- a/test/dataframe_tester_5.cc +++ b/test/dataframe_tester_5.cc @@ -1450,18 +1450,18 @@ static void test_ARIMAVisitor() { const auto result1 = ari.get_result(); assert(result1.size() == 3); - assert(std::fabs(result1[0] - 247.175) < 0.001); - assert(std::fabs(result1[1] - 197.294) < 0.001); - assert(std::fabs(result1[2] - 220.021) < 0.001); + assert(std::fabs(result1[0] - 245.745) < 0.001); + assert(std::fabs(result1[1] - 197.122) < 0.001); + assert(std::fabs(result1[2] - 219.314) < 0.001); df.single_act_visit("oscil", ari); const auto result2 = ari.get_result(); assert(result2.size() == 3); - assert(std::fabs(result2[0] - 1.77088) < 0.00001); - assert(std::fabs(result2[1] - 1.67015) < 0.00001); - assert(std::fabs(result2[2] - 1.74417) < 0.00001); + assert(std::fabs(result2[0] - 1.76961) < 0.00001); + assert(std::fabs(result2[1] - 1.66931) < 0.00001); + assert(std::fabs(result2[2] - 1.74256) < 0.00001); try { df.single_act_visit("constant", ari); @@ -1475,18 +1475,18 @@ static void test_ARIMAVisitor() { const auto result3 = ari.get_result(); assert(result3.size() == 3); - assert(std::fabs(result3[0] - 14.3335) < 0.0001); - assert(std::fabs(result3[1] - 13.09) < 0.0001); - assert(std::fabs(result3[2] - 14.6469) < 0.0001); + assert(std::fabs(result3[0] - 14.1842) < 0.0001); + assert(std::fabs(result3[1] - 13.0033) < 0.0001); + assert(std::fabs(result3[2] - 14.4171) < 0.0001); df.single_act_visit("decreasing", ari); const auto result4 = ari.get_result(); assert(result4.size() == 3); - assert(std::fabs(result4[0] - 7.42058) < 0.00001); - assert(std::fabs(result4[1] - 7.21897) < 0.00001); - assert(std::fabs(result4[2] - 7.11158) < 0.00001); + assert(std::fabs(result4[0] - 7.40899) < 0.00001); + assert(std::fabs(result4[1] - 7.2049) < 0.0001); + assert(std::fabs(result4[2] - 7.09123) < 0.00001); // Now some real data // @@ -1568,9 +1568,9 @@ static void test_HWESForecastVisitor() { const auto result2 = hwes2.get_result(); assert(result2.size() == 3); - assert(std::fabs(result2[0] - 1.73499) < 0.00001); - assert(std::fabs(result2[1] - 1.9216) < 0.00001); - assert(std::fabs(result2[2] - 1.76383) < 0.00001); + assert(std::fabs(result2[0] - 1.90718) < 0.00001); + assert(std::fabs(result2[1] - 1.74941) < 0.00001); + assert(std::fabs(result2[2] - 1.93602) < 0.00001); df.single_act_visit("constant", hwes); diff --git a/test/dataframe_tester_output.txt b/test/dataframe_tester_output.txt index 036964dc..3ed6c937 100644 --- a/test/dataframe_tester_output.txt +++ b/test/dataframe_tester_output.txt @@ -714,18 +714,18 @@ col_2:1::8, col_3:1::15, col_str:1::11, col_4:1::22, -INDEX:1::123450, -col_1:1::1, -col_2:1::8, -col_3:1::15, -col_str:1::11, -col_4:1::22, -INDEX:1::123450, -col_1:1::1, -col_2:1::8, -col_3:1::15, -col_str:1::11, -col_4:1::22, +INDEX:1::123451, +col_1:1::2, +col_2:1::9, +col_3:1::16, +col_str:1::22, +col_4:1::23, +INDEX:1::123452, +col_1:1::3, +col_2:1::10, +col_3:1::17, +col_str:1::33, +col_4:1::24, Testing write(json) ... Writing in JSON: @@ -768,28 +768,28 @@ Means of clusters are: 19.9820684374, 3.85685151881, 0.606906002983, 1.840001468 4.39684180804 | 1.12423086285, 5.32645061833 | 0.258879189718, 6.56770612254 | 1.88980716816, 4.97200415894 | 0.407118210797, 6.03564845001 | 0.751012545119, 2.94035843782 | 0.314775456319, 3.60858208918 | 2.8089375563, 2.85547872035 | 2.51220235177, 3.89964990999 | 0.856374176369, 3.02735723184 | 0.226780241358, 3.11375978646 | 1.37885008492, 3.23563684881 | 1.32453935853, 4.74247908133 | 2.89363607052, 3.96514946019 | 0.387975516108, 4.4204216991 | 0.657322986031, 5.4094749177 | 0.376668387939, 3.21957447237 | 2.63364446483, 3.3948416627 | 0.163593376023, 3.84539096788 | 0.297100044069, 5.92619382132 | 1.84933460179, 4.590537026 | 2.53226329642, 3.06375316017 | 1.05043548764, 5.07407491355 | 2.45226502214, 6.28652144612 | 0.975087587592, 6.58207129189 | 0.355147337604, 3.72469838853 | 1.61156994808, 6.52684709038 | 1.96726515223, 2.9631288797 | 0.857151949309, 2.97138456376 | 0.629994862143, 5.90181207884 | 1.652460374, 3.05207906639 | 0.495742043608, 2.89162903925 | 1.03099415637, 2.75375629284 | 1.33757933345, 4.2123939732 | 2.11995054879, 4.7521600173 | 1.26866382559, 9.20566066884 | 0.635108361464, 3.35937247059 | 0.237096980877, 6.7470372648 | 2.90489612946, 5.56776407333 | 1.45780412902, 2.97283115952 | 0.768662732119, 5.29033890732 | 1.46684775613, 3.33442909318 | 1.41208702577, 3.3456478028 | 3.15791207885, 5.34911249975 | 0.310107651423, 6.26122042314 | 0.228863156957, 5.58341660812 | 0.652454865698, 4.0839718193 | 0.433244844548, 3.4840279873 | 0.588712471498, 5.48604913149 | 1.29273094631, 2.69170609551 | 0.557907104774, 4.07594518419 | 0.0794660417283, 4.22286388466 | 2.27121671294, 2.78656421487 | 1.09789223697, 3.68130821292 | 0.42822664425, 5.99835163287 | 1.64461791151, 3.77323745304 | 2.3901193323, 4.28719972597 | 0.741702952331, 3.48535661343 | 0.547992783008, 2.91570692268 | 0.297979629089, 4.22482418509 | 1.60766663614, 2.87894041417 | 0.210360178912, 3.81005611699 | 2.09964844162, 5.79399808209 | 0.665058152339, 4.83290584104 | 0.917001805312, 4.54824347301 | 1.28630290164, 2.7904401791 | 0.303055664379, 4.08422388761 | 1.52920833284, 4.597008972 | 0.552201184549, 3.44381095773 | 0.166042909241, 3.50438824207 | 1.92655877373, 4.00278453337 | 0.724330335054, 4.13312479796 | 3.83162584598, 2.98368298402 | 0.248108338937, 2.99413350468 | 0.106464616772, 4.57993344663 | 2.78540985725, 7.41604024389 | 0.452783882865, 5.34126437112 | 0.222824052156, 3.28587796229 | 0.595307265723, 3.12354995598 | 0.167846498981, 3.3611648224 | 0.838456623147, 3.42429119576 | 0.789557072651, 5.8512550284 | 1.99281532804, 3.91105879269 | 0.452095791245, 3.69740704803 | 0.820050676208, 3.1442635652 | 2.21862308472, 2.92284337303 | 0.620210327906, 3.36458898403 | 0.714739708796, 3.19489739931 | 0.350048969031, 7.88833318096 | 0.365468571658, 5.00308578152 | 0.223766784791, 3.22840497378 | 0.539154032626, 4.32414826891 | 0.342917766201, 4.2372602031 | 2.62312171147, 4.03489667364 | 1.19724040465, 2.98468107894 | 1.70816624101, 3.70061721998 | 0.144003322043, 3.32759851079 | 3.22759915673, 3.5348883024 | 2.26666238419, 3.61986027727 | 1.42234783125, 3.13105489274 | 0.42684491591, 4.13584607426 | 0.38240113167, 4.35065180729 | 0.450662984793, 5.29386720306 | 1.13522472017, 4.92212961708 | 0.504918122593, 4.19696666531 | 0.542106242523, 4.70822553306 | 0.431753357932, 4.70496853934 | 0.585882562077, 2.82163939708 | 0.354211377662, 2.7749285673 | 0.773819372569, 4.20106078629 | 3.5468186015, 3.25806230422 | 1.16446490776, 2.97824153046 | 0.634719609118, 6.24049625894 | 0.604503799466, 7.18692448088 | 0.702479476078, 3.9767986394 | 2.46691462078, 4.90972585114 | 2.12411456206, 3.25824555983 | 1.42599273813, 3.1842857736 | 0.780409710085, 7.33108959719 | 0.948593648814, 6.88914480747 | 1.18840666953, 3.93008795772 | 1.03634566272, 4.86100091822 | 0.460655730352, 5.64255515254 | 0.727288715906, 4.05696325215 | 0.596443341098, 4.12078782301 | 0.766925538593, 6.46611201826 | 1.85524070679, 4.94726132676 | 1.97043281815, 6.83805063598 | 0.444160900675, 4.86890068342 | 0.625975194945, 6.76331265223 | 0.224926532447, 8.96799021152 | 1.90689894212, 3.16719895489 | 0.646323589286, 4.93474307859 | 1.07218770687, 4.37542292896 | 0.884930238015, 5.81515843154 | 5.06933221904, -3[-8.25876466996|-1.07032737808|11.5596438002] +3[-8.73774733197|-1.93743664834|11.0190319172] -3[10.3976182234|-9.24973794974|-2.29062664521] +3[5.71303266874|-11.1714075284|-5.73628993344] -3[-10.937153444|-1.56477415998|-10.4109668341] +3[-7.65801427656|6.70019604884|-10.3360783645] -3[7.80765358278|12.6157913257|-0.495282470215] +3[11.212359209|9.74359344371|4.1733220999] Testing affinity propagation visitor ... 84.698691042, 7.74053662821, 24.698691042, 54.698691042, 2.8343340698, -48.1376687354, 49.1902139452, 62.6673194908, 48.8946164157, 43.0265110822, 63.9863705151, 45.1236798139, 62.3261150813, 62.7901547432, 49.2432766496, 62.193377464, 44.8974351335, 63.4788238344, 61.1971951939, 49.2006581232, 47.0600229553, 62.6635007764, 61.7111017446, 47.0088268057, 49.6187295396, 46.9186390843, 48.9752505523, 43.6763809821, 46.0378652093, 42.7901547432, 49.1370436351, 46.0736577489, 61.7126417006, 43.7555409884, 64.1005086332, 63.0766002353, 62.1982671268, 44.0322445851, 47.5318004395, 63.6763809821, 63.2682485746, 49.0569935418, 44.5917265576, 64.0322445851, 43.9863705151, 43.1343657534, 63.4428973212, 42.7757641221, 61.7248938756, 45.2267926641, 44.1005086332, 62.7757641221, 43.0766002353, 43.4788238344, 47.24429173, 60.4663675459, 63.1343657534, 63.0265110822, 43.2682485746, 63.7555409884, 46.034928453, 46.0128310734, 61.387824026, 60.3607367105, 45.6928812986, 62.2520045807, 60.6189653073, 43.4428973212, +3.44289732124, 6.01283107341, 5.22679266414, 7.00882680571, 9.20065812318, 6.0378652093, 0.618965307272, 4.59172655755, 7.0600229553, 3.75554098836, 4.1005086332, 3.02651108225, 3.67638098212, 3.47882383437, 0.466367545914, 4.03224458506, 1.71264170064, 2.79015474321, 2.66350077638, 2.19337746397, 7.5318004395, 4.89743513354, 3.9863705151, 8.13766873538, 1.19719519395, 7.24429172995, 2.32611508127, 6.03492845303, 6.07365774893, 3.07660023531, 9.19021394518, 2.77576412214, 5.69288129863, 8.97525055228, 2.25200458072, 0.36073671047, 8.89461641573, 9.13704363512, 1.38782402596, 6.91863908428, 9.24327664962, 9.61872953961, 3.13436575339, 9.05699354185, 1.72489387559, 2.19826712681, 3.26824857456, 5.12367981394, 2.66731949081, 1.71110174459, -3.44289732124, 0.618965307272, 4.59172655755, 3.75554098836, 4.1005086332, 3.02651108225, 3.67638098212, 3.47882383437, 0.466367545914, 4.03224458506, 1.71264170064, 2.79015474321, 2.66350077638, 2.19337746397, 3.9863705151, 1.19719519395, 2.32611508127, 3.07660023531, 2.77576412214, 2.25200458072, 0.36073671047, 1.38782402596, 3.13436575339, 1.72489387559, 2.19826712681, 3.26824857456, 2.66731949081, 1.71110174459, +80.3607367105, 86.0736577489, 82.6635007764, 84.5917265576, 86.0378652093, 83.6763809821, 80.4663675459, 89.0569935418, 89.2006581232, 81.1971951939, 82.7901547432, 84.1005086332, 85.2267926641, 82.6673194908, 80.6189653073, 84.8974351335, 86.0128310734, 82.2520045807, 88.9752505523, 88.8946164157, 89.6187295396, 88.1376687354, 81.7111017446, 86.9186390843, 84.0322445851, 83.7555409884, 82.193377464, 83.4428973212, 83.2682485746, 82.3261150813, 89.1902139452, 81.7126417006, 81.7248938756, 83.0766002353, 86.034928453, 87.0600229553, 83.0265110822, 83.1343657534, 89.2432766496, 81.387824026, 89.1370436351, 82.1982671268, 83.4788238344, 85.6928812986, 87.5318004395, 85.1236798139, 83.9863705151, 87.0088268057, 82.7757641221, 87.24429173, -28.8946164157, 25.1236798139, 29.0569935418, 29.2432766496, 41.7248938756, 40.3607367105, 28.1376687354, 25.2267926641, 42.6635007764, 29.1370436351, 26.0128310734, 29.2006581232, 26.0736577489, 27.24429173, 26.0378652093, 41.7111017446, 23.4788238344, 23.0766002353, 42.193377464, 26.9186390843, 27.0600229553, 23.1343657534, 41.387824026, 23.6763809821, 42.6673194908, 40.4663675459, 29.6187295396, 40.6189653073, 23.4428973212, 27.0088268057, 42.3261150813, 42.1982671268, 41.7126417006, 23.2682485746, 26.034928453, 23.9863705151, 23.7555409884, 29.1902139452, 27.5318004395, 24.1005086332, 42.2520045807, 24.5917265576, 24.8974351335, 41.1971951939, 24.0322445851, 25.6928812986, 28.9752505523, +62.6673194908, 63.9863705151, 67.0600229553, 62.3261150813, 62.7901547432, 66.034928453, 62.193377464, 63.4788238344, 61.1971951939, 69.0569935418, 62.6635007764, 61.7111017446, 69.6187295396, 66.9186390843, 67.24429173, 65.1236798139, 64.5917265576, 61.7126417006, 64.1005086332, 63.0766002353, 62.1982671268, 63.6763809821, 69.1370436351, 65.6928812986, 63.2682485746, 66.0378652093, 68.1376687354, 64.0322445851, 66.0736577489, 65.2267926641, 67.5318004395, 69.2006581232, 63.4428973212, 68.8946164157, 61.7248938756, 62.7757641221, 60.4663675459, 63.1343657534, 63.0265110822, 66.0128310734, 63.7555409884, 61.387824026, 68.9752505523, 60.3607367105, 69.1902139452, 67.0088268057, 62.2520045807, 64.8974351335, 60.6189653073, 69.2432766496, -6.01283107341, 5.22679266414, 7.00882680571, 9.20065812318, 6.0378652093, 7.0600229553, 22.1982671268, 21.7248938756, 7.5318004395, 4.89743513354, 21.7126417006, 8.13766873538, 7.24429172995, 22.3261150813, 22.193377464, 22.2520045807, 6.03492845303, 20.6189653073, 6.07365774893, 9.19021394518, 5.69288129863, 23.0265110822, 8.97525055228, 8.89461641573, 9.13704363512, 6.91863908428, 9.24327664962, 22.7901547432, 20.4663675459, 22.7757641221, 21.7111017446, 9.61872953961, 22.6635007764, 21.387824026, 22.6673194908, 9.05699354185, 20.3607367105, 21.1971951939, 5.12367981394, +48.1376687354, 49.1902139452, 48.8946164157, 43.0265110822, 45.1236798139, 49.2432766496, 44.8974351335, 41.7248938756, 49.2006581232, 47.0600229553, 47.0088268057, 49.6187295396, 46.9186390843, 40.3607367105, 48.9752505523, 42.6635007764, 43.6763809821, 46.0378652093, 42.7901547432, 49.1370436351, 46.0736577489, 43.7555409884, 44.0322445851, 41.7111017446, 47.5318004395, 42.193377464, 49.0569935418, 41.387824026, 44.5917265576, 43.9863705151, 42.6673194908, 40.4663675459, 43.1343657534, 40.6189653073, 42.7757641221, 42.3261150813, 42.1982671268, 41.7126417006, 45.2267926641, 44.1005086332, 43.0766002353, 43.4788238344, 47.24429173, 43.2682485746, 46.034928453, 46.0128310734, 45.6928812986, 42.2520045807, 41.1971951939, 43.4428973212, -80.3607367105, 86.0736577489, 67.0600229553, 82.6635007764, 84.5917265576, 86.0378652093, 66.034928453, 83.6763809821, 80.4663675459, 89.0569935418, 89.2006581232, 81.1971951939, 82.7901547432, 84.1005086332, 69.0569935418, 85.2267926641, 82.6673194908, 69.6187295396, 66.9186390843, 80.6189653073, 84.8974351335, 86.0128310734, 82.2520045807, 67.24429173, 88.9752505523, 88.8946164157, 65.1236798139, 89.6187295396, 88.1376687354, 64.5917265576, 81.7111017446, 86.9186390843, 69.1370436351, 84.0322445851, 83.7555409884, 82.193377464, 65.6928812986, 83.4428973212, 83.2682485746, 82.3261150813, 66.0378652093, 89.1902139452, 81.7126417006, 81.7248938756, 83.0766002353, 68.1376687354, 86.034928453, 87.0600229553, 83.0265110822, 66.0736577489, 65.2267926641, 67.5318004395, 69.2006581232, 83.1343657534, 68.8946164157, 89.2432766496, 81.387824026, 89.1370436351, 82.1982671268, 83.4788238344, 85.6928812986, 66.0128310734, 87.5318004395, 68.9752505523, 85.1236798139, 83.9863705151, 69.1902139452, 67.0088268057, 87.0088268057, 64.8974351335, 82.7757641221, 87.24429173, 69.2432766496, +28.8946164157, 25.1236798139, 29.0569935418, 29.2432766496, 22.1982671268, 21.7248938756, 28.1376687354, 25.2267926641, 21.7126417006, 29.1370436351, 22.3261150813, 26.0128310734, 29.2006581232, 26.0736577489, 22.193377464, 27.24429173, 22.2520045807, 26.0378652093, 20.6189653073, 23.4788238344, 23.0766002353, 26.9186390843, 23.0265110822, 27.0600229553, 23.1343657534, 23.6763809821, 22.7901547432, 20.4663675459, 29.6187295396, 22.7757641221, 23.4428973212, 27.0088268057, 21.7111017446, 23.2682485746, 22.6635007764, 21.387824026, 26.034928453, 23.9863705151, 22.6673194908, 23.7555409884, 29.1902139452, 27.5318004395, 20.3607367105, 24.1005086332, 24.5917265576, 21.1971951939, 24.8974351335, 24.0322445851, 25.6928812986, 28.9752505523, Testing multi-column sort ... @@ -1671,10 +1671,9 @@ Index Type Size: 8 Testing get_view_by_idx(values) ... Testing groupby( ) ... -INDEX:1:,dbl_col:1:,dbl_col_2:1:,str_col:1:,int_col:1:,bool_col:1: -1,0,100,zz,1,0 -INDEX:1:,bool_col:1:,sum_dbl2:1:,cnt_dbl2:1: -1,0,100,1 +INDEX:1:,dbl_col:1:,dbl_col_2:1:,str_col:1:,int_col:1:, +1,0,100,zz,1, +INDEX:0:,bool_col:0:,sum_dbl2:0:,cnt_dbl2:0: Testing concat_view( ) ... @@ -1711,8 +1710,8 @@ STD,19.9912242114,19.91396027,20.1165782354,20.0127428788,11.1522066877,5068928 MIN,94.599998,97.739998,90.559998,94.769997,92.30722,1193000 MAX,198.050003,199.210007,195.880005,197.770004,155.360657,30490200 25%,137.6949995,138.720001,136.794998,137.794998,123.625126,3097700 -50%,149.899994,151,148.5,149.630005,129.661377,3893400 -75%,162,163,161.179993,162.070007,136.192581,5063500 +50%,149.875,151,148.5,149.630005,129.6561355,3892500 +75%,162,163,161.144997,162.0650025,136.16098,5063200 Testing T3MovingMeanVisitor{ } ... @@ -2001,57 +2000,57 @@ Testing duplication_mask( ) ... Testing get_top_n_data( ) ... INDEX:4::123453,123454,123456,123462, -col_1:4::4,5,7,13, -col_2:4::11,12,14,32, -col_3:4::18,19,21,19, +col_1:4::4.000000000000,5.000000000000,7.000000000000,13.000000000000, +col_2:4::11.000000000000,12.000000000000,14.000000000000,32.000000000000, +col_3:4::18.000000000000,19.000000000000,21.000000000000,19.000000000000, col_4:4::25,99,0,0, INDEX:4::123453,123454,123456,123462, -col_1:4::4,5,7,13, -col_2:4::11,12,14,32, -col_3:4::18,19,21,19, +col_1:4::4.000000000000,5.000000000000,7.000000000000,13.000000000000, +col_2:4::11.000000000000,12.000000000000,14.000000000000,32.000000000000, +col_3:4::18.000000000000,19.000000000000,21.000000000000,19.000000000000, col_4:4::25,99,0,0, Testing get_bottom_n_data( ) ... INDEX:4::123457,123458,123459,123461, -col_1:4::8,9,10,12, -col_2:4::20,22,23,31, -col_3:4::0.34,1.56,0.34,0.34, +col_1:4::8.000000000000,9.000000000000,10.000000000000,12.000000000000, +col_2:4::20.000000000000,22.000000000000,23.000000000000,31.000000000000, +col_3:4::0.340000000000,1.560000000000,0.340000000000,0.340000000000, col_4:4::0,0,0,0, INDEX:4::123457,123458,123459,123461, -col_1:4::8,9,10,12, -col_2:4::20,22,23,31, -col_3:4::0.34,1.56,0.34,0.34, +col_1:4::8.000000000000,9.000000000000,10.000000000000,12.000000000000, +col_2:4::20.000000000000,22.000000000000,23.000000000000,31.000000000000, +col_3:4::0.340000000000,1.560000000000,0.340000000000,0.340000000000, col_4:4::0,0,0,0, Testing get_above_quantile_data( ) ... INDEX:8::123450,123451,123452,123453,123454,123455,123456,123462, -col_1:8::1,2,3,4,5,6,7,13, -col_2:8::8,9,10,11,12,13,14,32, -col_3:8::15,16,15,18,19,16,21,19, +col_1:8::1.000000000000,2.000000000000,3.000000000000,4.000000000000,5.000000000000,6.000000000000,7.000000000000,13.000000000000, +col_2:8::8.000000000000,9.000000000000,10.000000000000,11.000000000000,12.000000000000,13.000000000000,14.000000000000,32.000000000000, +col_3:8::15.000000000000,16.000000000000,15.000000000000,18.000000000000,19.000000000000,16.000000000000,21.000000000000,19.000000000000, col_4:8::22,23,24,25,99,0,0,0, INDEX:8::123450,123451,123452,123453,123454,123455,123456,123462, -col_1:8::1,2,3,4,5,6,7,13, -col_2:8::8,9,10,11,12,13,14,32, -col_3:8::15,16,15,18,19,16,21,19, +col_1:8::1.000000000000,2.000000000000,3.000000000000,4.000000000000,5.000000000000,6.000000000000,7.000000000000,13.000000000000, +col_2:8::8.000000000000,9.000000000000,10.000000000000,11.000000000000,12.000000000000,13.000000000000,14.000000000000,32.000000000000, +col_3:8::15.000000000000,16.000000000000,15.000000000000,18.000000000000,19.000000000000,16.000000000000,21.000000000000,19.000000000000, col_4:8::22,23,24,25,99,0,0,0, Testing get_below_quantile_data( ) ... INDEX:6::123457,123458,123459,123460,123461,123466, -col_1:6::8,9,10,11,12,14, -col_2:6::20,22,23,30,31,1.89, -col_3:6::0.34,1.56,0.34,2.3,0.34,10, +col_1:6::8.000000000000,9.000000000000,10.000000000000,11.000000000000,12.000000000000,14.000000000000, +col_2:6::20.000000000000,22.000000000000,23.000000000000,30.000000000000,31.000000000000,1.890000000000, +col_3:6::0.340000000000,1.560000000000,0.340000000000,2.300000000000,0.340000000000,10.000000000000, col_4:6::0,0,0,0,0,0, INDEX:6::123457,123458,123459,123460,123461,123466, -col_1:6::8,9,10,11,12,14, -col_2:6::20,22,23,30,31,1.89, -col_3:6::0.34,1.56,0.34,2.3,0.34,10, +col_1:6::8.000000000000,9.000000000000,10.000000000000,11.000000000000,12.000000000000,14.000000000000, +col_2:6::20.000000000000,22.000000000000,23.000000000000,30.000000000000,31.000000000000,1.890000000000, +col_3:6::0.340000000000,1.560000000000,0.340000000000,2.300000000000,0.340000000000,10.000000000000, col_4:6::0,0,0,0,0,0, @@ -3031,3 +3030,53 @@ Testing HWESForecastVisitor{ } ... Testing LSTMForecastVisitor{ } ... Testing kshape_groups( ) ... + +Testing count( ) ... + +Testing class_count( ) ... + +Testing AnomalyDetectByKNNVisitor{ } ... + +Testing BIRCHVisitor{ } ... + +Testing get_data_by_birch( ) ... + +Testing get_md_stats( ) ... + +Testing KrigingVisitor{ } ... + +Testing asof_join( ) ... + +Testing crosstab( ) ... + +Testing pivot_table( ) ... + +Testing JarqueBeraTestVisitor{ } ... + +Testing LjungBoxTestVisitor{ } ... + +Testing DurbinWatsonVisitor{ } ... + +Testing SilhouetteScoreVisitor{ } ... + +Testing DaviesBouldinIndexVisitor{ } ... + +Testing CalinskiHarabaszVisitor{ } ... + +Testing AnomalyDetectByIsoForestVisitor{ } ... + +Testing DivergenceVisitor{ } ... + +Testing GradientVisitor{ } ... + +Testing JacobianVisitor{ } ... + +Testing LaplacianVisitor{ } ... + +Testing remove_data_by_isof( ) ... + +Testing read_chunked_data( ) ... + +Testing streamed_write( ) ... + +Testing StreamAppender{ } ...