diff --git a/include/sparseir/augment.hpp b/include/sparseir/augment.hpp
new file mode 100644
index 0000000..21c2799
--- /dev/null
+++ b/include/sparseir/augment.hpp
@@ -0,0 +1,166 @@
+#pragma once
+
+#include <Eigen/Dense>
+#include <cmath>
+#include <complex>
+#include <limits>
+#include <memory>
+#include <stdexcept>
+#include <vector>
+
+namespace sparseir {
+
+template <typename T>
+class AbstractAugmentation {
+public:
+    virtual T operator()(T tau) const = 0; // Ensure this is virtual
+    virtual T deriv(T tau, int n = 1) const = 0; // Add this line
+    virtual std::complex<T> hat(int n) const = 0; // Add this line
+    virtual ~AbstractAugmentation() = default;
+    virtual std::unique_ptr<AbstractAugmentation<T>> clone() const = 0;
+};
+
+template <typename S>
+class AugmentedBasis : public AbstractBasis<S> {
+public:
+    AugmentedBasis(std::shared_ptr<AbstractBasis<S>> basis,
+                   const std::vector<std::unique_ptr<AbstractAugmentation<S>>>& augmentations)
+        : _basis(basis), _augmentations(augmentations), _naug(augmentations.size()) {
+        //Error Handling: Check for null basis pointer
+        if (!_basis) {
+            throw std::invalid_argument("Basis cannot be null");
+        }
+        //Check for valid augmentations
+        for (const auto& aug : _augmentations) {
+            if (!aug) {
+                throw std::invalid_argument("Augmentation cannot be null");
+            }
+        }
+    }
+
+    size_t size() const { return _basis->size() + _naug; }
+
+    Eigen::VectorXd u(const Eigen::VectorXd& tau) const {
+        Eigen::VectorXd result(size());
+        for (size_t i = 0; i < _naug; ++i) {
+            result(i) = (*_augmentations[i])(tau(i));
+        }
+        for (size_t i = _naug; i < size(); ++i) {
+            result(i) = _basis->u(tau(i - _naug))(i - _naug);
+        }
+        return result;
+    }
+
+    Eigen::VectorXcf uhat(const Eigen::VectorXcf& wn) const {
+        Eigen::VectorXcf result(size());
+        for (size_t i = 0; i < _naug; ++i) {
+            result(i) = (*_augmentations[i]).hat(wn(i));
+        }
+        for (size_t i = _naug; i < size(); ++i) {
+            result(i) = _basis->uhat(wn(i - _naug))(i - _naug);
+        }
+        return result;
+    }
+
+    Eigen::VectorXd v(const Eigen::VectorXd& w) const {
+        return _basis->v(w);
+    }
+
+    Eigen::VectorXd s() const { return _basis->s(); }
+
+    double beta() const { return _basis->beta(); }
+
+    double wmax() const { return _basis->wmax(); }
+
+    std::shared_ptr<Statistics> statistics() const {
+        return _basis->statistics(); // Assuming _basis also returns a shared_ptr
+    }
+
+private:
+    std::shared_ptr<AbstractBasis<S>> _basis;
+    std::vector<std::unique_ptr<AbstractAugmentation<S>>> _augmentations;
+    size_t _naug;
+};
+
+
+template <typename T>
+class TauConst : public AbstractAugmentation<T> {
+public:
+    TauConst(T beta) : beta_(beta), norm_(1.0 / std::sqrt(beta)) {
+        if (beta_ <= 0) {
+            throw std::invalid_argument("beta must be positive");
+        }
+    }
+
+    T operator()(T tau) const override {
+        check_tau_range(tau, beta_);
+        return norm_;
+    }
+    T deriv(T tau, int n = 1) const  {
+        if (n == 0) return (*this)(tau);
+        return 0.0;
+    }
+    std::complex<T> hat(int n) const  {
+        return norm_ * std::sqrt(beta_);
+    }
+    std::unique_ptr<AbstractAugmentation<T>> clone() const override {
+        return std::make_unique<TauConst<T>>(*this);
+    }
+
+private:
+    T beta_;
+    T norm_;
+};
+
+template <typename T>
+class TauLinear : public AbstractAugmentation<T> {
+public:
+    TauLinear(T beta) : beta_(beta), norm_(std::sqrt(3.0 / beta)) {
+        if (beta_ <= 0) {
+            throw std::invalid_argument("beta must be positive");
+        }
+    }
+
+    T operator()(T tau) const override {
+        check_tau_range(tau, beta_);
+        return norm_ * (2.0 * tau / beta_ - 1.0);
+    }
+    T deriv(T tau, int n = 1) const override {
+        if (n == 1) return norm_ * 2.0 / beta_;
+        return 0.0;
+    }
+    std::complex<T> hat(int n) const override {
+        if (n == 0) return 0.0;
+        return norm_ * 2.0 * beta_ / (n * M_PI * std::complex<T>(0, 1));
+    }
+    std::unique_ptr<AbstractAugmentation<T>> clone() const override {
+        return std::make_unique<TauLinear<T>>(*this);
+    }
+
+private:
+    T beta_;
+    T norm_;
+};
+
+template <typename T>
+class MatsubaraConst : public AbstractAugmentation<T> {
+public:
+    MatsubaraConst(T beta) : beta_(beta) {
+        if (beta_ <= 0) {
+            throw std::invalid_argument("beta must be positive");
+        }
+    }
+
+    T operator()(T tau) const override { return std::numeric_limits<T>::quiet_NaN(); }
+    T deriv(T tau, int n = 1) const override { return std::numeric_limits<T>::quiet_NaN(); }
+    std::complex<T> hat(int n) const override { return 1.0; }
+    std::unique_ptr<AbstractAugmentation<T>> clone() const override {
+        return std::make_unique<MatsubaraConst<T>>(*this);
+    }
+
+private:
+    T beta_;
+};
+
+
+} // namespace sparseir
diff --git a/include/sparseir/basis.hpp b/include/sparseir/basis.hpp
index a840693..1148014 100644
--- a/include/sparseir/basis.hpp
+++ b/include/sparseir/basis.hpp
@@ -303,7 +303,7 @@ class FiniteTempBasis : public AbstractBasis<S> {
 
 private:
     // Placeholder statistics function
-    S statistics() const { return S(); }
+    std::shared_ptr<S> statistics() const { return std::make_shared<S>(); }
 
     // Default sampling points function
     Eigen::VectorXd
diff --git a/include/sparseir/poly.hpp b/include/sparseir/poly.hpp
index d88737f..7b8504d 100644
--- a/include/sparseir/poly.hpp
+++ b/include/sparseir/poly.hpp
@@ -872,6 +872,7 @@ inline PowerModel<double> power_model(const Statistics &stat,
 class BosonicStatistics : public Statistics {
 public:
     int zeta() const override { return 1; }
+    bool allowed(int n) const override { return n % 2 == 0 && n != 0; }
 };
 
 // PiecewiseLegendreFT class template
@@ -1000,6 +1001,7 @@ std::complex<double> PiecewiseLegendreFT<StatisticsType>::evalpoly(
 class FermionicStatistics : public Statistics {
 public:
     int zeta() const override { return -1; }
+    bool allowed(int n) const override { return n % 2 != 0; }
 };
 
 } // namespace sparseir
diff --git a/include/sparseir/sparseir-header-only.hpp b/include/sparseir/sparseir-header-only.hpp
index 9672402..3e21d2e 100644
--- a/include/sparseir/sparseir-header-only.hpp
+++ b/include/sparseir/sparseir-header-only.hpp
@@ -11,4 +11,4 @@
 #include "./kernel.hpp"
 #include "./sve.hpp"
 #include "./basis.hpp"
-
+#include "./augment.hpp"
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index c61c7e8..8de152c 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -25,6 +25,7 @@ add_executable(libsparseirtests
     svd.cxx
     sve.cxx
     basis.cxx
+    augment.cxx
 )
 
 target_link_libraries(libsparseirtests PRIVATE Catch2::Catch2WithMain)
diff --git a/test/augment.cxx b/test/augment.cxx
new file mode 100644
index 0000000..a65a1c5
--- /dev/null
+++ b/test/augment.cxx
@@ -0,0 +1,129 @@
+#include <Eigen/Dense>
+#include <algorithm>
+#include <cmath>
+#include <complex>
+#include <limits>
+#include <memory>
+#include <stdexcept>
+#include <vector>
+
+#include <catch2/catch_test_macros.hpp>
+
+#include <sparseir/sparseir-header-only.hpp>
+#include <xprec/ddouble-header-only.hpp>
+
+using namespace sparseir;
+using namespace std;
+
+TEST_CASE("Augmented bosonic basis") {
+    using T = double;
+    T beta = 1000.0;
+    T wmax = 2.0;
+
+    // Create bosonic basis
+    auto kernel = LogisticKernel(beta * wmax);
+    auto sve_result = compute_sve(kernel, 1e-6);
+    auto basis = make_shared<FiniteTempBasis<Bosonic>>(beta, wmax, 1e-6, kernel, sve_result);
+    // Create augmented basis with TauConst and TauLinear
+    /*
+    vector<unique_ptr<AbstractAugmentation<T>>> augmentations;
+    //augmentations.push_back(make_unique<TauConst<T>>(beta));
+    //augmentations.push_back(make_unique<TauLinear<T>>(beta));
+    AugmentedBasis<T> basis_aug(basis, augmentations);
+
+    REQUIRE(basis_aug.size() == basis->size() + 2);
+
+    // Define G(τ) = c - exp(-τ * pole) / (1 - exp(-β * pole))
+    T pole = 1.0;
+    T c = 1e-2;
+
+    // Create tau sampling points
+    auto tau_sampling = make_tau_sampling(basis_aug);
+    auto tau = tau_sampling.tau;
+    Eigen::VectorX<T> gtau(tau.size());
+    for (size_t i = 0; i < tau.size(); ++i) {
+        gtau(i) = c - exp(-tau(i) * pole) / (1 - exp(-beta * pole));
+    }
+    T magn = gtau.maxCoeff();
+
+    // This illustrates that "naive" fitting is a problem if the fitting matrix
+    // is not well-conditioned.
+    Eigen::MatrixXd tau_matrix = tau_sampling.matrix;
+    Eigen::VectorX<T> gl_fit_bad = tau_matrix.completeOrthogonalDecomposition().solve(gtau);
+    Eigen::VectorX<T> gtau_reconst_bad = tau_matrix * gl_fit_bad;
+    REQUIRE(!gtau_reconst_bad.isApprox(gtau, 1e-13 * magn));
+    REQUIRE(gtau_reconst_bad.isApprox(gtau, 5e-16 * tau_matrix.norm() * magn));
+    REQUIRE(tau_matrix.norm() > 1e7);
+    REQUIRE(tau_matrix.rows() == basis_aug.size());
+    REQUIRE(tau_matrix.cols() == tau.size());
+
+    // Now do the fit properly
+    Eigen::VectorX<T> gl_fit = tau_matrix.colPivHouseholderQr().solve(gtau);
+    Eigen::VectorX<T> gtau_reconst = tau_matrix * gl_fit;
+
+    REQUIRE(gtau_reconst.isApprox(gtau, 1e-14 * magn));
+    */
+}
+
+
+TEST_CASE("Vertex basis with stat = $stat", "[augment]") {
+     /*
+     for (const shared_ptr<Statistics>& stat : {make_shared<Fermionic>(), make_shared<Bosonic>()}) {
+        using T = double;
+        T beta = 1000.0;
+        T wmax = 2.0;
+        auto basis = make_shared<FiniteTempBasis<T>>(stat, beta, wmax, 1e-6);
+        vector<unique_ptr<AbstractAugmentation<T>>> augmentations;
+        augmentations.push_back(make_unique<MatsubaraConst<T>>(beta));
+        AugmentedBasis<T> basis_aug(basis, augmentations);
+        REQUIRE(!basis_aug.uhat.empty());
+
+        // G(iν) = c + 1 / (iν - pole)
+        T pole = 1.0;
+        T c = 1.0;
+        auto matsu_sampling = make_matsubara_sampling(basis_aug);
+        Eigen::VectorXcf gi_n(matsu_sampling.wn.size());
+        for (size_t i = 0; i < matsu_sampling.wn.size(); ++i) {
+            complex<T> iwn(0, matsu_sampling.wn(i));
+            gi_n(i) = c + 1.0 / (iwn - pole);
+        }
+        Eigen::VectorXcf gl = fit(matsu_sampling, gi_n);
+        Eigen::VectorXcf gi_n_reconst = evaluate(matsu_sampling, gl);
+        REQUIRE(gi_n_reconst.isApprox(gi_n, gi_n.maxCoeff() * 1e-7));
+    }
+}
+
+
+TEST_CASE("unit tests", "[augment]") {
+    using T = double;
+    T beta = 1000.0;
+    T wmax = 2.0;
+    auto basis = make_shared<FiniteTempBasis<T>>(Bosonic(), beta, wmax, 1e-6);
+    vector<unique_ptr<AbstractAugmentation<T>>> augmentations;
+    augmentations.push_back(make_unique<TauConst<T>>(beta));
+    augmentations.push_back(make_unique<TauLinear<T>>(beta));
+    AugmentedBasis<T> basis_aug(basis, augmentations);
+
+    SECTION("getindex") {
+        REQUIRE(basis_aug.u.size() == basis_aug.size());
+        REQUIRE(basis_aug.u[0]->operator()(0.0) == basis_aug[0](0.0));
+        REQUIRE(basis_aug.u[1]->operator()(0.0) == basis_aug[1](0.0));
+    }
+
+    size_t len_basis = basis->size();
+    size_t len_aug = len_basis + 2;
+
+    REQUIRE(basis_aug.size() == len_aug);
+    // REQUIRE(basis_aug.accuracy == basis->accuracy);
+    REQUIRE(basis_aug.Lambda() == beta * wmax);
+    REQUIRE(basis_aug.wmax() == wmax);
+
+    REQUIRE(basis_aug.u.size() == len_aug);
+    REQUIRE(basis_aug.u(0.8).size() == len_aug);
+    //REQUIRE(basis_aug.uhat(MatsubaraFreq<T>(4.0)).size() == len_aug);
+    REQUIRE(basis_aug.u.minCoeff() == 0.0);
+    REQUIRE(basis_aug.u.maxCoeff() == beta);
+
+    //Further tests omitted for brevity,  adapt as needed from Julia code.
+    */
+}
\ No newline at end of file