Just for fun, after reading Riley's answer, I decided to modify it into a simulated annealing optimization program. The basic idea here is: suppose you want a ball to find the lowest point in a terrain. In order to avoid finding the bottom of a shallow depression (local minimum), you allow the ball to bounce at random by a gradually decreasing amount. So, if the ball is bouncing in one of the shallow depressions with a deeper depression nearby, then eventually the "bounce amount" will nudge the ball out of the shallow depression, but to a point where it can no longer escape the deeper depression.
The exact formulation I'm using is roughly based on thermodynamics (where the term "simulated annealing" originated, I think, from a metallurgy technique "annealing" to create good alloys by a similar process of heating first, then gradually cooling). In my particular case, I was using "energy" as negative of a board's score, since simulated annealing minimizes an energy function whereas I wanted to find a maximum score.
Code:
#include <random>
#include <bitset>
#include <cmath>
#include <iostream>
#include <functional>
#include <unistd.h>
#include <time.h>
#include <boost/range/irange.hpp>
class Board {
public:
Board() = delete;
Board(std::mt19937& randGen) : m_randGen{randGen} { }
Board(const Board&) = default;
Board(Board&&) = default;
Board& operator=(const Board&) = default;
Board& operator=(Board&&) = default;
~Board() = default;
int score() const {
int result = 0;
for (int r : boost::irange(0, board_height)) {
for (int c : boost::irange(0, board_width)) {
if (cell(r, c))
continue;
for (int dr : boost::irange(-1, 2)) {
for (int dc : boost::irange(-1, 2)) {
if (checkedCell(r + dr, c + dc))
result++;
}
}
}
}
return result;
}
void mutate() {
std::uniform_int_distribution<> dist { 0, board_size - 1 };
int slot = dist(m_randGen.get());
m_board.flip(slot);
}
private:
static constexpr int board_height = 9;
static constexpr int board_width = 9;
static constexpr int board_size = board_height * board_width;
using internal_type = std::bitset<board_size>;
internal_type m_board;
std::reference_wrapper<std::mt19937> m_randGen;
bool cell(int r, int c) const {
return m_board[r * board_width + c];
}
bool checkedCell(int r, int c) const {
return (r >= 0 && r < board_height && c >= 0 && c < board_width) &&
cell(r, c);
}
friend std::ostream& operator<<(std::ostream& os, const Board& b) {
for (int r : boost::irange(0, board_height)) {
for (int c : boost::irange(0, board_width)) {
if (b.cell(r, c))
os << '*';
else {
int nbrCount = 0;
for (int dr : boost::irange(-1, 2)) {
for (int dc : boost::irange(-1, 2)) {
if (b.checkedCell(r + dr, c + dc))
nbrCount++;
}
}
os << char('0' + nbrCount);
}
}
os << '\n';
}
return os;
}
};
int main() {
std::mt19937 randGen;
std::seed_seq sseq{int(getpid()), int(time(nullptr))};
randGen.seed(sseq);
std::uniform_real_distribution<> probDist;
Board bestBoardSoFar{randGen}, currBoard{randGen};
int bestScoreSoFar = 0, currScore = 0;
[[maybe_unused]] constexpr double startTemp = 1000.0;
constexpr double stepTemp = 0.001;
constexpr int numTempSteps = 1000000;
for (int i : boost::irange(0, numTempSteps)) {
double temp = stepTemp * (numTempSteps - i);
Board newBoard = currBoard;
newBoard.mutate();
int newScore = newBoard.score();
if (newScore > bestScoreSoFar) {
bestBoardSoFar = newBoard;
bestScoreSoFar = newScore;
}
// Naively, we would do:
// double currBoardProb = std::exp(currScore / temp);
// double newBoardProb = std::exp(newScore / temp);
//
// However, at low temperatures, this could cause both to be
// inf values. To avoid this, we scale both down by the same
// factor of std::exp(currScore / temp)
double currBoardProb = 1.0;
double newBoardProb = std::exp((newScore - currScore) / temp);
if (probDist(randGen) > currBoardProb / (currBoardProb + newBoardProb)) {
currBoard = newBoard;
currScore = newScore;
}
}
std::cout << "Found solution with score " << bestScoreSoFar << '\n';
std::cout << bestBoardSoFar;
return 0;
}
Some sample outputs:
Found solution with score 183
*4*4*4*4*
*6*6*6*5*
*6*5**5*3
*4*45*6*3
355*5*6*3
****6*6*4
3545**5**
*4*5*6*53
*4*4*4*3*
Found solution with score 196
2*4*4*4*2
3*6*6*6*3
3*6*6*6*3
3*6*6*6*4
3*6*6*5**
3*6*6*654
3*6*5****
3*6*54664
2*4*3****
Found solution with score 200
*********
466666664
*********
466666664
*********
466666664
*********
466666664
*********
It's interesting that even when it doesn't find an optimal solution, the solutions sort of look like they have multiple "crystalization domains" which are locally of the form from the solutions with score 200.
.
So, given that it pretty consistently gets close to 200 (the worst run so far was 183) and hasn't yet beaten 200 in the several times I've run it, this would seem to be (weak) evidence that 200 could indeed be the global maximum.
