DAMYAAN - Editorial

PROBLEM LINK:

Practice
Contest: Division 1
Contest: Division 2
Contest: Division 3
Contest: Division 4

Author: wasd2401
Tester: wasd2401
Editorialist: iceknight1093

DIFFICULTY:

TBD

PREREQUISITES:

Dynamic programming on a DAG, probability

PROBLEM:

You’re given a DAG, each of whose edges has a probability.
When at a node u,

  • Edge E out of u will exist with probability p_E.
  • You can choose any existing edge and move along it; or wait for the next day for a (potentially) new set of edges to appear.

With optimal strategy, find the minimum expected time needed to travel from 1 to N.

EXPLANATION:

Let \mathbb{E}[u] denote the minimum expected time to reach node N if you’re currently at node u (or \infty, if you can’t reach N at all).

What should your optimal strategy from here look like?

Answer

An initial idea would be to pick a single edge, and keep waiting at u till it shows up and then use it.
This isn’t necessarily optimal though: for instance, if you have two edges that are similarly good, a better strategy might be to wait for at least one of them to show up, and then take whichever one shows up first.
Of course, the same logic dictates that choosing two edges isn’t necessarily optimal either: perhaps you might want a set of 3, or 4, or even more candidates.

More generally, there will exist some subset S of outgoing edges, such that you will choose the first one among them that appears (if multiple appear at the same time, of course you’ll need to have some tie-breaking rule as for which one to choose, i.e, a preferential order among these edges).

The questions we have to answer now are:

  • What will this subset S of candidate edges look like; and
  • What the tie-breaking rule should be.

The second question is easier to answer: notice that moving to vertex v (when the necessary edge exists) will result in an expected time of \mathbb{E}[v]+2 days to reach N.
So, if multiple edges exist, it’s optimal to take the one among them with minimum \mathbb{E}[v] value.

Now, we look at what exactly a subset of edges gets us to.
Suppose our candidate edges are to vertices v_1, v_2, \ldots, v_k, and exist with probabilities p_1, p_2, \ldots, p_k.
Let’s also suppose that this is our order of preference, i.e, \mathbb{E}[v_i] \leq \mathbb{E}[v_j] for i \lt j.

Then, the i-th edge will be chosen if and only if:

  1. It exists on this day, which has a probability of p_i; and
  2. Every previous edge doesn’t exist, which has a probability of (1-p_1)(1-p_2)\ldots (1-p_{i-1})

Further, if every edge doesn’t exist, we wait for the next day.
This gives us the formula

\mathbb{E}[u] = p_1\cdot (\mathbb{E}[v_1]+2) + (1-p_1)p_2\cdot (\mathbb{E}[v_2]+2) + \ldots + (1-p_1)(1-p_2)\ldots (1-p_k)\cdot (\mathbb{E}[u] + 1)

which is a linear equation in \mathbb{E}[u] and can be solved for it, if all the \mathbb{E}[v_i] values are known.
(Using this equation, you can also formally prove that preferring small \mathbb{E}[v_i] is better via an exchange argument).


The only question remaining is, which subset of edges S do we choose?

Let’s list out every edge out of u, sorted by their \mathbb{E}[v] values in ascending order.
Then, it can be proved that the optimal subset will be some prefix of these sorted edges.
This makes sense intuitively: if you’re willing to take some edge when it appears, you should also be willing to take some edge that results in less expected time if it appears.
A formal proof can, once again, be made algebraically by looking at the formula above, and seeing what happens when you insert a new term in the middle (with smaller \mathbb{E}[v] than some existing value).

This gives us a pretty simple implementation: sort the edges going out of u by expected time, and try every prefix of them as the optimal subset; then take the best answer.
Since the graph is a DAG, we don’t run into any cyclic dependencies, so this is enough.

TIME COMPLEXITY:

\mathcal{O}(N + M\log M) per test.

CODE:

Author's code (C++)
/*

*       *  *  ***       *       *****
 *     *   *  *  *     * *        *
  *   *    *  ***     *****       *
   * *     *  * *    *     *   *  *
    *      *  *  *  *       *   **

                                 *
                                * *
                               *****
                              *     *
        *****                *       *
      _*     *_
     | * * * * |                ***
     |_*  _  *_|               *   *
       *     *                 *  
        *****                  *  **
       *     *                  ***
  {===*       *===}
      *  IS   *                 ***
      *  IT   *                *   *
      * RATED?*                *  
      *       *                *  **
      *       *                 ***
       *     *
        *****                  *   *
                               *   *
                               *   *
                               *   *
                                ***   

*/

// 2 Years Tribute to Competitive Programming

#include <bits/stdc++.h>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/tree_policy.hpp>

using namespace __gnu_pbds;
using namespace std;

#define osl tree<ll, null_type, less<ll>, rb_tree_tag, tree_order_statistics_node_update>
#define ll long long
#define ld long double
#define forl(i, a, b) for(ll i = a; i < b; i++)
#define rofl(i, a, b) for(ll i = a; i > b; i--)
#define fors(i, a, b, c) for(ll i = a; i < b; i += c)
#define fora(x, v) for(auto x : v)
#define vl vector<ll>
#define vb vector<bool>
#define pub push_back
#define pob pop_back
#define fbo find_by_order
#define ook order_of_key
#define yesno(x) cout << ((x) ? "YES" : "NO")
#define all(v) v.begin(), v.end()

const ll N = 2e5 + 4;
const ll mod = 1e18 + 7;
// const ll mod = 998244353;

vector<vl> v(N);
vector<vector<ld>> prob(N);
vector<ld> dp(N);
ll modinverse(ll a) {
	ll m = mod, y = 0, x = 1;
	while (a > 1) {
		ll q = a / m;
		ll t = m;
		m = a % m;
		a = t;
		t = y;
		y = x - q * y;
		x = t;
	}
	if (x < 0) x += mod;
	return x;
}
ll gcd(ll a, ll b) {
	if (b == 0)
		return a;
	return gcd(b, a % b);
}
ll lcm(ll a, ll b) {
	return (a / gcd(a, b)) * b;
}
bool poweroftwo(ll n) {
	return !(n & (n - 1));
}
ll power(ll a, ll b, ll md = mod) {
	ll product = 1;
	a %= md;
	while (b) {
		if (b & 1) product = (product * a) % md;
		a = (a * a) % md;
		b /= 2;
	}
	return product % md;
}
void barfi(ll n){
	vector<vector<ld>> v1;
	forl(i,0,v[n].size()){
		ll x=v[n][i];
		if(dp[x]==mod) barfi(x);
		if(dp[x]!=-1){
			v1.pub({1+dp[x],prob[n][i]});
		}
	}
	if(v1.empty()){
		dp[n]=-1;
		return;
	}
	sort(all(v1));
	ld p1=1;
	ld sum=0;
	fora(x,v1){
		sum+=x[1]*p1*x[0];
		p1*=(1-x[1]);
		dp[n]=min(dp[n],(sum+1)/(1-p1));
		// if(n==3) cout<<x[0]<<' ';
	}
}
void panipuri() {
	ll n, m = 0, k = -1, c = 0, sum = 0, q = 0, ans = 0, p = 1;
	string s;
	bool ch = true;
	cin >> n>>m;
	forl(i,1,n+1){
		v[i].clear();
		prob[i].clear();
		dp[i]=mod;
	}
	dp[n]=0;
	forl(i,0,m){
		ll x,y;
		cin>>x>>y>>p>>q;
		v[x].pub(y);
		ld p1=p;
		prob[x].pub(p1/q);
	}
	// forl(i,1,10) cout<<dp[i]<<' ';
	barfi(1);
	cout<<fixed;
	cout<<setprecision(10);
	cout<<dp[1];
	return;
}
int main() {
	ios::sync_with_stdio(false);
	cin.tie(NULL);
	#ifndef ONLINE_JUDGE
	freopen("input.txt", "r", stdin);
	freopen("output.txt", "w", stdout);
	#endif
	int laddu = 1;
	cin >> laddu;
	forl(i, 1, laddu + 1) {
		// cout << "Case #" << i << ": ";
		panipuri();
		cout << '\n';
	}
}
Editorialist's code (C++)
// #pragma GCC optimize("O3,unroll-loops")
// #pragma GCC target("avx2,bmi,bmi2,lzcnt,popcnt")
#include "bits/stdc++.h"
using namespace std;
using ll = long long int;
mt19937_64 RNG(chrono::high_resolution_clock::now().time_since_epoch().count());

int main()
{
    ios::sync_with_stdio(false); cin.tie(0);

    int t; cin >> t;
    while (t--) {
        int n, m; cin >> n >> m;
        vector<vector<array<int, 4>>> adj(n);
        for (int i = 0; i < m; ++i) {
            int u, v, w, p, q; cin >> u >> v >> p >> q;
            w = 1;
            adj[--u].push_back({--v, w, p, q});
        }
        
        using ftype = double;
        vector<ftype> dp(n, -10);
        dp[n-1] = 0;
        auto calc = [&] (const auto &self, int u) -> void {
            if (u == n-1) return;
            if (dp[u] > -5) return;
            dp[u] = -1;
            vector<array<ftype, 2>> poss;
            for (auto &[v, w, p, q] : adj[u]) {
                self(self, v);
                if (dp[v] < -0.5) continue;

                // with probability p/q, w + dp[v]
                poss.push_back({1 + w + dp[v], ftype(p)/ftype(q)});
            }
            sort(begin(poss), end(poss)); // think about breaking ties?

            if (poss.empty()) return;
            dp[u] = 1e18;
            int k = poss.size();
            ftype mul = 1, sm = 0;
            for (int i = 0; i < k; ++i) {
                // take the prefix till here
                // (p1, x1), (p2, x2), ..., (pk, xk)
                // p1x2 + (1-p1)p2x2 + (1-p1)(1-p2)p3x3 + ... + prod(1-pi)(1 + dp[u]) = dp[u]
                auto [xi, pi] = poss[i];
                sm += mul*xi*pi;
                mul *= (1 - pi);

                dp[u] = min(dp[u], (sm + mul) / (1 - mul));
            }
        };
        calc(calc, 0);
        cout << fixed << setprecision(15) << dp[0] << '\n';
    }
}
3 Likes