sampling algorithms

Author: GaspareG

Chap.03/README.md

@@ -1,8 +1,9 @@
 # 3 Random Sampling
 
-- Drawing from all un-sampled positions
-- Dictionary of sampled positions
-- Sorting sampling
-- Scanning an selecting
-- Heap and random keys
-- Reservoir sampling
+- Drawing from all un-sampled positions ([drawing_sampling.hpp](drawing_sampling.hpp))
+- Dictionary of sampled positions ([dictionary_sampling.hpp](dictionary_sampling.hpp))
+- Sorting sampling ([sorting_sampling.hpp](sorting_sampling.hpp))
+- Scanning an selecting ([scanning_sampling.hpp](scanning_sampling.hpp))
+- Heap and random keys ([heap_sampling.hpp](heap_sampling.hpp))
+- Reservoir sampling ([reservoir_sampling.hpp](reservoir_sampling.hpp))
+

Chap.03/dictionary_sampling.hpp

@@ -0,0 +1,19 @@
+/*
+  License: GPL-3.0
+  Author: Gaspare Ferraro <[email protected]>
+*/
+
+#pragma once
+
+#include <vector>
+#include <set>
+#include <utility>
+#include <iterator>
+
+template <class T, class RandomNumberGenerator>
+std::vector<T> dictionary_sampling(const std::vector<T>& S, size_t m, RandomNumberGenerator& gen) {
+  std::set<T> D;
+  while( D.size() < m )
+    D.insert(S[gen() % S.size()]);
+  return std::vector<T>(D.begin(), D.end());
+}

Chap.03/drawing_sampling.hpp

@@ -9,7 +9,14 @@
 #include <utility>
 #include <iterator>
 
-template <class T>
-std::vector<T> drawing_sampling(const std::vector<T> S, int m) {
-  
+template <class T, class RandomNumberGenerator>
+std::vector<T> drawing_sampling(const std::vector<T>& S, size_t m, RandomNumberGenerator& gen){
+  std::vector<T> tmp(S.begin(), S.end());
+  std::vector<T> out(m, T());
+  for(size_t i=0; i<m; i++) {
+    int p = gen() % (tmp.size()-i);
+    out[i] = tmp[p];
+    std::swap(tmp[p], tmp[tmp.size()-i-1]);
+  }
+  return out;
 }

Chap.03/heap_sampling.hpp

@@ -0,0 +1,31 @@
+/*
+  License: GPL-3.0
+  Author: Gaspare Ferraro <[email protected]>
+*/
+
+#pragma once
+
+#include <vector>
+#include <utility>
+#include <iterator>
+#include <queue>
+#include <climits>
+
+template <class T, class RandomNumberGenerator>
+std::vector<T> heap_sampling(const std::vector<T>& S, size_t m, RandomNumberGenerator& gen) {
+  std::priority_queue< std::pair<double, T> > Q;
+  std::vector<T> out;
+  for(size_t i=0; i<m; i++) Q.push(std::pair<double, T>(1., S[0]));
+  for(T s : S) {
+    double r = (double) (gen() % INT_MAX)/INT_MAX;
+    if( r < Q.top().first ) {
+      Q.pop();
+      Q.push(std::make_pair(r, s));
+    } 
+  }
+  while( !Q.empty() ) {
+    out.push_back(Q.top().second);
+    Q.pop();
+  }
+  return out;
+}

Chap.03/reservoir_sampling.hpp

@@ -0,0 +1,23 @@
+/*
+  License: GPL-3.0
+  Author: Gaspare Ferraro <[email protected]>
+*/
+
+#pragma once
+
+#include <vector>
+#include <utility>
+#include <iterator>
+
+template <class T, class RandomNumberGenerator>
+std::vector<T> reservoir_sampling(const std::vector<T>& S, size_t m, RandomNumberGenerator& gen) {
+  std::vector<T> out(S.begin(), S.begin()+m);
+  auto start = S.begin()+m;
+  while( start != S.end() ) {
+    size_t h = gen() % (start-S.begin());
+    if( h < m )
+      out[h] = *start;
+    start++;
+  }
+  return out;
+}

Chap.03/scanning_sampling.hpp

@@ -0,0 +1,22 @@
+/*
+  License: GPL-3.0
+  Author: Gaspare Ferraro <[email protected]>
+*/
+
+#pragma once
+
+#include <vector>
+#include <utility>
+#include <iterator>
+
+template <class T, class RandomNumberGenerator>
+std::vector<T> scanning_sampling(const std::vector<T>& S, size_t m, RandomNumberGenerator& gen) {
+  std::vector<T> out(m, 0);
+  size_t s=0;
+  for(size_t i=0; i<S.size(); ++i) {
+    size_t p = gen() % (S.size()-i+1) ;
+    if( p < m-s )
+      out[s++] = S[i];
+  }
+  return out;
+}

Chap.03/sorting_sampling.hpp

@@ -0,0 +1,22 @@
+/*
+  License: GPL-3.0
+  Author: Gaspare Ferraro <[email protected]>
+*/
+
+#pragma once
+
+#include <vector>
+#include <set>
+#include <utility>
+#include <iterator>
+
+template <class T, class RandomNumberGenerator>
+std::vector<T> sorting_sampling(const std::vector<T>& S, size_t m, RandomNumberGenerator& gen) {
+  std::set<T> D;
+  while( D.size() < m ){
+    D.clear();
+    for(size_t i=0; i<m; ++i)
+    	D.insert(S[gen() % S.size()]);
+  }
+  return std::vector<T>(D.begin(), D.end());
+}

Chap.03/test.cpp

@@ -0,0 +1,47 @@
+#include <bits/stdc++.h>
+#include "dictionary_sampling.hpp"
+#include "drawing_sampling.hpp"
+#include "heap_sampling.hpp"
+#include "reservoir_sampling.hpp"
+#include "scanning_sampling.hpp"
+#include "sorting_sampling.hpp"
+
+#define NTEST 100000
+
+using namespace std;
+
+int main(){
+
+  random_device rd;
+  mt19937 rng(rd());
+
+  vector<int> S;
+  map<int, int> dictH;
+  map<int, int> drawH;
+  map<int, int> heapH;
+  map<int, int> reservoirH;
+  map<int, int> scanningH;
+  map<int, int> sortingH;
+
+  for(int i=0; i<100; i++) S.push_back(i);
+
+  for(int i=0; i<NTEST; i++) {
+    auto dict = dictionary_sampling<int>(S, 10, rng);
+    auto draw = drawing_sampling<int>(S, 10, rng);
+    auto heap = heap_sampling<int>(S, 10, rng);
+    auto reservoir = reservoir_sampling<int>(S, 10, rng);
+    auto scanning = scanning_sampling<int>(S, 10, rng);
+    auto sorting = sorting_sampling<int>(S, 10, rng);
+
+    for(int x: dict) dictH[x]++;
+    for(int x: draw) drawH[x]++;
+    for(int x: heap) heapH[x]++; 
+    for(int x: reservoir) reservoirH[x]++;
+    for(int x: scanning) scanningH[x]++;
+    for(int x: sorting) sortingH[x]++;
+  }
+
+  for(int i=0; i<100; i++)
+    printf("[%3d] %6d %6d %6d %6d %6d %6d\n", i, dictH[i],  drawH[i],  heapH[i],  reservoirH[i],  scanningH[i],  sortingH[i]);
+
+}

Chap.06/benchmark

@@ -0,0 +1,60 @@
+/*
+  License: GPL-3.0
+  Author: Gaspare Ferraro <[email protected]>
+*/
+
+#include <bits/stdc++.h>
+#include "binary_search_intersection.hpp"
+#include "merge_intersection.hpp"
+#include "doubling_search.hpp"
+#include "mutual_partitioning.hpp"
+
+using namespace std;
+
+vector<int> randomSet(int size, int max, int min=1) {
+
+  if( max < min ) return randomSet(size, min, max);
+  assert( max-min+1 > size );
+
+  random_device rd;
+  mt19937 rng(rd());
+
+  vector<int> out(size, 0);
+  for(int i=0; i<size; i++)
+    out[i] = min+i;
+  for(int i=size; i<max; i++)
+  {
+    int j = rng() % i;
+    if( j < size )
+      out[j] = i;
+  }
+  sort(out.begin(), out.end());
+  return out;    
+}
+
+int main()
+{
+
+  vector<int> size = {100, 1000, 10000, 100000, 1000000, 10000000};
+
+  vector<int> a = randomSet(size[4], 3*size[5]);
+  vector<int> b = randomSet(size[5], 3*size[5]);
+
+	for(int k=0; k<100; k++){
+		clock_t doubling = 0;
+		clock_t start = 0;
+
+		vector<int> c(b.size());
+		c.clear();
+		start = clock();
+		doubling_search_intersection<int>(a, b, c.begin());
+		doubling += clock()-start;
+		double doublingd = (double) doubling / CLOCKS_PER_SEC;
+
+		printf("| %8d | %8d | %.3fs |\n", size[4], size[5], doublingd);
+		fflush(stdout);  
+	}	
+
+  return 0;
+}
+

Chap.06/benchmark_doubling.cpp

@@ -24,41 +24,44 @@ int main()
 
   vector<int> a(sa.begin(), sa.end());
   vector<int> b(sb.begin(), sb.end());
-  vector<int> c(20);
+  vector<int> c(20), d(20), e(20), f(20), g(20);
 
   auto it = set_intersection(a.begin(), a.end(), b.begin(), b.end(), c.begin());
   c.resize(it-c.begin());
-
-  auto d = binary_search_intersection<int>(a, b);
-  auto e = merge_intersection<int>(a, b);
-  auto f = doubling_search_intersection<int>(a, b);
-  auto g = mutual_partitioning_intersection<int>(a, b);
-
   printf("stl set_intersection:\n\t");
   for(auto x : c)
     printf("%d ", x);
   printf("\n");
 
+  it = binary_search_intersection<int>(a, b, d.begin());
+  d.resize(it-d.begin());
   printf("binary_search_intersection:\n\t");
   for(auto x : d)
     printf("%d ", x);
   printf("\n");
 
+  it = merge_intersection<int>(a, b, e.begin());
+  e.resize(it-e.begin());
   printf("merge_intersection:\n\t");
   for(auto x : e)
     printf("%d ", x);
   printf("\n");
 
+  it = doubling_search_intersection<int>(a, b, f.begin());
+  f.resize(it-f.begin());
   printf("doubling_search_intersection:\n\t");
   for(auto x : f)
     printf("%d ", x);
   printf("\n");
 
+  it = mutual_partitioning_intersection<int>(a, b, g.begin());
+  g.resize(it-g.begin());
   printf("mutual_partitioning:\n\t");
   for(auto x : g)
     printf("%d ", x);
   printf("\n");
 
+
   return 0;
 }
 

Chap.06/test.cpp

@@ -9,12 +9,12 @@ Chapters:
 [3 Random Sampling](Chap.03/)
 ------------------
 
-- Drawing from all un-sampled positions
-- Dictionary of sampled positions
-- Sorting sampling
-- Scanning an selecting
-- Heap and random keys
-- Reservoir sampling
+- Drawing from all un-sampled positions ([drawing_sampling.hpp](Chap.03/drawing_sampling.hpp))
+- Dictionary of sampled positions ([dictionary_sampling.hpp](Chap.03/dictionary_sampling.hpp))
+- Sorting sampling ([sorting_sampling.hpp](Chap.03/sorting_sampling.hpp))
+- Scanning an selecting ([scanning_sampling.hpp](Chap.03/scanning_sampling.hpp))
+- Heap and random keys ([heap_sampling.hpp](Chap.03/heap_sampling.hpp))
+- Reservoir sampling ([reservoir_sampling.hpp](Chap.03/reservoir_sampling.hpp))
 
 [4 List Ranking](Chap.04/)
 ------------------