Capacitated Arc Routing Problem (CARP)

import sys
import hexaly.optimizer


def main(instance_file, output_file, str_time_limit):
    #
    # Read instance data
    #
    instance = CarpInstance(instance_file)

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

        # Sequence of edges visited and "serviced" by each truck
        edges_sequences_vars = [model.list(2 * instance.nb_required_edges)
                                for _ in range(instance.nb_trucks)]
        edges_sequences = model.array(edges_sequences_vars)

        # Create distance and cost arrays to be able to access it with an "at" operator
        demands = model.array(instance.demands_data)
        costs = model.array(instance.costs_data)
        dist_from_depot = model.array(instance.dist_from_depot_data)
        dist_to_depot = model.array(instance.dist_to_depot_data)
        edges_dist = model.array()
        for n in range(2 * instance.nb_required_edges):
            edges_dist.add_operand(model.array(instance.edges_dist_data[n]))

        # An edge must be serviced by at most one truck
        model.constraint(model.disjoint(edges_sequences))

        # An edge can be travelled in both directions but its demand must be 
        # satisfied only once
        for i in range(instance.nb_required_edges):
            model.constraint(
                model.contains(edges_sequences, 2 * i)
                + model.contains(edges_sequences, 2 * i + 1) 
                == 1)

        route_distances = [None] * instance.nb_trucks
        for k in range(instance.nb_trucks):
            sequence = edges_sequences_vars[k]
            c = model.count(sequence)

            # Quantity in each truck
            demand_lambda = model.lambda_function(lambda j: demands[j])
            route_quantity = model.sum(sequence, demand_lambda)
            # Capacity constraint : a truck must not exceed its capacity
            model.constraint(route_quantity <= instance.truck_capacity)

            # Distance travelled by each truck
            dist_lambda = model.lambda_function(
                lambda i:
                    costs[sequence[i]]
                    + model.at(
                        edges_dist,
                        sequence[i - 1],
                        sequence[i]))
            route_distances[k] = model.sum(model.range(1, c), dist_lambda) \
                + model.iif(
                    c > 0,
                    costs[sequence[0]] + dist_from_depot[sequence[0]] \
                        + dist_to_depot[sequence[c - 1]],
                    0)

        # Total distance travelled
        total_distance = model.sum(route_distances)

        # Objective: minimize the distance travelled
        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
        #  - number of routes
        #  - for each truck, the edges visited
        #
        if output_file:
            with open(output_file, 'w') as f:
                f.write("Objective function value : {}\nNumber of routes : {}\n".format(
                    total_distance.value, instance.nb_trucks))
                for k in range(instance.nb_trucks):
                    f.write("Sequence of truck {}: ".format(k + 1))
                    sequence = edges_sequences_vars[k].value
                    c = len(sequence)
                    for i in range(c):
                        f.write("({}, {})  ".format(instance.origins_data[sequence[i]],
                                instance.destinations_data[sequence[i]]))
                    f.write("\n")


class CarpInstance:

    def read_elem(self, filename):
        with open(filename) as f:
            return [str(elem) for elem in f.read().strip().split("\n")]

    # The input files follow the format of the DIMACS challenge
    def __init__(self, filename):
        file_it = iter(self.read_elem(filename))
        for _ in range(2):
            next(file_it)
        nb_nodes = int(next(file_it).strip().split(":")[1])
        self.nb_required_edges = int(next(file_it).strip().split(":")[1])
        nb_not_required_edges = int(next(file_it).strip().split(":")[1])
        self.nb_trucks = int(next(file_it).strip().split(":")[1])
        self.truck_capacity = int(next(file_it).strip().split(":")[1])
        for _ in range(3):
            next(file_it)
        self.demands_data = list()
        self.costs_data = list()
        self.origins_data = list()
        self.destinations_data = list()
        required_nodes = list()
        node_neighbors = [([0] * nb_nodes) for _ in range(nb_nodes)]
        for _ in range(self.nb_required_edges):
            elements = next(file_it)
            edge = tuple(map(int, elements.strip().split("   ")[0][2:-1].strip().split(",")))
            cost = int(elements.strip().split("   ")[1].strip().split()[1])
            demand = int(elements.strip().split("   ")[2].strip().split()[1])
            for _ in range(2):
                self.costs_data.append(cost)
                self.demands_data.append(demand)
            # even indices store direct edges, and odd indices store reverse edges
            self.origins_data.append(edge[0])
            self.destinations_data.append(edge[1])
            self.origins_data.append(edge[1])
            self.destinations_data.append(edge[0])
            if edge[0] not in required_nodes:
                required_nodes.append(edge[0])
            if edge[1] not in required_nodes:
                required_nodes.append(edge[1])
            node_neighbors[edge[0] - 1][edge[1] - 1] = cost
            node_neighbors[edge[1] - 1][edge[0] - 1] = cost
        if nb_not_required_edges > 0:
            next(file_it)
            for _ in range(nb_not_required_edges):
                elements = next(file_it)
                edge = tuple(map(int, elements.strip().split("   ")[0][2:-1].strip().split(",")))
                cost = int(elements.strip().split("   ")[1].strip().split()[1])
                node_neighbors[edge[0] - 1][edge[1] - 1] = cost
                node_neighbors[edge[1] - 1][edge[0] - 1] = cost
        depot_node = int(next(file_it).strip().split(":")[1])
        # Finds the shortest path from one "required node" to another
        nb_required_nodes = len(required_nodes)
        required_distances = list()
        for node in required_nodes:
            paths = self.shortest_path_finder(node, nb_nodes, node_neighbors)
            required_distances.append(paths)
        # Since we can explore the edges in both directions, we will represent all possible
        # edges with an index
        self.edges_dist_data = None
        self.find_distance_between_edges(nb_required_nodes, required_nodes, required_distances)
        self.dist_to_depot_data = None
        self.find_distance_to_depot(nb_required_nodes, depot_node,
                                    required_nodes, required_distances)
        self.dist_from_depot_data = None
        self.find_distance_from_depot(nb_required_nodes, nb_nodes, depot_node,
                                      required_nodes, required_distances, node_neighbors)

    # Finds the shortest path from one node "origin" to all the other nodes of the graph
    # thanks to the Dijkstra's algorithm
    def min_distance(self, nb_nodes, shortest_path, sptSet):
        min = sys.maxsize
        for i in range(nb_nodes):
            if shortest_path[i] < min and sptSet[i] == False:
                min = shortest_path[i]
                min_index = i
        return min_index

    def shortest_path_finder(self, origin, nb_nodes, node_neighbors):
        shortest_path = [sys.maxsize] * nb_nodes
        shortest_path[origin - 1] = 0
        sptSet = [False] * nb_nodes
        for _ in range(nb_nodes):
            current_node = self.min_distance(nb_nodes, shortest_path, sptSet)
            sptSet[current_node] = True
            current_neighbors = node_neighbors[current_node]
            for neighbor in range(nb_nodes):
                if current_neighbors[neighbor] != 0:
                    distance = current_neighbors[neighbor]
                    if ((sptSet[neighbor] == False) and
                            (shortest_path[current_node] + distance < shortest_path[neighbor])):
                        shortest_path[neighbor] = distance + shortest_path[current_node]
        return shortest_path

    def find_distance_between_edges(self, nb_required_nodes, required_nodes, required_distances):
        self.edges_dist_data = [[None] * (2 * self.nb_required_edges)
                                for _ in range(2 * self.nb_required_edges)]
        for i in range(2 * self.nb_required_edges):
            for j in range(2 * self.nb_required_edges):
                if self.destinations_data[i] == self.origins_data[j]:
                    self.edges_dist_data[i][j] = 0
                else:
                    for k in range(nb_required_nodes):
                        if required_nodes[k] == self.destinations_data[i]:
                            self.edges_dist_data[i][j] = required_distances[k][
                                self.origins_data[j] - 1]

    def find_distance_to_depot(
            self, nb_required_nodes, depot_node, required_nodes, required_distances):
        self.dist_to_depot_data = [None] * (2 * self.nb_required_edges)
        for i in range(2 * self.nb_required_edges):
            if self.destinations_data[i] == depot_node:
                self.dist_to_depot_data[i] = 0
            else:
                for k in range(nb_required_nodes):
                    if required_nodes[k] == self.destinations_data[i]:
                        self.dist_to_depot_data[i] = required_distances[k][depot_node-1]

    def find_distance_from_depot(self, nb_required_nodes, nb_nodes, depot_node,
                                 required_nodes, required_distances, node_neighbors):
        self.dist_from_depot_data = [None] * (2 * self.nb_required_edges)
        for i in range(2 * self.nb_required_edges):
            if depot_node == self.origins_data[i]:
                self.dist_from_depot_data[i] = 0
            else:
                depot_is_required_node = False
                for k in range(nb_required_nodes):
                    if required_nodes[k] == depot_node:
                        depot_is_required_node = True
                        self.dist_from_depot_data[i] = required_distances[k][
                            self.origins_data[i] - 1]
                if not depot_is_required_node:
                    shortest_paths_from_depot = self.shortest_path_finder(
                        depot_node, nb_nodes, node_neighbors)
                    self.dist_from_depot_data[i] = shortest_paths_from_depot[self.origins_data[i] - 1]


if __name__ == '__main__':
    if len(sys.argv) < 2:
        print("Usage: python capacitated_arc_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, output_file, str_time_limit)

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

using namespace hexaly;
using namespace std;

class CapacitatedArcRouting {

public:
    // Data from the problem
    int nbTrucks;
    int truckCapacity;
    int nbNodes;
    int nbRequiredEdges;
    int nbRequiredNodes;
    int nbNotRequiredNodes;
    int depotNode;
    vector<int> demandsData;
    vector<int> costsData;
    vector<int> originsData;
    vector<int> destinationsData;
    vector<int> requiredNodes;
    vector<vector<int>> nodesDistData;
    vector<vector<int>> edgesDistData;
    vector<int> distToDepotData;
    vector<int> distFromDepotData;
    vector<vector<int>> nodeNeighborsData;

    // Hexaly Optimizer
    HexalyOptimizer optimizer;

    // Decision variables
    vector<HxExpression> edgesSequencesVars;
    HxExpression totalDistance;

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

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

        // Sequence of edges visited and "serviced" by each truck
        edgesSequencesVars.resize(nbTrucks);
        for (int i = 0; i < nbTrucks; ++i) {
            // A sequence can contain an edge in one direction or its reverse
            edgesSequencesVars[i] = model.listVar(2 * nbRequiredEdges);
        }
        HxExpression edgesSequences = model.array(edgesSequencesVars.begin(), edgesSequencesVars.end());

        // Create distance and cost arrays to be able to access it with an "at" operator
        HxExpression demands = model.array(demandsData.begin(), demandsData.end());
        HxExpression costs = model.array(costsData.begin(), costsData.end());
        HxExpression distFromDepot = model.array(distFromDepotData.begin(), distFromDepotData.end());
        HxExpression distToDepot = model.array(distToDepotData.begin(), distToDepotData.end());
        HxExpression edgesDist = model.array();
        for (int n = 0; n < 2 * nbRequiredEdges; ++n) {
            edgesDist.addOperand(model.array(edgesDistData[n].begin(), edgesDistData[n].end()));
        }

        // An edge must be serviced by at most one truck
        model.constraint(model.disjoint(edgesSequences));

        // An edge can be travelled in both directions but its demand must be satisfied only once
        for (int i = 0; i < nbRequiredEdges; ++i) {
            model.constraint(model.contains(edgesSequences, 2 * i) + model.contains(edgesSequences, 2 * i + 1) == 1);
        }

        vector<HxExpression> routeDistances(nbTrucks);
        for (int k = 0; k < nbTrucks; ++k) {
            HxExpression sequence = edgesSequencesVars[k];
            HxExpression c = model.count(sequence);

            // Quantity in each truck
            HxExpression demandLambda =
                model.createLambdaFunction([&](HxExpression j) { return demands[j]; });
            HxExpression routeQuantity = model.sum(sequence, demandLambda);
            // Capacity constraint : a truck must not exceed its capacity
            model.constraint(routeQuantity <= truckCapacity);

            // Distance travelled by each truck
            HxExpression distLambda = model.createLambdaFunction(
                [&](HxExpression i) { return costs[sequence[i]] + model.at(edgesDist, sequence[i - 1], sequence[i]); });
            routeDistances[k] =
                model.sum(model.range(1, c), distLambda) +
                model.iif(c > 0, costs[sequence[0]] + distFromDepot[sequence[0]] + distToDepot[sequence[c - 1]], 0);
        }

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

        // Objective : minimize the distance travelled
        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
     * - number of routes
     * - for each truck, the edges visited */
    void writeSolution(const string& fileName) {
        ofstream outfile;
        outfile.exceptions(ofstream::failbit | ofstream::badbit);
        outfile.open(fileName.c_str());

        outfile << "Objective function value : " << totalDistance.getValue() << endl;
        outfile << "Number of routes : " << nbTrucks << endl;
        for (int k = 0; k < nbTrucks; ++k) {
            outfile << "Sequence of truck " << k + 1 << ": ";
            HxCollection sequence = edgesSequencesVars[k].getCollectionValue();
            int c = sequence.count();
            for (int i = 0; i < c; ++i) {
                outfile << "(" << originsData[sequence[i]] << ", " << destinationsData[sequence[i]] << ")  ";
            }
            outfile << endl;
        }
    }

private:
    // The input files follow the format of the DIMACS challenge
    void readCarp(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 first = -1;
        int second = -1;
        for (int i = 0; i < 2; ++i)
            getline(infile, str);
        getline(infile, str);
        line = strdup(str.c_str());
        pch = strtok(line, " : ");
        pch = strtok(NULL, " : ");
        nbNodes = atoi(pch);
        getline(infile, str);
        line = strdup(str.c_str());
        pch = strtok(line, " : ");
        pch = strtok(NULL, " : ");
        nbRequiredEdges = atoi(pch);
        getline(infile, str);
        line = strdup(str.c_str());
        pch = strtok(line, " : ");
        pch = strtok(NULL, " : ");
        nbNotRequiredNodes = atoi(pch);
        getline(infile, str);
        line = strdup(str.c_str());
        pch = strtok(line, " : ");
        pch = strtok(NULL, " : ");
        nbTrucks = atoi(pch);
        getline(infile, str);
        line = strdup(str.c_str());
        pch = strtok(line, " : ");
        pch = strtok(NULL, " : ");
        truckCapacity = atoi(pch);
        for (int i = 0; i < 3; ++i)
            getline(infile, str);
        nodeNeighborsData.resize(nbNodes);
        for (int i = 0; i < nbNodes; ++i) {
            nodeNeighborsData[i].resize(nbNodes);
        }
        for (int i = 0; i < nbNodes; ++i) {
            for (int j = 0; j < nbNodes; ++j) {
                nodeNeighborsData[i][j] = 0;
            }
        }
        getline(infile, str);
        for (int i = 0; i < nbRequiredEdges; ++i) {
            extractEdge(str, pch, line, first, second, true);
            extractCost(str, pch, line, first, second, true);
            extractDemand(str, pch, line);
            getline(infile, str);
        }
        if (nbNotRequiredNodes > 0) {
            getline(infile, str);
            for (int i = 0; i < nbNotRequiredNodes; ++i) {
                extractEdge(str, pch, line, first, second, false);
                extractCost(str, pch, line, first, second, false);
                getline(infile, str);
            }
        }
        line = strdup(str.c_str());
        pch = strtok(line, " : ");
        pch = strtok(NULL, " : ");
        depotNode = atoi(pch);
        nbRequiredNodes = requiredNodes.size();
        vector<int> shortestPath;
        shortestPath.resize(nbNodes);
        for (int i = 0; i < nbRequiredNodes; ++i) {
            shortestPathFinder(requiredNodes[i], shortestPath);
            for (int j = 0; j < nbNodes; ++j) {
                if (shortestPath[j] == std::numeric_limits<int>::max()) {
                }
            }
            nodesDistData.push_back(shortestPath);
        }
        findDistanceBetweenEdges();
        findDistanceToDepot();
        findDistanceFromDepot(shortestPath);
    }

    // Finds the shortest path from one node "origin" to all the other nodes of the graph
    // thanks to the Dijkstra's algorithm.
    int minDistance(vector<int> shortestPath, vector<bool> sptSet) {
        // Initializes min value
        int min = std::numeric_limits<int>::max(), minIndex = -1;
        for (int i = 0; i < nbNodes; ++i) {
            if (sptSet[i] == false && shortestPath[i] <= min) {
                min = shortestPath[i];
                minIndex = i;
            }
        }
        return minIndex;
    }

    void shortestPathFinder(int origin, vector<int>& shortestPath) {
        vector<int> currentNeighbors;
        int currentNode;
        int distance;
        vector<bool> sptSet;
        sptSet.resize(nbNodes);
        for (int i = 0; i < nbNodes; ++i) {
            shortestPath[i] = std::numeric_limits<int>::max();
            sptSet[i] = false;
        }
        shortestPath[origin - 1] = 0;
        for (int count = 0; count < nbNodes; ++count) {
            currentNode = minDistance(shortestPath, sptSet);
            sptSet[currentNode] = true;
            currentNeighbors = nodeNeighborsData[currentNode];
            for (int neighbor = 0; neighbor < nbNodes; ++neighbor) {
                if (currentNeighbors[neighbor] != 0) {
                    distance = currentNeighbors[neighbor];
                    if ((sptSet[neighbor] == false) &&
                        (shortestPath[currentNode] + distance < shortestPath[neighbor])) {
                        shortestPath[neighbor] = shortestPath[currentNode] + distance;
                    }
                }
            }
        }
    }

    void findDistanceBetweenEdges() {
        edgesDistData.resize(2 * nbRequiredEdges);
        for (int i = 0; i < 2 * nbRequiredEdges; ++i) {
            edgesDistData[i].resize(2 * nbRequiredEdges);
        }
        for (int i = 0; i < 2 * nbRequiredEdges; ++i) {
            for (int j = 0; j < 2 * nbRequiredEdges; ++j) {
                if (destinationsData[i] == originsData[j]) {
                    edgesDistData[i][j] = 0;
                } else {
                    for (int k = 0; k < nbRequiredNodes; ++k) {
                        if (requiredNodes[k] == destinationsData[i]) {
                            edgesDistData[i][j] = nodesDistData[k][originsData[j] - 1];
                        }
                    }
                }
            }
        }
    }

    void findDistanceToDepot() {
        distToDepotData.resize(2 * nbRequiredEdges);
        for (int i = 0; i < 2 * nbRequiredEdges; ++i) {
            if (destinationsData[i] == depotNode) {
                distToDepotData[i] = 0;
            } else {
                for (int k = 0; k < nbRequiredNodes; ++k) {
                    if (requiredNodes[k] == destinationsData[i]) {
                        distToDepotData[i] = nodesDistData[k][depotNode - 1];
                    }
                }
            }
        }
    }

    void findDistanceFromDepot(vector<int>& shortestPath) {
        distFromDepotData.resize(2 * nbRequiredEdges);
        for (int i = 0; i < 2 * nbRequiredEdges; ++i) {
            if (depotNode == originsData[i]) {
                distFromDepotData[i] = 0;
            } else {
                bool depotIsRequiredNode = false;
                for (int k = 0; k < nbRequiredNodes; ++k) {
                    if (requiredNodes[k] == depotNode) {
                        depotIsRequiredNode = true;
                        distFromDepotData[i] = nodesDistData[k][originsData[i] - 1];
                    }
                }
                if (!depotIsRequiredNode) {
                    shortestPathFinder(depotNode, shortestPath);
                    distFromDepotData[i] = shortestPath[originsData[i] - 1];
                }
            }
        }
    }

    void extractEdge(const string& str, char* pch, char* line, int& first, int& second, bool req) {
        line = strdup(str.c_str());
        pch = strtok(line, " ");
        if (strcmp(pch, "(") == 0) {
            pch = strtok(NULL, ",");
            first = atoi(pch);
            pch = strtok(NULL, ")");
            second = atoi(pch);
        }

        if (req == true) {
            // Even indices store direct edges, and odd indices store reverse edges
            originsData.push_back(first);
            destinationsData.push_back(second);
            originsData.push_back(second);
            destinationsData.push_back(first);
            if (find(requiredNodes.begin(), requiredNodes.end(), first) == requiredNodes.end()) {
                requiredNodes.push_back(first);
            }
            if (find(requiredNodes.begin(), requiredNodes.end(), second) == requiredNodes.end()) {
                requiredNodes.push_back(second);
            }
        }
    }

    void extractCost(const string& str, char* pch, char* line, int& first, int& second, bool req) {
        line = strdup(str.c_str());
        pch = strtok(line, "coste");
        int cost;
        pch = strtok(NULL, "coste");
        cost = atoi(pch);
        if (req == true) {
            for (int i = 0; i < 2; ++i)
                costsData.push_back(cost);
        }
        nodeNeighborsData[first - 1][second - 1] = cost;
        nodeNeighborsData[second - 1][first - 1] = cost;
    }

    void extractDemand(const string& str, char* pch, char* line) {
        line = strdup(str.c_str());
        pch = strtok(line, "demanda");
        int demand;
        pch = strtok(NULL, "demanda");
        demand = atoi(pch);
        for (int i = 0; i < 2; ++i)
            demandsData.push_back(demand);
    }
};

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

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

    try {
        CapacitatedArcRouting model;
        model.readInstance(instanceFile);
        model.solve(atoi(strTimeLimit));
        if (outputFile != NULL)
            model.writeSolution(outputFile);
        return 0;
    } catch (const exception& e) {
        cerr << "An error occurred: " << e.what() << endl;
        return 1;
    }
}

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

public class CapacitatedArcRouting : IDisposable
{
    // Data from the problem
    int nbTrucks;
    int truckCapacity;
    int nbNodes;
    int nbRequiredEdges;
    int nbRequiredNodes;
    int nbNotRequiredEdges;
    int depotNode;
    int[] demandsData;
    int[] costsData;
    List<int> requiredNodes;
    int[] originsData;
    int[] destinationsData;
    int[] distToDepotData;
    int[] distFromDepotData;
    int[][] nodesDistData;
    int[][] edgesDistData;
    int[][] nodeNeighborsData;

    // Hexaly Optimizer
    HexalyOptimizer optimizer;

    // Decision variables
    HxExpression[] edgesSequencesVars;
    HxExpression totalCost;

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

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

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

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

        // Sequence of edges visited and "serviced" by each truck
        edgesSequencesVars = new HxExpression[nbTrucks];
        for (int i = 0; i < nbTrucks; ++i)
        {
            // A sequence can contain an edge in one direction or its reverse
            edgesSequencesVars[i] = model.List(2 * nbRequiredEdges);
        }
        HxExpression edgesSequences = model.Array(edgesSequencesVars);

        // Creates distance and cost arrays to be able to access it with an "at" operator
        HxExpression demands = model.Array(demandsData);
        HxExpression costs = model.Array(costsData);
        HxExpression distFromDepot = model.Array(distFromDepotData);
        HxExpression distToDepot = model.Array(distToDepotData);
        HxExpression edgesDist = model.Array(edgesDistData);

        // An edge must be serviced by at most one truck
        model.Constraint(model.Disjoint(edgesSequences));

        // An edge can be travelled in both directions but its demand must be satisfied only once
        for (int i = 0; i < nbRequiredEdges; ++i)
        {
            model.Constraint(
                model.Contains(edgesSequences, 2 * i) + model.Contains(edgesSequences, 2 * i + 1)
                    == 1
            );
        }

        HxExpression[] routeDistances = new HxExpression[nbTrucks];
        for (int k = 0; k < nbTrucks; ++k)
        {
            HxExpression sequence = edgesSequencesVars[k];
            HxExpression c = model.Count(sequence);

            // Quantity in each truck
            HxExpression demandLambda = model.LambdaFunction(j => demands[j]);
            HxExpression routeQuantity = model.Sum(sequence, demandLambda);
            // Capacity constraint : a truck must not exceed its capacity
            model.Constraint(routeQuantity <= truckCapacity);

            // Distance travelled by each truck
            HxExpression distLambda = model.LambdaFunction(
                i => costs[sequence[i]] + edgesDist[sequence[i - 1], sequence[i]]
            );
            routeDistances[k] =
                model.Sum(model.Range(1, c), distLambda)
                + model.If(
                    c > 0,
                    costs[sequence[0]] + distFromDepot[sequence[0]] + distToDepot[sequence[c - 1]],
                    0
                );
        }

        // Total Distance travelled
        totalCost = model.Sum(routeDistances);

        // Objective : minimize the distance travelled
        model.Minimize(totalCost);

        model.Close();

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

        optimizer.Solve();
    }

    /* Write the solution in a file with the following format:
    - total distance
    - number of routes
    - for each truck, the edges visited */
    void WriteSolution(string fileName)
    {
        using (StreamWriter output = new StreamWriter(fileName))
        {
            output.WriteLine("Objective function value : " + totalCost.GetValue());
            output.WriteLine("Number of routes : " + nbTrucks);
            for (int k = 0; k < nbTrucks; ++k)
            {
                output.Write("Sequence of truck " + (k + 1) + ": ");
                HxCollection sequence = edgesSequencesVars[k].GetCollectionValue();
                int c = sequence.Count();
                for (int i = 0; i < c; ++i)
                {
                    output.Write(
                        "("
                            + originsData[sequence[i]]
                            + ", "
                            + destinationsData[sequence[i]]
                            + ")  "
                    );
                }
                output.WriteLine();
            }
        }
    }

    public static void Main(string[] args)
    {
        if (args.Length < 1)
        {
            Console.WriteLine("Usage: CapacitatedArcRouting 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] : "20";

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

    // The input files follow the format of the DIMACS challenge
    private void ReadInputCarp(string fileName)
    {
        using (StreamReader input = new StreamReader(fileName))
        {
            string[] splitted;
            for (int i = 0; i < 2; ++i)
                input.ReadLine();
            splitted = input.ReadLine().Split(':');
            nbNodes = int.Parse(splitted[1]);
            splitted = input.ReadLine().Split(':');
            nbRequiredEdges = int.Parse(splitted[1]);
            splitted = input.ReadLine().Split(':');
            nbNotRequiredEdges = int.Parse(splitted[1]);
            splitted = input.ReadLine().Split(':');
            nbTrucks = int.Parse(splitted[1]);
            splitted = input.ReadLine().Split(':');
            truckCapacity = int.Parse(splitted[1]);
            for (int i = 0; i < 2; ++i)
                input.ReadLine();
            // Even indices store direct edges, and odd indices store reverse edges
            originsData = new int[2 * nbRequiredEdges];
            destinationsData = new int[2 * nbRequiredEdges];
            demandsData = new int[2 * nbRequiredEdges];
            costsData = new int[2 * nbRequiredEdges];
            nodeNeighborsData = new int[nbNodes][];
            for (int i = 0; i < nbNodes; ++i)
                nodeNeighborsData[i] = new int[nbNodes];
            for (int i = 0; i < nbNodes; ++i)
            {
                for (int j = 0; j < nbNodes; ++j)
                    nodeNeighborsData[i][j] = 0;
            }
            requiredNodes = new List<int>();
            input.ReadLine();
            for (int i = 0; i < nbRequiredEdges; ++i)
            {
                splitted = input.ReadLine().Split(' ');
                int firstNode = int.Parse(splitted[2].Remove(splitted[2].Length - 1));
                int secondNode = int.Parse(splitted[3].Remove(splitted[3].Length - 1));
                int cost = int.Parse(splitted[7]);
                int demand = int.Parse(splitted[11]);
                costsData[2 * i] = cost;
                demandsData[2 * i] = demand;
                costsData[2 * i + 1] = cost;
                demandsData[2 * i + 1] = demand;
                originsData[2 * i] = firstNode;
                destinationsData[2 * i] = secondNode;
                originsData[2 * i + 1] = secondNode;
                destinationsData[2 * i + 1] = firstNode;
                if (!requiredNodes.Contains(firstNode))
                    requiredNodes.Add(firstNode);
                if (!requiredNodes.Contains(secondNode))
                    requiredNodes.Add(secondNode);
                nodeNeighborsData[firstNode - 1][secondNode - 1] = cost;
                nodeNeighborsData[secondNode - 1][firstNode - 1] = cost;
            }
            if (nbNotRequiredEdges > 0)
            {
                input.ReadLine();
                for (int i = 0; i < nbNotRequiredEdges; ++i)
                {
                    splitted = input.ReadLine().Split(' ');
                    int firstNode = int.Parse(splitted[2].Remove(splitted[2].Length - 1));
                    int secondNode = int.Parse(splitted[3].Remove(splitted[3].Length - 1));
                    int cost = int.Parse(splitted[7]);
                    nodeNeighborsData[firstNode - 1][secondNode - 1] = cost;
                    nodeNeighborsData[secondNode - 1][firstNode - 1] = cost;
                }
            }
            splitted = input.ReadLine().Split(':');
            depotNode = int.Parse(splitted[1]);
        }
        nbRequiredNodes = requiredNodes.Count;
        int[] shortestPath = new int[nbNodes];
        FindRequiredPaths(ref shortestPath);
        FindDistanceBetweenEdges();
        FindDistanceToDepot();
        FindDistanceFromDepot(ref shortestPath);
    }

    // Finds the shortest path from one node "origin" to all the other nodes of the graph
    // thanks to the Dijkstra's algorithm
    public int MinDistance(int[] shortestPath, bool[] sptSet)
    {
        // Initializes min value
        int min = int.MaxValue,
            minIndex = -1;
        for (int i = 0; i < nbNodes; ++i)
            if (sptSet[i] == false && shortestPath[i] <= min)
            {
                min = shortestPath[i];
                minIndex = i;
            }
        return minIndex;
    }

    public void ShortestPathFinder(int origin, ref int[] shortestPath)
    {
        int[] currentNeighbors;
        int currentNode,
            distance;
        bool[] sptSet = new bool[nbNodes];
        for (int i = 0; i < nbNodes; ++i)
        {
            shortestPath[i] = int.MaxValue;
            sptSet[i] = false;
        }
        shortestPath[origin - 1] = 0;
        for (int count = 0; count < nbNodes; ++count)
        {
            currentNode = MinDistance(shortestPath, sptSet);
            sptSet[currentNode] = true;
            currentNeighbors = nodeNeighborsData[currentNode];
            for (int neighbor = 0; neighbor < nbNodes; ++neighbor)
            {
                if (currentNeighbors[neighbor] != 0)
                {
                    distance = currentNeighbors[neighbor];
                    if (
                        !sptSet[neighbor]
                        && shortestPath[currentNode] + distance < shortestPath[neighbor]
                    )
                    {
                        shortestPath[neighbor] = shortestPath[currentNode] + distance;
                    }
                }
            }
        }
    }

    public void FindRequiredPaths(ref int[] shortestPath)
    {
        nodesDistData = new int[nbRequiredNodes][];
        for (int i = 0; i < nbRequiredNodes; ++i)
            nodesDistData[i] = new int[nbNodes];
        for (int i = 0; i < nbRequiredNodes; ++i)
        {
            int node = requiredNodes[i];
            ShortestPathFinder(node, ref shortestPath);
            for (int j = 0; j < nbNodes; ++j)
                nodesDistData[i][j] = shortestPath[j];
        }
    }

    public void FindDistanceBetweenEdges()
    {
        edgesDistData = new int[2 * nbRequiredEdges][];
        for (int i = 0; i < 2 * nbRequiredEdges; ++i)
            edgesDistData[i] = new int[2 * nbRequiredEdges];
        for (int i = 0; i < 2 * nbRequiredEdges; ++i)
        {
            for (int j = 0; j < 2 * nbRequiredEdges; ++j)
            {
                if (destinationsData[i] == originsData[j])
                    edgesDistData[i][j] = 0;
                else
                {
                    for (int k = 0; k < nbRequiredNodes; ++k)
                    {
                        if (requiredNodes[k] == destinationsData[i])
                            edgesDistData[i][j] = nodesDistData[k][originsData[j] - 1];
                    }
                }
            }
        }
    }

    public void FindDistanceToDepot()
    {
        distToDepotData = new int[2 * nbRequiredEdges];
        for (int i = 0; i < 2 * nbRequiredEdges; ++i)
        {
            if (destinationsData[i] == depotNode)
                distToDepotData[i] = 0;
            else
            {
                for (int k = 0; k < nbRequiredNodes; ++k)
                {
                    if (requiredNodes[k] == destinationsData[i])
                        distToDepotData[i] = nodesDistData[k][depotNode - 1];
                }
            }
        }
    }

    public void FindDistanceFromDepot(ref int[] shortestPath)
    {
        distFromDepotData = new int[2 * nbRequiredEdges];
        for (int i = 0; i < 2 * nbRequiredEdges; ++i)
        {
            if (depotNode == originsData[i])
                distFromDepotData[i] = 0;
            else
            {
                bool depotIsRequiredNode = false;
                for (int k = 0; k < nbRequiredNodes; ++k)
                {
                    if (requiredNodes[k] == depotNode)
                    {
                        depotIsRequiredNode = true;
                        distFromDepotData[i] = nodesDistData[k][originsData[i] - 1];
                    }
                }
                if (!depotIsRequiredNode)
                {
                    ShortestPathFinder(depotNode, ref shortestPath);
                    distFromDepotData[i] = shortestPath[originsData[i] - 1];
                }
            }
        }
    }
}

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

public class CapacitatedArcRouting {

    // Data from the problem
    int nbTrucks;
    int truckCapacity;
    int nbNodes;
    int nbRequiredEdges;
    int nbRequiredNodes;
    int nbNotRequiredEdges;
    int depotNode;
    int[] demandsData;
    int[] costsData;
    int[] originsData;
    int[] destinationsData;
    List<Integer> requiredNodes;
    int[] distToDepotData;
    int[] distFromDepotData;
    int[][] nodesDistData;
    int[][] edgesDistData;
    int[][] nodeNeighborsData;

    // Hexaly Optimizer
    private final HexalyOptimizer optimizer;

    // Decision variables
    HxExpression[] edgesSequencesVars;
    HxExpression totalDistance;

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

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

        // Sequence of edges visited and "serviced" by each truck
        edgesSequencesVars = new HxExpression[nbTrucks];
        for (int i = 0; i < nbTrucks; ++i) {
            edgesSequencesVars[i] = model.listVar(2 * nbRequiredEdges);
        }
        HxExpression edgesSequences = model.array(edgesSequencesVars);

        // Create distance and cost arrays to be able to access it with an "at" operator
        HxExpression demands = model.array(demandsData);
        HxExpression costs = model.array(costsData);
        HxExpression distFromDepot = model.array(distFromDepotData);
        HxExpression distToDepot = model.array(distToDepotData);
        HxExpression edgesDist = model.array(edgesDistData);

        // An edge must be serviced by at most one truck
        model.constraint(model.disjoint(edgesSequences));

        // An edge can be travelled in both directions but its demand must be satisfied
        // only once
        for (int i = 0; i < nbRequiredEdges; ++i) {
            model.constraint(model
                .eq(model.sum(model.contains(edgesSequences, 2 * i), model.contains(edgesSequences, (2 * i) + 1)), 1));
        }
        HxExpression[] routeDistances = new HxExpression[nbTrucks];
        for (int k = 0; k < nbTrucks; ++k) {
            HxExpression sequence = edgesSequencesVars[k];
            HxExpression c = model.count(sequence);

            // Quantity in each truck
            HxExpression demandLambda = model.lambdaFunction(j -> model.at(demands, j));
            HxExpression routeQuantity = model.sum(sequence, demandLambda);
            // Capacity constraint : a truck must not exceed its capacity
            model.constraint(model.leq(routeQuantity, truckCapacity));

            // Distance travelled by each truck
            HxExpression distLambda = model.lambdaFunction(i -> model.sum(model.at(costs, model.at(sequence, i)),
                model.at(edgesDist, model.at(sequence, model.sub(i, 1)), model.at(sequence, i))));
            routeDistances[k] = model.sum(model.sum(model.range(1, c), distLambda),
                model.iif(model.gt(c, 0),
                    model.sum(model.at(costs, model.at(sequence, 0)), model.at(distFromDepot, model.at(sequence, 0)),
                        model.at(distToDepot, model.at(sequence, model.sub(c, 1)))),
                    0));
        }

        // Total Distance travelled
        totalDistance = model.sum(routeDistances);

        // Objective : minimize the distance travelled
        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
    // - number of routes
    // - for each truck, the edges visited
    private void WriteSolution(String filename) throws IOException {
        try (PrintWriter output = new PrintWriter(filename)) {
            output.println("Objective function value : " + totalDistance.getValue());
            output.println("Number of routes : " + nbTrucks);
            for (int k = 0; k < nbTrucks; ++k) {
                output.print("Sequence of truck " + (k + 1) + ": ");
                HxCollection sequence = edgesSequencesVars[k].getCollectionValue();
                int c = sequence.count();
                for (int i = 0; i < c; ++i) {
                    output.print("(" + originsData[(int) sequence.get(i)] + ", "
                        + destinationsData[(int) sequence.get(i)] + ")  ");
                }
                output.println();
            }
        }
    }

    // The input files follow the format of the DIMACS challenge
    private void readInstance(String filename) throws IOException {
        try (Scanner input = new Scanner(new File(filename))) {
            String[] splitted;
            for (int i = 0; i < 2; ++i)
                input.nextLine();
            splitted = input.nextLine().split(":");
            nbNodes = Integer.parseInt(splitted[1].trim());
            splitted = input.nextLine().split(":");
            nbRequiredEdges = Integer.parseInt(splitted[1].trim());
            splitted = input.nextLine().split(":");
            nbNotRequiredEdges = Integer.parseInt(splitted[1].trim());
            splitted = input.nextLine().split(":");
            nbTrucks = Integer.parseInt(splitted[1].trim());
            splitted = input.nextLine().split(":");
            truckCapacity = Integer.parseInt(splitted[1].trim());
            for (int i = 0; i < 2; ++i)
                input.nextLine();
            requiredNodes = new ArrayList<Integer>();
            // Even indices store direct edges, and odd indices store reverse edges
            demandsData = new int[2 * nbRequiredEdges];
            costsData = new int[2 * nbRequiredEdges];
            originsData = new int[2 * nbRequiredEdges];
            destinationsData = new int[2 * nbRequiredEdges];
            nodeNeighborsData = new int[nbNodes][];
            for (int j = 0; j < nbNodes; ++j) {
                nodeNeighborsData[j] = new int[nbNodes];
            }
            input.nextLine();
            for (int i = 0; i < nbRequiredEdges; ++i) {
                splitted = input.nextLine().split(" ");
                int firstNode = Integer.parseInt(splitted[2].substring(0, splitted[2].length() - 1));
                int secondNode = Integer.parseInt(splitted[3].substring(0, splitted[3].length() - 1));
                int cost = Integer.parseInt(splitted[7]);
                int demand = Integer.parseInt(splitted[11]);
                costsData[2 * i] = cost;
                demandsData[2 * i] = demand;
                costsData[2 * i + 1] = cost;
                demandsData[2 * i + 1] = demand;
                originsData[2 * i] = firstNode;
                destinationsData[2 * i] = secondNode;
                originsData[2 * i + 1] = secondNode;
                destinationsData[2 * i + 1] = firstNode;
                if (!requiredNodes.contains(firstNode)) {
                    requiredNodes.add(firstNode);
                }
                if (!requiredNodes.contains(secondNode)) {
                    requiredNodes.add(secondNode);
                }
                nodeNeighborsData[firstNode - 1][secondNode - 1] = cost;
                nodeNeighborsData[secondNode - 1][firstNode - 1] = cost;
            }
            if (nbNotRequiredEdges > 0) {
                input.nextLine();
                for (int i = 0; i < nbNotRequiredEdges; ++i) {
                    splitted = input.nextLine().split(" ");
                    int firstNode = Integer.parseInt(splitted[2].substring(0, splitted[2].length() - 1));
                    int secondNode = Integer.parseInt(splitted[3].substring(0, splitted[3].length() - 1));
                    int cost = Integer.parseInt(splitted[7]);
                    nodeNeighborsData[firstNode - 1][secondNode - 1] = cost;
                    nodeNeighborsData[secondNode - 1][firstNode - 1] = cost;
                }
            }
            splitted = input.nextLine().split(":");
            depotNode = Integer.parseInt(splitted[1].trim());
            nbRequiredNodes = requiredNodes.size();
            int[] shortestPath = new int[nbNodes];
            findRequiredPaths(shortestPath);
            findDistanceBetweenEdges();
            findDistanceToDepot();
            findDistanceFromDepot(shortestPath);
        }
    }

    // Finds the shortest path from one node "origin" to all the other nodes of the
    // graph
    // thanks to the Dijkstra's algorithm
    public int minDistance(int[] shortestPath, boolean[] sptSet) {
        int min = Integer.MAX_VALUE, minIndex = -1;
        for (int i = 0; i < nbNodes; ++i)
            if (sptSet[i] == false && shortestPath[i] <= min) {
                min = shortestPath[i];
                minIndex = i;
            }
        return minIndex;
    }

    public void shortestPathFinder(int origin, int[] shortestPath) {
        int[] currentNeighbors;
        int currentNode, distance;
        boolean[] sptSet = new boolean[nbNodes];
        for (int i = 0; i < nbNodes; ++i) {
            shortestPath[i] = Integer.MAX_VALUE;
            sptSet[i] = false;
        }
        shortestPath[origin - 1] = 0;
        for (int count = 0; count < nbNodes; ++count) {
            currentNode = minDistance(shortestPath, sptSet);
            sptSet[currentNode] = true;
            currentNeighbors = nodeNeighborsData[currentNode];
            for (int neighbor = 0; neighbor < nbNodes; ++neighbor) {
                if (currentNeighbors[neighbor] != 0) {
                    distance = currentNeighbors[neighbor];
                    if (!sptSet[neighbor] && shortestPath[currentNode] + distance < shortestPath[neighbor]) {
                        shortestPath[neighbor] = shortestPath[currentNode] + distance;
                    }
                }
            }
        }
    }

    public void findRequiredPaths(int[] shortestPath) {
        nodesDistData = new int[nbRequiredNodes][];
        for (int i = 0; i < nbRequiredNodes; ++i) {
            nodesDistData[i] = new int[nbNodes];
        }
        for (int i = 0; i < nbRequiredNodes; ++i) {
            int node = requiredNodes.get(i);
            shortestPathFinder(node, shortestPath);
            for (int j = 0; j < nbNodes; ++j) {
                nodesDistData[i][j] = shortestPath[j];
            }
        }
    }

    public void findDistanceBetweenEdges() {
        edgesDistData = new int[2 * nbRequiredEdges][];
        for (int i = 0; i < 2 * nbRequiredEdges; ++i) {
            edgesDistData[i] = new int[2 * nbRequiredEdges];
        }
        for (int i = 0; i < 2 * nbRequiredEdges; ++i) {
            for (int j = 0; j < 2 * nbRequiredEdges; ++j) {
                if (destinationsData[i] == originsData[j]) {
                    edgesDistData[i][j] = 0;
                } else {
                    for (int k = 0; k < nbRequiredNodes; ++k) {
                        if (requiredNodes.get(k) == destinationsData[i]) {
                            edgesDistData[i][j] = nodesDistData[k][originsData[j] - 1];
                        }
                    }
                }
            }
        }
    }

    public void findDistanceToDepot() {
        distToDepotData = new int[2 * nbRequiredEdges];
        for (int i = 0; i < 2 * nbRequiredEdges; ++i) {
            if (destinationsData[i] == depotNode) {
                distToDepotData[i] = 0;
            } else {
                for (int k = 0; k < nbRequiredNodes; ++k) {
                    if (requiredNodes.get(k) == destinationsData[i]) {
                        distToDepotData[i] = nodesDistData[k][depotNode - 1];
                    }
                }
            }
        }
    }

    public void findDistanceFromDepot(int[] shortestPath) {
        distFromDepotData = new int[2 * nbRequiredEdges];
        for (int i = 0; i < 2 * nbRequiredEdges; ++i) {
            if (depotNode == originsData[i]) {
                distFromDepotData[i] = 0;
            } else {
                boolean depotIsRequiredNode = false;
                for (int k = 0; k < nbRequiredNodes; ++k) {
                    if (requiredNodes.get(k) == depotNode) {
                        depotIsRequiredNode = true;
                        distFromDepotData[i] = nodesDistData[k][originsData[i] - 1];
                    }
                }
                if (!depotIsRequiredNode) {
                    shortestPathFinder(depotNode, shortestPath);
                    distFromDepotData[i] = shortestPath[originsData[i] - 1];
                }
            }
        }
    }

    public static void main(String[] args) {
        if (args.length < 1) {
            System.err.println("Usage: java CapacitatedArcRouting 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";

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