#include <stdlib.h>
#include <stdio.h>

struct MatrixSize {
    int rows;
    int cols;
    int operator[](int dim) { return *(&rows + dim); }
};

// low overhead version of blitz::Array
// default assigment operator and copy constructor suffice in this example
// since we do not need deep copy here and the destructor does nothing
template <typename TYPE>
class TMatrix {
    MatrixSize size;
    TYPE* values;
public:
    TMatrix() : values(0) { size.rows = size.cols = 0; }
    TMatrix(MatrixSize size) : size(size),
			       values(new TYPE[size.rows * size.cols]) {}
    ~TMatrix() { /*delete[] values;*/ }
    TYPE  operator()(int i, int j) const { return values[size.cols * i + j]; }
    TYPE& operator()(int i, int j)       { return values[size.cols * i + j]; }
    void resize(MatrixSize size) {
	this->size = size;
	if (values != 0) delete values;
	values = new TYPE[size.rows * size.cols];
    }
    void reference(TMatrix& other) {
	*this = other;
    }
};
typedef TMatrix<double> Matrix;

#define ROWS 0
#define COLS 1

// But using references here saves time !!
void reaction_diffusion_texture(MatrixSize matrix_size,
				double d_a, double d_b, double s,
				Matrix& a, Matrix& b, Matrix& a_out, Matrix& b_out,
				Matrix& beta)
{
    int i, j;
    for (i = 1; i < matrix_size[ROWS] - 1; ++i) {
	for (j = 1; j < matrix_size[COLS] - 1; ++j) {
	    a_out(i,j) = a(i,j)
		+ s * (16 - a(i,j)*b(i,j))
		+ d_a * ( a(i+1,j) + a(i-1,j) + a(i,j+1) + a(i,j-1)
			  - 4 * a(i,j) );
	    b_out(i,j) = b(i,j)
		+ s * (a(i,j) * b(i,j) - b(i,j) - beta(i,j))
		+ d_b * ( b(i+1,j) + b(i-1,j) + b(i,j+1) + b(i,j-1)
			  - 4 * b(i,j) );
	    if (a_out(i,j) < 0.0) a_out(i,j) = 0.0;
	    if (b_out(i,j) < 0.0) b_out(i,j) = 0.0;
	}
    }
}

int main(int argc, char* argv[]) {
    double init_a = 4;
    double init_b = 4;
    double d_a    = .125;
    double d_b    = .03125;
    double s      = .005;
    double mean_beta = 12;
    double beta_variation = .05;
    int   beta_random_seed = 1;
    int   i, j, k;
    int   num_iterations = 30000;

    Matrix matrix_a;
    Matrix matrix_b;
    Matrix matrix_a_out;
    Matrix matrix_b_out;
    Matrix matrix_swap_tmp;
    Matrix matrix_beta;

    MatrixSize matrix_size = {64+2, 64+2}; //shape(64+2, 64+2);

    matrix_a.resize(matrix_size);
    matrix_b.resize(matrix_size);
    matrix_a_out.resize(matrix_size);
    matrix_b_out.resize(matrix_size);
    matrix_swap_tmp.resize(matrix_size);
    matrix_beta.resize(matrix_size);

    srand(beta_random_seed);
    for (i = 0; i < matrix_size[ROWS]; ++i) {
	for (j = 0; j < matrix_size[COLS]; ++j) {
	    matrix_a(i,j) = init_a;
	    matrix_b(i,j) = init_b;
	    matrix_beta(i,j) = mean_beta
		+2.0 * (((double)rand() / RAND_MAX) - 0.5) * beta_variation;
	}
    }

    for (k = 0; k < num_iterations; ++k) {
	reaction_diffusion_texture(matrix_size, d_a, d_b, s,
				   matrix_a, matrix_b,
				   matrix_a_out, matrix_b_out, matrix_beta);
	matrix_swap_tmp.reference(matrix_a); matrix_a.reference(matrix_a_out); matrix_a_out.reference(matrix_swap_tmp);
	matrix_swap_tmp.reference(matrix_b); matrix_b.reference(matrix_b_out); matrix_b_out.reference(matrix_swap_tmp);
	// slower than above:
// 	swap(matrix_a, matrix_a_out);
// 	swap(matrix_b, matrix_b_out);
	// slower than swap:
// 	matrix_swap_tmp = matrix_a; matrix_a = matrix_a_out; matrix_a_out = matrix_swap_tmp;
// 	matrix_swap_tmp = matrix_b; matrix_b = matrix_b_out; matrix_b_out = matrix_swap_tmp;
   }

    printf("P2\n%d %d\n255\n", matrix_size[COLS] - 2, matrix_size[ROWS] - 2);
    for (i = 1; i < matrix_size[ROWS] - 1; ++i) {
	for (j = 1; j < matrix_size[COLS] - 1; ++j) {
	    printf("%d ", (int) matrix_a(i,j));
	}
	printf("\n");
    }
#ifdef BZ_DEBUG
    fprintf(stderr, "WARNING BZ_DEBUG IS ENABLED!\n");
#endif
    return 0;
}