Dominance problems (red-blue dominance and colorless dominance)

dominance.cc

#include <iostream>
#include <vector>
#include <cstdio>
#include <cstdlib>
#include <map>
#include <cmath>
#include <cstring>
#include <functional>
#include <algorithm>
#include <unordered_map>
#include <unordered_set>

#include <ctime>
double tick() {
  static clock_t oldtick;
  clock_t newtick = clock();
  double diff = 1.0*(newtick - oldtick) / CLOCKS_PER_SEC;
  oldtick = newtick;
  return diff;
}

using namespace std;

#define fst first
#define snd second
#define all(c) ((c).begin()), ((c).end())

template <class T>
ostream &operator<<(ostream &os, const vector<T> &x) {
  for (auto &a: x) os << a << " ";
  return os;
}

typedef vector<double> point;

namespace filter_impl {
  vector<point> A, B;
  // assume: I and J are sorted in 0-th coordinate
  void filter(vector<int> &I, const vector<int> &J, int d) {
    if (d == 1) { 
      // for 2-dim, apply plane sweep
      double th = -1.0/0.0;
      for (int i = I.size()-1, j = J.size()-1; i >= 0 || j >= 0; ) {
        if (j < 0 || (i >= 0 && A[I[i]][0] > B[J[j]][0])) {
          if (A[I[i]][1] <= th) I[i] = -1;
          --i;
        } else th = max(th, B[J[j--]][1]);
      }
      I.erase(remove(all(I), -1), I.end());
    } else if (I.size() <= 1 || J.size() <= 1) { 
      // for constant I and J, apply naive
      auto dominated = [&](int i) {
        for (int j: J) {
          bool dom = true;
          for (int k = 0; dom && k <= d; ++k) 
            if (A[i][k] > B[j][k]) dom = false;
          if (dom) return true;
        }
        return false;
      };
      I.erase(remove_if(all(I), dominated), I.end());
    } else {
      // otherwise, apply divide-and-conquer
      double th = B[J[rand() % J.size()]][d];
      vector<int> IL, IH, JL, JH;
      for_each(all(I), [&](int i) { 
        if (A[i][d] <= th) IL.push_back(i);
        else               IH.push_back(i);
      });
      I.clear();
      for_each(all(J), [&](int j) { 
        if      (B[j][d] < th)          JL.push_back(j);
        else if (B[j][d] > th)          JH.push_back(j);
        else if (JL.size() < JH.size()) JL.push_back(j);
        else                            JH.push_back(j);
      });
      filter(IL, JL, d);
      filter(IL, JH, d-1);
      filter(IH, JH, d);
      for (int l = 0, h = 0; l < IL.size() || h < IH.size(); ) {
        if      (l == IL.size()) I.push_back(IH[h++]);
        else if (h == IH.size()) I.push_back(IL[l++]);
        else if (A[IL[l]][0] < A[IH[h]][0]) I.push_back(IL[l++]);
        else                                I.push_back(IH[h++]);
      }
    }
  }
  vector<point> filter() {
    vector<int> I(A.size()); iota(all(I), 0);
    vector<int> J(B.size()); iota(all(J), 0);
    sort(all(I), [&](int i, int j) { return A[i][0] < A[j][0]; });
    sort(all(J), [&](int i, int j) { return B[i][0] < B[j][0]; });
    if (!A.empty()) filter(I, J, A[0].size()-1);
    vector<point> out;
    for_each(all(I), [&](int i) { out.push_back(A[i]); });
    return out;
  }
};
vector<point> filter(vector<point> A, vector<point> B) {
  filter_impl::A.swap(A);
  filter_impl::B.swap(B);
  return filter_impl::filter();
}
namespace dominance_impl { 
  vector<point> A;
  vector<point> dominance(int i, int j) {
    if (i   == j) return {};
    if (i+1 == j) return {A[i]};
    int m = (i + j) / 2;
    vector<point> AL = dominance(i, m);
    vector<point> AH = dominance(m, j);
    AL = filter(AL, AH);
    AL.insert(AL.end(), all(AH));
    return AL;
  }
};
// assume: A has no same points
vector<point> dominance(vector<point> A) {
  sort(all(A), [](const point &a, const point &b) {
    for (int d = 0; d < a.size(); ++d) 
      if (a[d] != b[d]) return a[d] < b[d];
    return false;
  });
  dominance_impl::A.swap(A);
  return dominance_impl::dominance(0, dominance_impl::A.size());
}

// verifier
vector<point> filter_naive(vector<point> A, vector<point> B) {
  auto dominate = [&](int i, int j) {
    for (int k = 0; k < A[i].size(); ++k)
      if (A[i][k] > B[j][k]) return false;
    return true;
  };
  vector<point> out;
  for (int i = 0; i < A.size(); ++i) {
    bool OK = true;
    for (int j = 0; j < B.size(); ++j) {
      if (dominate(i, j)) {
        OK = false;
        break;
      }
    }
    if (OK) out.push_back(A[i]);
  }
  return out;
}

vector<point> dominance_naive(vector<point> A) {
  vector<point> out;
  for (int i = 0; i < A.size(); ++i) {
    bool dominated = false;
    for (int j = 0; j < A.size(); ++j) {
      if (i == j) continue;
      bool dom = true;
      for (int d = 0; dom && d < A[i].size(); ++d)
        if (A[i][d] > A[j][d]) dom = false;
      if (dom) {
        dominated = true;
        break;
      }
    }
    if (!dominated) {
      out.push_back(A[i]);
    }
  }
  return out;
}

int main() {
  int n = 10000, d = 5;
  vector<point> P;
  for (int i = 0; i < n; ++i) {
    point p;
    double norm = 0;
    for (int k = 0; k < d; ++k) {
      p.push_back( rand()%n );
      norm += p.back() * p.back();
    }
    norm = sqrt(norm);
    for (int k = 0; k < d; ++k) {
      p[k] /= norm;
    }
    P.push_back(p);
  }
//  cout << "P = " << endl;
//  for (auto p: P) cout << p << endl;

  tick();
  auto Q = dominance_naive(P);
  /*
  cout << "dom(P) = " << endl;
  for (auto p: Q) cout << p << endl;
  */
  cout << tick() << endl;

  auto R = dominance(P);
  /*
  cout << "dom(P) = " << endl;
  for (auto p: R) cout << p << endl;
  */
  cout << tick() << endl;

  if (Q.size() != R.size()) exit(-1);
}

int main2() {
  /*
  int s = 14; 
  cin >> s;
  srand( s );
  */
  vector<point> P, Q;
  int n = 10000, d = 5;
  for (int i = 0; i < n; ++i) {
    point p;
    double norm = 0;
    for (int k = 0; k < d; ++k) {
      p.push_back( rand()%n );
      norm += p.back() * p.back();
    }
    norm = sqrt(norm);
    for (int k = 0; k < d; ++k) {
      p[k] /= norm;
    }
    P.push_back(p);
  }
  for (int i = 0; i < n; ++i) {
    point p;
    double norm = 0;
    for (int k = 0; k < d; ++k) {
      p.push_back( rand()%n );
      norm += p.back() * p.back();
    }
    norm = sqrt(norm);
    for (int k = 0; k < d; ++k) {
      p[k] *= (1.05 / norm);
    }
    Q.push_back(p);
  }
  //D.report();

  tick();
  int count = 0;
  auto X = filter_naive(P, Q);
  for (auto p: X) {
    //cout << p << endl;
    ++count;
  }
  cout << count << " " << tick() << endl;

  //cout << endl;

  auto Y = filter(P, Q);
  for (auto p: Y) {
    //cout << p << endl;
    --count;
  }
  cout << count << " " << tick() << endl;
  if (count != 0) exit(-1);
  // main();
}