Finding quickly how many neighbors you have...

Back to articles

…inside matrices (gotcha!).

Created on 2024/11/22

Today I had a course of computer programming (it’s just C course actually). There was a strange exercise we had to do: from a matrix, create a new one telling for each cell how many neighbors which share the same value they have.

Easy at a quick glance, but harder when thinking about edge cases.

First, I took a paper, because it’s always easier for me to visualize by writing and drawing.

There are four edge cases relating to the analyzed cell :

  1. It is on the first line.
  2. It is on the first column.
  3. It is on the last line.
  4. It is on the last column.

Four more edge cases happen when you think about corners.

So I used the rough way : checking each edge case, summing the boolean results, then removing one, since the analyzed cell cannot be its own neighbor.

Using C

I defined the algorithm inside a function declared as this : size_t countNeighbors(int **matrix, size_t m, size_t n, size_t i, size_t j).

m corresponds to rows and n to columns. I could have created a struct, but I already had done this here.

size_t countNeighbors(int **matrix, size_t m, size_t n, size_t i, size_t j)
{
  const int a = matrix[i][j];

  const size_t up = (i != 0 && matrix[i-1][j] == a);
  const size_t down = (i != m-1 && matrix[i+1][j] == a);
  const size_t left = (j != 0 && matrix[i][j-1] == a);
  const size_t right = (j != n-1 && matrix[i][j+1] == a);

  const size_t up_left = (i != 0 && j != 0 && matrix[i-1][j-1] == a);
  const size_t up_right = (i != 0 && j != n-1 && matrix[i-1][j+1] == a);
  const size_t down_left = (i != m-1 && j != 0 && matrix[i+1][j-1] == a);
  const size_t down_right = (i != m-1 && j != n-1 && matrix[i+1][j+1] == a);

  return up + down + left + right + up_left + up_right + down_left + down_right;
}

It’s definitely ugly, but it works, and is somewhat readable.

Using C++

To do so I created a MatrixInt class, since I didn’t want to bother with templates to simplify the algorithm. Later we will see the generic version, making sure the Matrix only accepts numeric types.

Instead an int **, I’ll be using a std::vector to not worry about manual memory management.

class MatrixInt {
  std::vector<std::vector<int>> m_elems;
  size_t m_rows, m_columns;
};

I’ll be also using std::optional<int> to safely get elements from the matrix, which in my opinion will make the code more readable than try-catching exceptions or spamming if statements. This means that all the following pieces of code will work starting from C++17, and some of them starting from C++20.

class MatrixInt {
  MatrixInt(size_t m, size_t n) : m_elems(m, std::vector<int>(n)), m_rows{m}, m_columns{n} {}
  // Hiding operator methods...

  std::optional<int> at(int i, int j) const
  {
    if (0 <= i && i < m_rows && 0 <= j && j < m_columns) {
      return m_elems[i][j];
    }
    return std::nullopt;
  }

  size_t countNeighborsAt(int i, int j) const
  {
    size_t result = 0;
    const int a = m_elems[i][j];

    // I didn't know this syntax before coding this program.
    for (int k : {-1, 0, 1}) {
      for (int l : {-1, 0, 1}) {
        result += (this->at(i+k, j+l) == a);
      }
    }
    return result - 1;
  }

  MatrixInt neighbors() const
  {
    MatrixInt result(m_rows, m_columns);
    for (int i = 0; i < m_rows; ++i) {
      for (int j = 0; j < m_columns; ++j) {
        result.m_elems[i][j] = this->countNeighborsAt(i, j);
      }
    }
    return result;
  }
};

You can now use the double for loop going from -1 to 1. The only concession here is that I have to use ints instead of size_t, because the at method must accept negative parameters.

Why do you put your { in a separate line only when defining a function or a method?

I’m following the coding style of suckless (most of it at least). It emphasizes the fact that you cannot nest functions in standard C and C++.

But hey, wouldn’t it be fun if we were using a more modern version of C++ (like C++20) and making the last method more pythonic?

#include <ranges>

class MatrixInt {
  // ...
  MatrixInt neighbors() const
  {
    #define ifor(instructions) for (size_t instructions)
    #define in :

    auto range = std::views::iota;
    MatrixInt result(m_rows, m_columns);
    ifor (i in range(0ull, m_rows)) {
      ifor (j in range(0ull, m_columns)) {
        result.m_elems[i][j] = this->countNeighborsAt(i, j);
      }
    }
    return result;
  }
};

That’s cursed…

Yeah, definitely. Please don’t do this in production code.

Can’t we write range(0, m_rows)?

If you do so, you’ll get an unclear error about substitution failure.

matrix.cpp:58:16: error: no matching function for call to object of type '__iota::__fn'
   58 |     ifor (i in range(0, m_rows)) {
      |                ^~~~~
matrix.cpp:53:44: note: expanded from macro 'ifor'
   53 |     #define ifor(instructions) for (size_t instructions)
      |                                            ^~~~~~~~~~~~
// shortened for the sake of readability
1 error generated.

Because 0 is an int literal, while m_rows is a size_t.

Using C++ and templates

Genericity

Let’s make our Matrix generic. To do so, you add template<typename T> before the class you want to be generic, and add <T> to generic inner types (like std::vector<T> and std::optional<T>).

template<typename T>
class Matrix
{
  std::vector<std::vector<T>> m_elems;
  size_t m_rows, m_columns;

public:
  Matrix(size_t m, size_t n)
    : m_elems(m, std::vector<T>(n)), m_rows{m}, m_columns{n} {}
  // Hiding operator methods...

  std::optional<T> at(int i, int j) const
  {
    if (0 <= i && i < m_rows && 0 <= j && j < m_columns) {
      return m_elems[i][j];
    }
    return std::nullopt;
  }

  size_t countNeighborsAt(int i, int j) const
  {
    size_t result = 0;
    const T a = m_elems[i][j];
    for (int k : {-1, 0, 1}) {
      for (int l : {-1, 0, 1}) {
        result += (this->at(i+k, j+l) == a);
      }
    }
    return result - 1;
  }

  // This is the tricky part
  Matrix<size_t> neighbors() const
  {
    Matrix<size_t> result(m_rows, m_columns);
    for (size_t i = 0; i < m_rows; ++i) {
      for (size_t j = 0; j < m_columns; ++j) {
        // Oddly, you cannot access result.m_elems, so operators are important here
        result[i][j] = this->countNeighborsAt(i, j);
      }
    }
    return result;
  }
};

Fixing construction helper

We have to ensure that T is a numeric type. Currently, in our program, you’d get an error (or worse, only a warning) if you create a Matrix<std::string> or Matrix<char>. This is because I’m using C varargs to construct matrices, and the second parameter of va_arg will be interpreted as an int if it can. So no chars.

Here is the monster:

#include <cstdarg>

template<typename T>
Matrix<T> createMatrix(size_t m, size_t n, ...)
{
  va_list args;
  va_start(args, n);

  Matrix<T> result(m, n);
  for (size_t i = 0; i < m; ++i) {
    for (size_t j = 0; j < n; ++j) {
      result[i][j] = va_arg(args, T);
    }
  }

  va_end(args);
  return result;
}

At the beginning I made this function in C. It isn’t type-safe, and doesn’t check the number of parameters.

Let’s fix that by using std::initializer_list and its size() method. It doesn’t provide an index operator though, so we’ll use an iterator instead.

#include <cassert>
#include <initializer_list>

template<typename T>
Matrix<T> createMatrix(size_t m, size_t n, std::initializer_list<T> list)
{
  assert(list.size() == m * n);

  Matrix<T> result(m, n);
  auto it = list.begin();

  for (size_t i = 0; i < m; ++i) {
    for (size_t j = 0; j < n; ++j) {
      result[i][j] = *it++;
    }
  }

  return result;
}

Type-checking

The easiest way is to add the type-checking inside the constructor.

#include <type_traits>

template<typename T>
class Matrix {
  // ...
  Matrix(size_t m, size_t n)
    : m_elems(m, std::vector<T>(n)), m_rows{m}, m_columns{n}
  {
    static_assert(std::is_arithmetic<T>::value, "Type parameter must be numeric.");
  }
};

This works, but you can bypass it by declaring an uninitialized pointer, such as std::unique_ptr<Matrix<string>>, and the code will compile, which may lead to subtle bugs. This is the reason why most of the time you should use static_assert inside functions and not methods/constructors. Also, keep in mind that functions using this technique cannot be overloaded for types not respecting the assertion.

We want the type-checking to happen at compile-time, but we also want it to happen at the time of the type declaration, even before calling the constructor.

In C++20, you can use concepts, which allow a readable and self-explanatory code.

#include <concepts>

template<typename T>
concept arithmetic = std::integral<T> or std::floating_point<T>;

template<arithmetic T>
class Matrix {
  // ...
};

The error message will also be very clear.

matrix.cpp:112:19: error: constraints not satisfied for class template 'Matrix' [with T = std::string]
  112 |   std::unique_ptr<Matrix<std::string>> ptr;
      |                   ^~~~~~~~~~~~~~~~~~~
...
1 error generated.

Note that template<arithmetic T> is equivalent to template<typename T> requires arithmetic<T>.

This technique satisfies our needs, but what if for some reasons, you cannot use C++20?

Enters std::enable_if. It will be less readable, but will do its job.

#include <type_traits>

template<
  typename T,
  typename = typename std::enable_if<std::is_arithmetic<T>::value, T>::type
>
class Matrix {
  // ...
};

The shown error will be less clear, but will point to the correct line containing it.

matrix.cpp:11:38: error: failed requirement 'std::is_arithmetic<std::string>::value'; 'enable_if' cannot be used to disable this declaration
   11 |   typename = typename std::enable_if<std::is_arithmetic<T>::value, T>::type
      |                                      ^~~~~~~~~~~~~~~~~~~~~~~~~~~~
matrix.cpp:111:19: note: in instantiation of default argument for 'Matrix<std::string>' required here
  111 |   std::unique_ptr<Matrix<std::string>> ptr;
      |                   ^~~~~~~~~~~~~~~~~~~
1 error generated.

Using Zig

I wanted to try out Zig and see how it fares against C and C++. Being a Zig beginner doesn’t help, since I dealt with many errors, and the error messages aren’t always clear. To construct our Matrix type, we’ll create a function that returns a type.

DISCLAIMER: this code uses the version 0.15.2. It may not work on later versions. Remember that below 1.0, Zig will evolve quickly, especially its standard library.

const std = @import("std");

fn Matrix(comptime T: type) type {
  // Type-checking is really straightforward
  const info = @typeInfo(T);
  if (info != .int and info != .float) {
      @compileError(std.fmt.comptimePrint("expected numeric type; found '{any}'.", .{T}));
  }

  return struct {
      const Self = @This();

      elems: std.ArrayList(std.ArrayList(T)),
      rows: usize,
      columns: usize,

      // Define methods here...
  };
}

Getting type information is done at compile-time, so the if statement also occurs at compile-time. Note that in Zig, functions starting with an @ are builtins defined by the compiler.

Should I use the *const Self as parameter type when I’m not modifying the struct?

You can simply use the Self type, because when structs are passed as parameters, Zig will pass it by reference if it is faster to do so rather than copying it.

Tests

Tests are very easy to write in Zig.

test "equality" {
    const allocator = std.testing.allocator;
    const values = &[_]u8{
        1, 1,
        1, 1,
        1, 1,
    };

    const matrix = try Matrix(u8).init(allocator, 3, 2, values);
    defer matrix.deinit();
    std.debug.assert(matrix.equalsSlice(values));

    const neighbors = try matrix.neighbors(allocator);
    defer neighbors.deinit();

    const expected = &[_]usize{
        3, 3,
        5, 5,
        3, 3,
    };
    std.debug.assert(neighbors.equalsSlice(expected));
}

Then you call zig build test.

Coding the program in Zig was tough, not because of the algorithm (I love the fact that optional types are builtin and easy to use), but, as a pure beginner, because of my numerous syntax errors. Error messages aren’t always clear, which can be annoying. The most annoying feature when coming from C or C++ is int casting. There is no implicit type promotion, so you have to learn the appropriate builtin functions to convert int types, but this article summarizes them very well. All in all I think it’s a good thing, as it creates a more predictable code.

Conclusion

The C code is straightforward but unsafe and not totally readable. Meanwhile, the C++ code is more complex, with the language offering multiple techniques, so it took more time to write it, and it always depends on your compiler version. The Zig code is safe, easily testable, and offers a more easier way to type-check comptime parameters. However, since the syntax differs a bit from C and C++, you have to take time to learn the syntax and features (Ziglings is a good start).

You can find the source code in the three languages here.