/*
Copyright 2026 Manuele Turini
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/
#include <Eigen/Dense>
#include <Eigen/Eigenvalues>
#include <algorithm>
#include <array>
#include <cmath>
#include <complex>
#include <cctype>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <limits>
#include <numeric>
#include <random>
#include <sstream>
#include <string>
#include <vector>
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
struct Particle {
double x = 0.0;
double y = 0.0;
};
struct RGB {
unsigned char r, g, b;
};
struct Metrics {
double shannon = 0.0; // Shannon entropy on a grid
int K = 0; // number of low modes chosen
double eigengap = 0.0; // spectral leap
double roughness = 0.0; // roughness of the low mode field
double spectral_concentration = 0.0; // concentration of low modes
double metamorphic_stability = 0.0; // stability under small disturbances
double omega = 0.0; // raw spectral index
double Om = 0.0; // modular / spectral order (penalized by entropy)
double Oo4 = 0.0; // 4-fold orientational order
double Oo6 = 0.0; // 6-fold orientational order
double Oe = 0.0; // order of exclusion / uniformity of distances
double mean_q4 = 0.0;
double cv_q4 = 0.0;
double mean_q6 = 0.0;
double cv_q6 = 0.0;
double mean_nn = 0.0;
double cv_nn = 0.0;
double frac_short_nn = 0.0;
};
struct Stats {
double mean = 0.0;
double stddev = 0.0;
double min = 0.0;
double max = 0.0;
};
static double clamp01(double x) {
return std::max(0.0, std::min(1.0, x));
}
static double sqr(double x) { return x * x; }
static double distancePeriodic(const Particle& a, const Particle& b) {
double dx = std::abs(a.x - b.x);
double dy = std::abs(a.y - b.y);
dx = std::min(dx, 1.0 - dx);
dy = std::min(dy, 1.0 - dy);
return std::sqrt(dx * dx + dy * dy);
}
static Particle wrapParticle(Particle p) {
while (p.x < 0.0) p.x += 1.0;
while (p.x >= 1.0) p.x -= 1.0;
while (p.y < 0.0) p.y += 1.0;
while (p.y >= 1.0) p.y -= 1.0;
return p;
}
static std::string sanitize(const std::string& s) {
std::string out;
for (char c : s) {
if (std::isalnum(static_cast<unsigned char>(c))) out.push_back(static_cast<char>(std::tolower(static_cast<unsigned char>(c))));
else if (c == ' ' || c == '-' || c == '/') out.push_back('_');
}
return out;
}
static std::string fmt(double x, int prec = 4) {
std::ostringstream oss;
oss << std::fixed << std::setprecision(prec) << x;
return oss.str();
}
static std::vector<Particle> generateRandomUniform(int N, std::mt19937& rng) {
std::uniform_real_distribution<double> U(0.0, 1.0);
std::vector<Particle> pts(N);
for (auto& p : pts) {
p.x = U(rng);
p.y = U(rng);
}
return pts;
}
static std::vector<Particle> generateGaussianBlobs(int N, std::mt19937& rng, int blobs = 6) {
std::uniform_real_distribution<double> U(0.1, 0.9);
std::normal_distribution<double> G(0.0, 0.045);
std::vector<Particle> centers(blobs);
for (auto& c : centers) {
c.x = U(rng);
c.y = U(rng);
}
std::uniform_int_distribution<int> pick(0, blobs - 1);
std::vector<Particle> pts;
pts.reserve(N);
for (int i = 0; i < N; ++i) {
auto c = centers[pick(rng)];
Particle p { c.x + G(rng), c.y + G(rng) };
pts.push_back(wrapParticle(p));
}
return pts;
}
static std::vector<Particle> generateLatticeClustered(int N, std::mt19937& rng, bool hexagonal) {
// 6 cluster in 2 x 3
const int rows = 2;
const int cols = 3;
const int clusters = rows * cols;
const int perCluster = N / clusters;
std::vector<Particle> pts;
pts.reserve(N);
std::uniform_real_distribution<double> Uj(-0.006, 0.006);
for (int r = 0; r < rows; ++r) {
for (int c = 0; c < cols; ++c) {
double cx = (c + 0.5) / cols;
double cy = (r + 0.5) / rows;
int n = perCluster;
int nx = 8;
int ny = n / nx;
if (ny <= 0) ny = 1;
double sx = 0.14 / std::max(1, nx - 1);
double sy = 0.14 / std::max(1, ny - 1);
for (int iy = 0; iy < ny; ++iy) {
for (int ix = 0; ix < nx; ++ix) {
if ((int)pts.size() >= N) break;
double x = cx + (ix - 0.5 * (nx - 1)) * sx;
double y = cy + (iy - 0.5 * (ny - 1)) * sy;
if (hexagonal && (iy % 2 == 1))
x += 0.5 * sx;
x += Uj(rng);
y += Uj(rng);
pts.push_back(wrapParticle({x, y}));
}
}
}
}
while ((int)pts.size() < N) {
pts.push_back(wrapParticle({0.5 + Uj(rng), 0.5 + Uj(rng)}));
}
return pts;
}
static std::vector<Particle> generatePoissonDisk(int N, std::mt19937& rng) {
std::uniform_real_distribution<double> U(0.0, 1.0);
std::vector<Particle> pts;
pts.reserve(N);
// r_min consistent with ~400 points
double rmin = 0.03;
int trials = 0;
while ((int)pts.size() < N && trials < 2000000) {
Particle p{U(rng), U(rng)};
bool ok = true;
for (const auto& q : pts) {
if (distancePeriodic(p, q) < rmin) {
ok = false;
break;
}
}
if (ok) pts.push_back(p);
++trials;
}
// If that's not enough, complete evenly
while ((int)pts.size() < N) {
pts.push_back({U(rng), U(rng)});
}
return pts;
}
static std::vector<Particle> generateDensityGradient(int N, std::mt19937& rng) {
std::uniform_real_distribution<double> U(0.0, 1.0);
std::vector<Particle> pts;
pts.reserve(N);
while ((int)pts.size() < N) {
double x = U(rng);
double y = U(rng);
double density = 0.15 + 0.85 * x; // linear gradient from left to right
if (U(rng) < density) {
pts.push_back({x, y});
}
}
return pts;
}
static double computeShannonEntropyGrid(const std::vector<Particle>& pts, int G = 20) {
std::vector<int> occ(G * G, 0);
for (const auto& p : pts) {
int ix = std::min(G - 1, std::max(0, int(p.x * G)));
int iy = std::min(G - 1, std::max(0, int(p.y * G)));
occ[iy * G + ix]++;
}
double N = double(pts.size());
double H = 0.0;
for (int c : occ) {
if (c <= 0) continue;
double p = c / N;
H -= p * std::log(p);
}
return H / std::log(double(G * G)); // normalized in [0,1]
}
static Eigen::MatrixXd buildAffinity(const std::vector<Particle>& pts, double sigma = 0.08) {
int N = (int)pts.size();
Eigen::MatrixXd W = Eigen::MatrixXd::Zero(N, N);
double inv2s2 = 1.0 / (2.0 * sigma * sigma);
for (int i = 0; i < N; ++i) {
for (int j = i + 1; j < N; ++j) {
double d = distancePeriodic(pts[i], pts[j]);
double w = std::exp(-d * d * inv2s2);
W(i, j) = W(j, i) = w;
}
}
return W;
}
static Eigen::MatrixXd buildNormalizedLaplacian(const Eigen::MatrixXd& W) {
int N = (int)W.rows();
Eigen::VectorXd d = W.rowwise().sum();
Eigen::VectorXd invSqrt = d.array().max(1e-12).inverse().sqrt();
Eigen::MatrixXd L = Eigen::MatrixXd::Identity(N, N);
for (int i = 0; i < N; ++i) {
for (int j = 0; j < N; ++j) {
if (i == j) continue;
L(i, j) = -W(i, j) * invSqrt(i) * invSqrt(j);
}
L(i, i) = 1.0 - W(i, i) * invSqrt(i) * invSqrt(i);
}
return L;
}
static int chooseKByEigengap(const Eigen::VectorXd& evals, int kmin = 3, int kmax = 12) {
kmax = std::min(kmax, (int)evals.size() - 2);
int bestK = kmin;
double bestGap = -1.0;
for (int k = kmin; k <= kmax; ++k) {
double gap = evals(k) - evals(k - 1);
if (gap > bestGap) {
bestGap = gap;
bestK = k;
}
}
return bestK;
}
static double computeRoughness(const std::vector<Particle>& pts, const Eigen::MatrixXd& lowModes) {
int N = (int)pts.size();
int K = (int)lowModes.cols();
double acc = 0.0;
int count = 0;
// neighbors within r
const double r = 0.08;
for (int i = 0; i < N; ++i) {
Eigen::VectorXd vi = lowModes.row(i);
for (int j = i + 1; j < N; ++j) {
double d = distancePeriodic(pts[i], pts[j]);
if (d <= r) {
Eigen::VectorXd vj = lowModes.row(j);
acc += (vi - vj).squaredNorm();
count++;
}
}
}
if (count == 0) return 0.0;
return acc / (count * K);
}
static double computeSpectralConcentration(const Eigen::VectorXd& evals, int K) {
// high score if the K low modes are small and separated from the rest
double low = 0.0;
for (int i = 0; i < K; ++i) low += evals(i);
low /= std::max(1, K);
double high = 0.0;
int M = 0;
for (int i = K; i < evals.size(); ++i) {
high += evals(i);
M++;
}
high /= std::max(1, M);
double score = (high - low) / std::max(1e-12, high);
return clamp01(score);
}
static double computeMetamorphicStability(const std::vector<Particle>& pts,
int trials,
std::mt19937& rng,
int Kref) {
std::normal_distribution<double> G(0.0, 0.0025);
Eigen::MatrixXd W0 = buildAffinity(pts);
Eigen::MatrixXd L0 = buildNormalizedLaplacian(W0);
Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> es0(L0);
Eigen::MatrixXd evecs0 = es0.eigenvectors().leftCols(Kref);
double score = 0.0;
for (int t = 0; t < trials; ++t) {
std::vector<Particle> pert = pts;
for (auto& p : pert) {
p.x += G(rng);
p.y += G(rng);
p = wrapParticle(p);
}
Eigen::MatrixXd W = buildAffinity(pert);
Eigen::MatrixXd L = buildNormalizedLaplacian(W);
Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> es(L);
Eigen::MatrixXd evecs = es.eigenvectors().leftCols(Kref);
// subspace similarity: ||U^T V||_F^2 / K
Eigen::MatrixXd C = evecs0.transpose() * evecs;
double subspaceSimilarity = C.squaredNorm() / double(Kref);
score += clamp01(subspaceSimilarity);
}
return score / std::max(1, trials);
}
static void computeNearestNeighborDistances(const std::vector<Particle>& pts,
std::vector<double>& nnd,
std::vector<std::array<int,6>>& nnIndices) {
const int N = (int)pts.size();
nnd.assign(N, std::numeric_limits<double>::max());
nnIndices.assign(N, std::array<int,6>{-1,-1,-1,-1,-1,-1});
for (int i = 0; i < N; ++i) {
std::vector<std::pair<double,int>> ds;
ds.reserve(N - 1);
for (int j = 0; j < N; ++j) {
if (i == j) continue;
ds.push_back({distancePeriodic(pts[i], pts[j]), j});
}
std::nth_element(ds.begin(), ds.begin() + 6, ds.end(),
[](const auto& a, const auto& b){ return a.first < b.first; });
std::sort(ds.begin(), ds.begin() + 6,
[](const auto& a, const auto& b){ return a.first < b.first; });
nnd[i] = ds[0].first;
for (int k = 0; k < 6; ++k) nnIndices[i][k] = ds[k].second;
}
}
static std::pair<double,double> meanAndCV(const std::vector<double>& v) {
double m = std::accumulate(v.begin(), v.end(), 0.0) / std::max<size_t>(1, v.size());
double s2 = 0.0;
for (double x : v) s2 += (x - m) * (x - m);
s2 /= std::max<size_t>(1, v.size());
double sd = std::sqrt(s2);
double cv = (m > 1e-12 ? sd / m : 0.0);
return {m, cv};
}
static double localBondOrderN(const std::vector<Particle>& pts,
const std::vector<std::array<int,6>>& nnIndices,
int i, int n) {
std::complex<double> sum(0.0, 0.0);
for (int k = 0; k < 6; ++k) {
int j = nnIndices[i][k];
if (j < 0) continue;
double dx = pts[j].x - pts[i].x;
double dy = pts[j].y - pts[i].y;
if (dx > 0.5) dx -= 1.0;
if (dx < -0.5) dx += 1.0;
if (dy > 0.5) dy -= 1.0;
if (dy < -0.5) dy += 1.0;
double theta = std::atan2(dy, dx);
std::complex<double> z(std::cos(n * theta), std::sin(n * theta));
sum += z;
}
sum /= 6.0;
return std::abs(sum);
}
static void computeOrientationalOrder(const std::vector<Particle>& pts,
Metrics& m) {
std::vector<double> nnd;
std::vector<std::array<int,6>> nnIndices;
computeNearestNeighborDistances(pts, nnd, nnIndices);
std::vector<double> q4(pts.size()), q6(pts.size());
for (int i = 0; i < (int)pts.size(); ++i) {
q4[i] = localBondOrderN(pts, nnIndices, i, 4);
q6[i] = localBondOrderN(pts, nnIndices, i, 6);
}
auto [mq4, cvq4] = meanAndCV(q4);
auto [mq6, cvq6] = meanAndCV(q6);
m.mean_q4 = mq4; m.cv_q4 = cvq4;
m.mean_q6 = mq6; m.cv_q6 = cvq6;
// score: high average + low variability
m.Oo4 = clamp01(mq4 * std::exp(-cvq4));
m.Oo6 = clamp01(mq6 * std::exp(-cvq6));
auto [mnn, cvnn] = meanAndCV(nnd);
m.mean_nn = mnn; m.cv_nn = cvnn;
int shortCount = 0;
for (double d : nnd) {
if (d < 0.7 * mnn) shortCount++;
}
m.frac_short_nn = shortCount / double(pts.size());
// exclusion order: uniform neighborhood and few too short distances
m.Oe = clamp01((1.0 - clamp01(cvnn)) * (1.0 - clamp01(2.0 * m.frac_short_nn)));
}
static Metrics analyzeSystem(const std::vector<Particle>& pts,
std::mt19937& rng,
int stabilityTrials = 10,
bool autoSelectK = true) {
Metrics m;
m.shannon = computeShannonEntropyGrid(pts, 20);
Eigen::MatrixXd W = buildAffinity(pts, 0.08);
Eigen::MatrixXd L = buildNormalizedLaplacian(W);
Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> es(L);
Eigen::VectorXd evals = es.eigenvalues();
Eigen::MatrixXd evecs = es.eigenvectors();
m.K = autoSelectK ? chooseKByEigengap(evals, 3, 12) : 6;
m.eigengap = evals(m.K) - evals(m.K - 1);
Eigen::MatrixXd lowModes = evecs.leftCols(m.K);
m.roughness = computeRoughness(pts, lowModes);
m.spectral_concentration = computeSpectralConcentration(evals, m.K);
m.metamorphic_stability = computeMetamorphicStability(pts, stabilityTrials, rng, m.K);
// raw spectral index
m.omega = clamp01(m.spectral_concentration * m.metamorphic_stability * std::exp(-8.0 * m.roughness));
// entropy penalty: modular order distinct from simple concentration
// if H ~ 1 ("very sparse" distribution), Om is reduced.
m.Om = clamp01(m.omega * std::exp(-5.0 * m.shannon));
computeOrientationalOrder(pts, m);
return m;
}
static Stats computeStats(const std::vector<double>& v) {
Stats s;
if (v.empty()) return s;
s.mean = std::accumulate(v.begin(), v.end(), 0.0) / v.size();
s.min = *std::min_element(v.begin(), v.end());
s.max = *std::max_element(v.begin(), v.end());
double acc = 0.0;
for (double x : v) acc += sqr(x - s.mean);
s.stddev = std::sqrt(acc / v.size());
return s;
}
static void printStatsLine(const std::string& name, const Stats& s) {
std::cout << std::left << std::setw(22) << name
<< " mean=" << std::setw(10) << std::fixed << std::setprecision(6) << s.mean
<< " std=" << std::setw(10) << s.stddev
<< " min=" << std::setw(10) << s.min
<< " max=" << std::setw(10) << s.max << "\n";
}
static void writeCSV(const std::string& filename,
const std::vector<std::string>& header,
const std::vector<std::vector<std::string>>& rows) {
std::ofstream f(filename.c_str());
for (size_t i = 0; i < header.size(); ++i) {
if (i) f << ",";
f << header[i];
}
f << "\n";
for (const auto& row : rows) {
for (size_t i = 0; i < row.size(); ++i) {
if (i) f << ",";
f << row[i];
}
f << "\n";
}
}
static void savePPM(const std::string& filename, const std::vector<RGB>& img, int W, int H) {
std::ofstream out(filename.c_str(), std::ios::binary);
out << "P6\n" << W << " " << H << "\n255\n";
out.write(reinterpret_cast<const char*>(img.data()), std::streamsize(img.size() * sizeof(RGB)));
}
static void drawFilledCircle(std::vector<RGB>& img, int W, int H, int cx, int cy, int r, RGB color) {
for (int dy = -r; dy <= r; ++dy) {
for (int dx = -r; dx <= r; ++dx) {
if (dx * dx + dy * dy <= r * r) {
int x = cx + dx;
int y = cy + dy;
if (x >= 0 && x < W && y >= 0 && y < H) {
img[y * W + x] = color;
}
}
}
}
}
static void drawRect(std::vector<RGB>& img, int W, int H, int x0, int y0, int x1, int y1, RGB color) {
x0 = std::max(0, x0); y0 = std::max(0, y0);
x1 = std::min(W, x1); y1 = std::min(H, y1);
for (int y = y0; y < y1; ++y)
for (int x = x0; x < x1; ++x)
img[y * W + x] = color;
}
// tiny 5x7 font for digits and letters used in labels
static const int FONT_W = 5;
static const int FONT_H = 7;
static const unsigned char fontDigits[10][FONT_H] = {
{0x0E,0x11,0x13,0x15,0x19,0x11,0x0E}, //0
{0x04,0x0C,0x04,0x04,0x04,0x04,0x0E}, //1
{0x0E,0x11,0x01,0x02,0x04,0x08,0x1F}, //2
{0x0E,0x11,0x01,0x06,0x01,0x11,0x0E}, //3
{0x02,0x06,0x0A,0x12,0x1F,0x02,0x02}, //4
{0x1F,0x10,0x1E,0x01,0x01,0x11,0x0E}, //5
{0x06,0x08,0x10,0x1E,0x11,0x11,0x0E}, //6
{0x1F,0x01,0x02,0x04,0x08,0x08,0x08}, //7
{0x0E,0x11,0x11,0x0E,0x11,0x11,0x0E}, //8
{0x0E,0x11,0x11,0x0F,0x01,0x02,0x0C} //9
};
static const unsigned char fontPunct[][FONT_H] = {
{0x00,0x00,0x00,0x00,0x00,0x00,0x00}, // space
{0x00,0x00,0x00,0x00,0x00,0x06,0x06}, // .
{0x00,0x00,0x00,0x1F,0x00,0x00,0x00}, // -
{0x00,0x00,0x00,0x00,0x00,0x0C,0x0C}, // ,
{0x00,0x04,0x04,0x04,0x00,0x00,0x00}, // :
{0x00,0x0E,0x11,0x01,0x02,0x00,0x02}, // ?
{0x00,0x00,0x00,0x1F,0x00,0x00,0x00}, // =
{0x00,0x04,0x0E,0x15,0x04,0x04,0x04}, // *
{0x00,0x00,0x00,0x00,0x00,0x00,0x00} // /
};
static const unsigned char fontLetters[][FONT_H] = {
{0x0E,0x11,0x11,0x1F,0x11,0x11,0x11}, // A
{0x1E,0x11,0x11,0x1E,0x11,0x11,0x1E}, // B
{0x0E,0x11,0x10,0x10,0x10,0x11,0x0E}, // C
{0x1C,0x12,0x11,0x11,0x11,0x12,0x1C}, // D
{0x1F,0x10,0x10,0x1E,0x10,0x10,0x1F}, // E
{0x1F,0x10,0x10,0x1E,0x10,0x10,0x10}, // F
{0x0E,0x11,0x10,0x17,0x11,0x11,0x0F}, // G
{0x11,0x11,0x11,0x1F,0x11,0x11,0x11}, // H
{0x0E,0x04,0x04,0x04,0x04,0x04,0x0E}, // I
{0x07,0x02,0x02,0x02,0x12,0x12,0x0C}, // J
{0x11,0x12,0x14,0x18,0x14,0x12,0x11}, // K
{0x10,0x10,0x10,0x10,0x10,0x10,0x1F}, // L
{0x11,0x1B,0x15,0x15,0x11,0x11,0x11}, // M
{0x11,0x19,0x15,0x13,0x11,0x11,0x11}, // N
{0x0E,0x11,0x11,0x11,0x11,0x11,0x0E}, // O
{0x1E,0x11,0x11,0x1E,0x10,0x10,0x10}, // P
{0x0E,0x11,0x11,0x11,0x15,0x12,0x0D}, // Q
{0x1E,0x11,0x11,0x1E,0x14,0x12,0x11}, // R
{0x0F,0x10,0x10,0x0E,0x01,0x01,0x1E}, // S
{0x1F,0x04,0x04,0x04,0x04,0x04,0x04}, // T
{0x11,0x11,0x11,0x11,0x11,0x11,0x0E}, // U
{0x11,0x11,0x11,0x11,0x11,0x0A,0x04}, // V
{0x11,0x11,0x11,0x15,0x15,0x15,0x0A}, // W
{0x11,0x11,0x0A,0x04,0x0A,0x11,0x11}, // X
{0x11,0x11,0x0A,0x04,0x04,0x04,0x04}, // Y
{0x1F,0x01,0x02,0x04,0x08,0x10,0x1F} // Z
};
static const unsigned char* getCharBitmap(char c) {
if (c >= '0' && c <= '9') return fontDigits[c - '0'];
if (c >= 'A' && c <= 'Z') return fontLetters[c - 'A'];
if (c >= 'a' && c <= 'z') return fontLetters[c - 'a'];
if (c == ' ') return fontPunct[0];
if (c == '.') return fontPunct[1];
if (c == '-') return fontPunct[2];
if (c == ',') return fontPunct[3];
if (c == ':') return fontPunct[4];
if (c == '?') return fontPunct[5];
if (c == '=') return fontPunct[6];
if (c == '*') return fontPunct[7];
if (c == '/') return fontPunct[8];
return fontPunct[0];
}
static void drawChar(std::vector<RGB>& img, int W, int H, int x, int y, char c, RGB col, int scale=2) {
const unsigned char* bm = getCharBitmap(c);
for (int row = 0; row < FONT_H; ++row) {
for (int colx = 0; colx < FONT_W; ++colx) {
if ((bm[row] >> (FONT_W - 1 - colx)) & 1) {
for (int dy = 0; dy < scale; ++dy) {
for (int dx = 0; dx < scale; ++dx) {
int xx = x + colx * scale + dx;
int yy = y + row * scale + dy;
if (xx >= 0 && xx < W && yy >= 0 && yy < H)
img[yy * W + xx] = col;
}
}
}
}
}
}
static void drawText(std::vector<RGB>& img, int W, int H, int x, int y,
const std::string& s, RGB col, int scale=2) {
int cursor = x;
for (char c : s) {
drawChar(img, W, H, cursor, y, c, col, scale);
cursor += (FONT_W + 1) * scale;
}
}
static RGB heatColor(double t) {
t = clamp01(t);
// gradient: blue -> cyan -> green -> yellow -> red
double r=0,g=0,b=0;
if (t < 0.25) {
double u = t / 0.25;
r = 0; g = u; b = 1;
} else if (t < 0.5) {
double u = (t - 0.25) / 0.25;
r = 0; g = 1; b = 1 - u;
} else if (t < 0.75) {
double u = (t - 0.5) / 0.25;
r = u; g = 1; b = 0;
} else {
double u = (t - 0.75) / 0.25;
r = 1; g = 1 - u; b = 0;
}
return RGB{(unsigned char)(255*r),(unsigned char)(255*g),(unsigned char)(255*b)};
}
static void saveParticleImage(const std::string& filename,
const std::vector<Particle>& pts,
const Metrics& m,
const std::string& title) {
const int W = 1100;
const int H = 760;
const int left = 60;
const int top = 70;
const int size = 620;
std::vector<RGB> img(W * H, RGB{255,255,255});
RGB black{0,0,0};
RGB gray{230,230,230};
RGB blue{30,80,220};
RGB dark{40,40,40};
// frame + grid
for (int x = left; x <= left + size; ++x) {
img[top * W + x] = black;
img[(top + size) * W + x] = black;
}
for (int y = top; y <= top + size; ++y) {
img[y * W + left] = black;
img[y * W + (left + size)] = black;
}
for (int g = 1; g < 10; ++g) {
int x = left + g * size / 10;
int y = top + g * size / 10;
for (int yy = top; yy <= top + size; ++yy) img[yy * W + x] = gray;
for (int xx = left; xx <= left + size; ++xx) img[y * W + xx] = gray;
}
// particles
for (const auto& p : pts) {
int x = left + int(p.x * size);
int y = top + int((1.0 - p.y) * size);
drawFilledCircle(img, W, H, x, y, 4, blue);
}
// title
drawText(img, W, H, 60, 18, title, dark, 3);
// metric bars
int bx = 740, by = 110;
int bw = 260, bh = 22;
std::vector<std::pair<std::string,double>> vals = {
{"Shannon", m.shannon},
{"Om", m.Om},
{"Oo4", m.Oo4},
{"Oo6", m.Oo6},
{"Oe", m.Oe}
};
for (size_t i = 0; i < vals.size(); ++i) {
int y = by + int(i) * 90;
drawText(img, W, H, bx, y, vals[i].first + " = " + fmt(vals[i].second, 4), dark, 2);
int ybar = y + 28;
for (int xx = 0; xx < bw; ++xx)
for (int yy = 0; yy < bh; ++yy)
img[(ybar + yy) * W + (bx + xx)] = RGB{245,245,245};
int fill = int(vals[i].second * bw);
RGB col = heatColor(vals[i].second);
drawRect(img, W, H, bx, ybar, bx + fill, ybar + bh, col);
// border
for (int xx = bx; xx < bx + bw; ++xx) {
img[ybar * W + xx] = black;
img[(ybar + bh) * W + xx] = black;
}
for (int yy = ybar; yy <= ybar + bh; ++yy) {
img[yy * W + bx] = black;
img[yy * W + (bx + bw)] = black;
}
}
// small notes
drawText(img, W, H, 735, 610, "K = " + std::to_string(m.K), dark, 2);
drawText(img, W, H, 735, 645, "eigengap = " + fmt(m.eigengap,4), dark, 2);
drawText(img, W, H, 735, 680, "roughness = " + fmt(m.roughness,4), dark, 2);
drawText(img, W, H, 735, 715, "stability = " + fmt(m.metamorphic_stability,4), dark, 2);
savePPM(filename, img, W, H);
}
static void saveEntropyMapImage(const std::string& filename,
const std::vector<Particle>& pts,
const std::string& title,
int G = 40) {
const int W = 920;
const int H = 760;
const int left = 60;
const int top = 70;
const int size = 620;
const int cell = size / G;
std::vector<RGB> img(W * H, RGB{255,255,255});
RGB black{0,0,0};
RGB gray{235,235,235};
std::vector<int> occ(G * G, 0);
for (const auto& p : pts) {
int ix = std::min(G - 1, std::max(0, int(p.x * G)));
int iy = std::min(G - 1, std::max(0, int(p.y * G)));
occ[iy * G + ix]++;
}
int cmax = *std::max_element(occ.begin(), occ.end());
if (cmax == 0) cmax = 1;
// heatmap
for (int iy = 0; iy < G; ++iy) {
for (int ix = 0; ix < G; ++ix) {
int c = occ[(G - 1 - iy) * G + ix];
double t = double(c) / cmax;
RGB col = heatColor(t);
int x0 = left + ix * cell;
int y0 = top + iy * cell;
drawRect(img, W, H, x0, y0, x0 + cell, y0 + cell, col);
}
}
// grid
for (int g = 0; g <= G; ++g) {
int x = left + g * cell;
int y = top + g * cell;
if (x >= left && x <= left + size) {
for (int yy = top; yy <= top + size; ++yy) img[yy * W + x] = gray;
}
if (y >= top && y <= top + size) {
for (int xx = left; xx <= left + size; ++xx) img[y * W + xx] = gray;
}
}
for (int x = left; x <= left + size; ++x) {
img[top * W + x] = black;
img[(top + size) * W + x] = black;
}
for (int y = top; y <= top + size; ++y) {
img[y * W + left] = black;
img[y * W + (left + size)] = black;
}
drawText(img, W, H, 60, 18, title, black, 3);
drawText(img, W, H, 710, 120, "densita locale", black, 2);
drawText(img, W, H, 710, 150, "su griglia", black, 2);
// legend
for (int i = 0; i < 200; ++i) {
double t = 1.0 - double(i) / 199.0;
RGB col = heatColor(t);
drawRect(img, W, H, 760, 240 + i, 790, 241 + i, col);
}
drawText(img, W, H, 805, 235, "alta", black, 2);
drawText(img, W, H, 805, 340, "media", black, 2);
drawText(img, W, H, 805, 440, "bassa", black, 2);
savePPM(filename, img, W, H);
}
template<typename Generator>
static void runBatch(const std::string& name,
Generator generator,
int runs, int N, int stabilityTrials,
std::mt19937& rng,
std::vector<std::vector<std::string>>& csvRuns,
std::vector<std::vector<std::string>>& csvSummary) {
std::vector<double> Hs, Ks, gaps, roughs, concs, stabs, omegas, Oms, Oo4s, Oo6s, Oes;
std::vector<double> mq4s, cvq4s, mq6s, cvq6s, mnns, cvnns, shorts;
std::vector<Particle> examplePts;
Metrics exampleM;
bool exampleSet = false;
std::cout << "\n[" << name << "]\n";
for (int i = 0; i < runs; ++i) {
std::vector<Particle> pts = generator(N, rng);
Metrics m = analyzeSystem(pts, rng, stabilityTrials, true);
Hs.push_back(m.shannon);
Ks.push_back((double)m.K);
gaps.push_back(m.eigengap);
roughs.push_back(m.roughness);
concs.push_back(m.spectral_concentration);
stabs.push_back(m.metamorphic_stability);
omegas.push_back(m.omega);
Oms.push_back(m.Om);
Oo4s.push_back(m.Oo4);
Oo6s.push_back(m.Oo6);
Oes.push_back(m.Oe);
mq4s.push_back(m.mean_q4);
cvq4s.push_back(m.cv_q4);
mq6s.push_back(m.mean_q6);
cvq6s.push_back(m.cv_q6);
mnns.push_back(m.mean_nn);
cvnns.push_back(m.cv_nn);
shorts.push_back(m.frac_short_nn);
csvRuns.push_back({
name,
std::to_string(i + 1),
fmt(m.shannon, 6),
std::to_string(m.K),
fmt(m.eigengap, 6),
fmt(m.roughness, 6),
fmt(m.spectral_concentration, 6),
fmt(m.metamorphic_stability, 6),
fmt(m.omega, 6),
fmt(m.Om, 6),
fmt(m.Oo4, 6),
fmt(m.Oo6, 6),
fmt(m.Oe, 6),
fmt(m.mean_q4, 6),
fmt(m.cv_q4, 6),
fmt(m.mean_q6, 6),
fmt(m.cv_q6, 6),
fmt(m.mean_nn, 6),
fmt(m.cv_nn, 6),
fmt(m.frac_short_nn, 6)
});
if (!exampleSet) {
examplePts = pts;
exampleM = m;
exampleSet = true;
}
if ((i + 1) % 10 == 0 || i == 0 || i + 1 == runs) {
std::cout << " run " << std::setw(3) << (i + 1)
<< " / " << runs
<< " Om=" << fmt(m.Om, 4)
<< " Oo4=" << fmt(m.Oo4, 4)
<< " Oo6=" << fmt(m.Oo6, 4)
<< " Oe=" << fmt(m.Oe, 4) << "\n";
}
}
std::cout << "\n=== " << name << " : summary on " << runs << " run ===\n";
printStatsLine("Shannon", computeStats(Hs));
printStatsLine("chosen K", computeStats(Ks));
printStatsLine("eigengap", computeStats(gaps));
printStatsLine("roughness R", computeStats(roughs));
printStatsLine("spectral concentration",computeStats(concs));
printStatsLine("metamorphic stability", computeStats(stabs));
printStatsLine("omega", computeStats(omegas));
printStatsLine("Om", computeStats(Oms));
printStatsLine("Oo4", computeStats(Oo4s));
printStatsLine("Oo6", computeStats(Oo6s));
printStatsLine("Oe", computeStats(Oes));
printStatsLine("mean q4", computeStats(mq4s));
printStatsLine("cv q4", computeStats(cvq4s));
printStatsLine("mean q6", computeStats(mq6s));
printStatsLine("cv q6", computeStats(cvq6s));
printStatsLine("mean nn", computeStats(mnns));
printStatsLine("cv nn", computeStats(cvnns));
printStatsLine("fraction short nn", computeStats(shorts));
auto addSummary = [&](const std::string& metric, const Stats& s) {
csvSummary.push_back({name, metric, fmt(s.mean, 6), fmt(s.stddev, 6), fmt(s.min, 6), fmt(s.max, 6)});
};
addSummary("Shannon", computeStats(Hs));
addSummary("chosen K", computeStats(Ks));
addSummary("eigengap", computeStats(gaps));
addSummary("roughness R", computeStats(roughs));
addSummary("spectral concentration", computeStats(concs));
addSummary("metamorphic stability", computeStats(stabs));
addSummary("omega", computeStats(omegas));
addSummary("Om", computeStats(Oms));
addSummary("Oo4", computeStats(Oo4s));
addSummary("Oo6", computeStats(Oo6s));
addSummary("Oe", computeStats(Oes));
addSummary("mean q4", computeStats(mq4s));
addSummary("cv q4", computeStats(cvq4s));
addSummary("mean q6", computeStats(mq6s));
addSummary("cv q6", computeStats(cvq6s));
addSummary("mean nn", computeStats(mnns));
addSummary("cv nn", computeStats(cvnns));
addSummary("fraction short nn", computeStats(shorts));
// Save example images
std::string base = sanitize(name);
saveParticleImage(base + "_particles.ppm", examplePts, exampleM, name);
saveEntropyMapImage(base + "_entropy_map.ppm", examplePts, "Density map - " + name);
std::cout << "\nExample " << name << ":\n"
<< " Shannon = " << fmt(exampleM.shannon, 6) << "\n"
<< " K = " << exampleM.K << "\n"
<< " eigengap= " << fmt(exampleM.eigengap, 6) << "\n"
<< " Om = " << fmt(exampleM.Om, 6) << "\n"
<< " Oo4 = " << fmt(exampleM.Oo4, 6) << "\n"
<< " Oo6 = " << fmt(exampleM.Oo6, 6) << "\n"
<< " Oe = " << fmt(exampleM.Oe, 6) << "\n"
<< " images: " << base << "_particles.ppm"
<< " " << base << "_entropy_map.ppm\n";
}
int main() {
std::mt19937 rng(123456u);
const int runs = 500;
const int N = 400;
const int stabilityTrials = 20;
std::cout << "Batch analysis started.\n";
std::cout << "Cases = 6\n";
std::cout << "Runs per case = " << runs << "\n";
std::cout << "Particles per run = " << N << "\n";
std::cout << "Stability trials / run = " << stabilityTrials << "\n";
std::cout << "Auto-select low modes = yes\n";
std::vector<std::vector<std::string>> csvRuns;
std::vector<std::vector<std::string>> csvSummary;
runBatch("Random uniform",
[](int N, std::mt19937& r){ return generateRandomUniform(N, r); },
runs, N, stabilityTrials, rng, csvRuns, csvSummary);
runBatch("Gaussian blobs",
[](int N, std::mt19937& r){ return generateGaussianBlobs(N, r, 6); },
runs, N, stabilityTrials, rng, csvRuns, csvSummary);
runBatch("Hex clustered",
[](int N, std::mt19937& r){ return generateLatticeClustered(N, r, true); },
runs, N, stabilityTrials, rng, csvRuns, csvSummary);
runBatch("Square clustered",
[](int N, std::mt19937& r){ return generateLatticeClustered(N, r, false); },
runs, N, stabilityTrials, rng, csvRuns, csvSummary);
runBatch("Poisson disk",
[](int N, std::mt19937& r){ return generatePoissonDisk(N, r); },
runs, N, stabilityTrials, rng, csvRuns, csvSummary);
runBatch("Density gradient",
[](int N, std::mt19937& r){ return generateDensityGradient(N, r); },
runs, N, stabilityTrials, rng, csvRuns, csvSummary);
writeCSV("batch_runs_final.csv",
{"case","run","Shannon","K","eigengap","roughness_R",
"spectral_concentration","metamorphic_stability","omega","Om","Oo4","Oo6","Oe",
"mean_q4","cv_q4","mean_q6","cv_q6","mean_nn","cv_nn","fraction_short_nn"},
csvRuns);
writeCSV("batch_summary_final.csv",
{"case","metric","mean","std","min","max"},
csvSummary);
std::cout << "\nWritten files:\n";
std::cout << " batch_runs_final.csv\n";
std::cout << " batch_summary_final.csv\n";
std::cout << " *_particles.ppm\n";
std::cout << " *_entropy_map.ppm\n";
std::cout << "\nNote: Om penalizes simple concentration via Shannon entropy,\n"
<< "while Oo4, Oo6, and Oe distinguish 4-, 6-symmetries\n"
<< "and exclusion/uniform spacing, respectively.\n";
return 0;
}