remove stl-cpp submodule

contrib/CMakeLists.txt

@@ -77,7 +77,6 @@ add_contrib (minizip-ng-cmake minizip-ng)
 add_contrib (snappy-cmake snappy)
 add_contrib (rocksdb-cmake rocksdb)
 add_contrib (thrift-cmake thrift)
-add_contrib (stl-cpp-cmake stl-cpp)
 # parquet/arrow/orc
 add_contrib (arrow-cmake arrow) # requires: snappy, thrift, double-conversion
 add_contrib (avro-cmake avro) # requires: snappy

contrib/stl-cpp

@@ -1 +0,0 @@
-Subproject commit f3e1cd02b4f0341a7539ed1a171b21b14b37ed51

contrib/stl-cpp-cmake/CMakeLists.txt

@@ -1,9 +0,0 @@
-option (ENABLE_SEASONAL "Enable seasonal decompose (stl-cpp) library" ${ENABLE_LIBRARIES})
-
-if (NOT ENABLE_SEASONAL)
-    message(STATUS "Not using seasonal decompose (stl-cpp)")
-    return()
-endif()
-add_library(_stl-cpp INTERFACE)
-target_include_directories(_stl-cpp INTERFACE ${ClickHouse_SOURCE_DIR}/contrib/stl-cpp/include)
-add_library(ch_contrib::seasonal-decompose ALIAS _stl-cpp)

src/CMakeLists.txt

@@ -396,10 +396,6 @@ if (TARGET ch_contrib::nuraft)
     target_link_libraries (clickhouse_compression PUBLIC ch_contrib::nuraft)
 endif()
 
-if (TARGET ch_contrib::seasonal-decompose)
-    target_link_libraries(clickhouse_common_io PUBLIC ch_contrib::seasonal-decompose)
-endif ()
-
 dbms_target_link_libraries (
     PRIVATE
         boost::filesystem

src/Functions/CMakeLists.txt

@@ -75,10 +75,6 @@ if (TARGET ch_contrib::base64)
     list (APPEND PRIVATE_LIBS ch_contrib::base64)
 endif()
 
-if (TARGET ch_contrib::seasonal-decompose)
-    list (APPEND PRIVATE_LIBS ch_contrib::seasonal-decompose)
-endif()
-
 if (ENABLE_NLP)
     list (APPEND PRIVATE_LIBS ch_contrib::cld2)
 endif()

src/Functions/seriesDecomposeSTL.cpp

@@ -1,6 +1,5 @@
 #include "config.h"
 
-#if USE_SEASONAL
 #ifdef __clang__
 #pragma clang diagnostic push
 #pragma clang diagnostic ignored "-Wold-style-cast"
@@ -8,7 +7,7 @@
 #pragma clang diagnostic ignored "-Wimplicit-float-conversion"
 #endif
 
-#include <stl.hpp>
+#include <Functions/stl.hpp>
 
 #ifdef __clang__
 #pragma clang diagnostic pop
@@ -69,83 +68,92 @@ public:
         const ColumnArray * array = checkAndGetColumn<ColumnArray>(array_ptr.get());
 
         const IColumn & src_data = array->getData();
+        const ColumnArray::Offsets & src_offsets = array->getOffsets();
 
         Float64 period;
 
+        auto ret = ColumnFloat32::create();
+        auto & res_data = ret->getData();
 
-        auto period_ptr = arguments[1].column->convertToFullColumnIfConst();
-        if (checkAndGetColumn<ColumnUInt8>(period_ptr.get()) || checkAndGetColumn<ColumnUInt16>(period_ptr.get())
-            || checkAndGetColumn<ColumnUInt32>(period_ptr.get()) || checkAndGetColumn<ColumnUInt64>(period_ptr.get()))
-            period = period_ptr->getUInt(0);
-        else if (checkAndGetColumn<ColumnFloat32>(period_ptr.get()) || checkAndGetColumn<ColumnFloat64>(period_ptr.get()))
+        ColumnArray::ColumnOffsets::MutablePtr col_offsets = ColumnArray::ColumnOffsets::create();
+        auto & col_offsets_data = col_offsets->getData();
+
+        auto root_offsets = ColumnArray::ColumnOffsets::create();
+        auto & root_offsets_data = root_offsets->getData();
+
+        ColumnArray::Offset prev_src_offset = 0;
+
+        for (size_t i = 0; i < src_offsets.size(); ++i)
         {
-            period = period_ptr->getFloat64(0);
-            if (isNaN(period) || !std::isfinite(period) || period < 0)
+            auto period_ptr = arguments[1].column->convertToFullColumnIfConst();
+            if (checkAndGetColumn<ColumnUInt8>(period_ptr.get()) || checkAndGetColumn<ColumnUInt16>(period_ptr.get())
+                || checkAndGetColumn<ColumnUInt32>(period_ptr.get()) || checkAndGetColumn<ColumnUInt64>(period_ptr.get()))
+                period = period_ptr->getUInt(i);
+            else if (checkAndGetColumn<ColumnFloat32>(period_ptr.get()) || checkAndGetColumn<ColumnFloat64>(period_ptr.get()))
+            {
+                period = period_ptr->getFloat64(i);
+                if (isNaN(period) || !std::isfinite(period) || period < 0)
+                    throw Exception(
+                        ErrorCodes::ILLEGAL_COLUMN,
+                        "Illegal value {} for second argument of function {}. Should be a positive number",
+                        period,
+                        getName());
+            }
+            else
                 throw Exception(
                     ErrorCodes::ILLEGAL_COLUMN,
-                    "Illegal value {} for second argument of function {}. Should be a positive number",
-                    period,
+                    "Illegal column {} of second argument of function {}",
+                    arguments[1].column->getName(),
+                    getName());
+
+
+            std::vector<Float32> seasonal;
+            std::vector<Float32> trend;
+            std::vector<Float32> residue;
+
+            ColumnArray::Offset curr_offset = src_offsets[i];
+
+            if (executeNumber<UInt8>(src_data, period, prev_src_offset, curr_offset, seasonal, trend, residue)
+                || executeNumber<UInt16>(src_data, period, prev_src_offset, curr_offset, seasonal, trend, residue)
+                || executeNumber<UInt32>(src_data, period, prev_src_offset, curr_offset, seasonal, trend, residue)
+                || executeNumber<UInt64>(src_data, period, prev_src_offset, curr_offset, seasonal, trend, residue)
+                || executeNumber<Int8>(src_data, period, prev_src_offset, curr_offset, seasonal, trend, residue)
+                || executeNumber<Int16>(src_data, period, prev_src_offset, curr_offset, seasonal, trend, residue)
+                || executeNumber<Int32>(src_data, period, prev_src_offset, curr_offset, seasonal, trend, residue)
+                || executeNumber<Int64>(src_data, period, prev_src_offset, curr_offset, seasonal, trend, residue)
+                || executeNumber<Float32>(src_data, period, prev_src_offset, curr_offset, seasonal, trend, residue)
+                || executeNumber<Float64>(src_data, period, prev_src_offset, curr_offset, seasonal, trend, residue))
+            {
+                res_data.insert(res_data.end(), seasonal.begin(), seasonal.end());
+                col_offsets_data.push_back(res_data.size());
+
+                res_data.insert(res_data.end(), trend.begin(), trend.end());
+                col_offsets_data.push_back(res_data.size());
+
+                res_data.insert(res_data.end(), residue.begin(), residue.end());
+                col_offsets_data.push_back(res_data.size());
+
+                root_offsets_data.push_back(col_offsets->size());
+
+                prev_src_offset = curr_offset;
+            }
+            else
+                throw Exception(
+                    ErrorCodes::ILLEGAL_COLUMN,
+                    "Illegal column {} of first argument of function {}",
+                    arguments[0].column->getName(),
                     getName());
         }
-        else
-            throw Exception(
-                ErrorCodes::ILLEGAL_COLUMN,
-                "Illegal column {} of second argument of function {}",
-                arguments[1].column->getName(),
-                getName());
-
-
-        std::vector<Float32> seasonal;
-        std::vector<Float32> trend;
-        std::vector<Float32> residue;
-
-
-        if (executeNumber<UInt8>(src_data, period, seasonal, trend, residue)
-            || executeNumber<UInt16>(src_data, period, seasonal, trend, residue)
-            || executeNumber<UInt32>(src_data, period, seasonal, trend, residue)
-            || executeNumber<UInt64>(src_data, period, seasonal, trend, residue)
-            || executeNumber<Int8>(src_data, period, seasonal, trend, residue)
-            || executeNumber<Int16>(src_data, period, seasonal, trend, residue)
-            || executeNumber<Int32>(src_data, period, seasonal, trend, residue)
-            || executeNumber<Int64>(src_data, period, seasonal, trend, residue)
-            || executeNumber<Float32>(src_data, period, seasonal, trend, residue)
-            || executeNumber<Float64>(src_data, period, seasonal, trend, residue))
-        {
-            auto ret = ColumnFloat32::create();
-            auto & res_data = ret->getData();
-
-            ColumnArray::ColumnOffsets::MutablePtr col_offsets = ColumnArray::ColumnOffsets::create();
-            auto & col_offsets_data = col_offsets->getData();
-
-            auto root_offsets = ColumnArray::ColumnOffsets::create();
-            auto & root_offsets_data = root_offsets->getData();
-
-            res_data.insert(res_data.end(), seasonal.begin(), seasonal.end());
-            col_offsets_data.push_back(res_data.size());
-
-            res_data.insert(res_data.end(), trend.begin(), trend.end());
-            col_offsets_data.push_back(res_data.size());
-
-            res_data.insert(res_data.end(), residue.begin(), residue.end());
-            col_offsets_data.push_back(res_data.size());
-
-            root_offsets_data.push_back(col_offsets->size());
-
-            ColumnArray::MutablePtr nested_array_col = ColumnArray::create(std::move(ret), std::move(col_offsets));
-            return ColumnArray::create(std::move(nested_array_col), std::move(root_offsets));
-        }
-        else
-            throw Exception(
-                ErrorCodes::ILLEGAL_COLUMN,
-                "Illegal column {} of first argument of function {}",
-                arguments[0].column->getName(),
-                getName());
+        ColumnArray::MutablePtr nested_array_col = ColumnArray::create(std::move(ret), std::move(col_offsets));
+        return ColumnArray::create(std::move(nested_array_col), std::move(root_offsets));
     }
 
     template <typename T>
     bool executeNumber(
         const IColumn & src_data,
         Float64 period,
+        ColumnArray::Offset & start, 
+        ColumnArray::Offset & end,
         std::vector<Float32> & seasonal,
         std::vector<Float32> & trend,
         std::vector<Float32> & residue) const
@@ -156,14 +164,15 @@ public:
 
         const PaddedPODArray<T> & src_vec = src_data_concrete->getData();
 
-        size_t len = src_vec.size();
+        chassert(start <= end);
+        size_t len = end - start;
         if (len < 4)
             throw Exception(ErrorCodes::BAD_ARGUMENTS, "At least four data points are needed for function {}", getName());
         else if (period > (len / 2))
             throw Exception(
                 ErrorCodes::BAD_ARGUMENTS, "The series should have data of at least two period lengths for function {}", getName());
 
-        std::vector<float> src(src_vec.begin(), src_vec.end());
+        std::vector<float> src((src_vec.begin() + start), (src_vec.begin() + end));
 
         try
         {
@@ -191,4 +200,3 @@ Decompose time series data based on STL(Seasonal-Trend Decomposition Procedure B
         .categories{"Time series analysis"}});
 }
 }
-#endif

src/Functions/stl.hpp

@@ -0,0 +1,510 @@
+/*!
+ * STL C++ v0.1.3
+ * https://github.com/ankane/stl-cpp
+ * Unlicense OR MIT License
+ *
+ * Ported from https://www.netlib.org/a/stl
+ *
+ * Cleveland, R. B., Cleveland, W. S., McRae, J. E., & Terpenning, I. (1990).
+ * STL: A Seasonal-Trend Decomposition Procedure Based on Loess.
+ * Journal of Official Statistics, 6(1), 3-33.
+ */
+
+#pragma once
+
+#include <algorithm>
+#include <cmath>
+#include <numeric>
+#include <optional>
+#include <stdexcept>
+#include <vector>
+
+namespace stl {
+
+bool est(const float* y, size_t n, size_t len, int ideg, float xs, float* ys, size_t nleft, size_t nright, float* w, bool userw, const float* rw) {
+    auto range = ((float) n) - 1.0;
+    auto h = std::max(xs - ((float) nleft), ((float) nright) - xs);
+
+    if (len > n) {
+        h += (float) ((len - n) / 2);
+    }
+
+    auto h9 = 0.999 * h;
+    auto h1 = 0.001 * h;
+
+    // compute weights
+    auto a = 0.0;
+    for (auto j = nleft; j <= nright; j++) {
+        w[j - 1] = 0.0;
+        auto r = fabs(((float) j) - xs);
+        if (r <= h9) {
+            if (r <= h1) {
+                w[j - 1] = 1.0;
+            } else {
+                w[j - 1] = pow(1.0 - pow(r / h, 3), 3);
+            }
+            if (userw) {
+                w[j - 1] *= rw[j - 1];
+            }
+            a += w[j - 1];
+        }
+    }
+
+    if (a <= 0.0) {
+        return false;
+    } else { // weighted least squares
+        for (auto j = nleft; j <= nright; j++) { // make sum of w(j) == 1
+            w[j - 1] /= a;
+        }
+
+        if (h > 0.0 && ideg > 0) { // use linear fit
+            auto a = 0.0;
+            for (auto j = nleft; j <= nright; j++) { // weighted center of x values
+                a += w[j - 1] * ((float) j);
+            }
+            auto b = xs - a;
+            auto c = 0.0;
+            for (auto j = nleft; j <= nright; j++) {
+                c += w[j - 1] * pow(((float) j) - a, 2);
+            }
+            if (sqrt(c) > 0.001 * range) {
+                b /= c;
+
+                // points are spread out enough to compute slope
+                for (auto j = nleft; j <= nright; j++) {
+                    w[j - 1] *= b * (((float) j) - a) + 1.0;
+                }
+            }
+        }
+
+        *ys = 0.0;
+        for (auto j = nleft; j <= nright; j++) {
+            *ys += w[j - 1] * y[j - 1];
+        }
+
+        return true;
+    }
+}
+
+void ess(const float* y, size_t n, size_t len, int ideg, size_t njump, bool userw, const float* rw, float* ys, float* res) {
+    if (n < 2) {
+        ys[0] = y[0];
+        return;
+    }
+
+    size_t nleft = 0;
+    size_t nright = 0;
+
+    auto newnj = std::min(njump, n - 1);
+    if (len >= n) {
+        nleft = 1;
+        nright = n;
+        for (size_t i = 1; i <= n; i += newnj) {
+            auto ok = est(y, n, len, ideg, (float) i, &ys[i - 1], nleft, nright, res, userw, rw);
+            if (!ok) {
+                ys[i - 1] = y[i - 1];
+            }
+        }
+    } else if (newnj == 1) { // newnj equal to one, len less than n
+        auto nsh = (len + 1) / 2;
+        nleft = 1;
+        nright = len;
+        for (size_t i = 1; i <= n; i++) { // fitted value at i
+            if (i > nsh && nright != n) {
+                nleft += 1;
+                nright += 1;
+            }
+            auto ok = est(y, n, len, ideg, (float) i, &ys[i - 1], nleft, nright, res, userw, rw);
+            if (!ok) {
+                ys[i - 1] = y[i - 1];
+            }
+        }
+    } else { // newnj greater than one, len less than n
+        auto nsh = (len + 1) / 2;
+        for (size_t i = 1; i <= n; i += newnj) { // fitted value at i
+            if (i < nsh) {
+                nleft = 1;
+                nright = len;
+            } else if (i >= n - nsh + 1) {
+                nleft = n - len + 1;
+                nright = n;
+            } else {
+                nleft = i - nsh + 1;
+                nright = len + i - nsh;
+            }
+            auto ok = est(y, n, len, ideg, (float) i, &ys[i - 1], nleft, nright, res, userw, rw);
+            if (!ok) {
+                ys[i - 1] = y[i - 1];
+            }
+        }
+    }
+
+    if (newnj != 1) {
+        for (size_t i = 1; i <= n - newnj; i += newnj) {
+            auto delta = (ys[i + newnj - 1] - ys[i - 1]) / ((float) newnj);
+            for (auto j = i + 1; j <= i + newnj - 1; j++) {
+                ys[j - 1] = ys[i - 1] + delta * ((float) (j - i));
+            }
+        }
+        auto k = ((n - 1) / newnj) * newnj + 1;
+        if (k != n) {
+            auto ok = est(y, n, len, ideg, (float) n, &ys[n - 1], nleft, nright, res, userw, rw);
+            if (!ok) {
+                ys[n - 1] = y[n - 1];
+            }
+            if (k != n - 1) {
+                auto delta = (ys[n - 1] - ys[k - 1]) / ((float) (n - k));
+                for (auto j = k + 1; j <= n - 1; j++) {
+                    ys[j - 1] = ys[k - 1] + delta * ((float) (j - k));
+                }
+            }
+        }
+    }
+}
+
+void ma(const float* x, size_t n, size_t len, float* ave) {
+    auto newn = n - len + 1;
+    auto flen = (float) len;
+    auto v = 0.0;
+
+    // get the first average
+    for (size_t i = 0; i < len; i++) {
+        v += x[i];
+    }
+
+    ave[0] = v / flen;
+    if (newn > 1) {
+        auto k = len;
+        auto m = 0;
+        for (size_t j = 1; j < newn; j++) {
+            // window down the array
+            v = v - x[m] + x[k];
+            ave[j] = v / flen;
+            k += 1;
+            m += 1;
+        }
+    }
+}
+
+void fts(const float* x, size_t n, size_t np, float* trend, float* work) {
+    ma(x, n, np, trend);
+    ma(trend, n - np + 1, np, work);
+    ma(work, n - 2 * np + 2, 3, trend);
+}
+
+void rwts(const float* y, size_t n, const float* fit, float* rw) {
+    for (size_t i = 0; i < n; i++) {
+        rw[i] = fabs(y[i] - fit[i]);
+    }
+
+    auto mid1 = (n - 1) / 2;
+    auto mid2 = n / 2;
+
+    // sort
+    std::sort(rw, rw + n);
+
+    auto cmad = 3.0 * (rw[mid1] + rw[mid2]); // 6 * median abs resid
+    auto c9 = 0.999 * cmad;
+    auto c1 = 0.001 * cmad;
+
+    for (size_t i = 0; i < n; i++) {
+        auto r = fabs(y[i] - fit[i]);
+        if (r <= c1) {
+            rw[i] = 1.0;
+        } else if (r <= c9) {
+            rw[i] = pow(1.0 - pow(r / cmad, 2), 2);
+        } else {
+            rw[i] = 0.0;
+        }
+    }
+}
+
+void ss(const float* y, size_t n, size_t np, size_t ns, int isdeg, size_t nsjump, bool userw, float* rw, float* season, float* work1, float* work2, float* work3, float* work4) {
+    for (size_t j = 1; j <= np; j++) {
+        size_t k = (n - j) / np + 1;
+
+        for (size_t i = 1; i <= k; i++) {
+            work1[i - 1] = y[(i - 1) * np + j - 1];
+        }
+        if (userw) {
+            for (size_t i = 1; i <= k; i++) {
+                work3[i - 1] = rw[(i - 1) * np + j - 1];
+            }
+        }
+        ess(work1, k, ns, isdeg, nsjump, userw, work3, work2 + 1, work4);
+        auto xs = 0.0;
+        auto nright = std::min(ns, k);
+        auto ok = est(work1, k, ns, isdeg, xs, &work2[0], 1, nright, work4, userw, work3);
+        if (!ok) {
+            work2[0] = work2[1];
+        }
+        xs = k + 1;
+        size_t nleft = std::max(1, (int) k - (int) ns + 1);
+        ok = est(work1, k, ns, isdeg, xs, &work2[k + 1], nleft, k, work4, userw, work3);
+        if (!ok) {
+            work2[k + 1] = work2[k];
+        }
+        for (size_t m = 1; m <= k + 2; m++) {
+            season[(m - 1) * np + j - 1] = work2[m - 1];
+        }
+    }
+}
+
+void onestp(const float* y, size_t n, size_t np, size_t ns, size_t nt, size_t nl, int isdeg, int itdeg, int ildeg, size_t nsjump, size_t ntjump, size_t nljump, size_t ni, bool userw, float* rw, float* season, float* trend, float* work1, float* work2, float* work3, float* work4, float* work5) {
+    for (size_t j = 0; j < ni; j++) {
+        for (size_t i = 0; i < n; i++) {
+            work1[i] = y[i] - trend[i];
+        }
+
+        ss(work1, n, np, ns, isdeg, nsjump, userw, rw, work2, work3, work4, work5, season);
+        fts(work2, n + 2 * np, np, work3, work1);
+        ess(work3, n, nl, ildeg, nljump, false, work4, work1, work5);
+        for (size_t i = 0; i < n; i++) {
+            season[i] = work2[np + i] - work1[i];
+        }
+        for (size_t i = 0; i < n; i++) {
+            work1[i] = y[i] - season[i];
+        }
+        ess(work1, n, nt, itdeg, ntjump, userw, rw, trend, work3);
+    }
+}
+
+void stl(const float* y, size_t n, size_t np, size_t ns, size_t nt, size_t nl, int isdeg, int itdeg, int ildeg, size_t nsjump, size_t ntjump, size_t nljump, size_t ni, size_t no, float* rw, float* season, float* trend) {
+    if (ns < 3) {
+        throw std::invalid_argument("seasonal_length must be at least 3");
+    }
+    if (nt < 3) {
+        throw std::invalid_argument("trend_length must be at least 3");
+    }
+    if (nl < 3) {
+        throw std::invalid_argument("low_pass_length must be at least 3");
+    }
+    if (np < 2) {
+        throw std::invalid_argument("period must be at least 2");
+    }
+
+    if (isdeg != 0 && isdeg != 1) {
+        throw std::invalid_argument("seasonal_degree must be 0 or 1");
+    }
+    if (itdeg != 0 && itdeg != 1) {
+        throw std::invalid_argument("trend_degree must be 0 or 1");
+    }
+    if (ildeg != 0 && ildeg != 1) {
+        throw std::invalid_argument("low_pass_degree must be 0 or 1");
+    }
+
+    if (ns % 2 != 1) {
+        throw std::invalid_argument("seasonal_length must be odd");
+    }
+    if (nt % 2 != 1) {
+        throw std::invalid_argument("trend_length must be odd");
+    }
+    if (nl % 2 != 1) {
+        throw std::invalid_argument("low_pass_length must be odd");
+    }
+
+    auto work1 = std::vector<float>(n + 2 * np);
+    auto work2 = std::vector<float>(n + 2 * np);
+    auto work3 = std::vector<float>(n + 2 * np);
+    auto work4 = std::vector<float>(n + 2 * np);
+    auto work5 = std::vector<float>(n + 2 * np);
+
+    auto userw = false;
+    size_t k = 0;
+
+    while (true) {
+        onestp(y, n, np, ns, nt, nl, isdeg, itdeg, ildeg, nsjump, ntjump, nljump, ni, userw, rw, season, trend, work1.data(), work2.data(), work3.data(), work4.data(), work5.data());
+        k += 1;
+        if (k > no) {
+            break;
+        }
+        for (size_t i = 0; i < n; i++) {
+            work1[i] = trend[i] + season[i];
+        }
+        rwts(y, n, work1.data(), rw);
+        userw = true;
+    }
+
+    if (no <= 0) {
+        for (size_t i = 0; i < n; i++) {
+            rw[i] = 1.0;
+        }
+    }
+}
+
+float var(const std::vector<float>& series) {
+    auto mean = std::accumulate(series.begin(), series.end(), 0.0) / series.size();
+    std::vector<float> tmp;
+    tmp.reserve(series.size());
+    for (auto v : series) {
+        tmp.push_back(pow(v - mean, 2));
+    }
+    return std::accumulate(tmp.begin(), tmp.end(), 0.0) / (series.size() - 1);
+}
+
+float strength(const std::vector<float>& component, const std::vector<float>& remainder) {
+    std::vector<float> sr;
+    sr.reserve(remainder.size());
+    for (size_t i = 0; i < remainder.size(); i++) {
+        sr.push_back(component[i] + remainder[i]);
+    }
+    return std::max(0.0, 1.0 - var(remainder) / var(sr));
+}
+
+class StlResult {
+public:
+    std::vector<float> seasonal;
+    std::vector<float> trend;
+    std::vector<float> remainder;
+    std::vector<float> weights;
+
+    inline float seasonal_strength() {
+        return strength(seasonal, remainder);
+    }
+
+    inline float trend_strength() {
+        return strength(trend, remainder);
+    }
+};
+
+class StlParams {
+    std::optional<size_t> ns_ = std::nullopt;
+    std::optional<size_t> nt_ = std::nullopt;
+    std::optional<size_t> nl_ = std::nullopt;
+    int isdeg_ = 0;
+    int itdeg_ = 1;
+    std::optional<int> ildeg_ = std::nullopt;
+    std::optional<size_t> nsjump_ = std::nullopt;
+    std::optional<size_t> ntjump_ = std::nullopt;
+    std::optional<size_t> nljump_ = std::nullopt;
+    std::optional<size_t> ni_ = std::nullopt;
+    std::optional<size_t> no_ = std::nullopt;
+    bool robust_ = false;
+
+public:
+    inline StlParams seasonal_length(size_t ns) {
+        this->ns_ = ns;
+        return *this;
+    }
+
+    inline StlParams trend_length(size_t nt) {
+        this->nt_ = nt;
+        return *this;
+    }
+
+    inline StlParams low_pass_length(size_t nl) {
+        this->nl_ = nl;
+        return *this;
+    }
+
+    inline StlParams seasonal_degree(int isdeg) {
+        this->isdeg_ = isdeg;
+        return *this;
+    }
+
+    inline StlParams trend_degree(int itdeg) {
+        this->itdeg_ = itdeg;
+        return *this;
+    }
+
+    inline StlParams low_pass_degree(int ildeg) {
+        this->ildeg_ = ildeg;
+        return *this;
+    }
+
+    inline StlParams seasonal_jump(size_t nsjump) {
+        this->nsjump_ = nsjump;
+        return *this;
+    }
+
+    inline StlParams trend_jump(size_t ntjump) {
+        this->ntjump_ = ntjump;
+        return *this;
+    }
+
+    inline StlParams low_pass_jump(size_t nljump) {
+        this->nljump_ = nljump;
+        return *this;
+    }
+
+    inline StlParams inner_loops(bool ni) {
+        this->ni_ = ni;
+        return *this;
+    }
+
+    inline StlParams outer_loops(bool no) {
+        this->no_ = no;
+        return *this;
+    }
+
+    inline StlParams robust(bool robust) {
+        this->robust_ = robust;
+        return *this;
+    }
+
+    StlResult fit(const float* y, size_t n, size_t np);
+    StlResult fit(const std::vector<float>& y, size_t np);
+};
+
+StlParams params() {
+    return StlParams();
+}
+
+StlResult StlParams::fit(const float* y, size_t n, size_t np) {
+    if (n < 2 * np) {
+        throw std::invalid_argument("series has less than two periods");
+    }
+
+    auto ns = this->ns_.value_or(np);
+
+    auto isdeg = this->isdeg_;
+    auto itdeg = this->itdeg_;
+
+    auto res = StlResult {
+        std::vector<float>(n),
+        std::vector<float>(n),
+        std::vector<float>(),
+        std::vector<float>(n)
+    };
+
+    auto ildeg = this->ildeg_.value_or(itdeg);
+    auto newns = std::max(ns, (size_t) 3);
+    if (newns % 2 == 0) {
+        newns += 1;
+    }
+
+    auto newnp = std::max(np, (size_t) 2);
+    auto nt = (size_t) ceil((1.5 * newnp) / (1.0 - 1.5 / (float) newns));
+    nt = this->nt_.value_or(nt);
+    nt = std::max(nt, (size_t) 3);
+    if (nt % 2 == 0) {
+        nt += 1;
+    }
+
+    auto nl = this->nl_.value_or(newnp);
+    if (nl % 2 == 0 && !this->nl_.has_value()) {
+        nl += 1;
+    }
+
+    auto ni = this->ni_.value_or(this->robust_ ? 1 : 2);
+    auto no = this->no_.value_or(this->robust_ ? 15 : 0);
+
+    auto nsjump = this->nsjump_.value_or((size_t) ceil(((float) newns) / 10.0));
+    auto ntjump = this->ntjump_.value_or((size_t) ceil(((float) nt) / 10.0));
+    auto nljump = this->nljump_.value_or((size_t) ceil(((float) nl) / 10.0));
+
+    stl(y, n, newnp, newns, nt, nl, isdeg, itdeg, ildeg, nsjump, ntjump, nljump, ni, no, res.weights.data(), res.seasonal.data(), res.trend.data());
+
+    res.remainder.reserve(n);
+    for (size_t i = 0; i < n; i++) {
+        res.remainder.push_back(y[i] - res.seasonal[i] - res.trend[i]);
+    }
+
+    return res;
+}
+
+StlResult StlParams::fit(const std::vector<float>& y, size_t np) {
+    return StlParams::fit(y.data(), y.size(), np);
+}
+
+}

src/configure_config.cmake

@@ -161,9 +161,6 @@ endif ()
 if (ENABLE_OPENSSL)
     set(USE_OPENSSL_INTREE 1)
 endif ()
-if (TARGET ch_contrib::seasonal-decompose)
-    set(USE_SEASONAL 1)
-endif()
 if (TARGET ch_contrib::fiu)
     set(FIU_ENABLE 1)
 endif()