CAdatatypes.h

Go to the documentation of this file.

 
#pragma once // Ensures that this file is only included once
             // during compilation
 
#include "CAutils.h"
#include <array>
#include <unordered_map>
#include <iostream>
#include <algorithm> // max_element
#include <utility>   // pair
#include <cmath>     // pow
 
// File path of the output data log
const std::string FILE_PATH = "Data/data.csv";
 
using MajorityCounter = std::unordered_map<int, int>;
 
namespace CAEnums
{
    enum Neighborhood
    {
        VonNeumann,
        Moore
    };
 
    enum Boundary
    {
        Periodic,
        Walled,
        CutOff
    };
 
    enum Rule
    {
        Majority,
        Parity,
        Custom
    };
 
    enum ErrorCode
    {
        CellsAlreadyInitialized = -1,
        CellsAreNull = -2,
        CellsMalloc = -3,
        InvalidCellState = -4,
        InvalidCellStateCondition = -5,
        InvalidRadius = -6,
        InvalidNumStates = -7,
        NeighborhoodCellsMalloc = -8,
        CustomRuleIsNull = -9,
        RadiusLargerThanDimensions = -10
    };
}
 
class BaseCellularAutomata
{
public:
    CAEnums::Boundary boundary_type;         
    int boundary_radius;                     
    CAEnums::Neighborhood neighborhood_type; 
    int num_states;                          
    CAEnums::Rule rule_type;                 
    int axis1_dim;                           
    int axis2_dim;                           
    int axis3_dim;                           
 
    BaseCellularAutomata();
 
    virtual ~BaseCellularAutomata(){};
 
    int setup_neighborhood(CAEnums::Neighborhood neighborhood_type);
 
    int setup_cell_states(int num_states);
 
    int setup_rule(CAEnums::Rule rule_type);
 
    void print_error_status(CAEnums::ErrorCode error);
};
 
template <typename T>
class CellularAutomata : public BaseCellularAutomata
{
private:
    T *vector;        
    T *next_vector;   
    T **matrix;       
    T **next_matrix;  
    T ***tensor;      
    T ***next_tensor; 
    int steps_taken;  
 
    int set_new_cell_state(int *cell_index, int index_size,
                           T *neighborhood_cells, int neighborhood_size,
                           T &new_cell_state, void(custom_rule)(int *, int, T *, int, T &))
    {
        int sum = 0;                         // sum of cells within boundary_radius for Parity rule
        MajorityCounter state_votes_counter; // counter to keep track of votes for Majority rule
        initialize_majority_rule_counter(state_votes_counter, num_states);
 
        switch (rule_type)
        {
        case CAEnums::Custom:
            if (vector == nullptr && matrix == nullptr && tensor == nullptr)
            {
                new_cell_state.state = -1;
                return CAEnums::CellsAreNull;
            }
            if (custom_rule == nullptr)
            {
                new_cell_state.state = -1;
                return CAEnums::CustomRuleIsNull;
            }
            else
            {
                // custom_rule should set the new_cell_state
                custom_rule(cell_index, index_size, neighborhood_cells, neighborhood_size, new_cell_state);
            }
            break;
        case CAEnums::Parity:
            for (int i = 0; i < neighborhood_size; i++)
            {
                // update sum with current cell value
                sum += neighborhood_cells[i].state;
            }
            new_cell_state.state = sum % num_states; // store the parity state as the new state
            break;
        case CAEnums::Majority:
            for (int i = 0; i < neighborhood_size; i++)
            {
                // increment the cell state's number of votes
                auto it = state_votes_counter.find(neighborhood_cells[i].state);
                if (it != state_votes_counter.end())
                {
                    it->second++;
                }
            }
            // find the max_element in the counter based on the pairs' second variable (number of "votes")
            auto max_elem = max_element(state_votes_counter.begin(), state_votes_counter.end(),
                                        less_than_votes);
            new_cell_state.state = max_elem->first; // set the majority state as the new state
            break;
        }
        return 0;
    }
 
    int get_state_from_neighborhood_1d(int *cell_index, int index_size, T &new_cell_state,
                                       void(custom_rule)(int *, int, T *, int, T &))
    {
        int error_code = 0;         // store error code return by other methods
        int i = cell_index[0];      // get i-th index from array
        int neighborhood_index = 0; // keep track of neighborhood array as we iterate through CA cells
 
        // allocate memory and check if operation was successful
        T *neighborhood_cells = malloc_neighborhood_array(index_size);
        if (neighborhood_cells == nullptr)
        {
            return CAEnums::NeighborhoodCellsMalloc;
        }
 
        // VonNeumann and Moore do not differ for 1d (vector) case
        switch (boundary_type)
        {
        case CAEnums::Periodic:
            generate_periodic_neighborhood_1d(cell_index, neighborhood_cells, neighborhood_index);
            error_code = set_new_cell_state(cell_index, index_size, neighborhood_cells, neighborhood_index, new_cell_state, custom_rule);
            break;
        case CAEnums::Walled: // with walled boundaries the edge cells never change
            // check if i is a boundary cell
            if (i == 0 || i == axis1_dim - 1)
            {
                new_cell_state = vector[i]; // keep boundary cell state
                break;
            }
            // fall-through to CutOff case when cell isn't a boundary cell
        case CAEnums::CutOff:
            generate_cutoff_neighborhood_1d(cell_index, neighborhood_cells, neighborhood_index);
            error_code = set_new_cell_state(cell_index, index_size, neighborhood_cells, neighborhood_index, new_cell_state, custom_rule);
            break;
        }
        delete[] neighborhood_cells;
        return error_code;
    }
 
    int get_state_from_neighborhood_2d(int *cell_index, int index_size, T &new_cell_state,
                                       void(custom_rule)(int *, int, T *, int, T &))
    {
        int error_code = 0;         // store error code return by other methods
        int i = cell_index[0];      // get i-th index from array
        int j = cell_index[1];      // get j-th index from array
        int neighborhood_index = 0; // keep track of neighborhood array as we iterate through CA cells
 
        // allocate memory and check if operation was successful
        T *neighborhood_cells = malloc_neighborhood_array(index_size);
        if (neighborhood_cells == nullptr)
        {
            return CAEnums::NeighborhoodCellsMalloc;
        }
 
        switch (boundary_type)
        {
        case CAEnums::Periodic:
            generate_periodic_neighborhood_2d(cell_index, neighborhood_cells, neighborhood_index);
            error_code = set_new_cell_state(cell_index, index_size, neighborhood_cells, neighborhood_index, new_cell_state, custom_rule);
            break;
        case CAEnums::Walled: // with walled boundaries the edge cells never change
            // check if i,j is a boundary cell
            if ((i == 0 || i == axis1_dim - 1) || (j == 0 || j == axis2_dim - 1))
            {
                new_cell_state = matrix[i][j]; // keep boundary cell state
                break;
            }
            // fall-through to CutOff case when cell isn't a boundary cell
        case CAEnums::CutOff:
            generate_cutoff_neighborhood_2d(cell_index, neighborhood_cells, neighborhood_index);
            error_code = set_new_cell_state(cell_index, index_size, neighborhood_cells, neighborhood_index, new_cell_state, custom_rule);
            break;
        }
        delete[] neighborhood_cells;
        return error_code;
    }
 
    int get_state_from_neighborhood_3d(int *cell_index, int index_size, T &new_cell_state,
                                       void(custom_rule)(int *, int, T *, int, T &))
    {
        int error_code = 0;         // store error code return by other methods
        int i = cell_index[0];      // get i-th index from array
        int j = cell_index[1];      // get j-th index from array
        int k = cell_index[2];      // get k-th index from array
        int neighborhood_index = 0; // keep track of neighborhood array as we iterate through CA cells
 
        // allocate memory and check if operation was successful
        T *neighborhood_cells = malloc_neighborhood_array(index_size); // index_size == rank
        if (neighborhood_cells == nullptr)
        {
            return CAEnums::NeighborhoodCellsMalloc;
        }
 
        switch (boundary_type)
        {
        case CAEnums::Periodic:
            generate_periodic_neighborhood_3d(cell_index, neighborhood_cells, neighborhood_index);
            error_code = set_new_cell_state(cell_index, index_size, neighborhood_cells, neighborhood_index, new_cell_state, custom_rule);
            break;
        case CAEnums::Walled:
            // check if i,j,k is a boundary cell
            if ((i == 0 || i == axis1_dim - 1) || (j == 0 || j == axis2_dim - 1) || (k == 0 || j == axis3_dim - 1))
            {
                // with walled boundaries the edge cells never change
                new_cell_state = tensor[i][j][k]; // keep boundary cell state
                break;
            }
            // fall-through to CutOff case when cell isn't a boundary cell
        case CAEnums::CutOff:
            generate_cutoff_neighborhood_3d(cell_index, neighborhood_cells, neighborhood_index);
            error_code = set_new_cell_state(cell_index, index_size, neighborhood_cells, neighborhood_index, new_cell_state, custom_rule);
            break;
        }
        delete[] neighborhood_cells;
        return error_code;
    }
 
    void generate_periodic_neighborhood_1d(int *cell_index, T *neighborhood_cells, int &neighborhood_index)
    {
        int i = cell_index[0]; // get i-th index from array
        int periodic_i;        // axis1_dim index used by Periodic boundary type
 
        for (int di = -boundary_radius; di <= boundary_radius; di++)
        {
            periodic_i = get_periodic_index(i, di, axis1_dim);
            // add state to flattened array
            neighborhood_cells[neighborhood_index] = vector[periodic_i];
            neighborhood_index++;
        }
    }
 
    void generate_cutoff_neighborhood_1d(int *cell_index, T *neighborhood_cells, int &neighborhood_index)
    {
        int i = cell_index[0]; // get i-th index from array
        int neighbor_i;        // neighboring cell's i-th index
 
        for (int di = -boundary_radius; di <= boundary_radius; di++)
        {
            neighbor_i = di + i;
            // exclude cells that are out of bounds
            if (neighbor_i <= 0 || neighbor_i >= axis1_dim)
            {
                // outside bounds; don't include cell state in the sum/counter
                continue;
            }
            // add state to flattened array
            neighborhood_cells[neighborhood_index] = vector[di];
            neighborhood_index++;
        }
    }
 
    void generate_periodic_neighborhood_2d(int *cell_index, T *neighborhood_cells, int &neighborhood_index)
    {
        int i = cell_index[0]; // get i-th index from array
        int j = cell_index[1]; // get j-th index from array
        int periodic_i;        // axis1_dim index used by Periodic boundary type
        int periodic_j;        // axis2_dim index used by Periodic boundary type
 
        if (neighborhood_type == CAEnums::VonNeumann)
        {
            for (int di = -boundary_radius; di <= boundary_radius; di++)
            {
                periodic_i = get_periodic_index(i, di, axis1_dim);
                for (int dj = -boundary_radius; dj <= boundary_radius; dj++)
                {
                    // exclude diagonal cells from neighborhood when VonNeumann is selected
                    if (is_diagonal_neighboring_cell_2d(di, dj))
                    {
                        // current di,dj cell is a diagonal neighbor so exclude it
                        continue;
                    }
                    periodic_j = get_periodic_index(j, dj, axis2_dim);
                    // add state to flattened array
                    neighborhood_cells[neighborhood_index] = matrix[periodic_i][periodic_j];
                    neighborhood_index++;
                }
            }
        }
        else // Moore neighborhood
        {
            for (int di = -boundary_radius; di <= boundary_radius; di++)
            {
                periodic_i = get_periodic_index(i, di, axis1_dim);
                for (int dj = -boundary_radius; dj <= boundary_radius; dj++)
                {
                    periodic_j = get_periodic_index(j, dj, axis2_dim);
                    // add state to flattened array
                    neighborhood_cells[neighborhood_index] = matrix[periodic_i][periodic_j];
                    neighborhood_index++;
                }
            }
        }
    }
 
    void generate_cutoff_neighborhood_2d(int *cell_index, T *neighborhood_cells, int &neighborhood_index)
    {
        int i = cell_index[0]; // get i-th index from array
        int j = cell_index[1]; // get j-th index from array
        int neighbor_i;        // neighboring cell's i-th index
        int neighbor_j;        // neighboring cell's j-th index
 
        if (neighborhood_type == CAEnums::VonNeumann)
        {
            for (int di = -boundary_radius; di <= boundary_radius; di++)
            {
                neighbor_i = di + i;
                for (int dj = -boundary_radius; dj <= boundary_radius; dj++)
                {
                    neighbor_j = dj + j;
                    // exclude diagonal cells from neighborhood when VonNeumann is selected
                    if (is_diagonal_neighboring_cell_2d(di, dj))
                    {
                        // current di,dj cell is a diagonal neighbor so exclude it
                        continue;
                    }
                    // exclude cells that are out of bounds
                    if ((neighbor_i <= 0 || neighbor_i >= axis1_dim) || (neighbor_j <= 0 || neighbor_j >= axis2_dim))
                    {
                        // outside bounds; don't include cell state in the sum/counter
                        continue;
                    }
                    // add state to flattened array
                    neighborhood_cells[neighborhood_index] = matrix[neighbor_i][neighbor_j];
                    neighborhood_index++;
                }
            }
        }
        else // Moore neighborhood
        {
            for (int di = -boundary_radius; di <= boundary_radius; di++)
            {
                neighbor_i = di + i;
                for (int dj = -boundary_radius; dj <= boundary_radius; dj++)
                {
                    neighbor_j = dj + j;
                    // exclude cells that are out of bounds
                    if ((neighbor_i <= 0 || neighbor_i >= axis1_dim) || (neighbor_j <= 0 || neighbor_j >= axis2_dim))
                    {
                        // outside bounds; don't include cell state in the sum/counter
                        continue;
                    }
                    // add state to flattened array
                    neighborhood_cells[neighborhood_index] = matrix[neighbor_i][neighbor_j];
                    neighborhood_index++;
                }
            }
        }
    }
 
    void generate_periodic_neighborhood_3d(int *cell_index, T *neighborhood_cells, int &neighborhood_index)
    {
        int i = cell_index[0]; // get i-th index from array
        int j = cell_index[1]; // get j-th index from array
        int k = cell_index[2]; // get k-th index from array
        int periodic_i;        // axis1_dim index used by Periodic boundary type
        int periodic_j;        // axis2_dim index used by Periodic boundary type
        int periodic_k;        // axis3_dim index used by Periodic boundary type
 
        if (neighborhood_type == CAEnums::VonNeumann)
        {
            for (int di = -boundary_radius; di <= boundary_radius; di++)
            {
                periodic_i = get_periodic_index(i, di, axis1_dim);
                for (int dj = -boundary_radius; dj <= boundary_radius; dj++)
                {
                    periodic_j = get_periodic_index(j, dj, axis2_dim);
                    for (int dk = -boundary_radius; dk <= boundary_radius; dk++)
                    {
                        // exclude diagonal cells from neighborhood when VonNeumann is selected
                        if (is_diagonal_neighboring_cell_3d(di, dj, dk))
                        {
                            // current di,dj,dk cell is a diagonal neighbor so exclude it
                            continue;
                        }
                        periodic_k = get_periodic_index(k, dk, axis3_dim);
                        // add state to flattened array
                        neighborhood_cells[neighborhood_index] = tensor[periodic_i][periodic_j][periodic_k];
                        neighborhood_index++;
                    }
                }
            }
        }
        else // Moore neighborhood
        {
            for (int di = -boundary_radius; di <= boundary_radius; di++)
            {
                periodic_i = get_periodic_index(i, di, axis1_dim);
                for (int dj = -boundary_radius; dj <= boundary_radius; dj++)
                {
                    periodic_j = get_periodic_index(j, dj, axis2_dim);
                    for (int dk = -boundary_radius; dk <= boundary_radius; dk++)
                    {
                        periodic_k = get_periodic_index(k, dk, axis3_dim);
                        // add state to flattened array
                        neighborhood_cells[neighborhood_index] = tensor[periodic_i][periodic_j][periodic_k];
                        neighborhood_index++;
                    }
                }
            }
        }
    }
 
    void generate_cutoff_neighborhood_3d(int *cell_index, T *neighborhood_cells, int &neighborhood_index)
    {
        int i = cell_index[0]; // get i-th index from array
        int j = cell_index[1]; // get j-th index from array
        int k = cell_index[2]; // get k-th index from array
        int neighbor_i;        // neighboring cell's i-th index
        int neighbor_j;        // neighboring cell's j-th index
        int neighbor_k;        // neighboring cell's k-th index
 
        if (neighborhood_type == CAEnums::VonNeumann)
        {
            for (int di = -boundary_radius; di <= boundary_radius; di++)
            {
                neighbor_i = di + i;
                for (int dj = -boundary_radius; dj <= boundary_radius; dj++)
                {
                    neighbor_j = dj + j;
                    for (int dk = -boundary_radius; dk <= boundary_radius; dk++)
                    {
                        neighbor_k = dk + k;
                        // exclude diagonal cells from neighborhood when VonNeumann is selected
                        if (is_diagonal_neighboring_cell_3d(di, dj, dk))
                        {
                            // current di,dj,dk cell is a diagonal neighbor so exclude it
                            continue;
                        }
                        // exclude cells that are out of bounds
                        if ((neighbor_i <= 0 || neighbor_i >= axis1_dim) ||
                            (neighbor_j <= 0 || neighbor_j >= axis2_dim) ||
                            (neighbor_k <= 0 || neighbor_k >= axis3_dim))
                        {
                            // outside bounds; don't include cell state in the sum/counter
                            continue;
                        }
                        // add state to flattened array
                        neighborhood_cells[neighborhood_index] = tensor[neighbor_i][neighbor_j][neighbor_k];
                        neighborhood_index++;
                    }
                }
            }
        }
        else // Moore neighborhood
        {
            for (int di = -boundary_radius; di <= boundary_radius; di++)
            {
                neighbor_i = di + i;
                for (int dj = -boundary_radius; dj <= boundary_radius; dj++)
                {
                    neighbor_j = dj + j;
                    for (int dk = -boundary_radius; dk <= boundary_radius; dk++)
                    {
                        neighbor_k = dk + k;
                        // exclude cells that are out of bounds
                        if ((neighbor_i <= 0 || neighbor_i >= axis1_dim) ||
                            (neighbor_j <= 0 || neighbor_j >= axis2_dim) ||
                            (neighbor_k <= 0 || neighbor_k >= axis3_dim))
                        {
                            // outside bounds; don't include cell state in the sum/counter
                            continue;
                        }
                        // add state to flattened array
                        neighborhood_cells[neighborhood_index] = tensor[neighbor_i][neighbor_j][neighbor_k];
                        neighborhood_index++;
                    }
                }
            }
        }
    }
 
    T *malloc_neighborhood_array(int rank)
    {
        /*
         * Generate a flatten array of the cell's neighborhood.
         * The neighborhood array can then be utilized for Majority, Parity, or Custom rule
         */
        int max_neighborhood_size = get_neighborhood_size(rank, boundary_radius, neighborhood_type);
        // allocate memory and check if operation was successful
        T *neighborhood_cells = new (std::nothrow) T[max_neighborhood_size];
        return neighborhood_cells;
    }
 
    int append_log()
    {
        std::ofstream file;
        file.open(FILE_PATH, std::ios::app);
 
        if (vector != nullptr)
        {
            for (int i = 0; i < axis1_dim; i++)
            {
                file << vector[i].state << ",";
            }
        }
        else if (matrix != nullptr)
        {
            for (int j = 0; j < axis1_dim; j++)
            {
                for (int k = 0; k < axis2_dim; k++)
                {
                    file << matrix[j][k].state << ",";
                }
            }
        }
        else if (tensor != nullptr)
        {
            for (int i = 0; i < axis1_dim; i++)
            {
                for (int j = 0; j < axis2_dim; j++)
                {
                    for (int k = 0; k < axis3_dim; k++)
                    {
                        file << tensor[i][j][k].state << ",";
                    }
                }
            }
        }
        else
        {
            return CAEnums::CellsAreNull;
        }
        file << "\n";
        file.close();
        return (0);
    }
 
    int create_log(void)
    {
        std::ofstream file;
        file.open(FILE_PATH, std::ios::trunc | std::ios::out); // If it is pre-existing content, it should be erased
 
        file << this->num_states << ",\n"; // First Line: Number of the states.
        // Rule and function
        file << this->axis1_dim << ","
             << this->axis2_dim << ","
             << this->axis3_dim << ",\n"; // Second Line: Dimensions of each axis.
        file.close();
        return 0;
    }
 
public:
    T *get_vector()
    {
        return vector;
    }
 
    T *get_next_vector()
    {
        return next_vector;
    }
 
    T **get_matrix()
    {
        return matrix;
    }
 
    T **get_next_matrix()
    {
        return next_matrix;
    }
 
    T ***get_tensor()
    {
        return tensor;
    }
 
    T ***get_next_tensor()
    {
        return next_tensor;
    }
 
    int setup_boundary(CAEnums::Boundary bound_type, int radius)
    {
        if (radius <= 0)
        {
            return CAEnums::InvalidRadius;
        }
        if (vector != nullptr)
        {
            if (radius > axis1_dim / 2)
            {
                return CAEnums::RadiusLargerThanDimensions;
            }
        }
        else if (matrix != nullptr)
        {
            if (radius > axis1_dim / 2 || radius > axis2_dim / 2)
            {
                return CAEnums::RadiusLargerThanDimensions;
            }
        }
        else if (tensor != nullptr)
        {
            if (radius > axis2_dim / 2 || radius > axis3_dim / 2)
            {
                return CAEnums::RadiusLargerThanDimensions;
            }
        }
 
        this->boundary_type = bound_type;
        this->boundary_radius = radius;
        return 0;
    }
 
    CellularAutomata() : BaseCellularAutomata()
    {
        steps_taken = 0;
        vector = nullptr;
        next_vector = nullptr;
        matrix = nullptr;
        next_matrix = nullptr;
        tensor = nullptr;
        next_tensor = nullptr;
    }
 
    ~CellularAutomata()
    {
        // Free each sub-array
        if (vector != nullptr)
        {
            // Free the array of pointers
            delete[] vector;
            delete[] next_vector;
        }
 
        if (matrix != nullptr)
        {
            // Free each sub-array
            for (int i = 0; i < axis1_dim; ++i)
            {
                delete[] matrix[i];
                delete[] next_matrix[i];
            }
            // Free the array of pointers
            delete[] matrix;
            delete[] next_matrix;
        }
 
        if (tensor != nullptr)
        {
            // Free each sub-array
            for (int i = 0; i < axis1_dim; ++i)
            {
                for (int j = 0; j < axis2_dim; j++)
                {
                    delete[] tensor[i][j];
                    delete[] next_tensor[i][j];
                }
 
                delete[] tensor[i];
                delete[] next_tensor[i];
            }
 
            // Free the array of pointers
            delete[] tensor;
            delete[] next_tensor;
        }
    }
 
    int setup_dimensions_1d(int axis1_dim, int fill_value = 0)
    {
        if (vector != nullptr)
        {
            return CAEnums::CellsAlreadyInitialized;
        }
 
        this->axis1_dim = axis1_dim;
        vector = new (std::nothrow) T[axis1_dim];
        next_vector = new (std::nothrow) T[axis1_dim];
 
        if (vector == nullptr || next_vector == nullptr)
        {
            return CAEnums::CellsMalloc;
        }
 
        // initialize vector filled with zeros
        for (int j = 0; j < axis1_dim; j++)
        {
            vector[j].state = fill_value;
            next_vector[j].state = fill_value;
        }
 
        create_log();
        return 0;
    }
 
    int setup_dimensions_2d(int axis1_dim, int axis2_dim, int fill_value = 0)
    {
        if (matrix != nullptr)
        {
            return CAEnums::CellsAlreadyInitialized;
        }
 
        this->axis1_dim = axis1_dim;
        this->axis2_dim = axis2_dim;
        matrix = new (std::nothrow) T *[axis1_dim];
        next_matrix = new (std::nothrow) T *[axis1_dim];
 
        if (matrix == nullptr || next_matrix == nullptr)
        {
            return CAEnums::CellsMalloc;
        }
 
        for (int i = 0; i < axis1_dim; i++)
        {
            matrix[i] = new (std::nothrow) T[axis2_dim];
            next_matrix[i] = new (std::nothrow) T[axis2_dim];
        }
 
        // initialize matrix filled with zeros
        for (int j = 0; j < axis1_dim; j++)
        {
            for (int k = 0; k < axis2_dim; k++)
            {
                matrix[j][k].state = fill_value;
                next_matrix[j][k].state = fill_value;
            }
        }
 
        create_log();
        return 0;
    }
 
    int setup_dimensions_3d(int axis1_dim, int axis2_dim, int axis3_dim, int fill_value = 0)
    {
        if (tensor != nullptr)
        {
            return CAEnums::CellsAlreadyInitialized;
        }
 
        this->axis1_dim = axis1_dim;
        this->axis2_dim = axis2_dim;
        this->axis3_dim = axis3_dim;
        tensor = new (std::nothrow) T **[axis1_dim];
        next_tensor = new (std::nothrow) T **[axis1_dim];
 
        if (tensor == nullptr || next_tensor == nullptr)
        {
            return CAEnums::CellsMalloc;
        }
 
        for (int i = 0; i < axis1_dim; i++)
        {
            tensor[i] = new (std::nothrow) T *[axis2_dim];
            next_tensor[i] = new (std::nothrow) T *[axis2_dim];
 
            for (int j = 0; j < axis2_dim; j++)
            {
                tensor[i][j] = new (std::nothrow) T[axis3_dim];
                next_tensor[i][j] = new (std::nothrow) T[axis3_dim];
            }
        }
 
        // initialize matrix filled with zeros
        for (int i = 0; i < axis1_dim; i++)
        {
            for (int j = 0; j < axis2_dim; j++)
            {
                for (int k = 0; k < axis3_dim; k++)
                {
                    tensor[i][j][k].state = fill_value;
                    next_tensor[i][j][k].state = fill_value;
                }
            }
        }
 
        create_log();
        return 0;
    }
 
    int init_condition(int x_state, double prob)
    {
        if (!(x_state < num_states))
        {
            return CAEnums::InvalidCellStateCondition;
        }
 
        srand(time(NULL));
        double random_cell_state;
 
        if (vector != nullptr)
        {
            for (int i = 0; i < axis1_dim; i++)
            {
                random_cell_state = (double)rand() / RAND_MAX;
                if (random_cell_state < prob)
                {
                    vector[i].state = x_state;
                }
            }
        }
        else if (matrix != nullptr)
        {
            for (int j = 0; j < axis1_dim; j++)
            {
                for (int k = 0; k < axis2_dim; k++)
                {
                    random_cell_state = (double)rand() / RAND_MAX;
                    if (random_cell_state < prob)
                    {
                        matrix[j][k].state = x_state;
                    }
                }
            }
        }
        else if (tensor != nullptr)
        {
            for (int i = 0; i < axis1_dim; i++)
            {
                for (int j = 0; j < axis2_dim; j++)
                {
                    for (int k = 0; k < axis3_dim; k++)
                    {
                        random_cell_state = (double)rand() / RAND_MAX;
                        if (random_cell_state < prob)
                        {
                            tensor[i][j][k].state = x_state;
                        }
                    }
                }
            }
        }
        else
        {
            return CAEnums::CellsAreNull;
        }
        return 0;
    }
 
    int step(void(custom_rule)(int *, int, T *, int, T &))
    {
        int error_code = 0; // store error code return by other methods
        T new_cell_state;   // stores the cell's new state
        T empty_cell_state; // cell state for zeroing out old states
        int index_size;     // number of indices required to address the cell
 
        if (vector != nullptr)
        {
            // store the main cell's index in cell_index for custom rule type
            index_size = 1;
#ifdef ENABLE_OMP
#pragma omp parallel for firstprivate(error_code) private(new_cell_state)
#endif
            for (int i = 0; i < axis1_dim; i++)
            {
                int cell_index[index_size] = {i};
                new_cell_state = vector[i];
                error_code = get_state_from_neighborhood_1d(cell_index, index_size, new_cell_state, custom_rule);
                if (error_code < 0)
                {
#ifdef ENABLE_OMP
#pragma omp cancel for
#endif
#ifndef ENABLE_OMP
                    /*
                     * Return error code when omp is not enabled
                     * otherwise the above pragma will exit out of the for loop(s)
                     */
                    return error_code;
#endif
                }
                /*
                 * The update cell if new_cell_state is no empty_state.
                 * Avoids overwriting the motion of cells.
                 */
                if (new_cell_state != empty_cell_state)
                {
                    next_vector[cell_index[0]] = new_cell_state;
                }
            }
            // store next cell state to the current cell state for the next time step
            swap_states<T>(vector, next_vector, axis1_dim);
        }
        else if (matrix != nullptr)
        {
            // store the main cell's index in cell_index for custom rule type
            index_size = 2;
#ifdef ENABLE_OMP
#pragma omp parallel for firstprivate(error_code) private(new_cell_state)
#endif
            for (int i = 0; i < axis1_dim; i++)
            {
                for (int j = 0; j < axis2_dim; j++)
                {
                    int cell_index[index_size] = {i, j};
                    new_cell_state = matrix[i][j];
                    error_code = get_state_from_neighborhood_2d(cell_index, index_size, new_cell_state, custom_rule);
                    if (error_code < 0)
                    {
#ifdef ENABLE_OMP
#pragma omp cancel for
#endif
#ifndef ENABLE_OMP
                        /*
                         * Return error code when omp is not enabled
                         * otherwise the above pragma will exit out of the for loop(s)
                         */
                        return error_code;
#endif
                    }
                    /*
                     * The update cell if new_cell_state is no empty_state.
                     * Avoids overwriting the motion of cells.
                     */
                    if (new_cell_state != empty_cell_state)
                    {
                        next_matrix[cell_index[0]][cell_index[1]] = new_cell_state;
                    }
                }
            }
            // store next cell state to the current cell state for the next time step
            swap_states<T>(matrix, next_matrix, axis1_dim, axis2_dim);
        }
        else if (tensor != nullptr)
        {
            index_size = 3;
#ifdef ENABLE_OMP
#pragma omp parallel for firstprivate(error_code) private(new_cell_state)
#endif
            for (int i = 0; i < axis1_dim; i++)
            {
                for (int j = 0; j < axis2_dim; j++)
                {
                    for (int k = 0; k < axis3_dim; k++)
                    {
                        int cell_index[index_size] = {i, j, k};
                        new_cell_state = tensor[i][j][k];
                        error_code = get_state_from_neighborhood_3d(cell_index, index_size, new_cell_state, custom_rule);
                        if (error_code < 0)
                        {
#ifdef ENABLE_OMP
#pragma omp cancel for
#endif
#ifndef ENABLE_OMP
                            /*
                             * Return error code when omp is not enabled
                             * otherwise the above pragma will exit out of the for loop(s)
                             */
                            return error_code;
#endif
                        }
                        /*
                         * The update cell if new_cell_state is no empty_state.
                         * Avoids overwriting the motion of cells.
                         */
                        if (new_cell_state != empty_cell_state)
                        {
                            next_tensor[cell_index[0]][cell_index[1]][cell_index[2]] = new_cell_state;
                        }
                    }
                }
            }
            // store next cell state to the current cell state for the next time step
            swap_states<T>(tensor, next_tensor, axis1_dim, axis2_dim, axis3_dim);
        }
        else
        {
            return CAEnums::CellsAreNull;
        }
 
        steps_taken++;
        // Appending the step to the file log
        error_code = append_log();
        return error_code;
    }
 
    int step()
    {
        return step(nullptr); // return step(func) error code
    }
 
    int print_grid()
    {
        if (vector != nullptr)
        {
            for (int i = 0; i < axis1_dim; i++)
            {
                std::cout << vector[i].state << " ";
            }
            std::cout << std::endl;
        }
        else if (matrix != nullptr)
        {
            for (int j = 0; j < axis1_dim; j++)
            {
                for (int k = 0; k < axis2_dim; k++)
                {
                    std::cout << matrix[j][k].state << " ";
                }
                std::cout << std::endl;
            }
        }
        else if (tensor != nullptr)
        {
            for (int i = 0; i < axis1_dim; i++)
            {
                std::cout << "Printing " << i << "'th slice of Tensor" << std::endl;
                for (int j = 0; j < axis2_dim; j++)
                {
                    for (int k = 0; k < axis3_dim; k++)
                    {
                        std::cout << tensor[i][j][k].state << " ";
                    }
                    std::cout << std::endl;
                }
            }
        }
        else
        {
            return CAEnums::CellsAreNull;
        }
        return 0;
    }
 
    static int get_neighborhood_size(int rank, int radius, CAEnums::Neighborhood neighborhood_type)
    {
        if (neighborhood_type == CAEnums::VonNeumann)
        {
            return (2 * rank * radius) + 1; // +1 to include cell of interest
        }
        else // CAEnums::Moore neighborhood
        {
            return pow((2 * radius + 1), rank);
        }
    }
 
    static void initialize_majority_rule_counter(MajorityCounter &counter, int num_states)
    {
        // sets the counter for every cell state type to 0
        for (int j = 0; j < num_states; j++)
        {
            counter.insert(std::make_pair(j, 0));
        }
    }
};
 
// Declare specialized int methods. Methods are defined in cellularautomata.cpp
 
template <>
int CellularAutomata<int>::setup_dimensions_1d(int axis1_dim, int fill_value);
 
template <>
int CellularAutomata<int>::setup_dimensions_2d(int axis1_dim, int axis2_dim, int fill_value);
 
template <>
int CellularAutomata<int>::setup_dimensions_3d(int axis1_dim, int axis2_dim, int axis3_dim, int fill_value);
 
template <>
int CellularAutomata<int>::init_condition(int x_state, double prob);
 
template <>
int CellularAutomata<int>::set_new_cell_state(int *cell_index, int index_size,
                                              int *neighborhood_cells, int neighborhood_size,
                                              int &new_cell_state, void(custom_rule)(int *, int, int *, int, int &));
 
template <>
int CellularAutomata<int>::print_grid();
 
template <>
int CellularAutomata<int>::append_log();
 
template <>
int CellularAutomata<int>::step(void(custom_rule)(int *, int, int *, int, int &));