+void output_cycle_plot(const std::vector<pulse> &pulses, double calibration_factor)
+{
+ FILE *fp = fopen("cycles.plot", "w");
+ for (unsigned i = 0; i < pulses.size(); ++i) {
+ double cycles = pulses[i].len * calibration_factor * C64_FREQUENCY;
+ fprintf(fp, "%f %f\n", pulses[i].time, cycles);
+ }
+ fclose(fp);
+}
+
+std::pair<int, double> find_closest_point(double x, const std::vector<float> &points)
+{
+ int best_point = 0;
+ double best_dist = (x - points[0]) * (x - points[0]);
+ for (unsigned j = 1; j < train_snap_points.size(); ++j) {
+ double dist = (x - points[j]) * (x - points[j]);
+ if (dist < best_dist) {
+ best_point = j;
+ best_dist = dist;
+ }
+ }
+ return std::make_pair(best_point, best_dist);
+}
+
+float eval_badness(const std::vector<pulse>& pulses, double calibration_factor)
+{
+ double sum_badness = 0.0;
+ for (unsigned i = 0; i < pulses.size(); ++i) {
+ double cycles = pulses[i].len * calibration_factor * C64_FREQUENCY;
+ if (cycles > 2000.0) cycles = 2000.0; // Don't make pauses arbitrarily bad.
+ std::pair<int, double> selected_point_and_sq_dist = find_closest_point(cycles, train_snap_points);
+ sum_badness += selected_point_and_sq_dist.second;
+ }
+ return sqrt(sum_badness / (pulses.size() - 1));
+}
+
+void find_kmeans(const std::vector<pulse> &pulses, double calibration_factor, const std::vector<float> &initial_centers)
+{
+ std::vector<float> last_centers = initial_centers;
+ std::vector<float> sums;
+ std::vector<float> num;
+ sums.resize(initial_centers.size());
+ num.resize(initial_centers.size());
+ for ( ;; ) {
+ for (unsigned i = 0; i < initial_centers.size(); ++i) {
+ sums[i] = 0.0f;
+ num[i] = 0;
+ }
+ for (unsigned i = 0; i < pulses.size(); ++i) {
+ double cycles = pulses[i].len * calibration_factor * C64_FREQUENCY;
+ // Ignore heavy outliers, which are almost always long pauses.
+ if (cycles > 2000.0) {