+#if FIND_OPTIMAL_STREAM_ASSIGNMENT
+double find_best_assignment(const int *medoids, int *assignment)
+{
+ double current_score = 0.0;
+ for (int i = 0; i < 64; ++i) {
+ int best_medoid = medoids[0];
+ float best_medoid_score = kl_dist[i][medoids[0]];
+ for (int j = 1; j < NUM_CLUSTERS; ++j) {
+ if (kl_dist[i][medoids[j]] < best_medoid_score) {
+ best_medoid = medoids[j];
+ best_medoid_score = kl_dist[i][medoids[j]];
+ }
+ }
+ assignment[i] = best_medoid;
+ current_score += best_medoid_score;
+ }
+ return current_score;
+}
+
+double find_inv_sum(const SymbolStats &stats)
+{
+ double s = 0.0;
+ for (unsigned j = 0; j < NUM_SYMS; ++j) {
+ s += stats.freqs[j] + 0.5;
+ }
+ return 1.0 / s;
+}
+
+void find_optimal_stream_assignment(int base)
+{
+ // k-means init; make random assignments
+ std::random_device rd;
+ std::mt19937 g(rd());
+ std::uniform_int_distribution<> u(0, NUM_CLUSTERS - 1);
+ int assignment[64];
+ for (unsigned i = 0; i < 64; ++i) {
+ assignment[i] = u(g);
+ }
+ double inv_sum_coeffs[64];
+ for (unsigned i = 0; i < 64; ++i) {
+ inv_sum_coeffs[i] = find_inv_sum(stats[i + base]);
+ }
+
+ for (unsigned iter = 0; iter < 1000; ++iter) {
+ // make new clusters based on the current assignments
+ SymbolStats clusters[NUM_CLUSTERS];
+ for (unsigned i = 0; i < NUM_CLUSTERS; ++i) {
+ clusters[i].clear();
+ }
+ for (unsigned i = 0; i < 64; ++i) {
+ for (unsigned j = 0; j < NUM_SYMS; ++j) {
+ clusters[assignment[i]].freqs[j] += stats[i + base].freqs[j];
+ }
+ }
+
+ double inv_sum_clusters[NUM_CLUSTERS];
+ for (unsigned i = 0; i < NUM_CLUSTERS; ++i) {
+ inv_sum_clusters[i] = find_inv_sum(clusters[i]);
+ }
+
+ // find new assignments based on distance to the clusters
+ bool any_changed = false;
+ double total_d = 0.0;
+ for (unsigned i = 0; i < 64; ++i) {
+ int best_assignment = -1;
+ double best_distance = HUGE_VAL;
+ for (unsigned j = 0; j < NUM_CLUSTERS; ++j) {
+ double d = 0.0;
+ for (unsigned k = 0; k < NUM_SYMS; ++k) {
+ double p1 = (clusters[j].freqs[k] + 0.5) * inv_sum_clusters[j];
+ double p2 = (stats[i + base].freqs[k] + 0.5) * inv_sum_coeffs[i];
+
+ // K-L divergence is asymmetric; this is a hack.
+ d += p1 * log(p1 / p2);
+ d += p2 * log(p2 / p1);
+ }
+ if (d < best_distance) {
+ best_assignment = j;
+ best_distance = d;
+ }
+ }
+ if (assignment[i] != best_assignment) {
+ any_changed = true;
+ }
+ assignment[i] = best_assignment;
+ total_d += best_distance;
+ }
+ printf("iter %u: %.3f\n", iter, total_d);
+ if (!any_changed) break;
+ }
+ printf("\n");
+ std::unordered_map<int, int> rmap;
+ for (int i = 0; i < 64; ++i) {
+ if (i % 8 == 0) printf("\n");
+ if (!rmap.count(assignment[i])) {
+ rmap.emplace(assignment[i], rmap.size());
+ }
+ printf("%d, ", rmap[assignment[i]]);
+ }
+ printf("\n");
+}
+#endif
+