odetimet.h

Go to the documentation of this file.

/////////////////////////////////////////////////////////////////////////////
///
/// @file odetimet.h
///
/// A generic method for computing a time-t map for ODEs.
///
/// This file defines a generic class which uses the integration routines
/// of the CAPD library in order to compute a time-t map for ODEs.
///
/// @author Pawel Pilarczyk
///
/////////////////////////////////////////////////////////////////////////////
 
// Copyright (C) 1997-2014 by Pawel Pilarczyk.
//
// This file is part of my research software package. This is free software:
// you can redistribute it and/or modify it under the terms of the GNU
// General Public License as published by the Free Software Foundation,
// either version 3 of the License, or (at your option) any later version.
//
// This software is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this software; see the file "license.txt". If not,
// please, see <https://www.gnu.org/licenses/>.
 
// Started on October 18, 2008. Last revision: September 4, 2014.
 
 
#ifndef _CMGRAPHS_ODETIMET_H_
#define _CMGRAPHS_ODETIMET_H_
 
 
// include some standard C++ header files
#include <string>
#include <vector>
#include <iostream>
#include <algorithm>
#include <memory>
 
// include all the necessary header files from the CHomP library
#include "chomp/system/config.h"
#include "chomp/system/textfile.h"
#include "chomp/multiwork/mwdata.h"
 
// include local header files
#include "typeintv.h"
#include "typerect.h"
#include "maptight.h"
#include "datavector.h"
#include "odeintegr.h"
 
// include all the necessary header files from the CAPD library
#ifdef _CAPD29_
#include "capd/vectalg/lib.h"
#include "capd/dynset/lib.h"
#include "capd/dynsys/lib.h"
#include "capd/dynsys/Taylor.hpp"
#include "capd/dynsys/TaylorHOE.hpp"
//#include "capd/capdlib.h"
#else
//#include "capd/capdlib.h"
#define CAPD_USER_NAMESPACE capd
#define CAPD_DEFAULT_DIMENSION spaceDim
#include "capd/fdcapdlib.h"
#undef CAPD_USER_NAMESPACE
#undef CAPD_DEFAULT_DIMENSION
#endif
 
 
// --------------------------------------------------
// ----------------- the time-t map -----------------
// --------------------------------------------------
 
/// This class defines a map for the nonlinear density dependent
/// overcompensatory Leslie population model.
/// The interpretation of parameters depends on the dimension.
class MapOdeTimeT: public MapTight
{
public:
        /// The constructor which takes the formula for the ODE
        /// in the format used in the CAPD library.
        /// The names of parameters in the formula must be "pa", "pb", etc.
        /// The order of the Taylor method to use must be given (e.g. 2),
        /// the integration step size, and the number of steps.
        /// The total translation in the flow equals step * n.
        /// Note: The step should be a representable number
        /// (for example, a multiple of a power of 2, like 3.0/256).
        MapOdeTimeT (const std::string &formula, int order,
                double step, int n);
 
        /// The destructor.
        ~MapOdeTimeT ();
 
        /// Sets the parameters to the given intervals defined by their
        /// left and right coordinates.
        /// This function actually calls the original function from the
        /// master class, and also resets the Taylor class iterator.
        void setParam (const double *left, const double *right, int n);
 
        /// Computes a rectangular set that is an outer enclosure
        /// of the image of the given box.
        void compute (const IntervalVectorType &x,
                RectSetType &y, int dim,
                const spcCoord *coord, int subdiv) const;
 
        /// Sets a bounding box for the area of interest.
        /// Every time the provided number of itereates is made,
        /// the iterated set is verified whether it is contained
        /// in the bounding box. If it is completely outside
        /// then it is projected to the boundary of the box
        /// and the iterations are halted.
        void setBoundingBox (const double *left, const double *right,
                int dim, int stepsBetweenChecks);
 
        /// Sets the phase space for an automatic generation
        /// of bounding boxes.
        template <class Array1, class Array2>
        static void setPhaseSpace (const Array1 &offset, const Array2 &width,
                int dim, int subdivDepth, int stepsBetweenChecks = 2);
 
        /// Sets the bounding box to the default one, as defined
        /// when setting the phase space.
        void setDefaultBoundingBox ();
 
        /// Adjusts parameters if necessary after a successful computation
        /// or after an unsuccessful computation of the graph of the map
        /// at a specific subdivision depth.
        /// Returns true iff it is necessary to repeat the computation.
        bool adjust (bool successful, int subdiv);
 
        /// Gets a single line of information on optimization
        /// of the computation of the map, so that it can be used
        /// by maps created for the same and nearby parameter values.
        std::string getOptInfo () const;
 
        /// Uses a prevously saved line of information on optimization
        /// of the computation of the map for the same or nearby
        /// parameter values. The map object may be fed with several lines
        /// and should interprete them to set up its optimization parameters.
        void useOptInfo (const std::string &info);
 
        /// Saves the internal parameters and other data of the map
        /// to the mwData structure.
        void saveInternals (chomp::multiwork::mwData &data) const;
 
        /// Retrieves the internal paramters and other data of the map
        /// from the mwData structure.
        void loadInternals (chomp::multiwork::mwData &data);
 
        /// Returns an approximation of the currently recommended
        /// integration time t for the time-t map
        /// at the given subdivision depth.
        double getTimeStep (int subdiv) const;
 
protected:
        /// A bounding box for the area where the trajectories are computed.
        /// If an iterated set gets outside the given ranges of coordinates
        /// then it is projected to the boundary of this box
        /// and the iterating process is halted.
        /// The array contains left1, right1, left2, right2, etc.
        double *boundingBox;
 
        /// The dimension of the bounding box.
        int boundingBoxDim;
 
private:
        /// The interval matrix type.
        typedef capd::IMatrix IntervalMatrixType;
 
        /// The type of a map to be used for the Taylor method.
        typedef capd::IMap IntervalMapType;
 
#ifdef _CAPD29_
        /// The Taylor method type.
        typedef capd::dynsys::TaylorHOE<IntervalMapType> TaylorType;
//      typedef capd::dynsys::Taylor<IntervalMapType> TaylorType;
//      typedef capd::ITaylor TaylorType;
#else
        /// The type of fixed step control.
//      typedef capd::dynsys::FixedStepControl<IntervalType> StepControl;
 
        /// The Taylor method type.
//      typedef capd::dynsys::TaylorHOE<IntervalMapType,StepControl> TaylorType;
        typedef capd::dynsys::Solver<IntervalMapType> TaylorType;
//      typedef capd::dynsys::SolverHOE<IntervalMapType> TaylorType;
//      typedef capd::ITaylor TaylorType;
#endif
 
        /// The formula for the map in the text format used in CAPD.
        std::string mapFormula;
 
        /// An object for adjusting ODE integration parameters.
        mutable OdeIntegrParam odeIntegrParam;
 
        /// The map object to be used in the Taylor method object.
        IntervalMapType *theMap;
 
        /// The default bounding box to be used if no bounding box is set.
        static std::vector<double> defaultBoundingBox;
 
        /// The number of steps after which the box should be checked
        /// if it is still contained in the bounding box.
        int checkEvery;
 
        /// The default number of steps to check the bounding box.
        static int defaultCheckEvery;
 
        /// The copy constructor should not be used.
        MapOdeTimeT (const MapOdeTimeT &);
 
        /// The assignment operator should not be used.
        MapOdeTimeT &operator = (const MapOdeTimeT &);
 
        /// Verifies if the given point is outside the bounding box
        /// set with the method "setBoundingBox".
        /// Returns true if the bounding box was set, its dimension conforms
        /// to the dimension of the vector, and the point is outside the box
        /// (weak inequality); in this case the point is projected
        /// onto the hyperplane containing the corresponding face
        /// of the bounding box. Otherwise returns false.
        bool checkOutside (double *x, int dim) const;
 
        /// A version of the checkOutside procedure for interval vectors.
        bool checkOutside (IntervalVectorType &v, int dim) const;
 
        /// Verifies if the given interval vector is completely outside
        /// the bounding box set with the method "setBoundingBox".
        /// Returns true if the bounding box was set, its dimension conforms
        /// to the dimension of the vector, and the interval vector
        /// is outside the interior of the box; in this case the vector
        /// is projected onto the hyperplane containing the corresponding
        /// face of the bounding box. Otherwise returns false.
        bool checkOutside (double *left, double *right, int dim) const;
 
        /// Iterates the given Lohner set, using the Taylor method objects.
        void iterateLohnerSet (RectSetType &lohnerSet,
                TaylorType &theTaylorMethod, int steps, int dim) const;
 
        /// The numbers of successful integrations.
        mutable long integrations;
 
}; /* class MapOdeTimeT */
 
// --------------------------------------------------
 
inline void MapOdeTimeT::setDefaultBoundingBox ()
{
        // if there is no default bounding box then do nothing
        if (defaultBoundingBox. empty ())
                return;
 
        int dim = defaultBoundingBox. size () / 2;
        double *left = new double [dim];
        double *right = new double [dim];
        std::vector<double>::const_iterator boxIter =
                defaultBoundingBox. begin ();
        for (int d = 0; d < dim; ++ d)
        {
                left [d] = *(boxIter ++);
                right [d] = *(boxIter ++);
        }
        setBoundingBox (left, right, dim, defaultCheckEvery);
        delete [] right;
        delete [] left;
        return;
} /* MapOdeTimeT::setDefaultBoundingBox */
 
inline MapOdeTimeT::MapOdeTimeT (const std::string &formula,
        int order, double step, int n):
        boundingBox (0), boundingBoxDim (0),
        mapFormula (formula), odeIntegrParam (order, step, n),
        theMap (0), checkEvery (0), integrations (0)
{
        // test the interval arithmetic to make sure it workd properly
        testIntervals ();
 
        // set the bounding box to the default one if any
        setDefaultBoundingBox ();
 
#ifdef _CAPD29_
        // check if the map formula has parameters listed explicitly
        if (mapFormula. find ("par:") != std::string::npos)
        {
                chomp::homology::sout << "WARNING: The phrase "
                        "\"par:\" is not supported in CAPD 2.9.\n";
        }
#else
        // check if the map formula has parameters listed explicitly
        if (mapFormula. find ("par:") == std::string::npos)
        {
                chomp::homology::sout << "WARNING: Couldn't find the phrase "
                        "\"par:\" in the map formula.\n";
        }
#endif
 
        // swich the rounding direction to the neutral one
        resetRounding ();
        return;
} /* MapOdeTimeT::MapOdeTimeT */
 
inline void MapOdeTimeT::setParam (const double *left, const double *right,
        int n)
{
        // set the new values of parameters
        MapType::setParam (left, right, n);
        odeIntegrParam. setParam (left, right, n);
 
        // delete the previously used Map object (if any)
        if (theMap)
        {
                delete theMap;
                theMap = 0;
        }
 
        // create a new map object
        theMap = new IntervalMapType (mapFormula);
#ifdef _CAPD29_
        char paramName [3];
        paramName [0] = 'p';
        paramName [2] = '\0';
#endif
        for (int i = 0; i < n; ++ i)
        {
                IntervalType paramInterval (left [i], right [i]);
#ifdef _CAPD29_
                paramName [1] = static_cast<char> ('a' + i);
                theMap -> setParameter (paramName, paramInterval);
#else
                theMap -> setParameter (i, paramInterval);
#endif
        }
 
        // swich the rounding direction to the neutral one
        resetRounding ();
        return;
} /* MapOdeTimeT::setParam */
 
inline MapOdeTimeT::~MapOdeTimeT ()
{
        using chomp::homology::sbug;
 
        if (theMap)
                delete theMap;
        if (boundingBox)
                delete [] boundingBox;
 
        // show the number of integrations for each Taylor method
        if (integrations)
                sbug << "OdeTimeT integrations: " << integrations << ".\n";
 
        // show the final number of steps at each subdivision depth
/*      sbug << "Numbers of steps (from depth " << initialDepth << " up):";
        for (int subdiv = initialDepth; (subdiv <= finalDepth) &&
                (subdiv < static_cast<int> (stepsAdd. size ())); ++ subdiv)
        {
                sbug << " " << (numberOfSteps + stepsAdd [subdiv]);
        }
        sbug << ".\n";
*/
        // show the step multipliers
/*      sbug << "Step size multipliers (from depth " <<
                initialDepth << " up):";
        for (int subdiv = initialDepth; (subdiv <= finalDepth) &&
                (subdiv < static_cast<int> (stepMult. size ())); ++ subdiv)
        {
                sbug << " " << stepMult [subdiv];
        }
        sbug << ".\n";
*/
        // output C++ code to put in the map in order to define
        // the arrays stepsAdd and stepMult
/*      sbug << "C++ code for the map file to re-do the computation:\n";
        sbug << "\tstepsAddFactor = 0;\n";
        int stepsAddSize = stepsAdd. size ();
        if (stepsAddSize)
        {
                sbug << "\tconst int N = " << stepsAddSize << ";\n";
                for (int i = 0; i < stepsAddSize; ++ i)
                {
                        sbug << (i ? ", " : "\tint add [N] = {") <<
                                stepsAdd [i];
                }
                sbug << "};\n";
                sbug << "\tstepsAdd = std::vector<int> (N);\n";
                sbug << "\tstd::copy (add, add + N, stepsAdd. begin ());\n";
        }
        int stepMultSize = stepMult. size ();
        if (stepMultSize)
        {
                sbug << "\tconst int M = " << stepMultSize << ";\n";
                int precision = std::cout. precision ();
                int hiprec = 15;
                sbug << std::setprecision (hiprec);
                for (int i = 0; i < stepMultSize; ++ i)
                {
                        sbug << (i ? ", " : "\tdouble mult [M] = {") <<
                                stepMult [i];
                }
                sbug << std::setprecision (precision);
                sbug << "};\n";
                sbug << "\tstepMult = std::vector<double> (M);\n";
                sbug << "\tstd::copy (mult, mult + M, stepMult. begin ());\n";
        }
*/
        // swich the rounding direction to the neutral one
        resetRounding ();
        return;
} /* MapOdeTimeT::~MapOdeTimeT */
 
inline void MapOdeTimeT::iterateLohnerSet (RectSetType &lohnerSet,
        TaylorType &theTaylorMethod, int steps, int dim) const
{
        int counter = 0;
        for (int step = 0; step < steps; ++ step)
        {
                // move the set under the flow
                lohnerSet. move (theTaylorMethod);
 
                // check if the set went outside the given region
                if (++ counter >= checkEvery)
                {
                        IntervalVectorType enclosure (lohnerSet);
                        if (checkOutside (enclosure, dim))
                                break;
                        counter = 0;
                }
        }
 
        return;
} /* MapOdeTimeT::iterateLohnerSet */
 
inline void MapOdeTimeT::compute (const IntervalVectorType &x,
        RectSetType &y, int dim, const spcCoord *coord, int subdiv) const
{
        // create the Lohner set
        IntervalVectorType center (dim);
        for (int i = 0; i < dim; ++ i)
        {
                double middle = (x [i]. leftBound () +
                        x [i]. rightBound ()) / 2;
                center [i] = IntervalType (middle, middle);
        }
        RectSetType y0 (center, x - center);
 
        // create a new Taylor method object for the integration
        int order (0);
        double stepSize (0);
        int numberOfSteps (0);
        odeIntegrParam. get (order, stepSize, numberOfSteps, subdiv);
//      StepControl stepCtrl (IntervalType (stepSize, stepSize));
#ifdef _CAPD29_
        TaylorType theTaylorMethod (*theMap, order,
                IntervalType (stepSize, stepSize));
#else
        TaylorType theTaylorMethod (*theMap, order);
        theTaylorMethod. setStep (IntervalType (stepSize, stepSize));
#endif
 
        // try running the integration process
        try
        {
                y = y0;
                iterateLohnerSet (y, theTaylorMethod, numberOfSteps, dim);
                ++ integrations;
 
                // swich the rounding direction to the neutral one
                resetRounding ();
                return;
        }
        catch (const std::exception &e)
        {
                // show detailed information about the problem
                if (false)
                {
                        using chomp::homology::sbug;
                        sbug << "FAILURE:\n" << e. what () << "\n";
                }
                resetRounding ();
        }
 
        // show detailed debug information about the wrong box
        if (false)
        {
                using chomp::homology::sbug;
                sbug << "[Integration failure at subdiv depth " <<
                        subdiv << ".]\n";
                if (coord)
                {
                        sbug << "Trouble source cube: [";
                        for (int i = 0; i < dim; ++ i)
                                sbug << (i ? "," : "") << coord [i];
                        sbug << "].\n";
                        sbug << "Real coords: " << x << ".\n";
                }
                else
                {
                        sbug << "Trouble source box: " << x << ".\n";
                }
        }
 
        // swich the rounding direction to the neutral one
        resetRounding ();
 
        // throw an error message indicating the failure
        throw "Unable to integrate the ODE - the set became too large.";
        return;
} /* MapOdeTimeT::compute */
 
inline void MapOdeTimeT::setBoundingBox (const double *left,
        const double *right, int dim, int stepsBetweenChecks)
{
        using chomp::homology::sbug;
        if (boundingBox)
                delete [] boundingBox;
        boundingBoxDim = dim;
        boundingBox = new double [dim << 1];
        double *boxPtr (boundingBox);
        sbug << "* Setting the bounding box for the ODE solver:\n";
        for (int i = 0; i < dim; ++ i)
        {
                sbug << (i ? " x [" : "* [") << left [i] << "," <<
                        right [i] << "]";
                *(boxPtr ++) = left [i];
                *(boxPtr ++) = right [i];
        }
        sbug << "\n";
        checkEvery = stepsBetweenChecks;
        return;
} /* MapOdeTimeT::setBoundingBox */
 
template <class Array1, class Array2>
inline void MapOdeTimeT::setPhaseSpace (const Array1 &offset,
        const Array2 &width, int dim, int subdivDepth,
        int stepsBetweenChecks)
{
        // do nothing if the provided values seem to be wrong
        if ((dim <= 0) || (subdivDepth <= 0) ||
                (subdivDepth >= static_cast<int> (8 * sizeof (int) - 1)))
        {
                return;
        }
 
        // prepare for setting a new default bounding box
        defaultBoundingBox. clear ();
        int intWidth = 1 << subdivDepth;
 
        // calculate the left and right bounds
        for (int d = 0; d < dim; ++ d)
        {
                double widthFraction = (width [d] / intWidth) * 1.1;
                double left = offset [d] - widthFraction;
                double right = offset [d] + width [d] + widthFraction;
                defaultBoundingBox. push_back (left);
                defaultBoundingBox. push_back (right);
        }
 
        // set the default number of steps between checks
        defaultCheckEvery = stepsBetweenChecks;
 
        return;
} /* MapOdeTimeT::setPhaseSpace */
 
inline bool MapOdeTimeT::checkOutside (double *x, int dim) const
{
        if (!boundingBox || (boundingBoxDim != dim))
                return false;
        
        bool outside = false;
        const double *boxPtr (boundingBox);
        for (int i = 0; i < dim; ++ i)
        {
                if (x [i] <= *boxPtr)
                {
                        x [i] = *boxPtr;
                        outside = true;
                }
                ++ boxPtr;
                if (x [i] >= *boxPtr)
                {
                        x [i] = *boxPtr;
                        outside = true;
                }
                ++ boxPtr;
        }
 
        return outside;
} /* MapOdeTimeT::checkOutside */
 
inline bool MapOdeTimeT::checkOutside (double *left, double *right,
        int dim) const
{
        if (!boundingBox || (boundingBoxDim != dim))
                return false;
        
        bool outside = false;
        const double *boxPtr (boundingBox);
        for (int i = 0; i < dim; ++ i)
        {
                if (right [i] <= *boxPtr)
                {
                        left [i] = right [i] = *boxPtr;
                        outside = true;
                }
                ++ boxPtr;
                if (left [i] >= *boxPtr)
                {
                        left [i] = right [i] = *boxPtr;
                        outside = true;
                }
                ++ boxPtr;
        }
 
        return outside;
} /* MapOdeTimeT::checkOutside */
 
inline bool MapOdeTimeT::checkOutside (IntervalVectorType &v, int dim) const
{
        if (!boundingBox || (boundingBoxDim != dim))
                return false;
        
        bool outside = false;
        const double *boxPtr (boundingBox);
        for (int i = 0; i < dim; ++ i)
        {
                if (v [i]. rightBound () < *boxPtr)
                {
                        v [i] = *boxPtr;
                        outside = true;
                }
                ++ boxPtr;
                if (v [i]. leftBound () > *boxPtr)
                {
                        v [i] = *boxPtr;
                        outside = true;
                }
                ++ boxPtr;
        }
 
        return outside;
} /* MapOdeTimeT::checkOutside */
 
inline bool MapOdeTimeT::adjust (bool successful, int subdiv)
{
        return odeIntegrParam. adjust (successful, subdiv);
} /* MapOdeTimeT::adjust */
 
inline std::string MapOdeTimeT::getOptInfo () const
{
        return odeIntegrParam. getInfo ();
} /* MapOdeTimeT::getOptInfo */
 
inline void MapOdeTimeT::useOptInfo (const std::string &info)
{
        odeIntegrParam. useInfo (info);
        return;
} /* MapOdeTimeT::useOptInfo */
 
inline void MapOdeTimeT::saveInternals (chomp::multiwork::mwData &data) const
{
        // save the internals of the parent class object
        MapTight::saveInternals (data);
 
        // save the information about the map optimization parameters
        data << odeIntegrParam;
 
        return;
} /* MapOdeTimeT::saveInternals */
 
inline void MapOdeTimeT::loadInternals (chomp::multiwork::mwData &data)
{
        // load the internals of the parent class object
        MapTight::loadInternals (data);
 
        // load the information about the map optimization parameters
        data >> odeIntegrParam;
 
        return;
} /* MapOdeTimeT::loadInternals */
 
inline double MapOdeTimeT::getTimeStep (int subdiv) const
{
        int order (0);
        double step (0);
        int n (0);
        odeIntegrParam. get (order, step, n, subdiv);
        return (n * step);
} /* MapOdeTimeT::getTimeStep */
 
 
#endif // _CMGRAPHS_ODETIMET_H_