odeintegr.h

Go to the documentation of this file.

/////////////////////////////////////////////////////////////////////////////
///
/// @file odeintegr.h
///
/// ODE integration parameters: guessing, adjusting, saving.
///
/// @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 March 24, 2013. Last revision: November 27, 2013.
 
 
#ifndef _CMGRAPHS_ODEINTEGR_H_
#define _CMGRAPHS_ODEINTEGR_H_
 
 
// include some standard C++ header files
#include <string>
#include <sstream>
#include <vector>
#include <iostream>
#include <iomanip>
#include <algorithm>
#include <memory>
#include <utility>
 
// 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 "datavector.h"
#include "mapparam.h"
#include "datatext.h"
 
 
// --------------------------------------------------
// ----------- ODE integration parameters -----------
// --------------------------------------------------
 
/// A class whose objects are responsible for adjusting the integration
/// parameters such as step size and the number of steps for an ODE.
class OdeIntegrParam
{
public:
        /// The constructor which takes the suggested values for the order
        /// of the method, the step size, and the number of steps.
        OdeIntegrParam (int order, double step, int n);
 
        /// The destructor.
        ~OdeIntegrParam ();
 
        /// Sets the parameters to the given intervals defined by their
        /// left and right coordinates.
        void setParam (const double *left, const double *right, int n);
 
        /// Provides recommended values for integration:
        /// the order of the Taylor method, the size of the integration step,
        /// and the number of steps.
        void get (int &order, double &step, int &n, int subdiv);
 
        /// 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);
 
        /// Returns a line (or lines) with a summary of information
        /// about the knowledge gathered through the "adjust" function,
        /// for later use.
        std::string getInfo () const;
 
        /// Analyzes a previously saved line with a summary of information
        /// and adjusts the data accordingly.
        void useInfo (const std::string &info);
 
        /// Saves the data for serialization.
        void saveData (chomp::multiwork::mwData &data) const;
 
        /// Retrieves a serialized object.
        void loadData (chomp::multiwork::mwData &data);
 
private:
        /// The copy constructor should not be used.
        OdeIntegrParam (const OdeIntegrParam &);
 
        /// The assignment operator should not be used.
        OdeIntegrParam &operator = (const OdeIntegrParam &);
 
        /// The parameters of the map.
        MapParam param;
 
        /// Step sizes at subsequent subdivision depths.
        std::vector<double> stepSizes;
 
        /// Strategy for trying another step size at each subdivision depth.
        /// 0 = increase if successful (no attempt has yet been made);
        /// 1 = increase if successful (previous attempt was successful);
        /// 2 = an unsuccessful attempt encountered, don't try anything else
        ///     if successful, decrease if unsuccessful;
        /// 3 = don't increase even if successful,
        ///     it is known that this is a good guess already.
        std::vector<int> strategy;
 
        /// Extends the internal arrays to make sure that the entries
        /// for the given subdivision depth are valid.
        void extendVectors (int subdiv);
 
        /// The order of the Taylor method suggested to use.
        int defaultOrder;
 
        /// The size of the step suggested for itegrating the flow.
        double defaultStepSize;
 
        /// The number of steps suggested to use.
        int defaultNumberOfSteps;
 
        /// Was this information set up based on previous results?
        bool usingPrevious;
 
        /// Were the parameters changed, so that it is worth to share them?
        bool changed;
 
}; /* class OdeIntegrParam */
 
// --------------------------------------------------
 
inline OdeIntegrParam::OdeIntegrParam (int order, double step, int n):
        defaultOrder (order), defaultStepSize (step),
        defaultNumberOfSteps (n),
        usingPrevious (false), changed (false)
{
        return;
} /* OdeIntegrParam::OdeIntegrParam */
 
inline OdeIntegrParam::~OdeIntegrParam ()
{
        return;
} /* OdeIntegrParam::~OdeIntegrParam */
 
inline OdeIntegrParam::OdeIntegrParam (const OdeIntegrParam &)
{
        throw "Trying to use the copy constructor of OdeIntegrParam.";
        return;
} /* OdeIntegrParam::OdeIntegrParam */
 
inline OdeIntegrParam &OdeIntegrParam::operator = (const OdeIntegrParam &)
{
        throw "Trying to use the assignment operator of OdeIntegrParam.";
        return *this;
} /* OdeIntegrParam::OdeIntegrParam */
 
inline void OdeIntegrParam::setParam (const double *left, const double *right,
        int n)
{
        param. define (left, right, n);
        return;
} /* OdeIntegrParam::setParam */
 
inline void OdeIntegrParam::get (int &order, double &step, int &n,
        int subdiv)
{
        // make sure that the entries for the subdivision depth are defined
        extendVectors (subdiv);
 
        // prepare the order of the Taylor method
        order = defaultOrder;
 
        // prepare the step size
        step = stepSizes [subdiv];
 
        // prepare the number of steps
        n = defaultNumberOfSteps;
 
        // show info on the step size if necessary
        static int prevSubdiv = 0;
        if (subdiv > prevSubdiv)
        {
                using chomp::homology::sbug;
                sbug << "Using step size " << step << " at depth " <<
                        subdiv << ".\n";
                prevSubdiv = subdiv;
        }
 
        return;
} /* OdeIntegrParam::get */
 
inline bool OdeIntegrParam::adjust (bool successful, int subdiv)
{
        using chomp::homology::sbug;
 
        // extend the vectors to make sure that the entries for the given
        // subdivision depth are defined
        extendVectors (subdiv);
 
        // prepare coefficients for increasing and decreasing step size
        const double incCoef = 2.0;
        const double decCoef = 0.9;
 
        // if successful then increase the step size, unless not feasible
        if (successful)
        {
                if ((strategy [subdiv] == 2) || (strategy [subdiv] == 3))
                        return false;
                sbug << "Good!\nSuccessful integration at depth " <<
                        subdiv << " with step size = " <<
                        stepSizes [subdiv] << ". Trying ";
                stepSizes [subdiv] *= incCoef;
                sbug << stepSizes [subdiv] << ".\n";
                changed = true;
                return true;
        }
        changed = true;
 
        // make a note that an unsuccessful attempt took place
        strategy [subdiv] = 2;
 
        // show a message that the integration failed
        sbug << "Integration at depth " << subdiv <<
                " failed with step size " << stepSizes [subdiv] << ". ";
 
        // decrease the step size and suggest to try again
        // (in a more sophisticated version, the decrease can be different,
        // depending on the largest successful step size
        // encountered at previous computations)
        stepSizes [subdiv] *= decCoef;
        sbug << "Trying " << stepSizes [subdiv] << ".\n";
 
        return true;
} /* OdeIntegrParam::adjust */
 
inline void OdeIntegrParam::extendVectors (int subdiv)
{
        // extend the vector with step sizes
        while (static_cast<int> (stepSizes. size ()) <= subdiv)
        {
                if (static_cast<int> (stepSizes. size ()) <= initialDepth)
                {
                        stepSizes. push_back (defaultStepSize);
                }
                else if (static_cast<int> (stepSizes. size ()) <= finalDepth)
                {
                        stepSizes. push_back
                                (2 * stepSizes [stepSizes. size () - 1]);
                }
                else
                {
                        stepSizes. push_back
                                (stepSizes [stepSizes. size () - 1]);
                }
                changed = true;
        }
 
        // extend the vector with strategy numbers
        while (static_cast<int> (strategy. size ()) <= subdiv)
        {
                strategy. push_back (0);
        }
 
        return;
} /* OdeIntegrParam::extendVectors */
 
inline std::string OdeIntegrParam::getInfo () const
{
        // if there is no new information then don't return anything
        if (usingPrevious && !changed)
                return std::string ();
 
        // prepare a string stream to write the data to
        std::ostringstream s;
 
        // save the map parameters
        s << param;
 
#ifdef ODEPARAM_APPROXIMATE
        // save an identification code
        s << ' ' << 'S';
 
        // save the vector stepSizes
        s << ' ' << stepSizes. size ();
        int precision = s. precision ();
        int hiprec = 12;
        s << std::setprecision (hiprec);
        for (size_t n = 0; n < stepSizes. size (); ++ n)
                s << ' ' << stepSizes [n];
        s << std::setprecision (precision);
#else
        // save an identification code
        s << ' ' << 'B';
 
        // serialize the new information
        chomp::multiwork::mwData data;
        data << param;
        data << stepSizes;
 
        // encode the binary data to a simple string
        s << ' ' << data2text (data);
#endif
        // write an ending control character
        s << ' ' << 'Q';
 
        // return the string
        return s. str ();
} /* OdeIntegrParam::getInfo */
 
inline void OdeIntegrParam::useInfo (const std::string &info)
{
        // prepare a string stream to read the data from
        using chomp::homology::ignorecomments;
        std::istringstream s (info);
 
        // read the parameters of the map and check if they agree
        MapParam _param;
        s >> _param;
        if (!intersect (param, _param) || !shareInteriors (param, _param))
                return;
        ignorecomments (s);
 
#ifdef ODEPARAM_APPROXIMATE
        // read the identification code
        if (s. get () != 'S')
                return;
 
        // read the vector stepSizes
        size_t stepSizesSize = 0;
        s >> stepSizesSize;
        stepSizes. resize (stepSizesSize);
        for (size_t n = 0; n < stepSizesSize; ++ n)
                s >> stepSizes [n];
#else
        // read the identification code
        if (s. get () != 'B')
                return;
        ignorecomments (s);
 
        // decode the binary data from the input string
        std::string str;
        s >> str;
        chomp::multiwork::mwData data;
        text2data (str, data);
 
        // extract the parameters and check if they agree
        MapParam _paramExact;
        data >> _paramExact;
        if (!(param == _paramExact))
                return;
 
        // extract the step sizes
        data >> stepSizes;
#endif
        // make sure that there is the right control character at the end
        // of the input line
        ignorecomments (s);
        if (s. get () != 'Q')
                throw "Missing 'Q' at the end of ODE integr info.";
 
        // make a note that the previously tested settings will be used
        usingPrevious = true;
        strategy. clear ();
        strategy. resize (stepSizes. size (), 3);
 
        return;
} /* OdeIntegrParam::useInfo */
 
inline void OdeIntegrParam::saveData (chomp::multiwork::mwData &data) const
{
        data << stepSizes;
        data << strategy;
 
        return;
} /* OdeIntegrParam::saveData */
 
inline void OdeIntegrParam::loadData (chomp::multiwork::mwData &data)
{
        data >> stepSizes;
        data >> strategy;
 
        return;
} /* OdeIntegrParam::loadData */
 
// --------------------------------------------------
 
/// Serializes an object with ODE parameters.
inline chomp::multiwork::mwData &operator << (chomp::multiwork::mwData &data,
        const OdeIntegrParam &odeIntegrParam)
{
        odeIntegrParam. saveData (data);
        return data;
} /* operator << */
 
/// Deserializes an object with ODE parameters.
inline chomp::multiwork::mwData &operator >> (chomp::multiwork::mwData &data,
        OdeIntegrParam &odeIntegrParam)
{
        odeIntegrParam. loadData (data);
        return data;
} /* operator >> */
 
 
#endif // _CMGRAPHS_ODEINTEGR_H_