odeapprox.h

Go to the documentation of this file.

/////////////////////////////////////////////////////////////////////////////
///
/// @file odeapprox.h
///
/// A generic method for computing a time-t map for ODEs
/// in an approximate way.
///
/// This file defines a generic class which uses a non-rigorous integration
/// method based on Runge-Kutta applied to sample points in an interval box
/// in order to compute an approximation of 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 December 16, 2010. Last revision: May 30, 2011.
 
 
#ifndef _CMGRAPHS_ODEAPPROX_H_
#define _CMGRAPHS_ODEAPPROX_H_
 
 
// include some standard C++ header files
#include <string>
#include <vector>
#include <iostream>
 
// include local header files
#include "maptype.h"
#include "odetimet.h"
 
// include all the necessary header files from the CHomP library
#include "chomp/system/config.h"
#include "chomp/system/textfile.h"
 
 
// --------------------------------------------------
// ----------------- 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 MapOdeApprox: public MapOdeTimeT
{
public:
        /// The constructor in which one sets the step size
        /// and the number of steps for the non-rigorous integration,
        /// the number of sample points in each direction,
        /// and the margin to add around the box (as a fraction
        /// of its actual width).
        MapOdeApprox (double step, int nSteps, int nSamples,
                double margin = 0.01);
 
        /// The destructor.
        ~MapOdeApprox ();
 
        /// Sets the parameters to the given intervals defined by their
        /// left and right coordinates.
        /// This function calls the original function from the
        /// class MapType, and skips the function from MapOdeTimeT.
        void setParam (const double *left, const double *right, int n);
 
        /// Computes the image of a box whose left and right coordinates
        /// are given. Fills in the images with the left and right
        /// coordinates of the image box.
        void compute (const double *xleft, const double *xright,
                double *yleft, double *yright, int dim,
                const spcCoord *coord, int subdiv) const;
 
private:
        /// The size of the step for itegrating the flow.
        double step_;
 
        /// The number of steps to use.
        int nSteps_;
 
        /// The number of samples in each direction.
        int nSamples_;
 
        /// The margin to add around the box (a fraction of its width).
        double margin_;
 
        /// The vector field. Parameters (if any) should be determined
        /// using the methods of the class MapType: getLeftParam
        /// and getRightParam.
        virtual void vectorField (double t, const double *x, double *y,
                int dim) const = 0;
 
        /// The copy constructor should not be used.
        MapOdeApprox (const MapOdeApprox &);
 
        /// The assignment operator should not be used.
        MapOdeApprox &operator = (const MapOdeApprox &);
 
        /// The Runge-Kutta method for iterating a point,
        /// using the vector field defined in this class.
        /// The buffer should be preallocated for 5*dim entries.
        void RungeKutta (double t, double *x, double step, int dim,
                double *buffer) const;
 
        /// Computes the ranges of coordinates.
        static void computeRanges (const double *samples, int nSamples,
                int dim, double *minimal, double *maximal);
 
        /// Makes a correction to the vector field to make it fade
        /// gradually outside the bounding box if one is in defined
        /// and in use for the parent class MapOdeTimeT.
        void fadeVectorField (const double *x, double *y, int dim) const;
 
}; /* class MapOdeApprox */
 
// --------------------------------------------------
 
inline MapOdeApprox::MapOdeApprox (double step, int nSteps, int nSamples,
        double margin):
        MapOdeTimeT ("", 0, 0, 0),
        step_ (step), nSteps_ (nSteps), nSamples_ (nSamples),
        margin_ (margin)
{
//      chomp::homology::sbug << "DEBUG: MapOdeApprox constructor.\n";
        if (step_ <= 0)
                throw "Non-positive step requested in MapOdeApprox.";
        if (nSteps_ <= 0)
        {
                throw "Non-positive number of steps requested "
                        "in MapOdeApprox.";
        }
        if (nSamples_ <= 0)
        {
                throw "Non-positive number of samples requested "
                        "in MapOdeApprox.";
        }
        static bool displayed = false;
        if (!displayed)
        {
                chomp::homology::sbug << "ODE approx: " << nSteps <<
                        " steps of " << step << " (1/" << (1.0 / step) <<
                        ") each, " << nSamples <<
                        " samples in each dir.\n";
                displayed = true;
        }
        return;
} /* MapOdeApprox::MapOdeApprox */
 
inline MapOdeApprox::~MapOdeApprox ()
{
//      chomp::homology::sbug << "DEBUG: MapOdeApprox destructor.\n";
        return;
} /* MapOdeApprox::~MapOdeApprox */
 
inline void MapOdeApprox::setParam (const double *left, const double *right,
        int n)
{
        // set the new values of parameters (skip the MapOdeTimeT method)
        MapType::setParam (left, right, n);
 
        return;
} /* MapOdeApprox::setParam */
 
inline void MapOdeApprox::fadeVectorField (const double *x, double *y,
        int dim) const
{
        if (!boundingBox || (boundingBoxDim != dim))
                return;
 
        // a fraction of the width of the bounding box for the fading
        const double fadeBegin = 0.1;
        const double fadeEnd = 0.2;
 
        // the sides of the bounding box (left0, right0, left1, right1, etc.)
        const double *sides = boundingBox;
 
        // make the fading in each direction separately
        for (int i = 0; i < dim; ++ i)
        {
                // determine the left and right sides of the bounding box
                double left = * (sides ++);
                double right = * (sides ++);
 
                // compute the distance from the box if the point is outside
                double distance = 0;
                if (x [i] < left)
                        distance = left - x [i];
                else if (x [i] > right)
                        distance = x [i] - right;
                else
                        continue;
 
                // if the distance is too small then skip this case
                double distBegin = fadeBegin * (right - left);
                if (distance <= distBegin)
                        continue;
 
                // fade the vector field outside the box
                double distEnd = fadeEnd * (right - left);
                if (distance >= distEnd)
                        y [i] = 0;
                else
                {
                        y [i] *= (distEnd - distance) /
                                (distEnd - distBegin);
                }
        }
 
        return;
} /* MapOdeApprox::fadeVectorField */
 
inline void MapOdeApprox::RungeKutta (double t, double *x, double step,
        int dim, double *buffer) const
{
        double *k1 = buffer;
        double *k2 = k1 + dim;
        double *k3 = k2 + dim;
        double *k4 = k3 + dim;
        double *tmp = k4 + dim;
 
        vectorField (t, x, k1, dim);
        fadeVectorField (x, k1, dim);
 
        for (int i = 0; i < dim; ++ i)
                tmp [i] = x [i] + 0.5 * step * k1 [i];
        vectorField (t + 0.5 * step, tmp, k2, dim);
        fadeVectorField (tmp, k2, dim);
 
        for (int i = 0; i < dim; ++ i)
                tmp [i] = x [i] + 0.5 * step * k2 [i];
        vectorField (t + 0.5 * step, tmp, k3, dim);
        fadeVectorField (tmp, k3, dim);
 
        for (int i = 0; i < dim; ++ i)
                tmp [i] = x [i] + step * k3 [i];
        vectorField (t + step, tmp, k4, dim);
        fadeVectorField (tmp, k4, dim);
 
        for (int i = 0; i < dim; ++ i)
        {
                x [i] += step / 6 * (k1 [i] + 2 * k2 [i] +
                        2 * k3 [i] + k4 [i]);
        }
 
        return;
} /* MapOdeApprox::RungeKutta */
 
inline void MapOdeApprox::computeRanges (const double *samples, int nSamples,
        int dim, double *minimal, double *maximal)
{
        const double *sample = samples;
        for (int n = 0; n < nSamples; ++ n)
        {
                for (int i = 0; i < dim; ++ i)
                {
                        if (!n || (minimal [i] > sample [i]))
                                minimal [i] = sample [i];
                        if (!n || (maximal [i] < sample [i]))
                                maximal [i] = sample [i];
                }
                sample += dim;
        }
        return;
} /* MapOdeApprox::computeRanges */
 
inline void MapOdeApprox::compute (const double *xleft, const double *xright,
        double *yleft, double *yright, int dim, const spcCoord *, int) const
{
        // for debugging only
        using chomp::homology::sbug;
 
        // calculate the number of points to iterate
        int countPoints = nSamples_;
        int nSamplesInterior = (nSamples_ > 2) ? (nSamples_ - 2) : 0;
        int interiorPoints = nSamplesInterior;
        for (int i = 1; i < dim; ++ i)
        {
                countPoints *= nSamples_;
                interiorPoints *= nSamplesInterior;
        }
        countPoints -= interiorPoints;
 
        // allocate an array for the coordinates of sample points to iterate
        double *samples = new double [(countPoints + 5) * dim];
        double *buffer = samples + (countPoints * dim);
 
        // compute the coordinates of the sample points to iterate
        int *counters = new int [dim];
        for (int i = 0; i < dim; ++ i)
                counters [i] = nSamples_ - 1;
        int curPoint = 0;
        double *curSample = samples;
        while (curPoint < countPoints)
        {
                // increase the multi-dimensional counter
                for (int curIndex = 0; curIndex < dim; ++ curIndex)
                {
                        if (++ counters [curIndex] < nSamples_)
                                break;
                        counters [curIndex] = 0;
                }
 
                // skip this counter if this is not a boundary point
                bool boundary = false;
                for (int i = 0; i < dim; ++ i)
                {
                        if ((counters [i] == 0) ||
                                (counters [i] == nSamples_ - 1))
                        {
                                boundary = true;
                                break;
                        }
                }
                // debug begin
        //      sbug << "DEBUG: Counters " << curPoint << ":";
        //      for (int i = 0; i < dim; ++ i)
        //              sbug << ' ' << counters [i];
        //      sbug << (boundary ? " [bd]\n" : " [int]\n");
                // debug end
                if (!boundary)
                        continue;
 
                // compute the coordinates of the point
                for (int i = 0; i < dim; ++ i)
                {
                        if (counters [i] == 0)
                                *curSample = xleft [i];
                        else if (counters [i] == nSamples_ - 1)
                                *curSample = xright [i];
                        else
                        {
                                double fraction = static_cast<double>
                                        (counters [i]) / (nSamples_ - 1);
                                *curSample = xleft [i] + fraction *
                                        (xright [i] - xleft [i]);
                        }
                        ++ curSample;
                }
                ++ curPoint;
        }
        // debug check: the counter should be reset to the original one
        for (int i = 0; i < dim; ++ i)
        {
                if (counters [i] != nSamples_ - 1)
                        throw "Non-zero counter, bad programming, sorry!";
        }
        delete [] counters;
 
        // iterate the points
//      int counter = 0;
        bool interrupted = false;
        for (int step = 0; step < nSteps_; ++ step)
        {
                // move the set of sample points one step forward
                double t = step * step_;
                double *sample = samples;
                for (int n = 0; n < countPoints; ++ n)
                {
                        RungeKutta (t, sample, step_, dim, buffer);
                //      checkOutside (sample, dim);
                        sample += dim;
                }
 
                // check if the set went outside the given region
        /*      if (++ counter >= checkEvery)
                {
                        // compute the ranges of the coordinates
                        computeRanges (samples, countPoints, dim,
                                yleft, yright);
                        interrupted = checkOutside (yleft, yright, dim);
                        if (interrupted)
                                break;
                        counter = 0;
                }
        */
        }
 
        // compute the ranges of the coordinates of the set of points
        if (!interrupted)
                computeRanges (samples, countPoints, dim, yleft, yright);
 
        // add the margin around the resulting box
        if (margin_ > 0)
        {
                for (int i = 0; i < dim; ++ i)
                {
                        double collar = (yright [i] - yleft [i]) * margin_;
                        yleft [i] -= collar;
                        yright [i] += collar;
                }
        }
 
        delete [] samples;
        return;
} /* MapOdeApprox::compute */
 
 
#endif // _CMGRAPHS_ODEAPPROX_H_