/*
 * Travelling salesperson with time windows, using an open variable list.
 *
 * Companion to the post "Planning with open variable lists in Gecode".
 * Written for Gecode 6.5.0. The A* and CABS engines require a Gecode build
 * configured with GECODE_ENABLE_EXPERIMENTAL=ON (or --enable-experimental).
 *
 * The instances, instance reader, bound, and dominance test follow
 * examples/tsptw.cpp in Gecode 6.5.0, which uses fixed arrays for the tour.
 * Here the tour is an open list whose positions, visit times, and travel
 * times are created as the tour grows.
 *
 * Copyright: Mikael Zayenz Lagerkvist, 2026.
 * Distributed under the MIT license, as is Gecode.
 */

#include <gecode/driver.hh>
#include <gecode/int.hh>
#include <gecode/minimodel.hh>

#include <algorithm>

using namespace Gecode;

namespace {

  /// Built-in TSPTW instance
  class Builtin {
  public:
    /// Number of nodes, node 0 is the depot
    int n;
    /// Travel times
    const int* d;
    /// Ready and due time for each node
    const int* w;
  };

  /// Instance rbg010a (Ascheuer)
  const int PA_n = 11;
  const int PA_d[PA_n*PA_n] = {
    0,0,0,0,0,0,0,0,0,0,0,
    71,0,85,85,85,96,85,85,97,85,77,
    50,65,0,77,77,62,77,77,62,77,74,
    64,74,88,0,88,80,88,88,82,88,87,
    44,54,68,68,0,60,68,68,62,68,67,
    51,75,51,51,51,0,51,51,80,51,65,
    53,64,78,78,78,68,0,78,69,78,76,
    51,67,69,69,69,75,69,0,75,69,63,
    43,61,71,71,71,43,71,71,0,71,68,
    53,64,78,78,78,68,78,78,69,0,80,
    42,66,52,52,52,70,52,52,70,52,0
  };
  const int PA_w[2*PA_n] = {
    0,9396, 0,865, 399,1299, 550,1450, 685,1585, 819,1719,
    1208,2108, 1412,2312, 1701,2330, 2413,3313, 3798,4698
  };

  /// Instance n20w20.001 (Dumas et al.)
  const int PB_n = 21;
  const int PB_d[PB_n*PB_n] = {
    0,19,17,34,7,20,10,17,28,15,23,29,23,29,21,20,9,16,21,13,12,
    19,0,10,41,26,3,27,25,15,17,17,14,18,48,17,6,21,14,17,13,31,
    17,10,0,47,23,13,26,15,25,22,26,24,27,44,7,5,23,21,25,18,29,
    34,41,47,0,36,39,25,51,36,24,27,38,25,44,54,45,25,28,26,28,27,
    7,26,23,36,0,27,11,17,35,22,30,36,30,22,25,26,14,23,28,20,10,
    20,3,13,39,27,0,26,27,12,15,14,11,15,49,20,9,20,11,14,11,30,
    10,27,26,25,11,26,0,26,31,14,23,32,22,25,31,28,6,17,21,15,4,
    17,25,15,51,17,27,26,0,39,31,38,38,38,34,13,20,26,31,36,28,27,
    28,15,25,36,35,12,31,39,0,17,9,2,11,56,32,21,24,13,11,15,35,
    15,17,22,24,22,15,14,31,17,0,9,18,8,39,29,21,8,4,7,4,18,
    23,17,26,27,30,14,23,38,9,9,0,11,2,48,33,23,17,7,2,10,27,
    29,14,24,38,36,11,32,38,2,18,11,0,13,57,31,20,25,14,13,17,36,
    23,18,27,25,30,15,22,38,11,8,2,13,0,47,34,24,16,7,2,10,26,
    29,48,44,44,22,49,25,34,56,39,48,57,47,0,46,48,31,42,46,40,21,
    21,17,7,54,25,20,31,13,32,29,33,31,34,46,0,11,29,28,32,25,33,
    20,6,5,45,26,9,28,20,21,21,23,20,24,48,11,0,23,19,22,17,32,
    9,21,23,25,14,20,6,26,24,8,17,25,16,31,29,23,0,11,15,9,10,
    16,14,21,28,23,11,17,31,13,4,7,14,7,42,28,19,11,0,5,3,21,
    21,17,25,26,28,14,21,36,11,7,2,13,2,46,32,22,15,5,0,8,25,
    13,13,18,28,20,11,15,28,15,4,10,17,10,40,25,17,9,3,8,0,19,
    12,31,29,27,10,30,4,27,35,18,27,36,26,21,33,32,10,21,25,19,0
  };
  const int PB_w[2*PB_n] = {
    0,408, 62,68, 181,205, 306,324, 214,217, 51,61,
    102,129, 175,186, 250,263, 3,23, 21,49, 79,90,
    78,96, 140,154, 354,386, 42,63, 2,13, 24,42,
    20,33, 9,21, 275,300
  };

  /// Two customers, where the best tour waits at customer 1
  const int PC_n = 3;
  const int PC_d[PC_n*PC_n] = {
    0,2,9,
    4,0,1,
    2,8,0
  };
  const int PC_w[2*PC_n] = {
    0,20, 5,10, 0,12
  };

  /// The instances
  const Builtin instances[] = {
    {PA_n, PA_d, PA_w},
    {PB_n, PB_d, PB_w},
    {PC_n, PC_d, PC_w}
  };
  /// Number of instances
  const unsigned int n_instances = sizeof(instances)/sizeof(Builtin);

  /// Reader spec for files from López-Ibáñez and Blum's collection
  const char* instance_spec = R"(format tsptw

comment "#"
n: 2..1000
dist[n][n]: 0..gecode.int_max / 2 / n

type time = 0..gecode.int_max / 2
repeat n as nodes {
  each node
  ready: time
  due: time
}

check node.ready <= node.due "the time window of a node is invalid"
check nodes[0].ready == 0 "the depot (node 0) must open at time 0"
)";

  /// TSPTW instance
  class Problem {
  public:
    /// Number of nodes, node 0 is the depot (0 if invalid)
    int n;
    /// Travel times
    IntSharedArray d;
    /// Ready and due time for each node
    IntSharedArray w;
    /// Initialize from a built-in instance number or an instance file
    Problem(const InstanceOptions& opt) : n(0) {
      unsigned int b;
      if (opt.builtin(b)) {
        if (b >= n_instances)
          return;
        n = instances[b].n;
        d = IntSharedArray(IntArgs(n*n, instances[b].d));
        w = IntSharedArray(IntArgs(2*n, instances[b].w));
      } else if (opt.is_file()) {
        Instance::Data data = opt.data();
        int m = static_cast<int>(data.integer("n"));
        IntArgs ready = data.ints("ready"), due = data.ints("due");
        IntArgs win(2*m);
        for (int i=0; i<m; i++) {
          win[2*i] = ready[i]; win[2*i+1] = due[i];
        }
        d = IntSharedArray(data.ints("dist"));
        w = IntSharedArray(win);
        n = m;
      }
    }
    /// Whether an instance was found
    bool valid(void) const {
      return n > 0;
    }
  };

}

/**
 * \brief Travelling salesperson with time windows as an open list
 *
 * Starting and ending at the depot, visit every customer once such that
 * service at each customer starts within its time window (waiting is
 * allowed). Minimize the total travel time.
 *
 * The open list x holds the customer at each position. When a position is
 * created, its transition creates the start of service and the travel time
 * so far for that position, and connects them to the previous position.
 * The tour may end only after every customer has been visited.
 */
class OpenTSPTW : public IntMinimizeScript {
protected:
  /// The instance
  const Problem p;
  /// Travel time from the depot, ready time, and due time of each node
  IntSharedArray from, ready, due;
  /// Travel time back to the depot, directly and via other nodes
  IntSharedArray to, back;
  /// Travel time between each pair of distinct customers
  TupleSet legs;
  /// Customer at each position of the tour
  OpenIntVarList x;
  /// Start of service at each position
  OpenIntVarList t;
  /// Travel time from the depot to each position
  OpenIntVarList c;
  /// Total travel time, including the return to the depot
  IntVar total;
  /// Return a variable equal to a[y]
  IntVar lookup(const IntSharedArray& a, IntVar y) {
    IntVar v(*this, 0, Int::Limits::max);
    element(*this, a, y, v);
    return v;
  }
public:
  /// Search engines
  enum {
    SEARCH_BAB, SEARCH_ASTAR, SEARCH_CABS
  };
  /// The actual problem
  OpenTSPTW(const InstanceOptions& opt)
    : IntMinimizeScript(opt), p(opt),
      x(*this, IntSet(1,p.n-1),
        [](Space& home, OpenIntVarList x, int i) {
          static_cast<OpenTSPTW&>(home).visit(x,i);
        }, p.n-1),
      t(*this, IntSet(0,p.w[1]), p.n-1),
      c(*this, IntSet(0,p.w[1]), p.n-1),
      total(*this, 0, p.w[1]) {
    const int n = p.n;
    IntArgs f(n), r(n), u(n), b(n), s(n*n);
    for (int i=0; i<n*n; i++)
      s[i] = p.d[i];
    // Shortest travel times, as the travel times need not be metric
    for (int k=0; k<n; k++)
      for (int i=0; i<n; i++)
        for (int j=0; j<n; j++)
          s[i*n+j] = std::min(s[i*n+j], s[i*n+k]+s[k*n+j]);
    for (int i=0; i<n; i++) {
      f[i] = p.d[i]; b[i] = p.d[i*n];
      r[i] = p.w[2*i]; u[i] = p.w[2*i+1];
    }
    from = IntSharedArray(f); ready = IntSharedArray(r);
    due = IntSharedArray(u); to = IntSharedArray(b);
    IntArgs shortest(n);
    for (int i=0; i<n; i++)
      shortest[i] = s[i*n];
    back = IntSharedArray(shortest);
    legs = TupleSet(3);
    for (int j=1; j<n; j++)
      for (int k=1; k<n; k++)
        if (j != k)
          legs.add({j, k, p.d[j*n+k]});
    legs.finalize();

    rel(*this, t.length() == x.length());
    rel(*this, c.length() == x.length());
    distinct(*this, x);
    // Every tour visits at least one customer
    rel(*this, x.length() >= 1);
    // Assign the next customer on every alternative: CABS prefix dominance
    // does not account for value exclusions on the unassigned suffix.
    branch(*this, x, INT_VAR_NONE(), INT_VALUES_MIN());
  }
  /// Connect the new position \a i to the previous one
  void visit(OpenIntVarList x, int i) {
    t.materialize(*this, i+1);
    c.materialize(*this, i+1);
    if (failed())
      return;
    if (i == 0) {
      element(*this, from, x[0], c[0]);
      rel(*this, t[0] == max(lookup(ready,x[0]), c[0]));
    } else {
      IntVar leg(*this, 0, p.w[1]);
      extensional(*this, IntVarArgs({x[i-1], x[i], leg}), legs);
      rel(*this, c[i] == c[i-1] + leg);
      rel(*this, t[i] == max(lookup(ready,x[i]), t[i-1] + leg));
    }
    rel(*this, t[i] <= lookup(due,x[i]));
    // Every completion must still return to the depot in time
    IntVar home_trip = lookup(back,x[i]);
    rel(*this, t[i] + home_trip <= p.w[1]);
    rel(*this, total >= c[i] + home_trip);
    // Ending here requires every customer in the prefix, and a direct return
    BoolVar last = expr(*this, x.length() == i+1);
    IntVarArgs prefix(i+1);
    for (int j=0; j<=i; j++)
      prefix[j] = x[j];
    for (int k=1; k<p.n; k++)
      member(*this, prefix, IntVar(*this,k,k), imp(last));
    IntVar direct = lookup(to,x[i]);
    rel(*this, last >> (t[i] + direct <= p.w[1]));
    rel(*this, last >> (total == c[i] + direct));
  }
  /// Return cost
  virtual IntVar cost(void) const {
    return total;
  }
  /// Return the number of positions assigned from the start
  int prefix(void) const {
    int k = 0;
    while ((k < x.known_length()) && x[k].assigned() &&
           t[k].assigned() && c[k].assigned())
      k++;
    return k;
  }
  /// Store the customers in the first \a k positions in \a v
  void visited(int k, Support::BitSet<Region>& v) const {
    for (int i=0; i<k; i++)
      v.set(static_cast<unsigned int>(x[i].val()));
  }
  /// Whether the first \a k positions of this and \a s visit the same customers
  bool same_visited(const OpenTSPTW& s, int k) const {
    Region r;
    Support::BitSet<Region> v(r, static_cast<unsigned int>(p.n));
    visited(k, v);
    // Both prefixes visit k different customers, so inclusion is equality
    for (int i=0; i<k; i++)
      if (!v.get(static_cast<unsigned int>(s.x[i].val())))
        return false;
    return true;
  }
  /// Whether this partial tour dominates \a s
  bool dominates(const OpenTSPTW& s) const {
    int k = prefix();
    if ((k == 0) || (k != s.prefix()) || (x[k-1].val() != s.x[k-1].val()))
      return false;
    return (t[k-1].val() <= s.t[k-1].val()) &&
      (c[k-1].val() <= s.c[k-1].val()) && same_visited(s, k);
  }
  /// Return a key for dominance
  std::size_t key(void) const {
    int k = prefix();
    Region r;
    unsigned int n = static_cast<unsigned int>(p.n);
    Support::BitSet<Region> v(r, n);
    visited(k, v);
    // The customers in increasing order, so that equal sets hash equally
    std::size_t h = 0;
    for (unsigned int i = v.next(0); i < n; i = v.next(i+1))
      cmb_hash(h, i);
    cmb_hash(h, k > 0 ? x[k-1].val() : 0);
    return h;
  }
  /// Constructor for cloning \a s
  OpenTSPTW(OpenTSPTW& s)
    : IntMinimizeScript(s), p(s.p),
      from(s.from), ready(s.ready), due(s.due), to(s.to), back(s.back),
      legs(s.legs) {
    x.update(*this, s.x);
    t.update(*this, s.t);
    c.update(*this, s.c);
    total.update(*this, s.total);
  }
  /// Copy during cloning
  virtual Space* copy(void) {
    return new OpenTSPTW(*this);
  }
  /// Print solution
  virtual void print(std::ostream& os) const {
    os << "\tTour: 0";
    for (int i=0; i<x.known_length(); i++)
      os << " -> " << x[i];
    os << " -> 0" << std::endl
       << "\tTravel time: " << total << std::endl
       << "\tPositions created: " << x.known_length() << std::endl;
  }
};

int
main(int argc, char* argv[]) {
  InstanceOptions opt("OpenTSPTW", instance_spec);
  opt.instance("0");
  opt.solutions(0);
  opt.search(OpenTSPTW::SEARCH_BAB);
  opt.search(OpenTSPTW::SEARCH_BAB, "bab", "branch and bound");
#ifdef GECODE_HAS_EXPERIMENTAL
  opt.search(OpenTSPTW::SEARCH_ASTAR, "astar", "A*");
  opt.search(OpenTSPTW::SEARCH_CABS, "cabs", "complete anytime beam search");
#endif
  opt.parse(argc,argv);
  try {
    Problem p(opt);
    if (!p.valid()) {
      std::cerr << "Error: instance must be a number between 0 and "
                << n_instances-1 << " or a file" << std::endl;
      return 1;
    }
  } catch (const Instance::Error& e) {
    std::cerr << e.what() << std::endl;
    return 1;
  }
  switch (opt.search()) {
  case OpenTSPTW::SEARCH_BAB:
    IntMinimizeScript::run<OpenTSPTW,BAB,InstanceOptions>(opt);
    break;
#ifdef GECODE_HAS_EXPERIMENTAL
  case OpenTSPTW::SEARCH_ASTAR:
  case OpenTSPTW::SEARCH_CABS:
    {
      Search::HeuristicGuidance g(Search::HeuristicGuidance::MINIMIZE,
        [](const Space& s) -> std::int64_t {
          return static_cast<const OpenTSPTW&>(s).cost().min();
        });
      g.dominates = [](const Space& a, const Space& b) {
        return static_cast<const OpenTSPTW&>(a)
          .dominates(static_cast<const OpenTSPTW&>(b));
      };
      g.dominance_key = [](const Space& s) {
        return static_cast<const OpenTSPTW&>(s).key();
      };
      opt.guidance(&g);
      if (opt.search() == OpenTSPTW::SEARCH_ASTAR)
        IntMinimizeScript::run<OpenTSPTW,AStar,InstanceOptions>(opt);
      else
        IntMinimizeScript::run<OpenTSPTW,CABS,InstanceOptions>(opt);
    }
    break;
#endif
  }
  return 0;
}