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¶

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.

LSP
Python
C++
C#
Java

Execution:: localsolver smallest_circle.lsp inFileName=instances/10points.txt [lsTimeLimit=] [solFileName=]

/********** smallest_circle.lsp **********/

use io;

/* Reads instance data */
function input() {
     usage = "\nUsage: localsolver smallest_circle.lsp "
         + "inFileName=inputFile [solFileName=outputFile] [lsTimeLimit=timeLimit]\n";

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

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

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

/* Declares 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 1..nbPoints](pow(x - coordX[i], 2) + pow(y - coordY[i], 2)));
    minimize r;
}

/* Parameterizes the solver. */
function param() {
    if (lsTimeLimit == nil) lsTimeLimit = 6;
}

/* Writes the solution in a file */
function output() {
    if (solFileName != nil) // write solution file if needed
    {
        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=%LS_HOME%\bin\

python smallest_circle.py instances\10points.txt

Execution (Linux): export PYTHONPATH=/opt/localsolver_9_5/bin/

python smallest_circle.py instances/10points.txt

########## smallest_circle.py ##########

import localsolver
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 localsolver.LocalSolver() as ls:

    #
    # Reads 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]

    #
    # Declares the optimization model
    #
    model = ls.model

    # Numerical decisions
    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()

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

    ls.solve()

    #
    # Writes 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%LS_HOME%\include /link %LS_HOME%\bin\localsolver95.lib

smallest_circle instances\10points.txt

Compilation / Execution (Linux): g++ smallest_circle.cpp -I/opt/localsolver_9_5/include -llocalsolver95 -lpthread -o smallest_circle

./smallest_circle instances/10points.txt

//********* smallest_circle.cpp *********

#include <iostream>
#include <sstream>
#include <fstream>
#include <vector>
#include "localsolver.h"

using namespace localsolver;
using namespace std;

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

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

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

    // Solver. 
    LocalSolver localsolver;

    // LS Program variables 
    LSExpression x;
    LSExpression y;
    
    // Objective 
    LSExpression r;

    // Reads 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) {
        // Declares the optimization model. 
        LSModel model = localsolver.getModel();

        // Numerical decisions
        x = model.floatVar(minX, maxX);
        y = model.floatVar(minY, maxY);

        // Distance between the origin and the point i
        vector<LSExpression> 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();            

        // Parameterizes the solver. 
        localsolver.getParam().setTimeLimit(limit);

        localsolver.solve();
    }

    // Writes 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 %LS_HOME%\bin\localsolvernet.dll .

csc SmallestCircle.cs /reference:localsolvernet.dll

SmallestCircle instances\10points.txt

/********** SmallestCircle.cs **********/

using System;
using System.IO;
using localsolver;

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;

    // Solver
    LocalSolver localsolver;

    // LS Program variables
    LSExpression x;
    LSExpression y;

    // Objective
    LSExpression r;

    public SmallestCircle()
    {
        localsolver = new LocalSolver();
    }

    // Reads 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 (localsolver != null)
            localsolver.Dispose();
    }

    public void Solve(int limit)
    {
        // Declares the optimization model.
        LSModel model = localsolver.GetModel();

        // Numerical decisions
        x = model.Float(minX, maxX);
        y = model.Float(minY, maxY);

        // Distance between the origin and the point i
        LSExpression[] radius = new LSExpression[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();

        // Parameterizes the solver.
        localsolver.GetParam().SetTimeLimit(limit);

        localsolver.Solve();
    }

    // Writes 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 %LS_HOME%\bin\localsolver.jar

java -cp %LS_HOME%\bin\localsolver.jar;. SmallestCircle instances\10points.txt

Compilation / Execution (Linux): javac SmallestCircle.java -cp /opt/localsolver_9_5/bin/localsolver.jar

java -cp /opt/localsolver_9_5/bin/localsolver.jar:. SmallestCircle instances/10points.txt

/********** SmallesrCircle.java **********/

import java.util.*;
import java.io.*;
import localsolver.*;

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;

    // Solver.
    private final LocalSolver localsolver;

    // LS Program variables
    private LSExpression x;
    private LSExpression y;

    // Objective i
    private LSExpression r;

    private SmallestCircle(LocalSolver localsolver) {
        this.localsolver = localsolver;
    }

    // Reads 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) {
        // Declares the optimization model.
        LSModel model = localsolver.getModel();

        // Numerical decisions
        x = model.floatVar(minX, maxX);
        y = model.floatVar(minY, maxY);

        // Distance between the origin and the point i
        LSExpression[] radius = new LSExpression[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();

        // Parameterizes the solver.
        localsolver.getParam().setTimeLimit(limit);

        localsolver.solve();
    }

    // Writes 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 (LocalSolver localsolver = new LocalSolver()) {
            SmallestCircle model = new SmallestCircle(localsolver);
            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);
        }
    }
}