.. _program_listing_file_src_ompl_control_ODESolver.h:

Program Listing for File ODESolver.h
====================================

|exhale_lsh| :ref:`Return to documentation for file <file_src_ompl_control_ODESolver.h>` (``src/ompl/control/ODESolver.h``)

.. |exhale_lsh| unicode:: U+021B0 .. UPWARDS ARROW WITH TIP LEFTWARDS

.. code-block:: cpp

   /*********************************************************************
   * Software License Agreement (BSD License)
   *
   *  Copyright (c) 2011, Rice University
   *  All rights reserved.
   *
   *  Redistribution and use in source and binary forms, with or without
   *  modification, are permitted provided that the following conditions
   *  are met:
   *
   *   * Redistributions of source code must retain the above copyright
   *     notice, this list of conditions and the following disclaimer.
   *   * Redistributions in binary form must reproduce the above
   *     copyright notice, this list of conditions and the following
   *     disclaimer in the documentation and/or other materials provided
   *     with the distribution.
   *   * Neither the name of the Rice University nor the names of its
   *     contributors may be used to endorse or promote products derived
   *     from this software without specific prior written permission.
   *
   *  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
   *  "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
   *  LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
   *  FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
   *  COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
   *  INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
   *  BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
   *  LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
   *  CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
   *  LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
   *  ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
   *  POSSIBILITY OF SUCH DAMAGE.
   *********************************************************************/
   
   /* Author: Ryan Luna */
   
   #ifndef OMPL_CONTROL_ODESOLVER_
   #define OMPL_CONTROL_ODESOLVER_
   
   #include "ompl/control/Control.h"
   #include "ompl/control/SpaceInformation.h"
   #include "ompl/control/StatePropagator.h"
   #include "ompl/util/Console.h"
   #include "ompl/util/ClassForward.h"
   
   #include <boost/numeric/odeint/stepper/runge_kutta4.hpp>
   #include <boost/numeric/odeint/stepper/runge_kutta_cash_karp54.hpp>
   #include <boost/numeric/odeint/stepper/controlled_runge_kutta.hpp>
   #include <boost/numeric/odeint/stepper/generation/make_controlled.hpp>
   #include <boost/numeric/odeint/integrate/integrate_const.hpp>
   #include <boost/numeric/odeint/integrate/integrate_adaptive.hpp>
   namespace odeint = boost::numeric::odeint;
   #include <functional>
   #include <cassert>
   #include <utility>
   #include <vector>
   
   namespace ompl
   {
       namespace control
       {
           OMPL_CLASS_FORWARD(ODESolver);
   
   
           class ODESolver
           {
           public:
               using StateType = std::vector<double>;
   
               using ODE = std::function<void(const StateType &, const Control *, StateType &)>;
   
               using PostPropagationEvent =
                   std::function<void(const base::State *, const Control *, double, base::State *)>;
   
               ODESolver(SpaceInformationPtr si, ODE ode, double intStep)
                 : si_(std::move(si)), ode_(std::move(ode)), intStep_(intStep)
               {
               }
   
               virtual ~ODESolver() = default;
   
               void setODE(const ODE &ode)
               {
                   ode_ = ode;
               }
   
               double getIntegrationStepSize() const
               {
                   return intStep_;
               }
   
               void setIntegrationStepSize(double intStep)
               {
                   intStep_ = intStep;
               }
   
               const SpaceInformationPtr &getSpaceInformation() const
               {
                   return si_;
               }
   
               static StatePropagatorPtr getStatePropagator(ODESolverPtr solver,
                                                            const PostPropagationEvent &postEvent = nullptr)
               {
                   class ODESolverStatePropagator : public StatePropagator
                   {
                   public:
                       ODESolverStatePropagator(const ODESolverPtr& solver, PostPropagationEvent pe)
                         : StatePropagator(solver->si_), solver_(solver), postEvent_(std::move(pe))
                       {
                           if (solver == nullptr)
                               OMPL_ERROR("ODESolverPtr does not reference a valid ODESolver object");
                       }
   
                       void propagate(const base::State *state, const Control *control, double duration,
                                      base::State *result) const override
                       {
                           ODESolver::StateType reals;
                           si_->getStateSpace()->copyToReals(reals, state);
                           solver_->solve(reals, control, duration);
                           si_->getStateSpace()->copyFromReals(result, reals);
   
                           if (postEvent_)
                               postEvent_(state, control, duration, result);
                       }
   
                   protected:
                       ODESolverPtr solver_;
                       ODESolver::PostPropagationEvent postEvent_;
                   };
                   return std::make_shared<ODESolverStatePropagator>(std::move(solver), postEvent);
               }
   
           protected:
               virtual void solve(StateType &state, const Control *control, double duration) const = 0;
   
               const SpaceInformationPtr si_;
   
               ODE ode_;
   
               double intStep_;
   
               // Functor used by the boost::numeric::odeint stepper object
               struct ODEFunctor
               {
                   ODEFunctor(ODE o, const Control *ctrl) : ode(std::move(o)), control(ctrl)
                   {
                   }
   
                   // boost::numeric::odeint will callback to this method during integration to evaluate the system
                   void operator()(const StateType &current, StateType &output, double /*time*/)
                   {
                       ode(current, control, output);
                   }
   
                   ODE ode;
                   const Control *control;
               };
           };
   
           template <class Solver = odeint::runge_kutta4<ODESolver::StateType>>
           class ODEBasicSolver : public ODESolver
           {
           public:
               ODEBasicSolver(const SpaceInformationPtr &si, const ODESolver::ODE &ode, double intStep = 1e-2)
                 : ODESolver(si, ode, intStep)
               {
               }
   
           protected:
               void solve(StateType &state, const Control *control, double duration) const override
               {
                   Solver solver;
                   ODESolver::ODEFunctor odefunc(ode_, control);
                   odeint::integrate_const(solver, odefunc, state, 0.0, duration, intStep_);
               }
           };
   
           template <class Solver = odeint::runge_kutta_cash_karp54<ODESolver::StateType>>
           class ODEErrorSolver : public ODESolver
           {
           public:
               ODEErrorSolver(const SpaceInformationPtr &si, const ODESolver::ODE &ode, double intStep = 1e-2)
                 : ODESolver(si, ode, intStep)
               {
               }
   
               ODESolver::StateType getError()
               {
                   return error_;
               }
   
           protected:
               void solve(StateType &state, const Control *control, double duration) const override
               {
                   ODESolver::ODEFunctor odefunc(ode_, control);
   
                   if (error_.size() != state.size())
                       error_.assign(state.size(), 0.0);
   
                   Solver solver;
                   solver.adjust_size(state);
   
                   double time = 0.0;
                   while (time < duration + std::numeric_limits<float>::epsilon())
                   {
                       solver.do_step(odefunc, state, time, intStep_, error_);
                       time += intStep_;
                   }
               }
   
               mutable ODESolver::StateType error_;
           };
   
           template <class Solver = odeint::runge_kutta_cash_karp54<ODESolver::StateType>>
           class ODEAdaptiveSolver : public ODESolver
           {
           public:
               ODEAdaptiveSolver(const SpaceInformationPtr &si, const ODESolver::ODE &ode, double intStep = 1e-2)
                 : ODESolver(si, ode, intStep), maxError_(1e-6), maxEpsilonError_(1e-7)
               {
               }
   
               double getMaximumError() const
               {
                   return maxError_;
               }
   
               void setMaximumError(double error)
               {
                   maxError_ = error;
               }
   
               double getMaximumEpsilonError() const
               {
                   return maxEpsilonError_;
               }
   
               void setMaximumEpsilonError(double error)
               {
                   maxEpsilonError_ = error;
               }
   
           protected:
               void solve(StateType &state, const Control *control, double duration) const override
               {
                   ODESolver::ODEFunctor odefunc(ode_, control);
                   auto solver = make_controlled(1.0e-6, 1.0e-6, Solver());
                   odeint::integrate_adaptive(solver, odefunc, state, 0.0, duration, intStep_);
               }
   
               double maxError_;
   
               double maxEpsilonError_;
           };
       }
   }
   
   #endif