Can you beat the optimizer?

Question

This is an optimization/algorithms challenge. Given two sets of points X and Y, which are to be thought of as points on a circle of circumference T, the challenge is to find the rotation which minimizes the distance from the set X to the set Y. Here the distance is just the sum of the absolute distance from each point in X to its closest point in Y. Here closest is really just the minimum distance you need to travel along the circle to reach a point in Y (remembering that you can go clockwise or anti-clockwise). Note that this measure of distance is not symmetric.

In order to be able to create the test data, you will need to be able to run Python, but your answer can be in any language of your choosing that is available without cost to run on Linux. You can also use any libraries you choose as long as they are easy to install on Linux and available without cost.

The test data will be created using the following Python code. It outputs two lines to standard out. The first is circle X and the second is circle Y.

import numpy as np
#Set T to 1000, 10000, 100000, 1000000, 10000000, 1000000000
T = 1000
N = int(T/5)

#Set the seed to make data reproducible
np.random.seed(seed = 10)

timesY = np.sort(np.random.uniform(low = 0.0, high = T, size = N))
indices = np.sort(np.random.randint(0,N, N/10))
rotation =  np.random.randint(0,T, 1)[0]
timesX = [(t + rotation + np.random.random()) % T for t in timesY[indices]]

for tX in timesX:
    print tX,
print
for tY in timesY:
    print tY,

For those who don't run python, I have also uploaded some test files. Each file has three lines. The first line is the points for circle X, the second is the points for circle Y, and the last is the score you get from the generic optimizer in the python code below. This score will be in no way optimal and you should do much better.

http://wikisend.com/download/359660/1000.txt.lzma
http://wikisend.com/download/707350/10000.txt.lzma
http://wikisend.com/download/582950/100000.txt.lzma
http://wikisend.com/download/654002/1000000.txt.lzma
http://wikisend.com/download/416516/10000000.txt.lzma

Here is also a toy worked example of the distance between two sets of point to make the distance function clear.

Let timesX = [2, 4, 98] and timesY = [1, 7, 9, 99]. Let us assume the circumference of the circle is 100. As it stands the distance is (2-1)+(7-4)+(99-98) = 5. Now if we rotate timesX by 3 we get [5, 7, 1] which has total absolute distance 2 to timesY. In this case this appears to be the optimal solution.

Here is some sample code which assumes that the circles are called timesX and timesY, and outputs the minimum distance it finds. It doesn't, however, do a very good job, as you will see.

from bisect import bisect_left
from scipy.optimize import minimize_scalar
def takeClosest(myList, myNumber, T):
    """
    Assumes myList is sorted. Returns closest value to myNumber in a circle of circumference T.
    """
    pos = bisect_left(myList, myNumber)
    if (pos == 0 and myList[pos] != myNumber):
        before = myList[pos - 1] - T
        after = myList[0]
    elif (pos == len(myList)):
        before = myList[pos-1]
        after = myList[0] + T
    else:
        before = myList[pos - 1]
        after = myList[pos]
    if after - myNumber < myNumber - before:
        return after
    else:
        return before

def circle_dist(timesX, timesY):
    dist = 0
    for t in timesX:
        closest_number = takeClosest(timesY, t, T)
        dist += np.abs(closest_number - t)
    return dist

def dist(x):
    timesX_rotated = [(tX+x) % T for tX in timesX]
    return circle_dist(timesX_rotated, timesY)

#Now perform the optimization
res = minimize_scalar(dist, bounds=(0,T), method='bounded', options={ 'disp':True})
print res.fun

Scoring

You should run your code for T = 1000, 10000, 100000, 1000000, 10000000. Your code should output both the optimum rotation and the distance you get for it so that it can be checked. The scoring criterion is as follows.

• You should calculate the sum of the distances you find for the five values of T. Low is better and the best possible is 0. You aren't allowed to cheat by looking at what rotation was used to create the data in the first place.
• Then you should calculate 551431.326508571/(the sum you just calculated) and report that as your score.

---

Edit 1. Changed score so that high is better and removed time restriction so you can just run it on your own computer.

Edit 2 Note that the function dist can be used as a verifier. Just put the points for circle X in timesX and the points for circle Y in timesY and set x to be the rotation you want to test. It will return a score which you can compare to what you get.

Edit 3 Fixed the data creation code which wasn't exactly the code I had used to create the data.

Answer 1

C++, 10.60236

If we assume nothing about the data then it seems unlikely that we can do much better than brute force/Monte Carlo plus maybe some local minimization. The distance/rotation graphs for T=1,000 and T=10,000 are already pretty noisy:

However, assuming this is fair game, taking into account the fact that X is simply a rotated subset of Y with some small error, recovering the approximate rotation and minimizing locally is easy:

Compile with g++ mindist.cpp -omindist -std=c++11 -O3.

Run with ./mindist T < dataset (e.g. ./mindist 1000 < 1000.txt).

#include <iostream>
#include <vector>
#include <string>
#include <sstream>
#include <algorithm>
#include <iterator>
#include <tuple>
#include <cmath>

using namespace std;

/* Positive modulo */
template <typename T>
T mod(T x, T y) { x = fmod(x, y); return x >= 0 ? x : x + y; }

/* Less-than comparison with a margin. Can be used to search for elements
 * within +-epsilon. */
template <typename T>
struct epsilon_less {
    T epsilon;
    bool operator()(const T& x, const T& y) const { return x < y - epsilon; }
};

/* Calculates the distance between X and Y on a circle with circumference T,
 * where X is rotated by `rotation`. */
template <typename P, typename Ps>
P distance(P T, const Ps& X, const Ps& Y, P rotation) {
    P dist = 0;
    for (P x : X) {
        x = mod(x + rotation, T);
        auto l = lower_bound(Y.begin(), Y.end(), x);
        auto r = prev(l);
        dist += min(*l - x, x - *r);
    }
    return dist;
}

/* Minimizes a function f on the interval [l, r]. f is assumed to be "well-
 * behaved" within this range. This is a simple n-iterations binary search
 * based on the sign of the derivative of f. */
template <typename T, typename F>
T minimize(F f, T l, T r, int n = 32) {
    auto fl = f(l), fr = f(r);
    for (; n; --n) {
        T m = (l + r) / 2;
        auto fm = f(m);
        T dm = (r - l) / 1e5;
        auto fm2 = f(m + dm);
        auto df_dm = (fm2 - fm) / dm;
        (df_dm < 0 ? tie(l, fl) : tie(r, fr)) = make_tuple(m, fm);
    }
    return (l + r) / 2;
}

int main(int argc, const char* argv[]) {
    typedef double point;

    const point T = stod(argv[1]);
    const point error_margin = /*+-*/1;

    // Read input. Assumed to be sorted.
    vector<point> X, Y(1);
    for (auto* P : {&X, &Y}) {
        string line; getline(cin, line);
        stringstream s(line);
        copy(
            istream_iterator<point>(s), istream_iterator<point>(),
            back_inserter(*P)
        );
    }
    // Add wrap-around points to Y to simplify distance calculation
    Y[0] = Y.back() - T; Y.push_back(Y[1] + T);

    for (point y : Y) {
        // For each point y in Y we rotation X so that the last point in X
        // aligns with y.
        auto rotation = mod(y - X.back(), T);
        // See if the rest of X aligns within the error margin
        if (all_of(
            X.begin(), X.end(), [=, &Y] (point x) {
                return binary_search(
                    Y.begin(), Y.end(), mod(x + rotation, T),
                    epsilon_less<point>{error_margin}
                );
            }
        )) {
            // We found a rotation. Minimize the distance within the error
            // margin.
            rotation = minimize(
                [&] (point r) { return distance(T, X, Y, r); },
                rotation - error_margin, rotation + error_margin
            );
            // Print results and exit.
            cout.precision(12);
            cout << "Rotation: " << rotation << endl;
            cout << "Distance: " << distance(T, X, Y, rotation) << endl;
            return 0;
        }
    }
    // We didn't find a rotation. Something must have went wrong!
    return 1;
}

The program assumes X is a rotated subset of Y with an error margin of ±1 for each point in X. It tries to rotate X such that it aligns with Y within this error margin (that is, s.t. each point in X is at distance at most 1 from Y.) The reason that we can do this efficiently is that we don't have to try arbitrary rotations: It's enough to choose an arbitrary point x in X and rotate X such that x aligns with each of the points in Y. For one such rotation it's guaranteed that X and Y align within the error range; for every other rotation, it's very unlikely that even two points will align, so we can abandon nonoptimal rotations very quickly. It's extremely unlikely that two different rotations will both align X and Y, so once the program finds a rotation that does it doesn't bother to look any further. Once a rotation is found, it's adjusted within the error margin to minimize the overall distance.

Results:

T           Rotation         Distance
     1,000        494.64423       5.59189
    10,000      6,967.55630      45.03235
   100,000     24,454.48700     471.71127
 1,000,000    974,011.50572   4,671.75611
10,000,000  2,720,706.49829  46,816.15372

Total distance: 52,010.24535
Score: 10.60236

Answer 2

Java - 10.60236

As promised, Lembik:

import java.io.BufferedReader;
import java.io.FileReader;
import java.io.IOException;
import java.util.*;

public final class SpinDots {

public static void main( String[] sArgs ) {
    final int T = Integer.parseInt( sArgs[0] );
    List<Double> LDX = new LinkedList<>(),
        LDY = new LinkedList<>();

    double RS = 0.;

    try( BufferedReader BR = new BufferedReader( new FileReader( T + ".txt" ), 65536 ) ) {
        String SX = BR.readLine();
        StringTokenizer STX = new StringTokenizer( SX );

        while( STX.hasMoreTokens() )
            LDX.add( Double.parseDouble( STX.nextToken() ) );

        String SY = BR.readLine();
        StringTokenizer STY = new StringTokenizer( SY );

        while( STY.hasMoreTokens() )
            LDY.add( Double.parseDouble( STY.nextToken() ) );

        RS = Double.parseDouble( BR.readLine() );
        BR.close();
    } catch( IOException ioe ) {
        System.err.println( "Unable to read contents of file: " + T + ".txt" );
        System.exit( 1 );
    }

    final int NX = LDX.size(),
        NY = LDY.size(),
        RNY = 2*NY+3;

    final double[] DX = new double[NX],
        DY = new double[RNY];

    final int[] X = new int[NX],
        Y = new int[RNY];
    int[] YS = new int[RNY],
        YSP = new int[RNY];

    int dx = 0;
    for( Double dx_ : LDX ) {
        DX[dx] = dx_;
        X[dx] = (int)Math.round( dx_ );
        dx++;
    }
    Arrays.sort( X );
    Arrays.sort( DX );

    int dy = 1;
    for( Double dy_ : LDY ) {
        DY[dy] = dy_;
        DY[dy+NY] = dy_+T;
        Y[dy] = (int)Math.round( dy_ );
        Y[dy+NY] = Y[dy] + T;
        dy++;
    }
    DY[0] = DY[NY] - T;
    DY[2*NY+1] = DY[1] + 2*T;
    DY[2*NY+2] = DY[2] + 2*T;
    Y[0] = Y[NY] - T;
    Y[2*NY+1] = Y[1] + 2*T;
    Y[2*NY+2] = Y[2] + 2*T;

    System.out.println( "Computing optimum shift for: |X| = " + NX + ", |Y| = " + NY + ", T = " + T );

    int DO_ = scoreForShiftFirstPass( X, Y, YSP );
    List<Integer> LXO = new LinkedList<>(),
        LDXO = new LinkedList<>();

    LXO.add( 0 );
    LDXO.add( DO_ );
    for( int iShift = 1; iShift < T; iShift++ ) {
        int D_ = scoreForShift( X, Y, iShift, YSP, YS );
        int[] YT = YSP;
        YSP = YS;
        YS = YT;
        if( D_ < DO_ ) {
            DO_ = D_;
            Iterator<Integer> I = LXO.iterator(),
                DI = LDXO.iterator();

            while( DI.hasNext() ) {
                I.next();
                if( DI.next() > DO_ + NX ) {
                    DI.remove();
                    I.remove();
                }
            }
        }
        if( D_ <= DO_ + NX && LXO.size() < 100 ) {
            LXO.add( iShift );
            LDXO.add( D_ );
        }

        if( iShift %(T/100) == 0 )
            System.out.println( "Completed " + (iShift*100/T) + "% of coarse pass." );
    }

    System.out.println( "Completed 100% of coarse pass. Coarse optimum is: " + DO_ );
    System.out.println( "Identified " + LXO.size() + " fine filtering candidates." );

    double ddoo = 0.;
    double DDOO_ = Double.MAX_VALUE;
    for( int iShift : LXO ) {
        double DDO_ = doubleScoreForShiftFirstPass( DX, DY, iShift-1, YSP );
        double ddo = iShift;

        for( int iFineShift = -999; iFineShift < 1000; iFineShift++ ) {
            double dd = iShift + iFineShift/1000.0;
            double DD_ = doubleScoreForShift( DX, DY, dd, YSP, YS );
            int[] YT = YSP;
            YSP = YS;
            YS = YT;
            if( DD_ < DDO_ ) {
                DDO_ = DD_;
                ddo = dd;
            }
        }
        if( DDO_ < DDOO_ ) {
            DDOO_ = DDO_;
            ddoo = ddo;
        }
    }

    System.out.println( "Fine optimum is: " + DDOO_ );

    double ddooo = ddoo-1e-3;
    double DDOOO_ = doubleScoreForShiftFirstPass( DX, DY, ddoo-1e-3, YSP );
    for( int iHyperFineShift = -999; iHyperFineShift < 1000; iHyperFineShift++ ) {
        double dd = ddoo + iHyperFineShift/1.e6;
        double DD_ = doubleScoreForShift( DX, DY, dd, YSP, YS );
        int[] YT = YSP;
        YSP = YS;
        YS = YT;
        if( DD_ < DDOOO_ ) {
            DDOOO_ = DD_;
            ddooo = dd;
        }
    }
    System.out.println( "Optimal forward rotation in X: " + ddooo );
    System.out.println( "Optimal score: " + DDOOO_ );
    System.out.println( "Reference score: " + RS );
}

static int scoreForShift( int[] X, int[] Y, int iShift, int[] YSP, int[] YS ) {
    int D_ = 0;
    for( int x = 0; x < X.length; x++ ) {
        int x_ = X[x]+iShift;
        int y = YSP[x],
            y_ = Y[y];

        int D = x_ < y_ ? y_ - x_ : x_ - y_;
        int y_nx_ = Y[y+1];
        int D_nx = x_ < y_nx_ ? y_nx_ - x_ : x_ - y_nx_;
        while( D_nx <= D ) {
            D = D_nx;
            y_nx_ = Y[++y+1];
            D_nx = x_ < y_nx_ ? y_nx_ - x_ : x_ - y_nx_;
        }
        YS[x] = y;
        D_ += D;
    }
    return D_;
}

static int scoreForShiftFirstPass( int[] X, int[] Y, int[] YS ) {
    int D_ = 0;
    int y = 0;
    for( int x = 0; x < X.length; x++ ) {
        int x_ = X[x];
        int y_ = Y[y];

        int D = x_ < y_ ? y_ - x_ : x_ - y_;
        int y_nx_ = Y[y+1];
        int D_nx = x_ < y_nx_ ? y_nx_ - x_ : x_ - y_nx_;
        while( D_nx <= D ) {
            D = D_nx;
            y_nx_ = Y[++y+1];
            D_nx = x_ < y_nx_ ? y_nx_ - x_ : x_ - y_nx_;
        }
        YS[x] = y;
        D_ += D;
    }
    return D_;
}

static double doubleScoreForShift( double[] X, double[] Y, double d_Shift, int[] YSP, int[] YS ) {
    double D_ = 0;
    for( int x = 0; x < X.length; x++ ) {
        double x_ = X[x]+d_Shift;
        int y = YSP[x];
        double y_ = Y[y];

        double D = x_ < y_ ? y_ - x_ : x_ - y_;
        double y_nx_ = Y[y+1];
        double D_nx = x_ < y_nx_ ? y_nx_ - x_ : x_ - y_nx_;
        while( D_nx <= D ) {
            D = D_nx;
            y_nx_ = Y[++y+1];
            D_nx = x_ < y_nx_ ? y_nx_ - x_ : x_ - y_nx_;
        }
        YS[x] = y;
        D_ += D;
    }
    return D_;
}

static double doubleScoreForShiftFirstPass( double[] X, double[] Y, double d_Shift, int[] YS ) {
    double D_ = 0;
    int y = 0;
    for( int x = 0; x < X.length; x++ ) {
        double x_ = X[x]+d_Shift;
        double y_ = Y[y];

        double D = x_ < y_ ? y_ - x_ : x_ - y_;
        double y_nx_ = Y[y+1];
        double D_nx = x_ < y_nx_ ? y_nx_ - x_ : x_ - y_nx_;
        while( D_nx <= D ) {
            D = D_nx;
            y_nx_ = Y[++y+1];
            D_nx = x_ < y_nx_ ? y_nx_ - x_ : x_ - y_nx_;
        }
        YS[x] = y;
        D_ += D;
    }
    return D_;
}
}

My approach is similar to Ell's in that it does a grid search over all potentials (first at integer resolution, then at a resolution of 1/1,000ths, then at a resolution of 1/1,000,000ths). Unlike Ell's approach, mine doesn't require or make use of the fact that the error in X is tightly bounded. I literally evaluate all 10,000,000 possibilities during coarse filtering in the T = 10,000,000 case and store a list of candidates with distance small enough to be within "fine filtering" distance of the optimum.

Modular arithmetic is slow, and Java doesn't support template methods for native types. Hence my implementation "unwinds" Y to avoid the use of remainder operators, and the four working methods are near-identical clones of each other with differing types and 1-2 lines of code that would be handled in C/C++ using templates and macros.

I compute shifts sequentially, which allows the nearest neighbour search to incrementally tweak values from the previous shift. As a result, the routine completes in about 2 hours, which was acceptable since I slept through it. ;)

My outcomes are technically better than Ell's since I search to a higher resolution, but the improvement is utterly negligible, and probabilistic to boot. In particular, a higher resolution search only yields an improved score if the number of shifted X values falling in the set of intervals [Y-eps,Y] is different from the number of shifted X values falling in the set of intervals [Y,Y+eps], where eps is the coarser resolution (0.000001, in the case of Ell's implementation). I don't know if this is the case for any of the datasets in this problem. In any event, his code should "win" since it predates mine.

Output

Optimal forward rotation in X: 494.64500000000004
Optimal score: 5.591890690289915
Reference score: 33.0305238422

Optimal forward rotation in X: 6967.556309
Optimal score: 45.032350687407416
Reference score: 429.211871389

Optimal forward rotation in X: 24454.487003000002
Optimal score: 471.7112722703896
Reference score: 4676.23391624

Optimal forward rotation in X: 974011.505734
Optimal score: 4671.756109326729
Reference score: 49148.9676631

Optimal forward rotation in X: 2720706.498309
Optimal score: 46816.15372609813
Reference score: 497143.882534

Answer 3

Mathematica, 1.04834

Okay, maybe people just need to be shown how simple it is to perform better than the reference implementation.

So here is a random implementation. I'm simply trying out some 10,000 random rotations (basically, as many as I care to wait for) and pick the best one I find.

First@SortBy[(
  r = RandomReal[t];
  f = 
   Nearest[Flatten[{First[z = Sort@Mod[y + r, t]] + t, z, 
      Last@z - t}]];
  {r, Total[Abs[# - f[#][[1]]~Mod~t] & /@ x]}
  ) & /@ Range@100, Last@# &]

This assumes t to be the circumference, x to be the X list and y to be the Y list.

Results:

T          rotation   dist       Reference dist
1000       494.669    5.59189    33.0305
10000      6967.74    53.6250    429.212
100000     24452.5    2774.90    4676.23
1000000    974013.    27023.3    49149.0
10000000   9315070    496149.    497144.
                      =======
               Sum:   526006.

               Score: 1.04834

Unfortunately, T = 10⁷ takes too long to run more than 600 trials, such that I can hardly beat the reference implementation there (and even more unfortunately, this dominates the score). This still shows that there is a lot of improvement possible.

Fun fact: I actually wanted to submit a proper solution, by using one of Mathematica's minimisers on the distance function, but somehow I couldn't convince Mathematica to use my compound statement as a function to be minimised (and I wasn't patient enough to figure out why it didn't work), so I thought - well, then let's just Monte Carlo it.

Answer 4

C++, Optimal, O(|X|*|Y|*log|X|)

Where X = T / 50 and Y = T / 5. It may take about 2 days on my laptop to finish the T = 10000000 case.

Given two set of points X and Y. We denote the function f(x) to be:

f(x) = min(|x - y| for y in Y)

And then answer we want is:

dist = min(sum(f(x + t) for x in X) for t in [0, T))

However, it's impossible and unnecessary to try all the values of t. Instead we try the t values on "key points". My approach is a scanning from t = 0 to t = T with the help of a priority queue. There will have |X| * |Y| key points and trying each will need O(log|X|) operations. So the overall time complexity is O(|X|*|Y|*log|X|).

Score

               T            angle             dist         ref_dist
          1000.0     494.64422766       5.59189069      33.03052384
         10000.0    6967.55630010      45.03235069     429.21187139
        100000.0   24454.48700190     471.71127227    4676.23391624
       1000000.0  974011.50572340    4671.75610935   49148.96766310
      10000000.0  (still running...)

Code

#include <iostream>
#include <cstdio>
#include <cassert>
#include <string>
#include <sstream>
#include <cmath>
#include <vector>
#include <queue>
#include <algorithm>

using namespace std;

struct Answer {
    double angle;
    double dist;
};

double check(double T, const vector<double>& xs, const vector<double>& ys, Answer answer) {
    double dist = -answer.dist;
    double angle = answer.angle;
    for(double x : xs) {
        double x1 = x + angle;
        while(x1 > T) x1 -= T;
        int i = upper_bound(ys.begin(), ys.end(), x1) - ys.begin();
        double d;
        if(i == 0) {
            d = min(ys[0] - x1, x1 - (ys.back() - T));
        } else if(i == ys.size()) {
            d = min(ys[0] + T - x1, x1 - ys.back());
        } else {
            d = min(ys[i] - x1, x1 - ys[i - 1]);
        }
        dist += d;
    }
    return dist;
}

Answer solve(double T, const vector<double>& xs, vector<double> ys) {
    assert(ys.size() >= 1);

    sort(ys.begin(), ys.end());
    vector<double> ps;
    ps.push_back(ys.back() - T);
    for(int i = 0; i < ys.size(); i++) {
        ps.push_back((ps.back() + ys[i]) / 2);
        ps.push_back(ys[i]);
    }
    ps.push_back((ps.back() + ys.front() + T) / 2);
    ps.push_back(ys.front() + T);

    struct X {
        double npd;
        int sign;
        int j;
        bool operator<(const X& rhs)const {
            return npd > rhs.npd;
        }
    };
    priority_queue<X> que;
    int count = 0;
    Answer ans, cur;
    cur.angle = 0;
    cur.dist = 0;

    for(int i = 0; i < xs.size(); i++) {
        double x = xs[i];
        int j = upper_bound(ps.begin(), ps.end(), x) - ps.begin();
        assert(j > 0 && j < ps.size());
        X q;
        q.j = j;
        q.sign = j % 2 == 1 ? 1 : -1;
        if(q.sign == 1) {
            cur.dist += x - ps[j - 1];
        } else {
            cur.dist += ps[j] - x;
        }
        count += q.sign;
        q.npd = ps[j] - x;
        que.push(q);
    }

    ans = cur;
    while(cur.angle < T) {
        X q = que.top(); que.pop();
        double da = q.npd - cur.angle;
        cur.dist += da * count;
        cur.angle += da;
        assert(cur.dist >= 0);

        if(cur.dist < ans.dist) {
            ans = cur;
        }
        count -= q.sign * 2;
        q.sign *= -1;
        if(q.j == ps.size() - 1) {
            q.j = 3;
        } else {
            q.j ++;
        }
        int j = q.j;
        q.npd = cur.angle + ps[j] - ps[j - 1];
        que.push(q);
    }
    double d_dist = check(T, xs, ys, ans);
    cerr << "check: " << d_dist << endl;
    ans.dist += d_dist;

    return ans;
}

vector<double> read_nums() {
    stringstream ss;
    string line;
    while(line.empty()) getline(cin, line);
    ss << line;
    vector<double> xs;
    while(ss) {
        double x;
        if(ss >> x)
            xs.push_back(x);
    }
    return xs;
}

int main() {
#ifdef DEBUG
    freopen("1000.txt", "r", stdin);
#endif
    double T;
    cin >> T;
    vector<double> xs = read_nums();
    vector<double> ys = read_nums();
    double ref_dist;
    cin >> ref_dist;

    Answer ans = solve(T, xs, ys);

    printf("%16s %16s %16s %16s\n", "T", "angle", "dist", "ref_dist");
    printf("%16.1f %16.8f %16.8f %16.8f\n", T, ans.angle, ans.dist, ref_dist);
}

A bit of details

The function f(x) is something like:

So the graph of sum d(t) = sum(f(x + t) for x in X) is will be a set of connected segments. We need to try each t on the connecting points, with the help of a priority queue.