From bdd766c820d672b4ab0d7651974a5ef8c201b08a Mon Sep 17 00:00:00 2001 From: Niklas Goellmann Date: Mon, 1 Dec 2025 16:50:21 +0100 Subject: [PATCH 1/4] added possibility for population analysis based on grid --- src/trees/FunctionTree.cpp | 22 ++++++++++++++++++++++ src/trees/FunctionTree.h | 1 + src/utils/CompFunction.cpp | 9 +++++++++ src/utils/CompFunction.h | 1 + 4 files changed, 33 insertions(+) diff --git a/src/trees/FunctionTree.cpp b/src/trees/FunctionTree.cpp index 99557843..b33caa39 100644 --- a/src/trees/FunctionTree.cpp +++ b/src/trees/FunctionTree.cpp @@ -435,6 +435,28 @@ template T FunctionTree::integrate() const { return jacobian * result; } +/** @returns Integral of the function over parts of the computational domain */ +template T FunctionTree::integrate(int dim, bool greater) const { + + T result = 0.0; + for (int i = 0; i < this->rootBox.size(); i++) { + const FunctionNode &fNode = getRootFuncNode(i); + if (fNode.getUpperBounds()[dim] <= 0.0 && not greater) + result += fNode.integrate(); + if (fNode.getLowerBounds()[dim] >= 0.0 && greater) + result += fNode.integrate(); + } + + // Handle potential scaling + auto scaling_factor = this->getMRA().getWorldBox().getScalingFactors(); + auto jacobian = 1.0; + for (const auto &sf_i : scaling_factor) jacobian *= std::sqrt(sf_i); + // Square root of scaling factor in each diection. The seemingly missing + // multiplication by the square root of sf_i is included in the basis + + return jacobian * result; +} + /** @returns Integral of a representable function over the grid given by the tree */ template <> double FunctionTree<3, double>::integrateEndNodes(RepresentableFunction_M &f) { // traverse tree, and treat end nodes only diff --git a/src/trees/FunctionTree.h b/src/trees/FunctionTree.h index 9d976d6b..f016da4b 100644 --- a/src/trees/FunctionTree.h +++ b/src/trees/FunctionTree.h @@ -62,6 +62,7 @@ template class FunctionTree final : public MWTree, pub ~FunctionTree() override; T integrate() const; + T integrate(int dim, bool greater) const; double integrateEndNodes(RepresentableFunction_M &f); T evalf_precise(const Coord &r); T evalf(const Coord &r) const override; diff --git a/src/utils/CompFunction.cpp b/src/utils/CompFunction.cpp index 250796c2..55cfca44 100644 --- a/src/utils/CompFunction.cpp +++ b/src/utils/CompFunction.cpp @@ -187,6 +187,15 @@ template ComplexDouble CompFunction::integrate() const { return integral; } +template ComplexDouble CompFunction::integrate(int dim, bool greater) const { + ComplexDouble integral; + if (isreal()) + integral = CompD[0]->integrate(dim, greater); + else + integral = CompC[0]->integrate(dim, greater); + return integral; +} + template double CompFunction::norm() const { double norm = getSquareNorm(); if (norm > 0.0) norm = std::sqrt(norm); diff --git a/src/utils/CompFunction.h b/src/utils/CompFunction.h index 3b63a63b..0ca59fad 100644 --- a/src/utils/CompFunction.h +++ b/src/utils/CompFunction.h @@ -126,6 +126,7 @@ template class CompFunction { CompFunction paramCopy(bool alloc = false) const; ComplexDouble integrate() const; + ComplexDouble integrate(int dim, bool greater) const; double norm() const; double getSquareNorm() const; void alloc(int nalloc = 1, bool zero = true); From d0312d69aaba9352eaf31113176ce8470e6f39f5 Mon Sep 17 00:00:00 2001 From: Niklas Goellmann Date: Tue, 2 Dec 2025 10:18:18 +0100 Subject: [PATCH 2/4] minor cleanup --- src/trees/FunctionTree.cpp | 13 ++++++++----- src/trees/FunctionTree.h | 2 +- src/utils/CompFunction.cpp | 12 +++++++++--- src/utils/CompFunction.h | 2 +- 4 files changed, 19 insertions(+), 10 deletions(-) diff --git a/src/trees/FunctionTree.cpp b/src/trees/FunctionTree.cpp index b33caa39..7751f2eb 100644 --- a/src/trees/FunctionTree.cpp +++ b/src/trees/FunctionTree.cpp @@ -435,15 +435,18 @@ template T FunctionTree::integrate() const { return jacobian * result; } -/** @returns Integral of the function over parts of the computational domain */ -template T FunctionTree::integrate(int dim, bool greater) const { - +/** + * @brief Integrate over half of the space + * @param dim Dimension along which to split + * @param largerSide If true, integrate over the larger side (x>0 if dim=0) + * @return Integral of the function over parts of the computational domain */ +template T FunctionTree::integrate(int dim, bool largerSide) const { T result = 0.0; for (int i = 0; i < this->rootBox.size(); i++) { const FunctionNode &fNode = getRootFuncNode(i); - if (fNode.getUpperBounds()[dim] <= 0.0 && not greater) + if (fNode.getUpperBounds()[dim] <= 0.0 && not largerSide) result += fNode.integrate(); - if (fNode.getLowerBounds()[dim] >= 0.0 && greater) + if (fNode.getLowerBounds()[dim] >= 0.0 && largerSide) result += fNode.integrate(); } diff --git a/src/trees/FunctionTree.h b/src/trees/FunctionTree.h index f016da4b..6f8c9770 100644 --- a/src/trees/FunctionTree.h +++ b/src/trees/FunctionTree.h @@ -62,7 +62,7 @@ template class FunctionTree final : public MWTree, pub ~FunctionTree() override; T integrate() const; - T integrate(int dim, bool greater) const; + T integrate(int dim, bool largerSide) const; double integrateEndNodes(RepresentableFunction_M &f); T evalf_precise(const Coord &r); T evalf(const Coord &r) const override; diff --git a/src/utils/CompFunction.cpp b/src/utils/CompFunction.cpp index 55cfca44..eb2dcba2 100644 --- a/src/utils/CompFunction.cpp +++ b/src/utils/CompFunction.cpp @@ -187,12 +187,18 @@ template ComplexDouble CompFunction::integrate() const { return integral; } -template ComplexDouble CompFunction::integrate(int dim, bool greater) const { +/** + * @brief Integrate over half of the space + * @param dim Dimension along which to split + * @param largerSide If true, integrate over the larger side (x>0 if dim=0) + * @return ComplexDouble Resulting integral + */ +template ComplexDouble CompFunction::integrate(int dim, bool largerSide) const { ComplexDouble integral; if (isreal()) - integral = CompD[0]->integrate(dim, greater); + integral = CompD[0]->integrate(dim, largerSide); else - integral = CompC[0]->integrate(dim, greater); + integral = CompC[0]->integrate(dim, largerSide); return integral; } diff --git a/src/utils/CompFunction.h b/src/utils/CompFunction.h index 0ca59fad..70d6ade4 100644 --- a/src/utils/CompFunction.h +++ b/src/utils/CompFunction.h @@ -126,7 +126,7 @@ template class CompFunction { CompFunction paramCopy(bool alloc = false) const; ComplexDouble integrate() const; - ComplexDouble integrate(int dim, bool greater) const; + ComplexDouble integrate(int dim, bool largerSide) const; double norm() const; double getSquareNorm() const; void alloc(int nalloc = 1, bool zero = true); From 8bf70a80f4f4cfd30c92f20dc5e21a5e872e3ce9 Mon Sep 17 00:00:00 2001 From: Niklas Goellmann Date: Tue, 2 Dec 2025 11:00:49 +0100 Subject: [PATCH 3/4] also works if rootBox is not centered around 0,0,0 --- src/trees/FunctionTree.cpp | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/trees/FunctionTree.cpp b/src/trees/FunctionTree.cpp index 7751f2eb..f9a3d8a5 100644 --- a/src/trees/FunctionTree.cpp +++ b/src/trees/FunctionTree.cpp @@ -439,14 +439,15 @@ template T FunctionTree::integrate() const { * @brief Integrate over half of the space * @param dim Dimension along which to split * @param largerSide If true, integrate over the larger side (x>0 if dim=0) - * @return Integral of the function over parts of the computational domain */ + * @return Integral of the function over parts of the computational domain + */ template T FunctionTree::integrate(int dim, bool largerSide) const { T result = 0.0; for (int i = 0; i < this->rootBox.size(); i++) { const FunctionNode &fNode = getRootFuncNode(i); - if (fNode.getUpperBounds()[dim] <= 0.0 && not largerSide) + if (fNode.getLowerBounds()[dim] == this->rootBox.getLowerBound(dim) && not largerSide) result += fNode.integrate(); - if (fNode.getLowerBounds()[dim] >= 0.0 && largerSide) + if (fNode.getUpperBounds()[dim] == this->rootBox.getUpperBound(dim) && largerSide) result += fNode.integrate(); } From d6ee0a826fd836ff4100f76e3feee036d7331af0 Mon Sep 17 00:00:00 2001 From: Niklas Goellmann Date: Wed, 10 Dec 2025 10:27:30 +0100 Subject: [PATCH 4/4] moved documentation; added small unit test --- src/trees/FunctionTree.cpp | 6 ------ src/trees/FunctionTree.h | 8 ++++++++ src/utils/CompFunction.cpp | 6 ------ src/utils/CompFunction.h | 8 ++++++++ tests/trees/function_tree.cpp | 1 + 5 files changed, 17 insertions(+), 12 deletions(-) diff --git a/src/trees/FunctionTree.cpp b/src/trees/FunctionTree.cpp index f9a3d8a5..987cf416 100644 --- a/src/trees/FunctionTree.cpp +++ b/src/trees/FunctionTree.cpp @@ -435,12 +435,6 @@ template T FunctionTree::integrate() const { return jacobian * result; } -/** - * @brief Integrate over half of the space - * @param dim Dimension along which to split - * @param largerSide If true, integrate over the larger side (x>0 if dim=0) - * @return Integral of the function over parts of the computational domain - */ template T FunctionTree::integrate(int dim, bool largerSide) const { T result = 0.0; for (int i = 0; i < this->rootBox.size(); i++) { diff --git a/src/trees/FunctionTree.h b/src/trees/FunctionTree.h index 6f8c9770..ffc7cd45 100644 --- a/src/trees/FunctionTree.h +++ b/src/trees/FunctionTree.h @@ -62,7 +62,15 @@ template class FunctionTree final : public MWTree, pub ~FunctionTree() override; T integrate() const; + + /** + * @brief Integrate over half of the space + * @param dim Dimension along which to split + * @param largerSide If true, integrate over the larger side (x>0 if dim=0) + * @return Integral of the function over parts of the computational domain + */ T integrate(int dim, bool largerSide) const; + double integrateEndNodes(RepresentableFunction_M &f); T evalf_precise(const Coord &r); T evalf(const Coord &r) const override; diff --git a/src/utils/CompFunction.cpp b/src/utils/CompFunction.cpp index eb2dcba2..ab374608 100644 --- a/src/utils/CompFunction.cpp +++ b/src/utils/CompFunction.cpp @@ -187,12 +187,6 @@ template ComplexDouble CompFunction::integrate() const { return integral; } -/** - * @brief Integrate over half of the space - * @param dim Dimension along which to split - * @param largerSide If true, integrate over the larger side (x>0 if dim=0) - * @return ComplexDouble Resulting integral - */ template ComplexDouble CompFunction::integrate(int dim, bool largerSide) const { ComplexDouble integral; if (isreal()) diff --git a/src/utils/CompFunction.h b/src/utils/CompFunction.h index 70d6ade4..3979866a 100644 --- a/src/utils/CompFunction.h +++ b/src/utils/CompFunction.h @@ -126,7 +126,15 @@ template class CompFunction { CompFunction paramCopy(bool alloc = false) const; ComplexDouble integrate() const; + + /** + * @brief Integrate over half of the space + * @param dim Dimension along which to split + * @param largerSide If true, integrate over the larger side (x>0 if dim=0) + * @return Resulting integral + */ ComplexDouble integrate(int dim, bool largerSide) const; + double norm() const; double getSquareNorm() const; void alloc(int nalloc = 1, bool zero = true); diff --git a/tests/trees/function_tree.cpp b/tests/trees/function_tree.cpp index 99e5e709..ab657f8f 100644 --- a/tests/trees/function_tree.cpp +++ b/tests/trees/function_tree.cpp @@ -57,6 +57,7 @@ template void testZeroFunction() { } THEN("its squared norm is zero") { REQUIRE(tree.getSquareNorm() == Catch::Approx(0.0)); } THEN("it integrates to zero") { REQUIRE(tree.integrate() == Catch::Approx(0.0)); } + THEN("half of it integrates to zero") { REQUIRE(tree.integrate(0, true) == Catch::Approx(0.0)); } THEN("the dot product with itself is zero") { REQUIRE(dot(tree, tree) == Catch::Approx(0.0)); } } finalize(&mra);