Line data    Source code

       1             : #include "Static/algos/FrankWolfe.hpp"
       2             : 
       3             : #include <chrono>
       4             : #include <cmath>
       5             : #include <iomanip>
       6             : #include <iostream>
       7             : #include <memory>
       8             : #include <utility>
       9             : 
      10             : #include "Opt/QuadraticSolver.hpp"
      11             : #include "Opt/UnivariateSolver.hpp"
      12             : #include "Static/Solution.hpp"
      13             : #include "Static/algos/AllOrNothing.hpp"
      14             : #include "Static/algos/DijkstraAoN.hpp"
      15             : 
      16             : using namespace std;
      17             : using namespace Static;
      18             : 
      19             : typedef Network::Node Node;
      20             : typedef Network::Edge Edge;
      21             : 
      22             : typedef chrono::high_resolution_clock hrc;
      23             : 
      24           3 : FrankWolfe::FrankWolfe(
      25             :     AllOrNothing          &aon_,
      26             :     Opt::UnivariateSolver &solver_,
      27             :     Log::ProgressLogger   &logger_
      28           3 : ):
      29             :     aon(aon_),
      30             :     solver(solver_),
      31           3 :     logger(logger_) {}
      32             : 
      33           3 : void FrankWolfe::setStopCriteria(Time e) {
      34           3 :     epsilon = e;
      35           3 : }
      36             : 
      37           0 : void FrankWolfe::setIterations(unsigned it) {
      38           0 :     iterations = it;
      39           0 : }
      40             : 
      41           3 : Solution FrankWolfe::solve(
      42             :     const Network  &network,
      43             :     const Demand   &dem,
      44             :     const Solution &startingSolution
      45             : ) {
      46           3 :     supply = &network;
      47           3 :     demand = &dem;
      48             : 
      49           3 :     logger << std::fixed << std::setprecision(3);
      50             : 
      51           3 :     xn = startingSolution;
      52           3 :     zn = supply->evaluate(xn);
      53             : 
      54           3 :     xStar      = xn;
      55           3 :     alpha      = 0.0;
      56           3 :     lowerBound = 0;
      57             : 
      58           3 :     double expectedIterations = getExpectedIterations();
      59             : 
      60           3 :     double estimation1 = 0.176 * expectedIterations;
      61             : 
      62           3 :     const double ETA_DECAY = 0.1 / min(expectedIterations, (double)iterations);
      63             : 
      64           3 :     logger << Log::ProgressLogger::Elapsed(0)
      65           3 :            << Log::ProgressLogger::Progress(0)
      66           3 :            << Log::ProgressLogger::ETA(estimation1)
      67           3 :            << Log::ProgressLogger::StartText()
      68           3 :            << "it\talpha\tzn\tdelta\tlowerBound\tAbsGap\tRelGap\tt1\tt2"
      69           3 :            << Log::ProgressLogger::EndMessage();
      70             : 
      71           3 :     double t1 = 0, t2 = 0;
      72             : 
      73           3 :     Time initialAbsoluteGap = 0;
      74             : 
      75           3 :     const hrc::time_point tStart = hrc::now();
      76             : 
      77           6 :     Solution xnPrev = xn;
      78           3 :     Flow     znPrev = zn;
      79           8 :     for(size_t it = 0; it < iterations; ++it) {
      80           8 :         Time delta       = znPrev - zn;
      81           8 :         Time absoluteGap = zn - lowerBound;
      82           8 :         Time relativeGap = absoluteGap / zn;
      83             : 
      84           8 :         const hrc::time_point t = hrc::now();
      85             : 
      86           8 :         double elapsed = (double)chrono::duration_cast<chrono::nanoseconds>(t - tStart).count() * 1e-9;
      87             : 
      88             :         // Progress
      89           8 :         if(it == 0) initialAbsoluteGap = absoluteGap;
      90           8 :         Time progressEpsilon = pow(
      91           8 :             log(absoluteGap / initialAbsoluteGap) / log(epsilon / initialAbsoluteGap),
      92             :             12
      93           8 :         );
      94           8 :         Time progressIterations = Time(it + 1) / iterations;
      95           8 :         Time progress           = max(progressEpsilon, progressIterations);
      96          14 :         progress                = max(0.0, min(1.0, progress));
      97             : 
      98             :         // ETA
      99           8 :         double estimation = estimation1;
     100           8 :         if(it > 0 && progress > 0) {
     101           5 :             double estimation2 = elapsed / progress;
     102           5 :             estimation1        = (1 - ETA_DECAY) * estimation1 + ETA_DECAY * estimation2;
     103           5 :             estimation         = (1 - progress * progress) * estimation1 + progress * progress * estimation2;
     104             :         }
     105           8 :         double eta = estimation - elapsed;
     106             : 
     107           8 :         progress = elapsed / estimation;
     108             : 
     109           8 :         logger << Log::ProgressLogger::Elapsed(elapsed)
     110           8 :                << Log::ProgressLogger::Progress(progress)
     111           8 :                << Log::ProgressLogger::ETA(eta)
     112           8 :                << Log::ProgressLogger::StartText()
     113           8 :                << it
     114           8 :                << "\t" << alpha
     115           8 :                << "\t" << zn
     116           8 :                << "\t" << delta
     117           8 :                << "\t" << lowerBound
     118           8 :                << "\t" << absoluteGap
     119           8 :                << "\t" << relativeGap
     120           8 :                << "\t" << t1
     121           8 :                << "\t" << t2
     122           8 :                << Log::ProgressLogger::EndMessage();
     123             : 
     124           8 :         if(absoluteGap <= epsilon) {
     125           2 :             logger << Log::ProgressLogger::Elapsed(elapsed)
     126           2 :                    << Log::ProgressLogger::Progress(progress)
     127           2 :                    << Log::ProgressLogger::ETA(eta)
     128           2 :                    << Log::ProgressLogger::StartText()
     129           2 :                    << "Met relative gap criteria"
     130           2 :                    << Log::ProgressLogger::EndMessage();
     131           3 :             return xn;
     132             :         }
     133           6 :         if(znPrev < zn) {
     134           1 :             logger << Log::ProgressLogger::Elapsed(elapsed)
     135           1 :                    << Log::ProgressLogger::Progress(progress)
     136           1 :                    << Log::ProgressLogger::ETA(eta)
     137           1 :                    << Log::ProgressLogger::StartText()
     138           1 :                    << "Solution got worse"
     139           1 :                    << Log::ProgressLogger::EndMessage();
     140           1 :             return xnPrev;
     141             :         }
     142             : 
     143           5 :         xnPrev = xn;
     144           5 :         znPrev = zn;
     145             : 
     146           5 :         hrc::time_point a = hrc::now();
     147             : 
     148           5 :         step1();
     149             : 
     150           5 :         hrc::time_point b = hrc::now();
     151             : 
     152           5 :         xn = step2(xStar);
     153             : 
     154           5 :         hrc::time_point c = hrc::now();
     155             : 
     156           5 :         t1 = (double)chrono::duration_cast<chrono::nanoseconds>(b - a).count() * 1e-9;
     157           5 :         t2 = (double)chrono::duration_cast<chrono::nanoseconds>(c - b).count() * 1e-9;
     158             : 
     159           5 :         zn = supply->evaluate(xn);
     160             :     }
     161             : 
     162           0 :     return xn;
     163             : }
     164             : 
     165           3 : double FrankWolfe::getExpectedIterations() {
     166           3 :     double linearWithIterations = pow(-log10(epsilon / zn), 12);  // This variable has a linear relation with number of iterations
     167           3 :     double expectedIterations   = linearWithIterations / 92259.0869806;
     168           3 :     return min(expectedIterations, (double)iterations);
     169             : }
     170             : 
     171           5 : Solution FrankWolfe::step1() {
     172           5 :     SolutionBase xAoN = aon.solve(*supply, *demand, xn);
     173             : 
     174          10 :     unordered_set<Edge::ID>        edgeIDs;
     175          10 :     const unordered_set<Edge::ID> &xnEdges   = xn.getEdges();
     176          10 :     const unordered_set<Edge::ID> &xAoNEdges = xAoN.getEdges();
     177           5 :     edgeIDs.insert(xnEdges.begin(), xnEdges.end());
     178           5 :     edgeIDs.insert(xAoNEdges.begin(), xAoNEdges.end());
     179             : 
     180             :     // Update lower bound
     181           5 :     Time zApprox = zn;
     182          17 :     for(const Edge::ID &eid: edgeIDs) {
     183          12 :         Edge &e = supply->getEdge(eid);
     184             : 
     185          12 :         Flow xna   = xn.getFlowInEdge(eid);
     186          12 :         Flow xAoNa = xAoN.getFlowInEdge(eid);
     187             : 
     188          12 :         zApprox += e.calculateCost(xn) * (xAoNa - xna);
     189             :     }
     190           5 :     lowerBound = max(lowerBound, zApprox);
     191             : 
     192             :     // Get xStar
     193           5 :     xStar = xAoN;
     194             : 
     195          10 :     return xStar;
     196             : }
     197             : 
     198           5 : Solution FrankWolfe::step2(const Solution &xstar) {
     199             :     // clang-format off
     200         155 :     Opt::UnivariateSolver::Problem p = [
     201           5 :         &supply = as_const(supply),
     202           5 :         &xn     = as_const(xn),
     203           5 :         &xstar  = as_const(xstar)
     204         300 :     ](Opt::UnivariateSolver::Var a) -> Time {
     205         150 :         Solution      x = Solution::interpolate(xn, xstar, a);
     206         150 :         Time c = supply->evaluate(x);
     207         300 :         return c;
     208           5 :     };
     209             :     // clang-format on
     210             : 
     211           5 :     alpha = solver.solve(p);
     212          10 :     alpha = max(0.0, min(1.0, alpha));
     213             : 
     214           5 :     Solution x = Solution::interpolate(xn, xstar, alpha);
     215             : 
     216          10 :     return x;
     217             : }