Skip to content

Solvers: integrator step size #357 #818

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

Open
pietercollins wants to merge 2 commits into
base: master
Choose a base branch
from
Open

Solvers: integrator step size #357 #818

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
e3d775e
Fixed bugs in Bounder due to incorrect widening constants. Now wideni…
pietercollins Mar 20, 2025
dfda30b
Exported Bounders to Python interface. Added Box add and mul, and exp…
pietercollins Mar 20, 2025
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
10 changes: 8 additions & 2 deletions python/source/geometry_submodule.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -652,7 +652,8 @@ Void export_intervals(pybind11::module& module) {
module.def("cast_singleton", (FloatDPBounds(*)(Interval<FloatDPUpperBound> const&)) &cast_singleton);
module.def("cast_singleton", (FloatMPBounds(*)(Interval<FloatMPUpperBound> const&)) &cast_singleton);

module.def("image", (UpperIntervalType(*)(UpperIntervalType const&, ValidatedScalarUnivariateFunction const&)) &_image_);
module.def("widen", (FloatDPUpperInterval(*)(FloatDPUpperInterval const&, FloatDPUpperBound eps)) &widen);
module.def("image", (FloatDPUpperInterval(*)(FloatDPUpperInterval const&, ValidatedScalarUnivariateFunction const&)) &_image_);

template_<Interval> interval_template(module);
interval_template.instantiate<Dyadic>();
Expand Down Expand Up @@ -768,6 +769,12 @@ template<class BX> Void export_box(pybind11::module& module, std::string name=py
module.def("intersection", (BX(*)(const BX&,const BX&)) &intersection);
module.def("split", (Pair<BX,BX>(*)(BX const&)) &split);

if constexpr (Same<BX,FloatDPUpperBox>) {
box_class.def("__add__", &__add__<BX,BX>);
box_class.def("__radd__", &__radd__<BX,BX>);
box_class.def("__rmul__", &__rmul__<BX,IVL>);
}

if constexpr (Same<BX,BoxDomainType>) {
module.attr("BoxDomainType")=box_class;
} else if constexpr (Same<BX,BoxValidatedRangeType>) {
Expand Down Expand Up @@ -833,7 +840,6 @@ Void export_boxes(pybind11::module& module) {
pybind11::implicitly_convertible<FloatDPUpperBox,FloatDPApproximateBox>();
pybind11::implicitly_convertible<FloatDPLowerBox,FloatDPApproximateBox>();

module.def("widen", (FloatDPUpperBox(*)(FloatDPExactBox const&, FloatDP eps)) &widen);
module.def("image", (FloatDPUpperBox(*)(FloatDPUpperBox const&, ValidatedVectorMultivariateFunction const&)) &_image_);
module.def("image", (FloatDPUpperInterval(*)(FloatDPUpperBox const&, ValidatedScalarMultivariateFunction const&)) &_image_);
module.def("image", (FloatDPUpperBox(*)(FloatDPUpperInterval const&, ValidatedVectorUnivariateFunction const&)) &_image_);
Expand Down
28 changes: 28 additions & 0 deletions python/source/solver_submodule.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -78,6 +78,19 @@ class SolverWrapper
};


class BounderWrapper
: public pybind11::wrapper<BounderInterface>
{
public:
BounderInterface* clone() const {
return this->get_override("clone")(); }
Pair<StepSizeType,BoxDomainType> compute(const ValidatedVectorMultivariateFunction& vf, const ExactBoxType& D, StepSizeType h) const {
return this->get_override("compute")(vf,D,h); }
Void _write(OutputStream& os) const {
this->get_override("_write")(os); }
};


class IntegratorWrapper
: public pybind11::wrapper<IntegratorInterface>
{
Expand Down Expand Up @@ -120,6 +133,19 @@ Void export_solvers(pybind11::module& module)



Void export_bounders(pybind11::module& module)
{
pybind11::class_<Suggestion<StepSizeType>> suggested_step_size_class(module,"SuggestedStepSize");
module.def("suggest", [](StepSizeType const& h){return static_cast<Suggestion<StepSizeType>>(suggest(h));});

pybind11::class_<BounderInterface,BounderWrapper> bounder_interface_class(module,"BounderInterface");
bounder_interface_class.def("compute",(Pair<StepSizeType,UpperBoxType>(BounderInterface::*)(const ValidatedVectorMultivariateFunction&, const BoxDomainType&, const Suggestion<StepSizeType>&)const) &BounderInterface::compute);

pybind11::class_<EulerBounder,BounderInterface> euler_bounder_class(module,"EulerBounder");
euler_bounder_class.def(pybind11::init<>());
euler_bounder_class.def("compute",(Pair<StepSizeType,UpperBoxType>(EulerBounder::*)(const ValidatedVectorMultivariateFunction&, const BoxDomainType&, const Suggestion<StepSizeType>&)const) &EulerBounder::compute);
}

Void export_integrators(pybind11::module& module)
{
pybind11::class_<FlowStepModelType,pybind11::bases<ValidatedVectorMultivariateFunctionPatch>> flow_step_model_class(module,"FlowStepModelType");
Expand Down Expand Up @@ -185,6 +211,8 @@ Void export_integrators(pybind11::module& module)
Void solver_submodule(pybind11::module& module)
{
export_solvers(module);

export_bounders(module);
export_integrators(module);
}

Expand Down
7 changes: 7 additions & 0 deletions source/geometry/box.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -159,6 +159,13 @@ template<class I> List<typename Box<I>::VertexType> Box<I>::vertices() const {

template List<typename FloatDPExactBox::VertexType> Box<FloatDPExactInterval>::vertices() const;

FloatDPUpperBox operator+(const FloatDPUpperBox& bx1, FloatDPUpperBox bx2) {
return FloatDPUpperBox(bx1.dimension(),[&bx1,bx2](SizeType i){return bx1[i]+bx2[i];});
}
FloatDPUpperBox operator*(const FloatDPUpperInterval& ivl1, FloatDPUpperBox bx2) {
return FloatDPUpperBox(bx2.dimension(),[&ivl1,bx2](SizeType i){return ivl1*bx2[i];});
}

/*
FloatDPLowerBox under_approximation(const RealBox& rbx) {
FloatDPLowerBox bx(rbx.size());
Expand Down
4 changes: 4 additions & 0 deletions source/geometry/box.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -538,6 +538,10 @@ template<class I, class X> inline Vector<decltype(declval<I>()-declval<X>())> op
return cast_vector(bx1)-v2; }



FloatDPUpperBox operator+(const FloatDPUpperBox& bx1, FloatDPUpperBox bx2);
FloatDPUpperBox operator*(const FloatDPUpperInterval& ivl1, FloatDPUpperBox bx2);

inline FloatDPApproximateBox widen(const FloatDPApproximateBox& bx, FloatDPApproximation e) {
FloatDPApproximation eps(e);
FloatDPApproximateBox r(bx.dimension());
Expand Down
12 changes: 6 additions & 6 deletions source/solvers/bounder.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -56,9 +56,9 @@ Pair<StepSizeType,UpperBoxType> EulerBounder::compute(ValidatedVectorMultivariat
Pair<StepSizeType,UpperBoxType> EulerBounder::_compute(ValidatedVectorMultivariateFunction const& f, BoxDomainType const& D, StepSizeType const& t, BoxDomainType const& A, Suggestion<StepSizeType> const& hsug) const {
CONCLOG_SCOPE_CREATE;
const PositiveFloatDP BOX_RADIUS_WIDENING=cast_positive(0.25_exact);
const PositiveFloatDP NO_WIDENING=cast_positive(1.0_exact);
const PositiveFloatDP INITIAL_STARTING_WIDENING=cast_positive(2.0_exact);
const PositiveFloatDP INITIAL_REFINING_WIDENING=cast_positive(1.125_exact);
const PositiveFloatDP NO_WIDENING=cast_positive(0.0_exact);
const PositiveFloatDP INITIAL_STARTING_WIDENING=cast_positive(1.0_exact);
const PositiveFloatDP INITIAL_REFINING_WIDENING=cast_positive(0.125_exact);
const StepSizeType MINIMUM_STEP_SIZE(1,20u);
const CounterType EXPANSION_STEPS=4;
const CounterType REFINEMENT_STEPS=4;
Expand Down Expand Up @@ -118,7 +118,7 @@ Pair<StepSizeType,UpperBoxType> EulerBounder::_compute(ValidatedVectorMultivaria
}

UpperBoxType EulerBounder::_refinement(ValidatedVectorMultivariateFunction const& f, BoxDomainType const& D, IntervalDomainType const& T, BoxDomainType const& A, UpperBoxType const& B) const {
const PositiveFloatDP NO_WIDENING=cast_positive(1.0_exact);
const PositiveFloatDP NO_WIDENING=cast_positive(0.0_exact);
return _formula(f,D,T,A,B, NO_WIDENING,NO_WIDENING);
}

Expand All @@ -130,9 +130,9 @@ UpperBoxType EulerBounder::_formula(ValidatedVectorMultivariateFunction const& f
const bool is_autonomous = (f.argument_size() == D.dimension()+A.dimension());
UpperBoxType dom = is_autonomous ? product(B,rA) : product(B,rT,rA);

UpperBoxType wD = (INITIAL_BOX_WIDENING!=1) ? D : D + INITIAL_BOX_WIDENING*(D-D.midpoint());
UpperBoxType wD = (INITIAL_BOX_WIDENING==0) ? D : D + INITIAL_BOX_WIDENING*(D-D.midpoint());

return wD+(VECTOR_WIDENING*rH)*cast_vector(apply(f,dom));
return wD+(rH+VECTOR_WIDENING*rH)*cast_vector(apply(f,dom));
}

} // namespace Ariadne;
Loading