Smallest Circle

Principles learned

  • Add float decision variables

  • Define the float bounds from the data

  • Define the actual set of decision variables

  • Create a non-linear expression with operators “sqrt” and “pow”

Problem

../_images/smallest_circle.svg

Given a set of points in the plane, find the circle with minimal radius which contains all of them.

For more details, see: problem.html.

Download the example




Data

Each data file contains:

  • number of points

  • x and y coordinates of each point

Program

The decision variables in the model are x and y, respectively equal to the abscissa and the ordinate of the origin of the circle. The radius to minimize is deduced as the maximum distance between the origin and each point.

Execution:
hexaly smallest_circle.hxm inFileName=instances/10points.txt [hxTimeLimit=] [solFileName=]
use io;

/* Read instance data */
function input() {
    usage = "Usage: hexaly smallest_circle.hxm "
            + "inFileName=inputFile [solFileName=outputFile] [hxTimeLimit=timeLimit]";

    if (inFileName == nil) throw usage;
    local inFile = io.openRead(inFileName);

    nbPoints = inFile.readInt();
    for [i in 0...nbPoints] {
        coordX[i] = inFile.readInt();
        coordY[i] = inFile.readInt();
    }

    minX = min[i in 0...nbPoints](coordX[i]);
    minY = min[i in 0...nbPoints](coordY[i]);
    maxX = max[i in 0...nbPoints](coordX[i]);
    maxY = max[i in 0...nbPoints](coordY[i]);
}

/* Declare the optimization model */
function model() {

    // x, y are respectively the abscissa and the ordinate of the origin of the circle
    x <- float(minX, maxX);
    y <- float(minY, maxY);

    // Minimize the radius
    r <- sqrt(max[i in 0...nbPoints](pow(x - coordX[i], 2) + pow(y - coordY[i], 2)));
    minimize r;
}

/* Parametrize the solver */
function param() {
    if (hxTimeLimit == nil) hxTimeLimit = 6;
}

/* Write the solution in a file */
function output() {
    if (solFileName != nil) {
        println("Write solution into file '" + solFileName + "'");
        local solFile = io.openWrite(solFileName);
        solFile.println("x=", x.value);
        solFile.println("y=", y.value);
        solFile.println("r=", r.value);
    }
}
Execution (Windows)
set PYTHONPATH=%HX_HOME%\bin\python
python smallest_circle.py instances\10points.txt
Execution (Linux)
export PYTHONPATH=/opt/hexaly_13_0/bin/python
python smallest_circle.py instances/10points.txt
import hexaly.optimizer
import sys

if len(sys.argv) < 2:
    print("Usage: python smallest_circle.py inputFile [outputFile] [timeLimit]")
    sys.exit(1)


def read_integers(filename):
    with open(filename) as f:
        return [int(elem) for elem in f.read().split()]


with hexaly.optimizer.HexalyOptimizer() as optimizer:
    #
    # Read instance data
    #
    file_it = iter(read_integers(sys.argv[1]))
    # Number of points
    nb_points = next(file_it)

    # Point coordinates
    coord_x = [None] * nb_points
    coord_y = [None] * nb_points

    coord_x[0] = next(file_it)
    coord_y[0] = next(file_it)

    # Minimum and maximum value of the coordinates of the points
    min_x = coord_x[0]
    max_x = coord_x[0]
    min_y = coord_y[0]
    max_y = coord_y[0]

    for i in range(1, nb_points):
        coord_x[i] = next(file_it)
        coord_y[i] = next(file_it)
        if coord_x[i] < min_x:
            min_x = coord_x[i]
        else:
            if coord_x[i] > max_x:
                max_x = coord_x[i]
        if coord_y[i] < min_y:
            min_y = coord_y[i]
        else:
            if coord_y[i] > max_y:
                max_y = coord_y[i]

    #
    # Declare the optimization model
    #
    model = optimizer.model

    # x, y are respectively the abscissa and the ordinate of the origin of the circle
    x = model.float(min_x, max_x)
    y = model.float(min_y, max_y)

    # Distance between the origin and the point i
    radius = [(x - coord_x[i]) ** 2 + (y - coord_y[i]) ** 2 for i in range(nb_points)]

    # Minimize the radius r
    r = model.sqrt(model.max(radius))
    model.minimize(r)

    model.close()

    # Parameterize the optimizer
    if len(sys.argv) >= 4:
        optimizer.param.time_limit = int(sys.argv[3])
    else:
        optimizer.param.time_limit = 6

    optimizer.solve()

    #
    # Write the solution in a file
    #
    if len(sys.argv) >= 3:
        with open(sys.argv[2], 'w') as f:
            f.write("x=%f\n" % x.value)
            f.write("y=%f\n" % y.value)
            f.write("r=%f\n" % r.value)
Compilation / Execution (Windows)
cl /EHsc smallest_circle.cpp -I%HX_HOME%\include /link %HX_HOME%\bin\hexaly130.lib
smallest_circle instances\10points.txt
Compilation / Execution (Linux)
g++ smallest_circle.cpp -I/opt/hexaly_13_0/include -lhexaly130 -lpthread -o smallest_circle
./smallest_circle instances/10points.txt
#include "optimizer/hexalyoptimizer.h"
#include <fstream>
#include <iostream>
#include <sstream>
#include <vector>

using namespace hexaly;
using namespace std;

class SmallestCircle {
public:
    // Number of points
    int nbPoints;

    // Point coordinates
    vector<int> coordX;
    vector<int> coordY;

    // Minimum and maximum value of the coordinates of the points
    hxdouble minX;
    hxdouble minY;
    hxdouble maxX;
    hxdouble maxY;

    // Hexaly Optimizer
    HexalyOptimizer optimizer;

    // Hexaly Program variables
    HxExpression x;
    HxExpression y;

    // Objective
    HxExpression r;

    // Read instance data
    void readInstance(const string& fileName) {
        ifstream infile;
        infile.exceptions(ifstream::failbit | ifstream::badbit);
        infile.open(fileName.c_str());

        infile >> nbPoints;

        coordX.resize(nbPoints);
        coordY.resize(nbPoints);
        infile >> coordX[0];
        infile >> coordY[0];

        minX = coordX[0];
        maxX = coordX[0];
        minY = coordY[0];
        maxY = coordY[0];

        for (int i = 1; i < nbPoints; ++i) {
            infile >> coordX[i];
            infile >> coordY[i];
            if (coordX[i] < minX)
                minX = coordX[i];
            else if (coordX[i] > maxX)
                maxX = coordX[i];
            if (coordY[i] < minY)
                minY = coordY[i];
            else if (coordY[i] > maxY)
                maxY = coordY[i];
        }
    }

    void solve(int limit) {
        // Declare the optimization model
        HxModel model = optimizer.getModel();

        // x, y are respectively the abscissa and the ordinate of the origin of the circle
        x = model.floatVar(minX, maxX);
        y = model.floatVar(minY, maxY);

        // Distance between the origin and the point i
        vector<HxExpression> radius(nbPoints);
        for (int i = 0; i < nbPoints; ++i) {
            radius[i] = model.pow(x - coordX[i], 2) + model.pow(y - coordY[i], 2);
        }

        // Minimize the radius r
        r = model.sqrt(model.max(radius.begin(), radius.end()));

        model.minimize(r);
        model.close();

        // Parametrize the optimizer
        optimizer.getParam().setTimeLimit(limit);

        optimizer.solve();
    }

    /* Write the solution in a file */
    void writeSolution(const string& fileName) {
        ofstream outfile;
        outfile.exceptions(ofstream::failbit | ofstream::badbit);
        outfile.open(fileName.c_str());

        outfile << "x=" << x.getDoubleValue() << endl;
        outfile << "y=" << y.getDoubleValue() << endl;
        outfile << "r=" << r.getDoubleValue() << endl;
    }
};

int main(int argc, char** argv) {
    if (argc < 2) {
        cerr << "Usage: smallest_circle inputFile [outputFile] [timeLimit]" << endl;
        return 1;
    }

    const char* instanceFile = argv[1];
    const char* solFile = argc > 2 ? argv[2] : NULL;
    const char* strTimeLimit = argc > 3 ? argv[3] : "6";

    try {
        SmallestCircle model;
        model.readInstance(instanceFile);
        model.solve(atoi(strTimeLimit));
        if (solFile != NULL)
            model.writeSolution(solFile);
        return 0;
    } catch (const exception& e) {
        cerr << "An error occurred: " << e.what() << endl;
        return 1;
    }
}
Compilation / Execution (Windows)
copy %HX_HOME%\bin\Hexaly.NET.dll .
csc SmallestCircle.cs /reference:Hexaly.NET.dll
SmallestCircle instances\10points.txt
using System;
using System.IO;
using Hexaly.Optimizer;

public class SmallestCircle : IDisposable
{
    // Number of points
    int nbPoints;

    // Point coordinates
    double[] coordX;
    double[] coordY;

    // Minimum and maximum value of the coordinates of the points
    double minX;
    double minY;
    double maxX;
    double maxY;

    // Hexaly Optimizer
    HexalyOptimizer optimizer;

    // Hexaly Program variables
    HxExpression x;
    HxExpression y;

    // Objective
    HxExpression r;

    public SmallestCircle()
    {
        optimizer = new HexalyOptimizer();
    }

    // Read instance data
    public void ReadInstance(string fileName)
    {
        using (StreamReader input = new StreamReader(fileName))
        {
            nbPoints = int.Parse(input.ReadLine());
            coordX = new double[nbPoints];
            coordY = new double[nbPoints];

            string[] splittedCoord = input.ReadLine().Split(' ');
            coordX[0] = int.Parse(splittedCoord[0]);
            coordY[0] = int.Parse(splittedCoord[1]);

            minX = coordX[0];
            maxX = coordX[0];
            minY = coordY[0];
            maxY = coordY[0];

            for (int i = 1; i < nbPoints; ++i)
            {
                splittedCoord = input.ReadLine().Split(' ');
                coordX[i] = int.Parse(splittedCoord[0]);
                coordY[i] = int.Parse(splittedCoord[1]);

                minX = Math.Min(coordX[i], minX);
                maxX = Math.Max(coordX[i], maxX);
                minY = Math.Min(coordY[i], minY);
                maxY = Math.Max(coordY[i], maxY);
            }
        }
    }

    public void Dispose()
    {
        if (optimizer != null)
            optimizer.Dispose();
    }

    public void Solve(int limit)
    {
        // Declare the optimization model
        HxModel model = optimizer.GetModel();

        // x, y are respectively the abscissa and the ordinate of the origin of the circle
        x = model.Float(minX, maxX);
        y = model.Float(minY, maxY);

        // Distance between the origin and the point i
        HxExpression[] radius = new HxExpression[nbPoints];
        for (int i = 0; i < nbPoints; ++i)
            radius[i] = model.Pow(x - coordX[i], 2) + model.Pow(y - coordY[i], 2);

        // Minimize the radius r
        r = model.Sqrt(model.Max(radius));

        model.Minimize(r);
        model.Close();

        // Parametrize the optimizer
        optimizer.GetParam().SetTimeLimit(limit);

        optimizer.Solve();
    }

    /* Write the solution in a file */
    public void WriteSolution(string fileName)
    {
        using (StreamWriter output = new StreamWriter(fileName))
        {
            output.WriteLine("x=" + x.GetDoubleValue());
            output.WriteLine("y=" + y.GetDoubleValue());
            output.WriteLine("r=" + r.GetDoubleValue());
        }
    }

    public static void Main(string[] args)
    {
        if (args.Length < 1)
        {
            Console.WriteLine("Usage: SmallestCircle inputFile [outputFile] [timeLimit]");
            Environment.Exit(1);
        }

        string instanceFile = args[0];
        string outputFile = args.Length > 1 ? args[1] : null;
        string strTimeLimit = args.Length > 2 ? args[2] : "6";

        using (SmallestCircle model = new SmallestCircle())
        {
            model.ReadInstance(instanceFile);
            model.Solve(int.Parse(strTimeLimit));
            if (outputFile != null)
                model.WriteSolution(outputFile);
        }
    }
}
Compilation / Execution (Windows)
javac SmallestCircle.java -cp %HX_HOME%\bin\hexaly.jar
java -cp %HX_HOME%\bin\hexaly.jar;. SmallestCircle instances\10points.txt
Compilation / Execution (Linux)
javac SmallestCircle.java -cp /opt/hexaly_13_0/bin/hexaly.jar
java -cp /opt/hexaly_13_0/bin/hexaly.jar:. SmallestCircle instances/10points.txt
import java.util.*;
import java.io.*;
import com.hexaly.optimizer.*;

public class SmallestCircle {
    // Number of points
    private int nbPoints;

    // Point coordinates
    private int[] coordX;
    private int[] coordY;

    // Minimum and maximum value of the coordinates of the points
    private int minX;
    private int minY;
    private int maxX;
    private int maxY;

    // Hexaly Optimizer
    private final HexalyOptimizer optimizer;

    // Hexaly Program variables
    private HxExpression x;
    private HxExpression y;

    // Objective i
    private HxExpression r;

    private SmallestCircle(HexalyOptimizer optimizer) {
        this.optimizer = optimizer;
    }

    // Read instance data
    private void readInstance(String fileName) throws IOException {
        try (Scanner input = new Scanner(new File(fileName))) {
            nbPoints = input.nextInt();

            coordX = new int[nbPoints];
            coordY = new int[nbPoints];

            coordX[0] = input.nextInt();
            coordY[0] = input.nextInt();
            minX = coordX[0];
            maxX = coordX[0];
            minY = coordY[0];
            maxY = coordY[0];

            for (int i = 1; i < nbPoints; ++i) {
                coordX[i] = input.nextInt();
                coordY[i] = input.nextInt();
                minX = Math.min(coordX[i], minX);
                maxX = Math.max(coordX[i], maxX);
                minY = Math.min(coordY[i], minY);
                maxY = Math.max(coordY[i], maxY);
            }
        }
    }

    private void solve(int limit) {
        // Declare the optimization model
        HxModel model = optimizer.getModel();

        // x, y are respectively the abscissa and the ordinate of the origin of the circle
        x = model.floatVar(minX, maxX);
        y = model.floatVar(minY, maxY);

        // Distance between the origin and the point i
        HxExpression[] radius = new HxExpression[nbPoints];
        for (int i = 0; i < nbPoints; ++i) {
            radius[i] = model.sum();
            radius[i].addOperand(model.pow(model.sub(x, coordX[i]), 2));
            radius[i].addOperand(model.pow(model.sub(y, coordY[i]), 2));
        }

        // Minimize the radius r
        r = model.sqrt(model.max(radius));

        model.minimize(r);
        model.close();

        // Parametrize the optimizer
        optimizer.getParam().setTimeLimit(limit);

        optimizer.solve();
    }

    /* Write the solution in a file */
    private void writeSolution(String fileName) throws IOException {
        try (PrintWriter output = new PrintWriter(fileName)) {
            output.println("x=" + x.getDoubleValue());
            output.println("y=" + y.getDoubleValue());
            output.println("r=" + r.getDoubleValue());
        }
    }

    public static void main(String[] args) {
        if (args.length < 1) {
            System.err.println("Usage: java SmallestCircle inputFile [outputFile] [timeLimit]");
            System.exit(1);
        }

        String instanceFile = args[0];
        String outputFile = args.length > 1 ? args[1] : null;
        String strTimeLimit = args.length > 2 ? args[2] : "6";
        try (HexalyOptimizer optimizer = new HexalyOptimizer()) {
            SmallestCircle model = new SmallestCircle(optimizer);
            model.readInstance(instanceFile);
            model.solve(Integer.parseInt(strTimeLimit));
            if (outputFile != null) {
                model.writeSolution(outputFile);
            }
        } catch (Exception ex) {
            System.err.println(ex);
            ex.printStackTrace();
            System.exit(1);
        }
    }
}