attractor.h

Go to the documentation of this file.

/////////////////////////////////////////////////////////////////////////////
///
/// @file attractor.h
///
/// Construction of a cover of an attractor.
/// This file contains the definition of a function that constructs
/// a cover of an attractor for a given dynamical system, and computes
/// a directed graph representation of the dynamics on this cover.
///
/// @author Pawel Pilarczyk
///
/////////////////////////////////////////////////////////////////////////////

// Copyright (C) 1997-2010 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 2 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, write to the
// Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
// MA 02111-1307, USA.

// Started on March 13, 2009. Last revision: June 18, 2009.


#ifndef _FINRESDYN_ATTRACTOR_H_
#define _FINRESDYN_ATTRACTOR_H_


// include some standard C++ header files
#include <iostream>
#include <memory>
#include <cmath>

// include local header files
#include "graph.h"
#include "streams.h"
#include "rounding.h"


// --------------------------------------------------
// ------------------- edge adder -------------------
// --------------------------------------------------

/// Auxiliary class for adding edges to a graph.
template <class GraphType>
class EdgeAdder
{
public:
        /// The only allowed constructor of an edge adder object.
        EdgeAdder (GraphType &graph): g (graph) {}

        /// Adds an edge to the graph.
        void add (int edge) {g. addEdge (edge); return;}

private:
        /// A reference to the graph to which the vertices
        /// are to be added.
        GraphType &g;

}; /* class EdgeAdder */


// --------------------------------------------------
// ------------------- diameter2 --------------------
// --------------------------------------------------

/// Computes an upper bound for the square of the diameter of a given box.
template <class NumType>
inline NumType diameter2 (const int dim,
        const NumType *leftBounds, const NumType *rightBounds)
{
        NumType squares = 0;
        for (int i = 0; i < dim; ++ i)
        {
                NumType diff = tRounding<NumType>::sub_up
                        (rightBounds [i], leftBounds [i]);
                NumType diff2 = tRounding<NumType>::mul_up (diff, diff);
                squares = tRounding<NumType>::add_up (squares, diff2);
        }
        return squares;
} /* diameter2 */


// --------------------------------------------------
// ------------------- attractor --------------------
// --------------------------------------------------

/// Constructs a cover of an attractor for a given dynamical system.
/// Returns an upper bound for the square of the diameter of map images.
/// Throws an exception in case of failure.
template <class MapType, class CoverType, class GraphType, class NumArray>
inline double constructAttractor (const MapType &function, CoverType &cover,
        GraphType &graph, NumArray initialPoint, int iterCount,
        int maxCoverSize)
{
        // determine the type of numbers and the space dimension
        typedef typename CoverType::NumberType NumType;
        int dim = cover. dim ();

        // create an initial point close to the attractor and cover it
        std::auto_ptr<NumType> xPtr (new NumType [dim]);
        NumType *x = xPtr. get ();
        for (int i = 0; i < dim; ++ i)
                x [i] = initialPoint [i];
        function. initialPoint (x, iterCount);
        cover. cover (x);

        // prepare various objects to be used during the computations
        EdgeAdder <GraphType> edgeAdder (graph);
        std::auto_ptr<NumType> yPtr (new NumType [dim << 1]);
        NumType *leftBounds = yPtr. get ();
        NumType *rightBounds = leftBounds + dim;
        NumType *leftImage = leftBounds;
        NumType *rightImage = rightBounds;
        NumType maxDiameter2 = 0;

        // compute the images of the points until the entire attractor
        // has been covered
        int current = 0;
        while (current < cover. size ())
        {
                // add a vertex corresponding to this cover element
                graph. addVertex ();

                // get a box which covers the given element of the cover
                cover. get (current, leftBounds, rightBounds);

                // compute a rectangular cover of the image of this box
                function. image (leftBounds, rightBounds,
                        leftImage, rightImage);

                // cover the computed box with elements of the cover
                NumType diam2;
                cover. cover (leftImage, rightImage, edgeAdder, &diam2);

                // update the maximal square of the diameter of images
                if (maxDiameter2 < diam2)
                        maxDiameter2 = diam2;

                // interrupt if the cover has become too large
                if ((maxCoverSize > 0) && (cover. size () > maxCoverSize))
                        throw "Maximal allowed cover size exceeded.";

                // take the next element of the cover into consideration
                ++ current;

                // show progress indicator
                if (!(current % 1371))
                {
                        scon << std::setw (10) << cover. size () <<
                                "\b\b\b\b\b\b\b\b\b\b";
                }
        }

        return static_cast<double> (maxDiameter2);
} /* constructAttractor */


#endif // _FINRESDYN_ATTRACTOR_H_