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

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

}

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);
    }

    // TODO : Use stfu

    /*
        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) {

}

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() {
    // attractorKernel->iterate()
}

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;
}