awesome-attractor/Attractor.cpp

#include <cmath>
#include <iostream>
#include <fstream>
#include <cstdlib>
#include <iomanip>
#include <cassert>
#include <string>
using namespace std;

#include "Attractor.hpp"


/*
    Constructors & initialisers
*/
Attractor::Attractor():
        dim(3), par(4), formula(LORENZ), param(NULL), point(NULL), new_point(NULL) {
    // Default attractor: 3D Lorenz attrractor
    init_vector();
    init_param();

    param[0] = 0.001;       //dt
    param[1] = 10;          //sigma
    param[2] = 28;          //rho
    param[3] = 8.0/3.0;     //beta

}

Attractor::Attractor(unsigned int dimensions, FormulaChoice formula, unsigned int orde):
        dim(dimensions), formula(formula), orde(orde), param(NULL), point(NULL), new_point(NULL) {

    init(dimensions, formula, orde);
}

Attractor::Attractor(const char* const filename) {
    ifstream file(filename);

    cout << "Reading file " << filename << "..." << endl;

    if ( !file ) {
        cout << " Error reading file '" << filename << "' dying now..." << endl;
        exit(2);
    }

    /*
        file.attr:
            lorenz
            3
            0

            3.24454
            1.25
            ....
    */

    string fileFormula;
    file >> fileFormula;

    for ( unsigned int i = 0; fileFormula[i] != '\0'; i++ ) {
        fileFormula[i] = tolower(fileFormula[i]);
    }

    unsigned int fileDim;
    file >> fileDim;

    unsigned int fileOrde;
    file >> fileOrde;

    cout << " Formula: " << fileFormula << endl;
    cout << " Dimensions: " << fileDim << endl;
    cout << " Orde: " << fileOrde << endl;

    if ( fileFormula == "lorenz" )
        init(fileDim, LORENZ, fileOrde);
    else if ( fileFormula == "poly_n" )
        init(fileDim, POLY_N, fileOrde);
    else if ( fileFormula == "poly_a" )
        init(fileDim, POLY_A, fileOrde);
    else if ( fileFormula == "logistic" )
        init(fileDim, LOGISTIC, fileOrde);
    else if ( fileFormula == "unravel" )
        init(fileDim, UNRAVEL, fileOrde);
    else {
        cout << " Formula not (yet) supported" << endl;
        exit(3);
    }

    for ( unsigned int i = 0; i < par; i++ ) {
        file >> param[i];
        cout << " Parameter " << i << " set to " << param[i] << ",  ";
    }

    cout << endl << " Reading file complete" << endl;
}

void Attractor::init(unsigned int dim_in, FormulaChoice formula_in, unsigned int orde_in) {

    dim = dim_in;
    formula = formula_in;
    orde = orde_in;

    switch (formula) {
    case POLY_N: {
        double n_coef = orde + 1;
        for (unsigned int i = 2; i <= dim; i++) {
            n_coef = n_coef*(orde + i)/(i - 1);
        }

        par = (unsigned int) n_coef;
        break;
    }
    case LORENZ: {
        assert(dim == 3 || dim == 4);
        par = dim + 1;
        break;
    }
    case POLY_A: {
        par = dim;
        break;
    }
    case LOGISTIC: {
        par = dim;
        break;
    }
    case UNRAVEL: {
        assert(dim == 3);
        par = 7;
        break;
    }
    default: {
        cout << "Formula not recognized" << endl;
        exit(1);
        break;
    }
    }

    init_vector();
    init_param();
}

void Attractor::init_vector() {
    point = new double[dim];
    new_point = new double[dim];

    assert(point != NULL);
    assert(new_point != NULL);

    for ( unsigned int i = 0; i < dim; i++ ) {
        point[i] = new_point[i] = 0.9;
    }
}

void Attractor::init_param() {
    param = new double[par];

    assert(param != NULL);

    for ( unsigned int i = 0; i < dim; i++ ) {
        param[i] = 0.0;
    }
}

void Attractor::init_range() {

    // stabilize attractor
    for ( unsigned int i = 0; i < 100000; i++ ) {
        iterate();
        if ( !is_chaos() ) {
            cout << "Attractor died after " << i << " iterations" << endl;
            exit(0);
        }
    }

    // initialize ranges
    for ( vector<Projector*>::iterator it = projectors.begin(); it != projectors.end(); it++ ) {
        (*it)->extern_dim = dim;
        (*it)->intern_dim = 2;
        (*it)->init(point);
    }

    // update ranges
    for ( unsigned int i = 0; i < 100000; i++ ) {
        iterate();
        if ( !is_chaos() ) {
            cout << "Attractor died after " << i << " iterations" << endl;
            exit(0);
        }

        for ( vector<Projector*>::iterator it = projectors.begin(); it != projectors.end(); it++ ) {
            (*it)->update_range(point);
        }

    }
    for ( vector<Projector*>::iterator it = projectors.begin(); it != projectors.end(); it++ ) {
        (*it)->finish_range();
    }
}

bool Attractor::is_chaos() {
    /*
        check existence of attractor:
            Escaping
            Single point attractor
            Lyapunov exponent
    */

    double sum = 0;
    for ( unsigned int i = 0; i < dim; i++ ) {
        const double dist = new_point[i] - point[i];
        sum += dist*dist;
    }
    if ( sum >= 1.0e7 ) {
        // big change => Escaping
        return false;
    }
    if ( sum <= 1.0e-7 ) {
        // small change => singularity
        return false;
    }
    return true;
}

/*
    Iteration & Math
*/
void Attractor::iterate() {
    // pointer swap
    double * temp = point;
    point = new_point;
    new_point = temp;

    // calculations
    switch (formula) {
    case POLY_N:
        polynome();
        break;

    case POLY_A:
        poly_A();
        break;

    case POLY_2:
        poly_2();
        break;

    case LOGISTIC:
        logistic();
        break;

    default:
        cout << "Formula not recognized" << endl;
        exit(1);
        break;
    }

#ifdef HARDDEBUG
    cout << "Formula " << formula << " is done" << endl;
    cout << "Old Vector: ";
    for ( unsigned int i = 0; i < dim; i++ ) {
        cout << point[i] << " ";
    }
    cout << endl << "Fresh Vector: ";
    for ( unsigned int i = 0; i < dim; i++ ) {
        cout << new_point[i] << " ";
    }
    cout << endl;
#endif
}

void Attractor::polynome() {
    unsigned int m = 0;
    for ( unsigned int i = 0; i < dim; i++ ) {

#ifdef HARDDEBUG
        cout << "Entering new dimension: " << i << " With m = " << m << endl;
#endif

        new_point[i] = param[m];
        m++;
        recur(i, 0, 1, m);
    }

}

void Attractor::recur(unsigned int curr_dimension, unsigned int prev_i, unsigned int n, unsigned int& m, double prev_product) {
    double product;
    for (unsigned int i = prev_i; i < dim; i++) {

#ifdef HARDDEBUG
        for ( unsigned int j = 0; j < n; j++ )
            cout << "  ";
        cout << "Calculation in dimension: " << i << " With m = " << m << " And depth = " << n << endl;
#endif

        product = prev_product * point[i];
        new_point[curr_dimension] += param[m] * product;
        m++;
        if (n < orde) {
            recur(curr_dimension, i, n+1, m, product);
        }
    }
}

void Attractor::poly_A() {
    switch (dim) {
    case 3:
        new_point[0] = param[0] + point[1] - point[1]*point[2];
        new_point[1] = param[1] + point[2] - point[2]*point[0];
        new_point[2] = param[2] + point[0] - point[0]*point[1];
        break;
    default:
        for ( unsigned int i = 0; i < dim; i++ ) {
            new_point[i] = param[i] + point[(i+1) % dim] - point[(i+1) % dim]*point[(i+2) % dim];
        }
        break;
    }
}

void Attractor::poly_2() {
    exit(2);
}

void Attractor::logistic() {
    for ( unsigned int i = 0; i < dim; i++ ) {
        new_point[i] = param[i]*point[i]*(1.0 - point[i]);
    }
}


void Attractor::plot() {
    for ( vector<Projector *>::iterator it = projectors.begin(); it != projectors.end(); it++ ) {
        (*it)->plot(point);
    }
}


/*
    IO & control
*/
void Attractor::output() {
    for ( unsigned int i = 0; i < dim; i++ ) {
        cout << point[i] << " ";
    }
    cout << endl;
}