Clustered Vehicle Routing (cluVRP)

import hexaly.optimizer
import sys
import math

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

def main(instance_file, str_time_limit, output_file):
    # Read instance data
    nb_customers, nb_trucks, nb_clusters, truck_capacity, dist_matrix_data, dist_depot_data, \
        demands_data, clusters_data = read_input_cvrp(instance_file)


    with hexaly.optimizer.HexalyOptimizer() as optimizer:
        # Declare the optimization model
        model = optimizer.model

        # Create Hexaly arrays to be able to access them with an "at" operator
        demands = model.array(demands_data)
        dist_matrix = model.array(dist_matrix_data)
        clusters = model.array(clusters_data)
        dist_depot = model.array(dist_depot_data)
        
        # A list is created for each cluster, to determine the order within the cluster
        clusters_sequences = []
        for k in range(nb_clusters):
            clusters_sequences.append(model.list(len(clusters_data[k])))
            # All customers in the cluster must be visited 
            model.constraint(model.count(clusters_sequences[k]) == len(clusters_data[k]))

        clustersDistances = model.array()
        initialNodes = model.array()
        endNodes = model.array()
        for k in range(nb_clusters):
            sequence = clusters_sequences[k]
            c = model.count(sequence)
            # Distance traveled within cluster k
            clustersDistances_lambda = model.lambda_function(lambda i:
                    model.at(dist_matrix, clusters[k][sequence[i - 1]],
                    clusters[k][sequence[i]]))
            clustersDistances.add_operand(model.sum(model.range(1,c), 
                    clustersDistances_lambda))
            
            # First and last point when visiting cluster k
            initialNodes.add_operand(clusters[k][sequence[0]]) 
            endNodes.add_operand(clusters[k][sequence[c - 1]])
        
        # Sequences of clusters visited by each truck
        truckSequences = [model.list(nb_clusters) for _ in range(nb_trucks)]

        # Each cluster must be visited by exactly one truck
        model.constraint(model.partition(truckSequences))

        routeDistances = [None] * nb_trucks
        for k in range(nb_trucks):
            sequence = truckSequences[k]
            c = model.count(sequence)

            # The quantity needed in each route must not exceed the truck capacity
            demand_lambda = model.lambda_function(lambda j: demands[j])
            route_quantity = model.sum(sequence, demand_lambda)
            model.constraint(route_quantity <= truck_capacity)

            # Distance traveled by each truck
            # = distance in each cluster + distance between clusters + distance with depot at 
            # the beginning end at the end of a route
            routeDistances_lambda = model.lambda_function(lambda i:
                    model.at(clustersDistances, sequence[i]) + model.at(dist_matrix,
                    endNodes[sequence[i - 1]], initialNodes[sequence[i]]))
            routeDistances[k] = model.sum(model.range(1, c), routeDistances_lambda) \
                    + model.iif(c > 0, model.at(clustersDistances, sequence[0]) 
                    + dist_depot[initialNodes[sequence[0]]]
                    + dist_depot[endNodes[sequence[c - 1]]], 0)

        # Total distance traveled
        total_distance = model.sum(routeDistances)

        # Objective:  minimize the distance traveled
        model.minimize(total_distance)

        model.close()

        # Parameterize the optimizer
        optimizer.param.time_limit = int(str_time_limit)

        optimizer.solve() 
        # Write the solution in a file with the following format:
        # - total distance
        # - for each truck the customers visited (omitting the start/end at the depot)
        if output_file is not None:
            with open(output_file, 'w') as f:
                f.write("%d\n" % (total_distance.value))
                for k in range(nb_trucks):
                    # Values in sequence are in [0..nbCustomers - 1]. +2 is to put it back
                    # in [2..nbCustomers+1] as in the data files (1 being the depot)
                    for cluster in truckSequences[k].value:
                        for customer in clusters_sequences[cluster].value:
                            f.write("%d " % (clusters_data[cluster][customer] + 2))
                    f.write("\n")

# The input files follow the "Augerat" format
def read_input_cvrp(filename):
    file_it = iter(read_elem(filename))

    nb_nodes = 0
    while True:
        token = next(file_it)
        if token == "DIMENSION:":
            nb_nodes = int(next(file_it))
            nb_customers = nb_nodes - 1
        if token == "VEHICLES:":
            nb_trucks = int(next(file_it))
        elif token == "GVRP_SETS:":
            nb_clusters = int(next(file_it))
        elif token == "CAPACITY:":
            truck_capacity = int(next(file_it))
        elif token == "NODE_COORD_SECTION":
            break

    customers_x = [None] * nb_customers
    customers_y = [None] * nb_customers
    depot_x = 0
    depot_y = 0
    for n in range(nb_nodes):
        node_id = int(next(file_it))
        if node_id != n + 1:
            print("Unexpected index")
            sys.exit(1)
        if node_id == 1:
            depot_x = int(float(next(file_it)))
            depot_y = int(float(next(file_it)))
        else:
            # -2 because original customer indices are in 2..nbNodes
            customers_x[node_id - 2] = int(float(next(file_it)))
            customers_y[node_id - 2] = int(float(next(file_it)))

    distance_matrix = compute_distance_matrix(customers_x, customers_y)
    distance_depots = compute_distance_depots(depot_x, depot_y, customers_x, customers_y)
    
    token = next(file_it)
    if token != "GVRP_SET_SECTION":
        print("Expected token GVRP_SET_SECTION")
        sys.exit(1)
    clusters_data = [None]*nb_clusters
    for n in range(nb_clusters):
        node_id = int(next(file_it))
        if node_id != n + 1:
            print("Unexpected index")
            sys.exit(1)
        cluster = []
        value = int(next(file_it))
        while value != -1:
            # -2 because original customer indices are in 2..nbNodes
            cluster.append(value-2)
            value = int(next(file_it))
        clusters_data[n] = cluster
    token = next(file_it)
    if token != "DEMAND_SECTION":
        print("Expected token DEMAND_SECTION")
        sys.exit(1)

    demands = [None] * nb_clusters
    for n in range(nb_clusters):
        node_id = int(next(file_it))
        if node_id != n + 1:
            print("Unexpected index")
            sys.exit(1)
        demands[n] = int(next(file_it))
    return nb_customers, nb_trucks, nb_clusters, truck_capacity, distance_matrix, \
        distance_depots, demands, clusters_data

# Compute the distance matrix
def compute_distance_matrix(customers_x, customers_y):
    nb_customers = len(customers_x)
    distance_matrix = [[None for i in range(nb_customers)] for j in range(nb_customers)]
    for i in range(nb_customers):
        distance_matrix[i][i] = 0
        for j in range(nb_customers):
            dist = compute_dist(customers_x[i], customers_x[j], customers_y[i], customers_y[j])
            distance_matrix[i][j] = dist
            distance_matrix[j][i] = dist
    return distance_matrix

# Compute the distances to depot
def compute_distance_depots(depot_x, depot_y, customers_x, customers_y):
    nb_customers = len(customers_x)
    distance_depots = [None] * nb_customers
    for i in range(nb_customers):
        dist = compute_dist(depot_x, customers_x[i], depot_y, customers_y[i])
        distance_depots[i] = dist
    return distance_depots


def compute_dist(xi, xj, yi, yj):
    exact_dist = math.sqrt(math.pow(xi - xj, 2) + math.pow(yi - yj, 2))
    return int(math.floor(exact_dist + 0.5))


if __name__ == '__main__':
    if len(sys.argv) < 2:
        print("Usage: python clustered-vehicle-routing.py input_file [output_file] [time_limit]")
        sys.exit(1)

    instance_file = sys.argv[1]
    output_file = sys.argv[2] if len(sys.argv) > 2 else None
    str_time_limit = sys.argv[3] if len(sys.argv) > 3 else "20"

    main(instance_file, str_time_limit, output_file)

#include "optimizer/hexalyoptimizer.h"
#include <cmath>
#include <cstring>
#include <fstream>
#include <iostream>
#include <vector>

using namespace hexaly;
using namespace std;

class ClusteredCvrp {
public:
    // Hexaly Optimizer
    HexalyOptimizer optimizer;

    // Number of customers
    int nbCustomers;

    // Capacity of the trucks
    int truckCapacity;

    // Demand on each cluster
    vector<int> demandsData;

    // Customers in each cluster;
    vector<vector<int>> clustersData;

    // Distance matrix between customers
    vector<vector<int>> distMatrixData;

    // Distances between customers and depot
    vector<int> distDepotData;

    // Number of trucks
    int nbTrucks;

    // Number of clusters
    int nbClusters;

    // Decision variables
    vector<HxExpression> truckSequences;
    vector<HxExpression> clustersSequences;

    // Are the trucks actually used
    vector<HxExpression> trucksUsed;

    // Number of trucks used in the solution
    HxExpression nbTrucksUsed;

    // Distance traveled by all the trucks
    HxExpression totalDistance;

    /* Read instance data */
    void readInstance(const string& fileName) {
        readInputCvrp(fileName);
    }

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

        // Create HexalyOptimizer arrays to be able to access them with an "at" operator
        HxExpression demands = model.array(demandsData.begin(), demandsData.end());
        HxExpression distMatrix = model.array();
        for (int n = 0; n < nbCustomers; ++n) {
            distMatrix.addOperand(model.array(distMatrixData[n].begin(), 
                    distMatrixData[n].end()));
        }
        HxExpression clusters = model.array();
        for (int n = 0; n < clustersData.size(); ++n) {
            clusters.addOperand(model.array(clustersData[n].begin(), clustersData[n].end()));
        }
        HxExpression distDepot = model.array(distDepotData.begin(), distDepotData.end());
        
        // A list is created for each cluster, to determine the order within the cluster
        clustersSequences.resize(nbClusters);
        for (int k = 0; k < nbClusters; ++k) {
            int c = (int) clustersData[k].size();
            clustersSequences[k] = model.listVar(c);
            // All customers in the cluster must be visited 
            model.constraint(model.count(clustersSequences[k]) == c);
        }

        HxExpression clustersDistances =  model.array();
        HxExpression initialNodes = model.array();
        HxExpression endNodes = model.array();
        for (int k = 0; k < nbClusters; ++k) {
            HxExpression sequence = clustersSequences[k];
            HxExpression c = model.count(sequence);

            // Distance traveled within clsuter k
            HxExpression clustersDistances_lambda = model.createLambdaFunction(
                    [&](HxExpression i) { return model.at(distMatrix,
                    clusters[k][sequence[i - 1]], clusters[k][sequence[i]]); });
            clustersDistances.addOperand(model.sum(model.range(1, c), clustersDistances_lambda));

            // First and last point when visiting cluster k
            initialNodes.addOperand(clusters[k][sequence[0]]);
            endNodes.addOperand(clusters[k][sequence[c - 1]]);
        }

        // Sequence of clusters visited by each truck
        truckSequences.resize(nbTrucks);
        for (int k = 0; k < nbTrucks; ++k) {
            truckSequences[k] = model.listVar(nbClusters);
        }
        // All clusters must be visited by the trucks
        model.constraint(model.partition(truckSequences.begin(), truckSequences.end()));
        
        vector<HxExpression> routeDistances(nbTrucks);
        for (int k = 0; k < nbTrucks; ++k) {
            HxExpression sequence = truckSequences[k];
            HxExpression c = model.count(sequence);

            // The quantity needed in each route must not exceed the truck capacity
            HxExpression demandLambda =
                    model.createLambdaFunction([&](HxExpression j) { return demands[j]; });
            HxExpression routeQuantity = model.sum(sequence, demandLambda);
            model.constraint(routeQuantity <= truckCapacity);

            // Distance traveled by truck k
            // = distance in each cluster + distance between clusters + distance with depot 
            // at the beginning end at the end of a route
            HxExpression routeDistances_lambda = model.createLambdaFunction(
                    [&](HxExpression i) { return model.at(clustersDistances, sequence[i]) 
                    + model.at(distMatrix, endNodes[sequence[i - 1]], initialNodes[sequence[i]]);}
            );
            
            routeDistances[k] = model.sum(model.range(1, c), routeDistances_lambda)
                    + model.iif(c > 0, model.at(clustersDistances, sequence[0]) 
                    + distDepot[initialNodes[sequence[0]]] 
                    + distDepot[endNodes[sequence[c - 1]]], 0);
        }

        // Total distance traveled
        totalDistance = model.sum(routeDistances.begin(), routeDistances.end());

        // Objective: minimize the distance traveled
        model.minimize(totalDistance);

        model.close();

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

    /* Write the solution in a file with the following format:
     * - total distance
     * - for each truck the customers visited (omitting the start/end at the depot) */
    void writeSolution(const string& fileName) {
        ofstream outfile;
        outfile.exceptions(ofstream::failbit | ofstream::badbit);
        outfile.open(fileName.c_str());

        outfile << totalDistance.getValue() << endl;
        for (int k = 0; k < nbTrucks; ++k) {
            // Values in sequence are in [0..nbCustomers-1]. +2 is to put it back in 
            // [2..nbCustomers+1] as in the data files (1 being the depot)
            HxCollection customersCollection = truckSequences[k].getCollectionValue();
            for (int i = 0; i < customersCollection.count(); ++i) {
                int cluster = customersCollection[i];
                HxCollection clustersCollection = 
                        clustersSequences[cluster].getCollectionValue();
                for (int j = 0; j < clustersCollection.count(); ++j)
                    outfile << clustersData[cluster][clustersCollection[j]] + 2 << " ";
            }
            outfile << endl;
        }
    }

private:
    // The input files follow the "Augerat" format
    void readInputCvrp(const string& fileName) {
        ifstream infile(fileName.c_str());
        if (!infile.is_open()) {
            throw std::runtime_error("File cannot be opened.");
        }
        string str;
        char* pch;
        char* line;
        int nbNodes;
        while (true) {
            getline(infile, str);
            line = strdup(str.c_str());
            pch = strtok(line, " ");
            if (strcmp(pch, "DIMENSION:") == 0) {
                pch = strtok(NULL, " ");
                nbNodes = atoi(pch);
                nbCustomers = nbNodes - 1;
            } else if (strcmp(pch, "VEHICLES:") == 0) {
                pch = strtok(NULL, " ");
                nbTrucks = atoi(pch);
            } else if (strcmp(pch, "GVRP_SETS:") == 0) {
                pch = strtok(NULL, " ");
                nbClusters = atoi(pch);
            } else if (strcmp(pch, "CAPACITY:") == 0) {
                pch = strtok(NULL, "");
                truckCapacity = atoi(pch);
            } else if (strcmp(pch, "NODE_COORD_SECTION") == 0) {
                break;
            }
        }

        vector<int> customersX(nbCustomers);
        vector<int> customersY(nbCustomers);
        int depotX, depotY;
        for (int n = 1; n <= nbNodes; ++n) {
            int id;
            infile >> id;
            if (id != n) {
                throw std::runtime_error("Unexpected index");
            }
            if (n == 1) {
                infile >> depotX;
                infile >> pch;
                infile >> depotY;
                infile >> pch;
            } else {
                // -2 because original customer indices are in 2..nbNodes
                infile >> customersX[n - 2];
                infile >> pch;
                infile >> customersY[n - 2];
                infile >> pch;
            }
        }

        computeDistanceMatrix(depotX, depotY, customersX, customersY);

        getline(infile, str); // End the last line
        getline(infile, str);
        line = strdup(str.c_str());
        pch = strtok(line, " ");
        if (strcmp(pch, "GVRP_SET_SECTION") != 0) {
            throw std::runtime_error("Expected keyword GVRP_SET_SECTION");
        }
        for (int n = 1; n <= nbClusters; ++n) {
            vector<int> cluster;
            int id;
            infile >> id;
            if (id != n) {
                throw std::runtime_error("Unexpected index");
            }
            int data;
            infile >> data;
            while (data != -1) {
                // -2 because original customer indices are in 2..nbNodes
                cluster.push_back(data - 2);
                infile >> data;
            }
            clustersData.push_back(cluster);

        };

        getline(infile, str); // End the last line
        getline(infile, str);
        line = strdup(str.c_str());
        pch = strtok(line, " ");
        if (strcmp(pch, "DEMAND_SECTION") != 0) {
            throw std::runtime_error("Expected keyword DEMAND_SECTION");
        }
        demandsData.resize(nbClusters);
        for (int n = 1; n <= nbClusters; ++n) {
            int id;
            infile >> id;
            if (id != n) {
                throw std::runtime_error("Unexpected index");
            }
            int demand;
            infile >> demand;
            demandsData[n - 1] = demand;
        }
        infile.close();

    }

    // Compute the distance matrix
    void computeDistanceMatrix(int depotX, int depotY, const vector<int>& customersX, const vector<int>& customersY) {
        distMatrixData.resize(nbCustomers);
        for (int i = 0; i < nbCustomers; ++i) {
            distMatrixData[i].resize(nbCustomers);
        }
        for (int i = 0; i < nbCustomers; ++i) {
            distMatrixData[i][i] = 0;
            for (int j = i + 1; j < nbCustomers; ++j) {
                int distance = computeDist(customersX[i], customersX[j], customersY[i], customersY[j]);
                distMatrixData[i][j] = distance;
                distMatrixData[j][i] = distance;
            }
        }

        distDepotData.resize(nbCustomers);
        for (int i = 0; i < nbCustomers; ++i) {
            distDepotData[i] = computeDist(depotX, customersX[i], depotY, customersY[i]);
        }
    }

    int computeDist(int xi, int xj, int yi, int yj) {
        double exactDist = sqrt(pow((double)xi - xj, 2) + pow((double)yi - yj, 2));
        return floor(exactDist + 0.5);
    }

};

int main(int argc, char** argv) {
    if (argc < 1) {
        cerr << "Usage: clustered-vehicle-routing inputFile [outputFile] [timeLimit] " << endl;
        return 1;
    }
    const char* instanceFile = argc > 1 ? argv[1] : "instances/A-n32-k5-C11-V2.gvrp";
    const char* solFile = argc > 2 ? argv[2] : NULL;
    const char* strTimeLimit = argc > 3 ? argv[3] : "5";
    try {
        ClusteredCvrp 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;
    }
}

using System;
using System.IO;
using System.Collections.Generic;
using Hexaly.Optimizer;

public class cluCvrp : IDisposable
{
    // Hexaly Optimizer
    HexalyOptimizer optimizer;

    // Number of customers
    int nbCustomers;

    // Capacity of the trucks
    int truckCapacity;

    // Demand on each customer
    long[] demandsData;

    // Customers in each cluster 
    long[][] clustersData;

    // Distance matrix between customers
    long[][] distMatrixData;

    // Distances between customers and depot
    long[] distDepotData;

    // Number of trucks
    int nbTrucks;

    // Number of Clusters
    int nbClusters;

    // Decision variables
    HxExpression[] truckSequences;
    HxExpression[] clustersSequences;


    // Distance traveled by all the trucks
    HxExpression totalDistance;

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

    /* Read instance data */
    void ReadInstance(string fileName)
    {
        ReadInputcluCvrp(fileName);
    }

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

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

        HxExpression[] clustersDistances = new HxExpression[nbClusters];
        HxExpression[] routeDistances = new HxExpression[nbTrucks];
        HxExpression[] initialNodes = new HxExpression[nbClusters];
        HxExpression[] endNodes = new HxExpression[nbClusters];

        clustersSequences = new HxExpression[nbClusters];
        // A list is created for each cluster, to determine the order within the cluster
        for ( int k = 0; k < nbClusters; ++k)
        {
            int c = (int) clustersData[k].Length;
            clustersSequences[k] = model.List(c);
            // All customers in the cluster must be visited 
            model.Constraint(model.Count(clustersSequences[k]) == c);
        }

        // Create HexalyOptimizer arrays to be able to access them with an "at" operator
        HxExpression demands = model.Array(demandsData);
        HxExpression distDepot = model.Array(distDepotData);
        HxExpression distMatrix = model.Array(distMatrixData);
        HxExpression clusters = model.Array(clustersData);
        
        for (int k = 0; k < nbClusters; ++k)
        {
            HxExpression sequence = clustersSequences[k];
            HxExpression c = model.Count(sequence);

            HxExpression routeDistances_lambda = model.LambdaFunction(
                i => model.At(distMatrix,
                     clusters[k][sequence[i - 1]], clusters[k][sequence[i]]) );
            
            // Distance traveled within cluster k
            clustersDistances[k] = model.Sum(model.Range(1,c), routeDistances_lambda);

            // first and last point visited when traveled into cluster k
            initialNodes[k] = clusters[k][sequence[0]];
            endNodes[k] = clusters[k][sequence[c - 1]];
        }

        truckSequences = new HxExpression[nbTrucks];
        // Sequence of clusters visited by each truck
        for (int k = 0; k < nbTrucks; ++k)
            truckSequences[k] = model.List(nbClusters);

        // All customers must be visited by the trucks
        model.Constraint(model.Partition(truckSequences));

        HxExpression valueDistClusters = model.Array(clustersDistances);
        HxExpression initials = model.Array(initialNodes);
        HxExpression ends = model.Array(endNodes);
        for (int k = 0; k < nbTrucks; ++k)
        {
            HxExpression sequence = truckSequences[k];
            HxExpression c = model.Count(sequence);

            // The quantity needed in each route must not exceed the truck capacity
            HxExpression demandLambda = model.LambdaFunction(j => demands[j]);
            HxExpression routeQuantity = model.Sum(sequence, demandLambda);
            model.Constraint(routeQuantity <= truckCapacity);

            // Distance traveled by truck k
            // = distance in each cluster + distance between clusters + distance with depot 
            // at the beginning end at the end of a route
            HxExpression routeDistances_lambda = model.LambdaFunction(
                    i => model.Sum( model.At(valueDistClusters, sequence[i]),
                    model.At(distMatrix, ends[sequence[i - 1]], initials[sequence[i]])));
            routeDistances[k] = model.Sum(model.Range(1, c), routeDistances_lambda)
                    + model.If(c > 0, valueDistClusters[ model.At(sequence,0)]
                    + distDepot[initials[sequence[0]]] + distDepot[ends[sequence[c - 1]]], 0);
        }

        totalDistance = model.Sum(routeDistances);

        // Objective:  minimize the distance traveled
        model.Minimize(totalDistance);

        model.Close();

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

        optimizer.Solve();
    }

    /* Write the solution in a file with the following format:
     * -  total distance
     * - for each truck the customers visited (omitting the start/end at the depot) */
    void WriteSolution(string fileName)
    {
        using (StreamWriter output = new StreamWriter(fileName))
        {
            output.WriteLine(totalDistance.GetValue());
            for (int k = 0; k < nbTrucks; ++k)
            {
                // Values in sequence are in [0..nbCustomers-1]. +2 is to put it back in 
                // [2..nbCustomers+1] as in the data files (1 being the depot)
                HxCollection customersCollection = truckSequences[k].GetCollectionValue();
                for (int i = 0; i < customersCollection.Count(); ++i)
                {
                    long cluster = customersCollection[i];
                    HxCollection clustersCollection = 
                        clustersSequences[cluster].GetCollectionValue();
                    for (int j = 0; j < clustersCollection.Count(); ++j)
                        output.Write((clustersData[cluster][clustersCollection[j]] + 2) + " ");
                }
                output.WriteLine();
            }
        }
    }

    public static void Main(string[] args)
    {
        if (args.Length < 1)
        {
            Console.WriteLine("Usage: clustered-vehicle-routing.cs inputFile [solFile] [timeLimit]");
            Environment.Exit(1);
        }
        string instanceFile = args[0];
        string outputFile = args.Length > 1 ? args[1] : null;
        string strTimeLimit = args.Length > 2 ? args[2] : "5";
        string strNbTrucks = args.Length > 3 ? args[3] : "0";

        using (cluCvrp model = new cluCvrp())
        {
            model.ReadInstance(instanceFile);
            model.Solve(int.Parse(strTimeLimit));
            if (outputFile != null)
                model.WriteSolution(outputFile);
        }
    }

    // The input files follow the "Augerat" format
    private void ReadInputcluCvrp(string fileName)
    {
        using (StreamReader input = new StreamReader(fileName))
        {   
            int nbNodes = 0;
            string[] splitted;
            while (true) 
            {
                splitted = input.ReadLine().Split(':');
                if (splitted[0].Contains("DIMENSION"))
                {
                    nbNodes = int.Parse(splitted[1]);
                    nbCustomers = nbNodes - 1;
                }
                else if (splitted[0].Contains("VEHICLES"))
                    nbTrucks = int.Parse(splitted[1]);
                else if (splitted[0].Contains("GVRP_SETS"))
                    nbClusters = int.Parse(splitted[1]);
                else if (splitted[0].Contains("CAPACITY"))
                    truckCapacity = int.Parse(splitted[1]);
                else if (splitted[0].Contains("NODE_COORD_SECTION"))
                    break;
            }
            int[] customersX = new int[nbCustomers];
            int[] customersY = new int[nbCustomers];
            int depotX = 0,
                depotY = 0;
            char[] delimiterChars = { ' ', '.' };
            for (int n = 1; n <= nbNodes; ++n)
            {
                splitted = input.ReadLine().Split(delimiterChars);
                int id = int.Parse(splitted[0]);
                if ( id != n)
                    throw new Exception("Unexpected index");
                if (n == 1)
                {
                    depotX = int.Parse(splitted[1]);
                    depotY = int.Parse(splitted[3]);
                } 
                else
                {
                    // -2 because original customer indices are in 2..nbNodes
                    customersX[n - 2] = int.Parse(splitted[1]);
                    customersY[n - 2] = int.Parse(splitted[3]);
                }
            }

            ComputeDistanceMatrix(depotX, depotY, customersX, customersY);

            splitted = input.ReadLine().Split(' ');
            if (!splitted[0].Contains("GVRP_SET_SECTION"))
                throw new Exception("Expected keyword GVRP_SET_SECTION");

            clustersData = new long[nbClusters][];
            for (int n = 1; n <= nbClusters; ++n)
            {
                splitted = input
                    .ReadLine()
                    .Split((char[])null, StringSplitOptions.RemoveEmptyEntries);
                if (int.Parse(splitted[0]) != n)
                    throw new Exception("Unexpected index");
                List<long> cluster = new List<long>();
                int i = 1;
                var customer = int.Parse(splitted[1]);
                while (customer != -1)
                {
                    // -2 because original customer indices are in 2..nbNodes
                    cluster.Add(customer - 2);
                    i++;
                    customer = int.Parse(splitted[i]);
                }
                clustersData[n - 1] = cluster.ToArray();
            }    
            
            splitted = input.ReadLine().Split(' ');
            if (!splitted[1].Contains("DEMAND_SECTION"))
                throw new Exception("Expected keyword DEMAND_SECTION");

            demandsData = new long[nbCustomers];
            for (int n = 1; n <= nbClusters; ++n)
            {
                splitted = input
                    .ReadLine()
                    .Split((char[])null, StringSplitOptions.RemoveEmptyEntries);
                if (int.Parse(splitted[0]) != n)
                    throw new Exception("Unexpected index");
                var demand = int.Parse(splitted[1]);
                demandsData[n - 1] = demand;
            }

        }
    }

    // Compute the distance matrix
    private void ComputeDistanceMatrix(int depotX, int depotY, int[] customersX, int[] customersY)
    {
        distMatrixData = new long[nbCustomers][];
        for (int i = 0; i < nbCustomers; ++i)
            distMatrixData[i] = new long[nbCustomers];

        for (int i = 0; i < nbCustomers; ++i)
        {
            distMatrixData[i][i] = 0;
            for (int j = i + 1; j < nbCustomers; ++j)
            {
                long dist = ComputeDist(customersX[i], customersX[j], customersY[i], 
                    customersY[j]);
                distMatrixData[i][j] = dist;
                distMatrixData[j][i] = dist;
            }
        }

        distDepotData = new long[nbCustomers];
        for (int i = 0; i < nbCustomers; ++i)
            distDepotData[i] = ComputeDist(depotX, customersX[i], depotY, customersY[i]);
    }

    private long ComputeDist(int xi, int xj, int yi, int yj)
    {
        double exactDist = Math.Sqrt(Math.Pow(xi - xj, 2) + Math.Pow(yi - yj, 2));
        return Convert.ToInt64(Math.Round(exactDist));
    }
}

import java.util.*;
import java.io.*;
import com.hexaly.optimizer.*;

public class clustered_vehicle_routing {
    // Hexaly Optimizer
    private final HexalyOptimizer optimizer;

    // Number of customers
    int nbCustomers;

    // Capacity of the trucks
    private int truckCapacity;

    // Demand on each customer
    private long[] demandsData;

    // Customers in each cluster 
    private long[][] clustersData;

    // Distance matrix between customers
    private long[][] distMatrixData;

    // Distances between customers and depot
    private long[] distDepotData;

    // Number of trucks
    private int nbTrucks;

    // Number of clusters
    private int nbClusters;

    // Decision variables
    private HxExpression[] truckSequences;
    private HxExpression[] clustersSequences;

    // Distance traveled by all the trucks
    private HxExpression totalDistance;

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

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

        // Create HexalyOptimizer arrays to be able to access them with an "at" operator
        HxExpression demands = model.array(demandsData);
        HxExpression distDepot = model.array(distDepotData);
        HxExpression distMatrix = model.array(distMatrixData);
        HxExpression[] distClusters = new HxExpression[nbClusters];
        HxExpression[] distRoutes = new HxExpression[nbTrucks];
        HxExpression[] initialNodes = new HxExpression[nbClusters];
        HxExpression[] endNodes = new HxExpression[nbClusters];

        clustersSequences  = new HxExpression[nbClusters];
        // A list is created for each cluster, to determine the order within the cluster
        for (int k = 0; k < nbClusters; ++k) {
            clustersSequences[k] = model.listVar(clustersData[k].length);
            // All customers in the cluster must be visited 
            model.constraint(model.eq(model.count(clustersSequences[k]), 
                clustersData[k].length));
        }

        for (int k = 0; k < nbClusters; ++k) {
            HxExpression sequence = clustersSequences[k];
            HxExpression sequenceCluster = model.array(clustersData[k]);
            HxExpression c = model.count(sequence);

            // Distance traveled within cluster k
            HxExpression distLambda = model
                .lambdaFunction(i -> model.at(distMatrix, 
                    model.at(sequenceCluster, model.at(sequence, model.sub(i, 1))), 
                    model.at(sequenceCluster, model.at(sequence, i) ))
                );
            distClusters[k] = model.sum(model.range(1,c), distLambda);

            // First and last point when visiting cluster k
            initialNodes[k] = model.at(sequenceCluster, model.at(sequence, 0));
            endNodes[k]     = model.at(sequenceCluster, model.at(sequence, model.sub(c,1)));
        }

        truckSequences = new HxExpression[nbTrucks];
        // Sequence of clusters visited by each truck.
        for (int k = 0; k < nbTrucks; ++k)
            truckSequences[k] = model.listVar(nbClusters);
        // Each cluster must be visited by exactly one truck
        model.constraint(model.partition(truckSequences));

        HxExpression valueDistCluster = model.array(distClusters);
        HxExpression initials         = model.array(initialNodes);
        HxExpression ends             = model.array(endNodes);
        
        for (int k = 0; k < nbTrucks; ++k) {
            HxExpression sequence = truckSequences[k];
            HxExpression c = model.count(sequence);

            // The quantity needed in each route must not exceed the truck capacity
            HxExpression demandLambda = model.lambdaFunction(j -> model.at(demands, j));
            HxExpression routeQuantity = model.sum(sequence, demandLambda);
            model.constraint(model.leq(routeQuantity, truckCapacity));

            // Distance traveled by truck k
            // = distance in each cluster + distance between clusters + distance with depot 
            // at the beginning end at the end of a route
            HxExpression distLambda = model.lambdaFunction(i -> model.sum(
                    model.at(valueDistCluster, model.at(sequence,i)),  
                    model.at(distMatrix, model.at(ends, model.at(sequence, model.sub(i, 1))), 
                    model.at(initials, model.at(sequence, i)))) );
            distRoutes[k] = model.sum(
                    model.sum(model.range(1, c), distLambda),
                    model.iif(model.gt(c, 0), model.sum( 
                    model.at(valueDistCluster, model.at(sequence,0)),
                    model.at(distDepot, model.at(initials, model.at(sequence, 0))),
                    model.at(distDepot, model.at(ends, model.at(sequence, model.sub(c, 1))))),
                    0) );
        }

        totalDistance = model.sum(distRoutes);

        // Objective: minimize the distance traveled
        model.minimize(totalDistance);

        model.close();

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

        optimizer.solve();
    }

    /*
     * Write the solution in a file with the following format:
     * - number of trucks used and total distance
     * - for each truck the customers visited (omitting the start/end at the depot)
     */
    private void writeSolution(String fileName) throws IOException {
        try (PrintWriter output = new PrintWriter(fileName)) {
            output.println( totalDistance.getValue());
            for (int k = 0; k < nbTrucks; ++k) {

                // Values in sequence are in [0..nbCustomers-1]. +2 is to put it back in
                // [2..nbCustomers+1] as in the data files (1 being the depot)
                HxCollection customersCollection = truckSequences[k].getCollectionValue();
                for (int i = 0; i < customersCollection.count(); ++i) {
                    int cluster = (int) customersCollection.get(i);
                    HxCollection clustersCollection = 
                            clustersSequences[cluster].getCollectionValue();
                    for (int j = 0; j < clustersCollection.count(); ++j)
                        output.print((clustersData[cluster][(int) clustersCollection.get(j)] + 2)
                            + " ");
                }
                output.println();
            }
        }
    }

    // The input files follow the "Augerat" format
    private void readInstance(String fileName) throws IOException {
        try (Scanner input = new Scanner(new File(fileName))) {
            int nbNodes = 0;
            String[] splitted;
            while (true) {
                splitted = input.nextLine().split(":");
                if (splitted[0].contains("DIMENSION")) {
                    nbNodes = Integer.parseInt(splitted[1].trim());
                    nbCustomers = nbNodes - 1;
                } else if (splitted[0].contains("VEHICLES")) {
                    nbTrucks = Integer.parseInt(splitted[1].trim());
                } else if (splitted[0].contains("GVRP_SETS")) {
                    nbClusters = Integer.parseInt(splitted[1].trim());
                } else if (splitted[0].contains("CAPACITY")) {
                    truckCapacity = Integer.parseInt(splitted[1].trim());
                } else if (splitted[0].contains("NODE_COORD_SECTION")) {
                    break;
                }
            }
            int[] customersX = new int[nbCustomers];
            int[] customersY = new int[nbCustomers];
            int depotX = 0, depotY = 0;
            for (int n = 1; n <= nbNodes; ++n) {
                int id = input.nextInt();
                if (id != n)
                    throw new IOException("Unexpected index");
                if (n == 1) {
                    depotX = (int) input.nextFloat();
                    depotY = (int) input.nextFloat();
                } else {
                    // -2 because original customer indices are in 2..nbNodes
                    customersX[n - 2] = (int) input.nextFloat();
                    customersY[n - 2] = (int) input.nextFloat();
                }
            }
            computeDistanceMatrix(depotX, depotY, customersX, customersY);

            input.nextLine().split(""); // End the last line
            String name = input.nextLine();
            if (!name.contains("GVRP_SET_SECTION")) {
                throw new RuntimeException("Expected keyword GVRP_SET_SECTION");
            }
            clustersData = new long[nbClusters][];
            for (int n = 1; n <= nbClusters; ++n) {
                List<Long> cluster = new ArrayList<>();
                int id = input.nextInt();
                if (id != n)
                    throw new IOException("Unexpected index");
                long customer = input.nextInt();
                while (customer != -1) {
                    // -2 because original customer indices are in 2..nbNodes
                    cluster.add(customer - 2);
                    customer = input.nextInt();
                }
                long[] clusterArray = cluster.stream()
                        .mapToLong(Long::intValue).toArray();
                clustersData[n - 1] = clusterArray;
            }

            input.nextLine().split(""); 
            name = input.nextLine();
            if (!name.contains("DEMAND_SECTION")) {
                throw new RuntimeException("Expected keyword DEMAND_SECTION");
            }
            demandsData = new long[nbCustomers];
            for (int n = 1; n <= nbClusters; ++n) {
                int id = input.nextInt();
                if (id != n) throw new IOException("Unexpected index");
                int demand = input.nextInt();
                demandsData[n - 1] = demand;
            }
        }
    }

    // Compute the distance matrix
    private void computeDistanceMatrix(int depotX, int depotY, int[] customersX, int[] customersY) {
        distMatrixData = new long[nbCustomers][nbCustomers];
        for (int i = 0; i < nbCustomers; ++i) {
            distMatrixData[i][i] = 0;
            for (int j = i + 1; j < nbCustomers; ++j) {
                long dist = computeDist(customersX[i], customersX[j], customersY[i], customersY[j]);
                distMatrixData[i][j] = dist;
                distMatrixData[j][i] = dist;
            }
        }

        distDepotData = new long[nbCustomers];
        for (int i = 0; i < nbCustomers; ++i) {
            distDepotData[i] = computeDist(depotX, customersX[i], depotY, customersY[i]);
        }
    }

    private long computeDist(int xi, int xj, int yi, int yj) {
        double exactDist = Math.sqrt(Math.pow(xi - xj, 2) + Math.pow(yi - yj, 2));
        return Math.round(exactDist);
    }

    public static void main(String[] args) {
        if (args.length < 1) {
            System.err.println("Usage: java clustered-vehicle-routing inputFile [outputFile] [timeLimit]");
            System.exit(1);
        }
        try (HexalyOptimizer optimizer = new HexalyOptimizer()) {
            String instanceFile = args[0];
            String outputFile = args.length > 1 ? args[1] : null;
            String strTimeLimit = args.length > 2 ? args[2] : "20";

            clustered_vehicle_routing model = new clustered_vehicle_routing(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);
        }
    }
}