#ifndef ATTRACTORKERNEL_HPP
#define ATTRACTORKERNEL_HPP

#include "Logger.hpp"
#include <algorithm>
#include <cstdlib>

#include "stfu/stf.hpp"

class AttractorKernel {
private:

#pragma mark -
#pragma mark memory
	void allocate() {
		parameters = new double[numberOfParameters];
		vectorNew = new double[dimension];
		vectorOld = new double[dimension];
	}

	void dealloc() {
		delete[] vectorOld;
		vectorOld = NULL;
		delete[] vectorNew;
		vectorNew = NULL;
		delete[] parameters;
		parameters = NULL;
	}

protected:
	double* parameters;
	double* vectorNew;
	double* vectorOld;

	unsigned int numberOfParameters;
	unsigned int dimension;

	AttractorKernel(const unsigned int dimension, const unsigned int numberOfParameters) :
		numberOfParameters(numberOfParameters), dimension(dimension) {

		try {
			allocate();
		} catch(std::exception& e) {
			LogError("Couldn't construct Attractorkernel: %s\n", e.what());
			dealloc();
		}

		std::fill_n(parameters, numberOfParameters, 0.0);
		std::fill_n(vectorNew, dimension, 0.0);
		std::fill_n(vectorOld, dimension, 0.0);
	}

public:
	virtual ~AttractorKernel() {
		dealloc();
	}

#pragma mark -
#pragma mark parameters
	double& operator[](const unsigned int index) {
		return parameters[index];
	}
	double const& operator[](const unsigned int index) const {
		return parameters[index];
	}
	unsigned int getNumberOfParameters() const {
		return numberOfParameters;
	}

	virtual void generate_random_parameters() {
		for(unsigned int i = 0; i < numberOfParameters; ++i)
			parameters[i] = 2.0 * rand() / double(RAND_MAX) - 1.0;
		for(unsigned int i = 0; i < dimension; ++i)
			vectorNew[i] = vectorOld[i] = 6.0 * rand() / double(RAND_MAX) - 3.0;
	}

#pragma mark -
#pragma mark tests
	bool convergent() {
		double sum = 0.0;
		for(unsigned int i = 0; i < dimension; ++i) {
			sum += (vectorNew[i] - vectorOld[i])*(vectorNew[i] - vectorOld[i]);
		}
		return sum < 1.0e-6;
	}

	bool divergent() {
		double sum = 0.0;
		for(unsigned int i = 0; i < dimension; ++i) {
			sum += (vectorNew[i] - vectorOld[i])*(vectorNew[i] - vectorOld[i]);
		}
		return sum > 1.0e3;
	}

#pragma mark -
#pragma mark main
	// iterate his formula, implemented by subclasses
	virtual void operator()() = 0;

#pragma mark -
#pragma mark vector
	double const* vector() const {
		return vectorNew;
	}

	double const* previousVector() const {
		return vectorOld;
	}

	unsigned int getDimension() const {
		return dimension;
	}

#pragma mark -
#pragma mark factory function
	static AttractorKernel* createAttractorKernel(stfu::node& attractorKernel);

};


#endif // ATTRACTORKERNEL_HPP