/*
 * simulate a cellular automaton with periodic boundaries (torus-like)
 * OpenMP version
 *
 * (c) 2024 Dorian Stoll (C++ / Eigen port)
 * (c) 2016 Steffen Christgau (C99 port, modularization, parallelization)
 * (c) 1996,1997 Peter Sanders, Ingo Boesnach (original source)
 *
 * command line arguments:
 * #1: Number of lines
 * #2: Number of iterations to be simulated
 *
 */

/*
 * Eigen defaults to using column-major order for its container types.
 *
 * However, using it like that would make the lines parameter correspond to Eigens .cols() function
 * instead of the more fitting .rows(). To make the code more obvious, we change the default.
 */
#define EIGEN_DEFAULT_TO_ROW_MAJOR

#include <Eigen/Eigen>

#include <cassert>
#include <cstdint>
#include <cstdlib>
#include <iostream>
#include <stdexcept>

class Automat {
private:
	/* "ADT" State */
	using cell_state_t = uint8_t;

	using Field = Eigen::Matrix<cell_state_t, Eigen::Dynamic, Eigen::Dynamic>;

private:
	/* horizontal size of the configuration */
	static constexpr Eigen::Index XSIZE = 1024;
	static constexpr Eigen::Index LINE_SIZE = XSIZE + 2;

	/* annealing rule from ChoDro96 page 34
	 * the table is used to map the number of nonzero
	 * states in the neighborhood to the new state
	 */
	static constexpr std::array<cell_state_t, 10> ANNEAL {0, 0, 0, 0, 1, 0, 1, 1, 1, 1};

private:
	int m_lines;

	Field m_from {};
	Field m_to {};

public:
	Automat(const int lines) : m_lines {lines}
	{
		if (lines <= 0)
			throw std::invalid_argument {"There must be at least one line"};

		m_from.resize(lines + 2, LINE_SIZE);
		m_to.resize(lines + 2, LINE_SIZE);
	}

	void run(const int its)
	{
		this->init_config();

		for (int i = 0; i < its; i++) {
			this->boundary();
			this->simulate();

			std::swap(m_from, m_to);
		}
	}

	void report()
	{
		auto block = m_from(Eigen::seqN(1, m_lines), Eigen::all);

		for (Eigen::Index y = 0; y < block.rows(); y++) {
			block(y, 0) = 0;
			block(y, block.cols() - 1) = 0;
		}

		const auto start = block.data();
		const auto end = start + (block.size() * sizeof(*start));

		std::copy(start, end, std::ostreambuf_iterator(std::cout));
	}

private:
	/* random starting configuration */
	void init_config()
	{
		srand(424243);

		for (Eigen::Index y = 1; y <= m_to.rows() - 2; y++) {
			for (Eigen::Index x = 1; x <= m_from.cols() - 2; x++)
				m_from(y, x) = rand() % 2;
		}
	}

	/* treat torus like boundary conditions */
	void boundary()
	{
#pragma omp sections
		{
#pragma omp section
#pragma omp parallel for
			for (Eigen::Index y = 0; y < m_from.rows(); y++) {
				/* copy rightmost column to the buffer column 0 */
				m_from(y, 0) = m_from(y, m_from.cols() - 2);

				/* copy leftmost column to the buffer column XSIZE + 1 */
				m_from(y, m_from.cols() - 1) = m_from(y, 1);
			}

#pragma omp section
#pragma omp parallel for
			for (Eigen::Index x = 0; x < m_from.cols(); x++) {
				/* copy bottommost row to buffer row 0 */
				m_from(0, x) = m_from(m_from.rows() - 2, x);

				/* copy topmost row to buffer row lines + 1 */
				m_from(m_from.rows() - 1, x) = m_from(1, x);
			}
		}
	}

	/* make one simulation iteration with lines lines.
	 * old configuration is in from, new one is written to to.
	 */
	void simulate()
	{
#pragma omp parallel for
		for (Eigen::Index y = 1; y <= m_from.rows() - 2; y++) {
			for (Eigen::Index x = 1; x <= m_from.cols() - 2; x++)
				m_to(y, x) = transition(m_from, x, y);
		}
	}

	/* a: reference to matrix; x,y: coordinates; result: n-th element of anneal,
	   where n is the number of neighbors */
	cell_state_t transition(const Field &a, const Eigen::Index x, const Eigen::Index y) const
	{
		return ANNEAL[a(Eigen::seq(y - 1, y + 1), Eigen::seq(x - 1, x + 1)).sum()];
	}
};

int main(int argc, char **argv)
{
	assert(argc == 3);

	int lines = atoi(argv[1]);
	int its = atoi(argv[2]);

	Automat ca {lines};

	ca.run(its);
	ca.report();

	return EXIT_SUCCESS;
}