Commit c4195c3b authored by Dennis Gläser's avatar Dennis Gläser
Browse files

Merge branch 'feature/constraints-violation-info' into 'master'

Feature/constraints violation info

Closes #16 and #12

See merge request tools/frackit!197
parents 3b0ef154 3856cc1b
Pipeline #2728 passed with stages
in 11 minutes and 16 seconds
frackit_add_test(NAME example1
SOURCES example1.cc
LABELS example)
LABELS example
CMD_ARGS 5)
frackit_symlink_or_copy(FILES example1.py)
......@@ -15,7 +15,7 @@
#include <frackit/entitynetwork/networkbuilder.hh>
//! create a network of 3d quadrilaterals
int main()
int main(int argc, char** argv)
{
using namespace Frackit;
......@@ -85,8 +85,8 @@ int main()
// use the status class to define when to stop sampling
SamplingStatus status;
status.setTargetCount(idSet1, 15); // we want 15 entities in set 1
status.setTargetCount(idSet2, 15); // we want 15 entities in set 2
status.setTargetCount(idSet1, argc > 1 ? std::stoi(argv[1]) : 15); // per default, we want 15 entities in set 1
status.setTargetCount(idSet2, argc > 1 ? std::stoi(argv[1]) : 15); // per default, we want 15 entities in set 2
// start sampling into set 1 and keep alternating
bool sampleIntoSet1 = true;
......
frackit_add_test(NAME example2
SOURCES example2.cc
LABELS example)
LABELS example
CMD_ARGS 5)
frackit_symlink_or_copy(FILES example2.py)
......@@ -17,7 +17,7 @@
#include <frackit/entitynetwork/networkbuilder.hh>
//! create a network of 3d quadrilaterals contained in a cylinder
int main()
int main(int argc, char** argv)
{
using namespace Frackit;
......@@ -71,8 +71,8 @@ int main()
const Id idSet2(2);
SamplingStatus status;
status.setTargetCount(idSet1, 12);
status.setTargetCount(idSet2, 12);
status.setTargetCount(idSet1, argc > 1 ? std::stoi(argv[1]) : 12);
status.setTargetCount(idSet2, argc > 1 ? std::stoi(argv[1]) : 12);
bool sampleIntoSet1 = true;
......
......@@ -177,20 +177,20 @@ int main(int argc, char** argv)
// If the set this geometry belongs to is finished, skip the rest
if (status.finished(id))
{ status.increaseRejectedCounter(); continue; }
{ status.increaseRejectedCounter("set finished"); continue; }
// Moreover, we want to avoid small fragments (< 250 m²)
const auto containedArea = computeContainedMagnitude(geom, networkDomain);
if (containedArea < 350.0)
{ status.increaseRejectedCounter(); continue; }
{ status.increaseRejectedCounter("minimum contained area violation"); continue; }
// enforce constraints w.r.t. to the other entities
if (!constraintsMatrix.evaluate(entitySets, geom, id))
{ status.increaseRejectedCounter(); continue; }
if (const auto res = constraintsMatrix.evaluate(entitySets, geom, id); res.violationDetected())
{ status.increaseRejectedCounter(res.violationLabel()); continue; }
// enforce constraints w.r.t. the domain boundaries
if (!constraintsOnDomain.evaluate(domainBoundaryFaces, geom))
{ status.increaseRejectedCounter(); continue; }
if (const auto res = constraintsOnDomain.evaluate(domainBoundaryFaces, geom); res.violationDetected())
{ status.increaseRejectedCounter(res.violationLabel() + " (boundary)"); continue; }
// the geometry is admissible
entitySets.addEntity(geom, id);
......@@ -205,7 +205,9 @@ int main(int argc, char** argv)
const auto density = containedNetworkArea/domainVolume;
std::cout << "\nEntity density of the contained network: " << density << " m²/m³" << std::endl;
// print info on rejection events
status.printRejectionData();
std::cout << std::endl;
//////////////////////////////////////////////////////////////////////////
......
......@@ -143,24 +143,26 @@ while not status.finished():
# If the set this geometry belongs to is finished, skip the rest
if status.finished(id):
sampleIntoSet1 = not sampleIntoSet1
status.increaseRejectedCounter()
status.increaseRejectedCounter("set finished")
continue
# Moreover, we want to avoid small fragments (< 250 m²)
from frackit.geometry import computeContainedMagnitude
containedArea = computeContainedMagnitude(geom, networkDomain);
if containedArea < 350.0:
status.increaseRejectedCounter()
status.increaseRejectedCounter("minimum contained area violation")
continue
# enforce constraints w.r.t. to the other entities
if not constraintsMatrix.evaluate(entitySets, geom, id):
status.increaseRejectedCounter()
eval = constraintsMatrix.evaluate(entitySets, geom, id)
if eval.violationDetected():
status.increaseRejectedCounter(eval.violationLabel())
continue
# enforce constraints w.r.t. the domain boundaries
if not constraintsOnDomain.evaluate(domainBoundaryFaces, geom):
status.increaseRejectedCounter()
eval = constraintsOnDomain.evaluate(domainBoundaryFaces, geom)
if eval.violationDetected():
status.increaseRejectedCounter(eval.violationLabel() + " (boundary)")
continue
# the geometry is admissible
......@@ -178,6 +180,9 @@ while not status.finished():
density = containedNetworkArea/domainVolume;
print("\nEntity density of the contained network: {:f} m²/m³\n".format(density))
# print info on rejection events
status.printRejectionData();
print("")
##########################################################################
## 5. The entities of the network have been created. We can now ##
......
......@@ -47,6 +47,65 @@
namespace Frackit {
/*!
* \brief Structure to store the result
* of a constraints evaluation.
*/
struct ConstraintsEvaluation
{
enum class Violation
{
none,
distance,
intersectionAngle,
intersectionMagnitude,
intersectionDimension,
intersectionDistance
};
//! constructor
ConstraintsEvaluation(Violation v)
: violation_(v)
{}
//! return the type of violation
Violation violation() const
{ return violation_; }
//! return the violation label
std::string violationLabel() const
{ return violationLabel(violation_); }
//! this evaluates to true if no violation was detected
explicit operator bool() const
{ return violation_ == Violation::none; }
//! returns true if a violation was detected
bool violationDetected() const
{ return violation_ != Violation::none; }
//! return a label for the given violation
static std::string violationLabel(Violation v)
{
if (v == Violation::none) return "no violation";
else if (v == Violation::distance)
return "minimum distance violation";
else if (v == Violation::intersectionAngle)
return "minimum intersection angle violation";
else if (v == Violation::intersectionMagnitude)
return "minimum intersection magnitude violation";
else if (v == Violation::intersectionDimension)
return "maximum intersection dimensionality violation";
else if (v == Violation::intersectionDistance)
return "minimum intersection distance violation";
else
throw std::runtime_error("Unsupported violation flag");
}
private:
Violation violation_;
};
/*!
* \brief Forward declaration of the constraints class.
*/
......@@ -87,6 +146,9 @@ EntityNetworkConstraints<ST> makeDefaultConstraints()
template<class ST, class AE>
class EntityNetworkConstraints
{
using Result = ConstraintsEvaluation;
using Violation = typename Result::Violation;
public:
//! Export type used for constraints
using Scalar = ST;
......@@ -94,6 +156,9 @@ public:
//! Exort the engine used for angle computations
using AngleComputationEngine = AE;
//! Export the result of an evaluation
using ResultType = Result;
/*!
* \brief Default constructor.
* \note This is only available if the engines are default constructible
......@@ -183,10 +248,12 @@ public:
* \brief Check if a pair of geometries fulfills the constraints
* \param geo1 The first geometry
* \param geo2 The second geometry
* \returns True if all constraints are fulfilled, false otherwise
* \return ConstraintsEvaluation type, which emulates true if all
* constraints are fulfilled, false otherwise. Moreover,
* it stores the information which violation was detected.
*/
template<class Geo1, class Geo2>
bool evaluate(const Geo1& geo1, const Geo2& geo2) const
ResultType evaluate(const Geo1& geo1, const Geo2& geo2) const
{
// if the provided geometries are pointers (we support shared_ptr here)
// to the abstract base class, we have to first cast them into the derived types
......@@ -224,14 +291,18 @@ public:
* \tparam Geo2 The geometry type of the entity to be checked
* \param entitySet An entity set
* \param entity The geometry of the entity to be checked
* \returns True if all constraints are fulfilled, false otherwise
* \return ConstraintsEvaluation type, which emulates true if all
* constraints are fulfilled, false otherwise. Moreover,
* it stores the information which violation was detected.
*/
template<class Geo1, class Geo2>
bool evaluate(const std::vector<Geo1>& entitySet, const Geo2& entity) const
ResultType evaluate(const std::vector<Geo1>& entitySet, const Geo2& entity) const
{
return std::all_of(entitySet.begin(),
entitySet.end(),
[&] (const auto& e) { return evaluate(e, entity); });
for (const auto& e : entitySet)
if (const auto res = evaluate(e, entity); !res)
return res;
return {Violation::none};
}
/*!
......@@ -241,10 +312,12 @@ public:
* \tparam Geo2 The geometry type of the entities in the set
* \param entity The geometry of the entity to be checked
* \param entitySet An entity set
* \returns True if all constraints are fulfilled, false otherwise
* \return ConstraintsEvaluation type, which emulates true if all
* constraints are fulfilled, false otherwise. Moreover,
* it stores the information which violation was detected.
*/
template<class Geo1, class Geo2>
bool evaluate(const Geo1& entity, const std::vector<Geo2>& entitySet) const
ResultType evaluate(const Geo1& entity, const std::vector<Geo2>& entitySet) const
{ return evaluate(entitySet, entity); }
/*!
......@@ -254,14 +327,18 @@ public:
* \tparam Geo2 The geometry type of the entities in the second set
* \param entitySet1 The first entity set
* \param entitySet2 The second entity set
* \returns True if all constraints are fulfilled, false otherwise
* \return ConstraintsEvaluation type, which emulates true if all
* constraints are fulfilled, false otherwise. Moreover,
* it stores the information which violation was detected.
*/
template<class Geo1, class Geo2>
bool evaluate(const std::vector<Geo1>& entitySet1, const std::vector<Geo2>& entitySet2) const
ResultType evaluate(const std::vector<Geo1>& entitySet1, const std::vector<Geo2>& entitySet2) const
{
return std::all_of(entitySet1.begin(),
entitySet1.end(),
[&] (const auto& entity1) { return evaluate(entity1, entitySet2); });
for (const auto& e1 : entitySet1)
if (const auto res = evaluate(e1, entitySet2); !res)
return res;
return {Violation::none};
}
private:
......@@ -272,7 +349,7 @@ private:
* \returns True if all constraints are fulfilled, false otherwise
*/
template<class Geo1, class Geo2>
bool evaluate_(const Geo1& geo1, const Geo2& geo2) const
ResultType evaluate_(const Geo1& geo1, const Geo2& geo2) const
{
static constexpr int dim = DimensionalityTraits<Geo1>::geometryDimension();
static_assert(dim == DimensionalityTraits<Geo2>::geometryDimension(),
......@@ -280,40 +357,44 @@ private:
const bool checkIs = useMinIsMagnitude_ || useMinIsAngle_ || useMinIsDistance_;
if (!checkIs && !useMinDistance_)
return true;
return {Violation::none};
const auto isection = !useIntersectionEps_ ? intersect(geo1, geo2)
: intersect(geo1, geo2, intersectionEps_);
if ( !isEmptyIntersection(isection) )
{
// check if dimensionality constraint is violated
// check if dimensionality constraint is violated
if (!allowEquiDimIS_ && !ConstraintImpl::isAdmissibleDimension(isection, dim-1))
return false;
return {Violation::intersectionDimension};
// magnitude constraint
if ( useMinIsMagnitude_ && !ConstraintImpl::isAdmissibleMagnitude(isection, minIsMagnitude_) )
return false;
return {Violation::intersectionMagnitude};
// angle constraint
if ( useMinIsAngle_ && angleEngine_(geo1, geo2, isection) < minIsAngle_ )
return false;
return {Violation::intersectionAngle};
// constraint on distance of intersection to geometry boundaries
if ( !useMinIsDistance_ )
return true;
return {Violation::none};
using namespace ConstraintImpl;
if (!isAdmissibleDistanceToBoundary(isection, geo1, minIsDistance_))
return false;
return isAdmissibleDistanceToBoundary(isection, geo2, minIsDistance_);
return {Violation::intersectionDistance};
if (!isAdmissibleDistanceToBoundary(isection, geo2, minIsDistance_))
return {Violation::intersectionDistance};
}
// no intersection - check distance constraint
else if (useMinDistance_)
return computeDistance(geo1, geo2) >= minDistance_;
{
if (computeDistance(geo1, geo2) < minDistance_)
return {Violation::distance};
}
return true;
return {Violation::none};
}
AngleComputationEngine angleEngine_; //! Algorithms to compute angles between intersecting geometries
......
......@@ -57,6 +57,9 @@ public:
//! export underlying constraints type
using Constraints = C;
//! export evaluation result type
using ResultType = typename C::ResultType;
/*!
* \brief Add a constraint that should be fulfilled
* between the two entity sets defined in the
......@@ -88,10 +91,12 @@ public:
* \param id The id of the entity set to which the entity belongs.
*/
template<class EntitySets, class Entity>
bool evaluate(const EntitySets& entitySets,
const Entity& entity,
const Id& id)
ResultType evaluate(const EntitySets& entitySets,
const Entity& entity,
const Id& id) const
{
using Violation = typename ResultType::Violation;
auto it = constraintsIndexMap_.find(id.get());
if (it == constraintsIndexMap_.end())
throw std::runtime_error("evaluate: no constraints set for given id");
......@@ -102,19 +107,19 @@ public:
const auto& constraints = constraints_[idxPair.second];
// lambda for the evaluation of the constraints
bool isAdmissible;
ResultType res{Violation::none};
auto evalConstraintOnSet = [&] (const auto& entitySet) -> void
{ isAdmissible = constraints.evaluate(entitySet, entity); };
{ res = constraints.evaluate(entitySet, entity); };
// applyOnSet() returns false if lambda was not applied
// This is the case e.g. when the set is empty. Here,
// return false only if applyOnSet() was successful.
if (entitySets.applyOnSet(otherId, evalConstraintOnSet))
if (!isAdmissible)
return false;
if (!res)
return res;
}
return true;
return {Violation::none};
}
private:
......
......@@ -27,6 +27,7 @@
#include <frackit/geometry/disk.hh>
#include <frackit/geometry/quadrilateral.hh>
#include <frackit/geometry/cylindersurface.hh>
#include <frackit/entitynetwork/constraints.hh>
#include <frackit/python/geometry/brepwrapper.hh>
......@@ -43,11 +44,12 @@ namespace Detail {
: public EntityNetworkConstraints<ST>
{
using ParentType = EntityNetworkConstraints<ST>;
using Result = typename ParentType::ResultType;
public:
// evaluate the constraint between two geometries
template<class Geo1, class Geo2>
bool evaluateBinary(const Geo1& geo1, const Geo2& geo2) const
Result evaluateBinary(const Geo1& geo1, const Geo2& geo2) const
{ return ParentType::evaluate(getUnwrappedShape(geo1), getUnwrappedShape(geo2)); }
};
......@@ -93,27 +95,37 @@ namespace Detail {
using Disk = Frackit::Disk<ctype>;
using Quad_3 = Frackit::Quadrilateral<ctype, 3>;
using Poly_3 = Frackit::Polygon<ctype, 3>;
using CylSurface = Frackit::CylinderSurface<ctype>;
using Face = FaceWrapper;
registerBinaryEvaluator<Disk, Disk>(cls);
registerBinaryEvaluator<Disk, Quad_3>(cls);
registerBinaryEvaluator<Disk, Poly_3>(cls);
registerBinaryEvaluator<Disk, Face>(cls);
registerBinaryEvaluator<Disk, CylSurface>(cls);
registerBinaryEvaluator<Quad_3, Quad_3>(cls);
registerBinaryEvaluator<Quad_3, Poly_3>(cls);
registerBinaryEvaluator<Quad_3, Disk>(cls);
registerBinaryEvaluator<Quad_3, Face>(cls);
registerBinaryEvaluator<Quad_3, CylSurface>(cls);
registerBinaryEvaluator<Poly_3, Poly_3>(cls);
registerBinaryEvaluator<Poly_3, Quad_3>(cls);
registerBinaryEvaluator<Poly_3, Disk>(cls);
registerBinaryEvaluator<Poly_3, Face>(cls);
registerBinaryEvaluator<Poly_3, CylSurface>(cls);
registerBinaryEvaluator<Face, Face>(cls);
registerBinaryEvaluator<Face, Disk>(cls);
registerBinaryEvaluator<Face, Quad_3>(cls);
registerBinaryEvaluator<Face, Poly_3>(cls);
registerBinaryEvaluator<Face, CylSurface>(cls);
registerBinaryEvaluator<CylSurface, Face>(cls);
registerBinaryEvaluator<CylSurface, Disk>(cls);
registerBinaryEvaluator<CylSurface, Quad_3>(cls);
registerBinaryEvaluator<CylSurface, Poly_3>(cls);
}
} // end namespace Detail
......@@ -121,6 +133,29 @@ namespace Detail {
template<class ctype>
void registerConstraints(py::module& module)
{
// register evaluation result type
using Result = ConstraintsEvaluation;
py::class_<Result> evalResult(module, "ConstraintsEvaluation");
// register constraints violation
py::enum_<Result::Violation>(evalResult, "Violation")
.value("none", Result::Violation::none)
.value("distance", Result::Violation::distance)
.value("intersectionAngle", Result::Violation::intersectionAngle)
.value("intersectionMagnitude", Result::Violation::intersectionMagnitude)
.value("intersectionDimension", Result::Violation::intersectionDimension)
.value("intersectionDistance", Result::Violation::intersectionDistance);
// define members of ConstraintsEvaluation
evalResult.def(py::init<Result::Violation>());
evalResult.def("__bool__", &Result::operator bool, "returns true if no violation was detected");
evalResult.def("violation", &Result::violation, "return the type of detected violation");
evalResult.def("violationDetected", &Result::violationDetected, "return true if a violation was detected");
evalResult.def("violationLabel", py::overload_cast<>(&Result::violationLabel, py::const_),
"return a label for the type of detected violation");
evalResult.def("violationLabel", py::overload_cast<Result::Violation>(&Result::violationLabel),
"return a label for the given violation");
using Constraints = Detail::EntityNetworkConstraintsWrapper<ctype>;
py::class_<Constraints> cls(module, "_EntityNetworkConstraints");
cls.def(py::init());
......
......@@ -20,6 +20,7 @@
#define FRACKIT_PYTHON_SAMPLING_STATUS_HH
#include <pybind11/pybind11.h>
#include <pybind11/stl.h>
#include <frackit/common/id.hh>
#include <frackit/sampling/status.hh>
......@@ -54,21 +55,38 @@ void registerSamplingStatus(py::module& module)
&SamplingStatus::reset,
"Reset everything");
cls.def("increaseRejectedCounter",
&SamplingStatus::increaseRejectedCounter,
py::overload_cast<>(&SamplingStatus::increaseRejectedCounter),
"Register that a sample has been rejected");
cls.def("increaseRejectedCounter",
py::overload_cast<const std::string&>(&SamplingStatus::increaseRejectedCounter),
"Register that a sample has been rejected passing a specific rejection label");
cls.def("getCount",
py::overload_cast<>(&SamplingStatus::getCount),
py::overload_cast<>(&SamplingStatus::getCount, py::const_),
"Returns the overall number of entities");
cls.def("getCount",
py::overload_cast<const Id&>(&SamplingStatus::getCount),
py::overload_cast<const Id&>(&SamplingStatus::getCount, py::const_),
"Returns the number of entities for the given id");
cls.def("getRejectedCount",
py::overload_cast<>(&SamplingStatus::getRejectedCount, py::const_),
"Returns the overall number of rejected entities");
cls.def("getRejectedCount",
py::overload_cast<const std::string&>(&SamplingStatus::getRejectedCount, py::const_),
"Returns the number of rejected entities for the given rejection label");
cls.def("getRejectionData",
&SamplingStatus::getRejectionData,
"Returns the registered rejection event data");
using namespace py::literals;
cls.def("print",
&SamplingStatus::print,
"forceHeaderPrint"_a=false,
"Print the current status");
cls.def("printRejectionData",
&SamplingStatus::printRejectionData,
"Print information on registered rejection events");
}
} // end namespace Frackit::Python
......
......@@ -25,6 +25,7 @@
#ifndef FRACKIT_SAMPLING_STATUS_HH
#define FRACKIT_SAMPLING_STATUS_HH
#include <string>
#include <iomanip>
#include <iostream>
#include <algorithm>
......@@ -67,9 +68,11 @@ public:
void reset()
{
headerPrinted_ = false;
rejectedCount_ = 0;
unspecifiedRejectedCount_ = 0;
count_.clear();
targetCount_.clear();
rejectedCount_.clear();
}
/*!
......@@ -77,7 +80,9 @@ public:
*/
void resetCounters()
{
for (auto& count : count_) count.second = 0;
std::for_each(count_.begin(),
count_.end(),
[] (auto& idCountPair) { idCountPair.second = 0.0; });
}
/*!
......@@ -127,19 +132,31 @@ public:
*/
void increaseRejectedCounter()
{
rejectedCount_++;
unspecifiedRejectedCount_++;
}
/*!
* \brief Increase counter of rejected samples
* with a specific rejection label.
*/
void increaseRejectedCounter(const std::string& label)
{
rejectedCount_[label]++;
}