Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Resolve roots #50

Merged
merged 6 commits into from
Dec 13, 2024
Merged

Resolve roots #50

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
8103ed1
Resolve roots
terasakisatoshi Dec 12, 2024
5426fee
Suppress output
terasakisatoshi Dec 12, 2024
c5ad724
Fix
terasakisatoshi Dec 12, 2024
0596cce
use legval in _specfuncs.hpp
terasakisatoshi Dec 12, 2024
2b1bbfb
Fix sve
terasakisatoshi Dec 12, 2024
1c695f6
Fix bug in refine_grid
shinaoka Dec 13, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 2 additions & 1 deletion CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,7 @@ project(

# Eigen3
find_package (Eigen3 3.4 REQUIRED NO_MODULE)
find_package(Threads REQUIRED)

# libxprec
include(FetchContent)
Expand Down Expand Up @@ -52,7 +53,7 @@ set_target_properties(sparseir PROPERTIES
CXX_STANDARD_REQUIRED ON
)

target_link_libraries(sparseir Eigen3::Eigen XPrec::xprec)
target_link_libraries(sparseir Eigen3::Eigen XPrec::xprec Threads::Threads)

# Use library convention.
add_library(SparseIR::sparseir ALIAS sparseir)
Expand Down
2 changes: 1 addition & 1 deletion include/sparseir/_root.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -235,7 +235,7 @@ std::vector<double> discrete_extrema(F f, const std::vector<double> &xgrid)
std::vector<double> res;
for (size_t i = 0; i < a.size(); ++i) {
// abs ∘ f
auto abf = [f](double x) { return std::abs(f(x)); };
auto abf = [f](double x) { return std::fabs(f(x)); };
res.push_back(
bisect_discr_extremum(abf, a[i], b[i], absf_a[i], absf_b[i]));
}
Expand Down
85 changes: 14 additions & 71 deletions include/sparseir/poly.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -146,7 +146,9 @@ class PiecewiseLegendrePoly {
double x_tilde;
std::tie(i, x_tilde) = split(x);
Eigen::VectorXd coeffs = data.col(i);
double value = legval(x_tilde, coeffs) * norms[i];
// convert coeffs to std::vector<double>
std::vector<double> coeffs_vec(coeffs.data(), coeffs.data() + coeffs.size());
double value = legval<double>(x_tilde, coeffs_vec) * norms[i];
return value;
}

Expand Down Expand Up @@ -215,7 +217,7 @@ class PiecewiseLegendrePoly {


Eigen::VectorXd refine_grid(const Eigen::VectorXd& grid, int alpha) const {
Eigen::VectorXd refined(grid.size() * alpha);
Eigen::VectorXd refined((grid.size() - 1) * alpha + 1);

for (size_t i = 0; i < grid.size() - 1; ++i) {
double start = grid[i];
Expand Down Expand Up @@ -245,80 +247,19 @@ class PiecewiseLegendrePoly {
}
}

// Equivalent to find_all function in julia above
Eigen::VectorXd find_all(const Eigen::VectorXd& xgrid) const {
Eigen::VectorXd fx(xgrid.size());
for (size_t i = 0; i < xgrid.size(); ++i) {
fx(i) = static_cast<double>((*this)(xgrid[i]));
}
// Find direct hits (zeros)
std::vector<bool> hit(fx.size(), false);
std::vector<double> x_hit;
for (Eigen::Index i = 0; i < fx.size(); ++i) {
hit[i] = (fx(i) == 0.0);
if (hit[i]) {
x_hit.push_back(xgrid(i));
}
}

// Check for sign changes
std::vector<bool> sign_change(fx.size() - 1, false);
bool found_sign_change = false;

for (Eigen::Index i = 0; i < fx.size() - 1; ++i) {
bool different_signs = std::signbit(fx(i)) != std::signbit(fx(i + 1));
bool neither_zero = fx(i) != 0.0 && fx(i + 1) != 0.0;
sign_change[i] = different_signs && neither_zero;
if (sign_change[i]) {
found_sign_change = true;
}
}

if (!found_sign_change) {
// convert to Eigen::VectorXd
return Eigen::Map<Eigen::VectorXd>(x_hit.data(), x_hit.size());
}

// Collect points for bisection
std::vector<double> a, b;
std::vector<double> fa;

for (size_t i = 0; i < sign_change.size(); ++i) {
if (sign_change[i]) {
a.push_back(xgrid(i));
b.push_back(xgrid(i + 1));
fa.push_back(fx(i));
}
}

// Calculate epsilon for floating point types
double eps_x;
if (std::is_floating_point<double>::value) {
eps_x = std::numeric_limits<double>::epsilon() * xgrid.cwiseAbs().maxCoeff();
} else {
eps_x = 0;
}

// Perform bisection for each interval
std::vector<double> x_bisect;
for (size_t i = 0; i < a.size(); ++i) {
x_bisect.push_back(bisect(a[i], b[i], fa[i], eps_x));
}

// Combine and sort results
x_hit.insert(x_hit.end(), x_bisect.begin(), x_bisect.end());
std::sort(x_hit.begin(), x_hit.end());
// convert to Eigen::VectorXd
return Eigen::Map<Eigen::VectorXd>(x_hit.data(), x_hit.size());
}

// Roots function
Eigen::VectorXd roots(double tol = 1e-10) const
{
Eigen::VectorXd grid = this->knots;

Eigen::VectorXd refined_grid = refine_grid(grid, 2);
std::cout << "refined_grid: " << refined_grid.size() << "\n";
return find_all(refined_grid);
auto f = [this](double x) { return this->operator()(x); };
// convert to std::vector<double>
std::vector<double> refined_grid_vec(refined_grid.data(),
refined_grid.data() +
refined_grid.size());
std::vector<double> roots = find_all(f, refined_grid_vec);
return Eigen::Map<Eigen::VectorXd>(roots.data(), roots.size());
}

// Overloaded operators
Expand Down Expand Up @@ -366,6 +307,7 @@ class PiecewiseLegendrePoly {
int get_polyorder() const { return polyorder; }

private:
/*
// Helper function to compute legval
static double legval(double x, const Eigen::VectorXd &coeffs)
{
Expand All @@ -385,6 +327,7 @@ class PiecewiseLegendrePoly {
}
return result;
}
*/

// Helper function to split x into segment index i and x_tilde
std::pair<int, double> split(double x) const
Expand Down
2 changes: 1 addition & 1 deletion include/sparseir/sve.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -532,7 +532,7 @@ template <typename K>
auto compute_sve(K kernel,
double epsilon = std::numeric_limits<double>::quiet_NaN(),
double cutoff = std::numeric_limits<double>::quiet_NaN(),
std::string Twork = "Floatt64",
std::string Twork = "Float64",
int lmax = std::numeric_limits<int>::max(), int n_gauss = -1,
const std::string &svd_strat = "auto")
{
Expand Down
26 changes: 14 additions & 12 deletions test/poly.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -299,8 +299,8 @@ TEST_CASE("Overlap")
TEST_CASE("Roots")
{
// Initialize data and knots (from Julia code)
Eigen::MatrixXd data(16, 2);
data << 0.16774734206553019, 0.49223680914312595, -0.8276728567928646,
Eigen::VectorXd v(32);
v << 0.16774734206553019, 0.49223680914312595, -0.8276728567928646,
0.16912891046582143, -0.0016231275318572044, 0.00018381683946452256,
-9.699355027805034e-7, 7.60144228530804e-8, -2.8518324490258146e-10,
1.7090590205708293e-11, -5.0081401126025e-14, 2.1244236198427895e-15,
Expand All @@ -311,32 +311,34 @@ TEST_CASE("Roots")
2.8518324490258146e-10, 1.7090590205708293e-11, 5.0081401126025e-14,
2.1244236198427895e-15, -2.0478095258000225e-16, -2.676573801530628e-16,
-2.338165820094204e-16, -1.2050663212312096e-16;
//data.resize(16, 2);

Eigen::Map<Eigen::MatrixXd> data(v.data(), 16, 2);

Eigen::VectorXd knots(3);
knots << 0.0, 0.5, 1.0;
int l = 3;

sparseir::PiecewiseLegendrePoly pwlp(data, knots, l);
/*

// Find roots
Eigen::VectorXd roots = pwlp.roots();

// Print roots for debugging
std::cout << "Found roots: " << roots.transpose() << "\n";

/*
// Expected roots from Julia
Eigen::VectorXd expected_roots(3);
expected_roots << 0.1118633448586015, 0.4999999999999998, 0.8881366551413985;

// REQUIRE(roots.size() == expected_roots.size());
REQUIRE(roots.size() == expected_roots.size());
for(Eigen::Index i = 0; i < roots.size(); ++i) {
//REQUIRE(std::abs(roots[i] - expected_roots[i]) < 1e-10);
REQUIRE(std::fabs(roots[i] - expected_roots[i]) < 1e-10);
// Verifying the polynomial evaluates to zero
//at the roots with tight tolerance
REQUIRE(std::fabs(pwlp(roots[i])) < 1e-10);
// Verify roots are in domain
//REQUIRE(roots[i] >= knots[0]);
//REQUIRE(roots[i] <= knots[knots.size()-1]);
REQUIRE(roots[i] >= knots[0]);
REQUIRE(roots[i] <= knots[knots.size() - 1]);
// Verify these are actually roots
//REQUIRE(std::abs(pwlp(roots[i])) < 1e-10);
REQUIRE(std::abs(pwlp(roots[i])) < 1e-10);
}
*/
}
26 changes: 12 additions & 14 deletions test/sve.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -114,22 +114,20 @@ TEST_CASE("CentrosymmSVE", "[CentrosymmSVE]")

TEST_CASE("sve.cpp", "[compute_sve]")
{
// auto sve =
// sparseir::compute_sve<sparseir::LogisticKernel>(sparseir::LogisticKernel(10.0));

// Define a map to store SVEResult objects
// auto sve_logistic = std::map < int,
// sparseir::SVEResult<sparseir::LogisticKernel>>{
// {10,
// sparseir::compute_sve<sparseir::LogisticKernel>(sparseir::LogisticKernel(10.0))},
// {42,
// sparseir::compute_sve<sparseir::LogisticKernel>(sparseir::LogisticKernel(42.0))},
// {10000,
// sparseir::compute_sve<sparseir::LogisticKernel>(sparseir::LogisticKernel(10000.0))},
// //{100000000,
// sparseir::compute_sve<sparseir::LogisticKernel>(sparseir::LogisticKernel(10000.0),
// 1e-12)},
// };
auto sve_logistic = std::map < int,
sparseir::SVEResult<sparseir::LogisticKernel>>{
{10,
sparseir::compute_sve<sparseir::LogisticKernel>(sparseir::LogisticKernel(10.0))},
{42,
sparseir::compute_sve<sparseir::LogisticKernel>(sparseir::LogisticKernel(42.0))},
{10000,
sparseir::compute_sve<sparseir::LogisticKernel>(sparseir::LogisticKernel(10000.0))},
//{100000000,
// sparseir::compute_sve<sparseir::LogisticKernel>(sparseir::LogisticKernel(10000.0),
// 1e-12)},
};

SECTION("smooth with Λ =")
{
Expand Down
Loading