Smallest Circle Problem
Problem
The Smallest Circle Problem (also known as the Minimum Covering Circle Problem or Smallest Enclosing Circle Problem) is a computational geometry problem that consists in computing the smallest circle containing a given set of points in the Euclidian plane.
This problem is an example of a Facility Location Problem (the 1-Center Problem) in which the location of a new facility must be chosen to provide service to a number of customers, minimizing the farthest distance that any customer must travel to reach the new facility.
Principles learned
- Add float decision variables to model the coordinates of the circle’s origin
- Use the non-linear operators ‘sqrt’ and ‘pow’ to compute the circle’s radius
- Learn about Hexaly Optimizer’s modeling style: distinguish decision variables from intermediate expressions
Data
The format of the data files is as follows:
- First line: number of points
- Next lines: for each point, its x and y coordinates.
Program
The Hexaly model for the Smallest Circle Problem uses two float decision variables, representing the abscissa and ordinate of the circle’s origin. Using the nonlinear operators ‘sqrt’, ‘min’, and ‘pow’, we can compute the circle’s radius as the minimum of all distances from its origin to one of the points. Indeed, the radius does not need to be a decision variable. Since its value can be computed from that of the other decisions, it is simply an intermediate expression.
The objective function is to minimize the circle’s radius.
- 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\pythonpython smallest_circle.py instances\10points.txt
- Execution (Linux)
-
export PYTHONPATH=/opt/hexaly_13_0/bin/pythonpython 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.libsmallest_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.dllSmallestCircle 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.jarjava -cp %HX_HOME%\bin\hexaly.jar;. SmallestCircle instances\10points.txt
- Compilation / Execution (Linux)
-
javac SmallestCircle.java -cp /opt/hexaly_13_0/bin/hexaly.jarjava -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);
}
}
}