Heap Sampler, including implementation, repro script, and log parser

Makefile

@@ -9,7 +9,8 @@ LDLIBS=-lbsd -lboost_system -lboost_context -lboost_chrono -lboost_thread -lpthr
 
 BIN=prac
 SRCS=prac.cpp mpcio.cpp preproc.cpp online.cpp mpcops.cpp rdpf.cpp \
-    cdpf.cpp duoram.cpp cell.cpp bst.cpp avl.cpp heap.cpp
+    cdpf.cpp duoram.cpp cell.cpp bst.cpp avl.cpp heap.cpp \
+    heapsampler.cpp
 OBJS=$(SRCS:.cpp=.o)
 ASMS=$(SRCS:.cpp=.s)
 
@@ -43,6 +44,7 @@ online.o: online.hpp mpcio.hpp types.hpp bitutils.hpp corotypes.hpp mpcio.tcc
 online.o: options.hpp mpcops.hpp coroutine.hpp mpcops.tcc rdpf.hpp dpf.hpp
 online.o: prg.hpp aes.hpp rdpf.tcc duoram.hpp duoram.tcc cdpf.hpp cdpf.tcc
 online.o: cell.hpp heap.hpp shapes.hpp shapes.tcc bst.hpp avl.hpp
+online.o: heapsampler.hpp
 mpcops.o: mpcops.hpp types.hpp bitutils.hpp mpcio.hpp corotypes.hpp mpcio.tcc
 mpcops.o: coroutine.hpp mpcops.tcc
 rdpf.o: rdpf.hpp mpcio.hpp types.hpp bitutils.hpp corotypes.hpp mpcio.tcc
@@ -65,3 +67,7 @@ heap.o: types.hpp bitutils.hpp duoram.hpp mpcio.hpp corotypes.hpp mpcio.tcc
 heap.o: coroutine.hpp rdpf.hpp dpf.hpp prg.hpp aes.hpp rdpf.tcc mpcops.hpp
 heap.o: mpcops.tcc duoram.tcc cdpf.hpp cdpf.tcc cell.hpp options.hpp
 heap.o: shapes.hpp shapes.tcc heap.hpp
+heapsampler.o: heapsampler.hpp mpcio.hpp types.hpp bitutils.hpp corotypes.hpp
+heapsampler.o: mpcio.tcc coroutine.hpp heap.hpp options.hpp mpcops.hpp
+heapsampler.o: mpcops.tcc duoram.hpp rdpf.hpp dpf.hpp prg.hpp aes.hpp
+heapsampler.o: rdpf.tcc duoram.tcc cdpf.hpp cdpf.tcc

heapsampler.cpp

@@ -0,0 +1,222 @@
+#include "heapsampler.hpp"
+
+// In each 64-bit RegAS in the heap, the top bit (of the reconstructed
+// value) is 0, the next 42 bits are the random tag, and the low 21 bits
+// are the element value.
+
+#define HEAPSAMPLE_TAG_BITS 42
+#define HEAPSAMPLE_ELT_BITS 21
+
+// Make the next random tag
+void HeapSampler::make_randtag(MPCTIO &tio, yield_t &yield)
+{
+    // Make a uniformly random HEAPSAMPLE_TAG_BITS-bit tag. This needs
+    // to be RegXS in order for the sum (XOR) of P0 and P1's independent
+    // values to be uniform.
+    RegXS tagx;
+    tagx.randomize(HEAPSAMPLE_TAG_BITS);
+    mpc_xs_to_as(tio, yield, randtag, tagx);
+}
+
+// Compute the heap size (the smallest power of two strictly greater
+// than k) needed to store k elements
+static size_t heapsize(size_t k)
+{
+    size_t ret = 1;
+    while (ret <= k) {
+        ret <<= 1;
+    }
+    return ret;
+}
+
+// Return a random bit that reconstructs to 1 with probability k/m
+static RegBS weighted_coin(MPCTIO &tio, yield_t &yield, size_t k,
+    size_t m)
+{
+    RegAS limit;
+    limit.ashare = size_t((__uint128_t(k)<<63)/m) * !tio.player();
+    RegXS randxs;
+    randxs.randomize(63);
+    RegAS randas;
+    mpc_xs_to_as(tio, yield, randas, randxs);
+    CDPF cdpf = tio.cdpf(yield);
+    auto[lt, eq, gt] = cdpf.compare(tio, yield, randas-limit, tio.aes_ops());
+
+    return lt;
+}
+
+// Constructor for a HeapSampler that samples k items from a stream
+// of abritrary and unknown size, using O(k) memory
+HeapSampler::HeapSampler(MPCTIO &tio, yield_t &yield, size_t k)
+    : k(k), m(0), heap(tio.player(), heapsize(k))
+{
+    run_coroutines(tio, [&tio, this](yield_t &yield) {
+        heap.init(tio, yield);
+    }, [&tio, this](yield_t &yield) {
+        make_randtag(tio, yield);
+    });
+}
+
+// An element has arrived
+void HeapSampler::ingest(MPCTIO &tio, yield_t &yield, RegAS elt)
+{
+    ++m;
+    RegAS tagged_elt = (randtag << HEAPSAMPLE_ELT_BITS) + elt;
+    RegAS elt_to_insert = tagged_elt;
+
+    if (m > k) {
+        RegAS extracted_elt;
+        RegBS selection_bit;
+        run_coroutines(tio, [&tio, this, &extracted_elt](yield_t &yield) {
+            extracted_elt = heap.extract_min(tio, yield);
+        }, [&tio, this, &selection_bit](yield_t &yield) {
+            selection_bit = weighted_coin(tio, yield, k, m);
+        });
+        mpc_select(tio, yield, elt_to_insert, selection_bit,
+            extracted_elt, tagged_elt);
+    }
+
+    run_coroutines(tio, [&tio, this, elt_to_insert](yield_t &yield) {
+        heap.insert_optimized(tio, yield, elt_to_insert);
+    }, [&tio, this](yield_t &yield) {
+        make_randtag(tio, yield);
+    });
+}
+
+// The stream has ended; output min(k,m) randomly sampled elements.
+// After calling this function, the HeapSampler is reset to its
+// initial m=0 state.
+std::vector<RegAS> HeapSampler::close(MPCTIO &tio, yield_t &yield)
+{
+    size_t retsize = k;
+    if (m < k) {
+        retsize = m;
+    }
+    std::vector<RegAS> ret(retsize);
+    for (size_t i=0; i<retsize; ++i) {
+        ret[i] = heap.extract_min(tio, yield);
+        ret[i] &= ((size_t(1)<<HEAPSAMPLE_ELT_BITS)-1);
+    }
+    // Compare each output to (size_t(1)<<HEAPSAMPLE_ELT_BITS), since
+    // there may be a carry; if the output is greater than or equal to
+    // that value, fix the carry
+    RegAS limit;
+    limit.ashare = (size_t(1)<<HEAPSAMPLE_ELT_BITS) * !tio.player();
+
+    std::vector<coro_t> coroutines;
+    for (size_t i=0; i<retsize; ++i) {
+        coroutines.emplace_back([&tio, &ret, i, limit](yield_t &yield) {
+            CDPF cdpf = tio.cdpf(yield);
+            auto[lt, eq, gt] = cdpf.compare(tio, yield,
+                ret[i]-limit, tio.aes_ops());
+            RegAS fix, zero;
+            mpc_select(tio, yield, fix, gt^eq, zero, limit);
+            ret[i] -= fix;
+        });
+    }
+    run_coroutines(tio, coroutines);
+    heap.init(tio, yield);
+    return ret;
+}
+
+void heapsampler_test(MPCIO &mpcio, const PRACOptions &opts, char **args)
+{
+    size_t n = 100;
+    size_t k = 10;
+
+    // The number of elements to stream
+    if (*args) {
+        n = atoi(*args);
+        ++args;
+    }
+    // The size of the random sample
+    if (*args) {
+        k = atoi(*args);
+        ++args;
+    }
+
+    MPCTIO tio(mpcio, 0, opts.num_threads);
+    run_coroutines(tio, [&mpcio, &tio, n, k] (yield_t &yield) {
+
+        std::cout << "\n===== STREAMING =====\n";
+
+        HeapSampler sampler(tio, yield, k);
+
+        for (size_t i=0; i<n; ++i) {
+            // For ease of checking, just have the elements be in a
+            // simple sequence
+            RegAS elt;
+            elt.ashare = (i+1) * (1 + 0xfff*tio.player());
+            sampler.ingest(tio, yield, elt);
+        }
+
+        std::vector<RegAS> sample = sampler.close(tio, yield);
+
+        tio.sync_lamport();
+        mpcio.dump_stats(std::cout);
+        mpcio.reset_stats();
+        tio.reset_lamport();
+
+        std::cout << "\n===== CHECKING =====\n";
+        size_t expected_size = k;
+        if (n < k) {
+            expected_size = n;
+        }
+        assert(sample.size() == expected_size);
+        std::vector<value_t> reconstructed_sample(expected_size);
+
+        std::vector<coro_t> coroutines;
+        for (size_t i=0; i<expected_size; ++i) {
+            coroutines.emplace_back(
+                [&tio, &sample, i, &reconstructed_sample](yield_t &yield) {
+                    reconstructed_sample[i] = mpc_reconstruct(
+                        tio, yield, sample[i]);
+                });
+        }
+        run_coroutines(tio, coroutines);
+        if (tio.player() == 0) {
+            for (size_t i=0; i<expected_size; ++i) {
+                printf("%06lx\n", reconstructed_sample[i]);
+            }
+        }
+    });
+}
+
+void weighted_coin_test(MPCIO &mpcio, const PRACOptions &opts, char **args)
+{
+    size_t iters = 100;
+    size_t m = 100;
+    size_t k = 10;
+
+    // The number of iterations
+    if (*args) {
+        iters = atoi(*args);
+        ++args;
+    }
+    // The denominator
+    if (*args) {
+        m = atoi(*args);
+        ++args;
+    }
+    // The numerator
+    if (*args) {
+        k = atoi(*args);
+        ++args;
+    }
+
+    MPCTIO tio(mpcio, 0, opts.num_threads);
+    run_coroutines(tio, [&mpcio, &tio, iters, m, k] (yield_t &yield) {
+
+        size_t heads = 0, tails = 0;
+        for (size_t i=0; i<iters; ++i) {
+            RegBS coin = weighted_coin(tio, yield, k, m);
+            bool coin_rec = mpc_reconstruct(tio, yield, coin);
+            if (coin_rec) {
+                heads++;
+            } else {
+                tails++;
+            }
+            printf("%lu flips %lu heads %lu tails\n", i+1, heads, tails);
+        }
+    });
+}

heapsampler.hpp

@@ -0,0 +1,65 @@
+#ifndef __HEAPSAMPLER_HPP__
+#define __HEAPSAMPLER_HPP__
+
+#include <vector>
+#include "mpcio.hpp"
+#include "coroutine.hpp"
+#include "heap.hpp"
+
+// Implement a stream sampler using a MinHeap.  A stream sampler will
+// sample a random subset of k elements from an arbitrality long stream
+// of elements, while using only O(k) memory.  The number of elements in
+// the stream need not be known in advance.  Importantly, no party knows
+// which elements ended up in the sample.
+
+// We use the technique from "Path Oblivious Heap" by Elaine Shi
+// (IEEE Symsposium on Security and Privacy 2020):
+// https://eprint.iacr.org/2019/274.pdf; in particular, Section 5.2
+// "Oblivious Streaming Sampler with Applications to Distributed
+// Differential Privacy".  We correct the typo that the probability to
+// keep the mth item (for m>k) is listed as 1/m, but it should be k/m.
+// Note that the Shi paper is in the client-server setting, where the
+// client _is_ allowed to know which elements did and did not end up in
+// the sample (but the server, who stores the heap, is not). In our MPC
+// setting, no party may learn which elements ended up in the sample.
+
+class HeapSampler {
+    // The number of items to sample
+    size_t k;
+
+    // The number of items that have arrived so far
+    size_t m;
+
+    // The MinHeap with O(k) storage that when m>=k stores a
+    // uniformly random sample of size k of the m items seen so far
+    MinHeap heap;
+
+    // The next random tag to use
+    RegAS randtag;
+
+    // Make the next random tag
+    void make_randtag(MPCTIO &tio, yield_t &yield);
+
+public:
+    // Constructor for a HeapSampler that samples k items from a stream
+    // of abritrary and unknown size, using O(k) memory
+    HeapSampler(MPCTIO &tio, yield_t &yield, size_t k);
+
+    // An element has arrived
+    void ingest(MPCTIO &tio, yield_t &yield, RegAS elt);
+
+    // The stream has ended; output min(k,m) randomly sampled elements.
+    // After calling this function, the HeapSampler is reset to its
+    // initial m=0 state.
+    std::vector<RegAS> close(MPCTIO &tio, yield_t &yield);
+};
+
+// A unit test for the HeapSampler
+
+void heapsampler_test(MPCIO &mpcio, const PRACOptions &opts, char **args);
+
+// A unit test for the weighted_coin internal function
+
+void weighted_coin_test(MPCIO &mpcio, const PRACOptions &opts, char **args);
+
+#endif

online.cpp

@@ -10,6 +10,7 @@
 #include "shapes.hpp"
 #include "bst.hpp"
 #include "avl.hpp"
+#include "heapsampler.hpp"
 
 static void online_test(MPCIO &mpcio,
     const PRACOptions &opts, char **args)
@@ -1691,6 +1692,12 @@ void online_main(MPCIO &mpcio, const PRACOptions &opts, char **args)
     } else if (!strcmp(*args, "heap")) {
         ++args;
         Heap(mpcio, opts, args);
+    } else if (!strcmp(*args, "heapsampler")) {
+        ++args;
+        heapsampler_test(mpcio, opts, args);
+    } else if (!strcmp(*args, "weightedcoin")) {
+        ++args;
+        weighted_coin_test(mpcio, opts, args);
     } else {
         std::cerr << "Unknown mode " << *args << "\n";
     }

repro/parse_logs

@@ -46,6 +46,8 @@ while(<>) {
             &parse_read($cmdline);
         } elsif ($cmdline =~ /^-- b?bsearch/) {
             &parse_bsearch($cmdline);
+        } elsif ($cmdline =~ /^-- heapsampler/) {
+            &parse_heapsampler($cmdline);
         } elsif ($cmdline =~ /^-- heap/) {
             &parse_heap($cmdline);
         } elsif ($cmdline =~ /^-- avl/) {
@@ -197,7 +199,9 @@ sub parse_read {
             next;
         }
         last if /===== End/;
-        # The log was malformed
+        # Try to recover from a malformed log
+        last if /^Max MB:/;
+        # It was too malformed
         die "Malformed log" if /===== Running/;
         if (/^(\d+) message bytes sent/) {
             $online_kib[$who] = $1 / 1024;
@@ -261,7 +265,9 @@ sub parse_bsearch {
             next;
         }
         last if /===== End/;
-        # The log was malformed
+        # Try to recover from a malformed log
+        last if /^Max MB:/;
+        # It was too malformed
         die "Malformed log" if /===== Running/;
         if ($section eq "BINARY SEARCH") {
             if (/^(\d+) message bytes sent/) {
@@ -335,7 +341,9 @@ sub parse_heap {
             next;
         }
         last if /===== End/;
-        # The log was malformed
+        # Try to recover from a malformed log
+        last if /^Max MB:/;
+        # It was too malformed
         die "Malformed log" if /===== Running/;
         my $rightsection = 0;
         if ($section eq "Insert" && $oper eq "Ins") {
@@ -417,7 +425,9 @@ sub parse_avl {
             next;
         }
         last if /===== End/;
-        # The log was malformed
+        # Try to recover from a malformed log
+        last if /^Max MB:/;
+        # It was too malformed
         die "Malformed log" if /===== Running/;
         my $rightsection = 0;
         if ($section eq "INSERTS" && $oper eq "Ins") {
@@ -464,6 +474,71 @@ sub parse_avl {
     &accum_data(\%online_P2mem_mib_data, $label, $online_mem_mib[2]);
 }
 
+sub parse_heapsampler {
+    my $cmdline = $_[0];
+    my $who = 0;
+    my $section = '';
+    my @online_seconds = (0, 0, 0);
+    my @online_kib = (0, 0, 0);
+    my @online_latencies = (0, 0, 0);
+    my @online_mem_mib = (0, 0, 0);
+
+    unless ($cmdline =~ /heapsampler (\d+) (\d+)/) {
+        die "Cannot parse heapsampler cmdline: $cmdline";
+    }
+    my ($num, $size) = ($1, $2);
+    while(<>) {
+        if (/===== P([012]) output/) {
+            $who = $1;
+            next;
+        }
+        if (/===== ([A-Z ]+) =====/) {
+            $section = $1;
+            next;
+        }
+        last if /===== End/;
+        # Try to recover from a malformed log
+        last if /^Max MB:/;
+        # It was too malformed
+        die "Malformed log" if /===== Running/;
+        if ($section eq "STREAMING") {
+            if (/^(\d+) message bytes sent/) {
+                $online_kib[$who] = $1 / 1024;
+            } elsif (/^(\d+) Lamport clock/) {
+                $online_latencies[$who] = $1;
+            } elsif (/^(\d+) milliseconds wall clock/) {
+                $online_seconds[$who] = $1 / 1000;
+            } elsif (/^Mem: (\d+) KiB/) {
+                $online_mem_mib[$who] = $1 / 1024;
+            }
+        }
+        if ($who == 0 && /^Precomputed values used: (.*)/) {
+            my %used_resources = ();
+            &parse_resources(\%used_resources, $1);
+            my $preproc_resources_str = &serialize_resources(\%preproc_resources);
+            my $used_resources_str = &serialize_resources(\%used_resources);
+            if ($preproc_resources_str ne $used_resources_str) {
+                warn "Resource usage does not match preprocessing:\n" .
+                    "Preproc: $preproc_resources_str\n" .
+                    "Used: $used_resources_str\n ";
+            }
+        }
+    }
+    my $label = "PRAC heapsampler $netsetup$size $num";
+    &accum_data(\%preproc_s_data, $label, &maxarray(@preproc_seconds));
+    &accum_data(\%preproc_kib_data, $label, &avgarray(@preproc_kib));
+    &accum_data(\%preproc_latencies_data, $label, &maxarray(@preproc_latencies));
+    &accum_data(\%preproc_P0mem_mib_data, $label, $preproc_mem_mib[0]);
+    &accum_data(\%preproc_P1mem_mib_data, $label, $preproc_mem_mib[1]);
+    &accum_data(\%preproc_P2mem_mib_data, $label, $preproc_mem_mib[2]);
+    &accum_data(\%online_s_data, $label, &maxarray(@online_seconds));
+    &accum_data(\%online_kib_data, $label, &avgarray(@online_kib));
+    &accum_data(\%online_latencies_data, $label, &maxarray(@online_latencies));
+    &accum_data(\%online_P0mem_mib_data, $label, $online_mem_mib[0]);
+    &accum_data(\%online_P1mem_mib_data, $label, $online_mem_mib[1]);
+    &accum_data(\%online_P2mem_mib_data, $label, $online_mem_mib[2]);
+}
+
 sub maxarray {
     my $max = $_[0];
     foreach (@_) {

repro/repro

@@ -280,6 +280,54 @@ if [ "$whichexps" = "tab4" -o "$whichexps" = "all" ]; then
     done
 fi
 
+if [ "$whichexps" = "heapsampler" -o "$whichexps" = "all" ]; then
+    echo "Running heap sampler experiments..."
+    for iter in $(seq 1 $numiters); do
+        # Table 4
+        logname='heapsampler'
+        run -p m:1033 h:7497 a:3 s:118 r1:1 r2:6 r3:57 i1:6 i2:57 i1.3:2 i2.3:4 i3.3:57 c:610
+        run heapsampler 64 10
+        run -p m:2121 h:15561 a:6 s:246 r1:1 r2:6 r3:121 i1:6 i2:121 i1.3:2 i2.3:4 i3.3:121 c:1250
+        run heapsampler 128 10
+        run -p m:4297 h:31689 a:12 s:502 r1:1 r2:6 r3:249 i1:6 i2:249 i1.3:2 i2.3:4 i3.3:249 c:2530
+        run heapsampler 256 10
+        run -p m:8649 h:63945 a:24 s:1014 r1:1 r2:6 r3:505 i1:6 i2:505 i1.3:2 i2.3:4 i3.3:505 c:5090
+        run heapsampler 512 10
+        run -p m:17353 h:128457 a:48 s:2038 r1:1 r2:6 r3:1017 i1:6 i2:1017 i1.3:2 i2.3:4 i3.3:1017 c:10210
+        run heapsampler 1024 10
+        run -p m:34761 h:257481 a:96 s:4086 r1:1 r2:6 r3:2041 i1:6 i2:2041 i1.3:2 i2.3:4 i3.3:2041 c:20450
+        run heapsampler 2048 10
+        run -p m:1278 h:6237 a:4 s:167 r1:1 r2:6 r3:57 i1:6 i2:57 i1.3:2 i2.3:4 i3.3:8 i4.3:49 c:708
+        run heapsampler 64 30
+        run -p m:2686 h:14301 a:8 s:359 r1:1 r2:6 r3:121 i1:6 i2:121 i1.3:2 i2.3:4 i3.3:8 i4.3:113 c:1476
+        run heapsampler 128 30
+        run -p m:5502 h:30429 a:16 s:743 r1:1 r2:6 r3:249 i1:6 i2:249 i1.3:2 i2.3:4 i3.3:8 i4.3:241 c:3012
+        run heapsampler 256 30
+        run -p m:11134 h:62685 a:32 s:1511 r1:1 r2:6 r3:505 i1:6 i2:505 i1.3:2 i2.3:4 i3.3:8 i4.3:497 c:6084
+        run heapsampler 512 30
+        run -p m:22398 h:127197 a:64 s:3047 r1:1 r2:6 r3:1017 i1:6 i2:1017 i1.3:2 i2.3:4 i3.3:8 i4.3:1009 c:12228
+        run heapsampler 1024 30
+        run -p m:44926 h:256221 a:128 s:6119 r1:1 r2:6 r3:2041 i1:6 i2:2041 i1.3:2 i2.3:4 i3.3:8 i4.3:2033 c:24516
+        run heapsampler 2048 30
+        run -p m:3496 h:9891 a:11 s:521 r1:1 r2:6 r3:121 i1:6 i2:121 i1.3:2 i2.3:4 i3.3:8 i4.3:16 i5.3:32 i6.3:65 c:1800
+        run heapsampler 128 100
+        run -p m:7592 h:26019 a:23 s:1161 r1:1 r2:6 r3:249 i1:6 i2:249 i1.3:2 i2.3:4 i3.3:8 i4.3:16 i5.3:32 i6.3:193 c:3848
+        run heapsampler 256 100
+        run -p m:15784 h:58275 a:47 s:2441 r1:1 r2:6 r3:505 i1:6 i2:505 i1.3:2 i2.3:4 i3.3:8 i4.3:16 i5.3:32 i6.3:449 c:7944
+        run heapsampler 512 100
+        run -p m:32168 h:122787 a:95 s:5001 r1:1 r2:6 r3:1017 i1:6 i2:1017 i1.3:2 i2.3:4 i3.3:8 i4.3:16 i5.3:32 i6.3:961 c:16136
+        run heapsampler 1024 100
+        run -p m:64936 h:251811 a:191 s:10121 r1:1 r2:6 r3:2041 i1:6 i2:2041 i1.3:2 i2.3:4 i3.3:8 i4.3:16 i5.3:32 i6.3:1985 c:32520
+        run heapsampler 2048 100
+        run -p m:18994 h:45675 a:57 s:3083 r1:1 r2:6 r3:120 r4:385 i1:6 i2:120 i3:385 i1.3:2 i2.3:4 i3.3:8 i4.3:16 i5.3:32 i6.3:64 i7.3:128 i8.3:257 c:9613
+        run heapsampler 512 300
+        run -p m:40498 h:110187 a:121 s:6667 r1:1 r2:6 r3:120 r4:897 i1:6 i2:120 i3:897 i1.3:2 i2.3:4 i3.3:8 i4.3:16 i5.3:32 i6.3:64 i7.3:128 i8.3:769 c:20365
+        run heapsampler 1024 300
+        run -p m:83506 h:239211 a:249 s:13835 r1:1 r2:6 r3:120 r4:1921 i1:6 i2:120 i3:1921 i1.3:2 i2.3:4 i3.3:8 i4.3:16 i5.3:32 i6.3:64 i7.3:128 i8.3:1793 c:41869
+        run heapsampler 2048 300
+    done
+fi
+
 now=`date`
 echo "$now: Experiments complete"
 
@@ -409,6 +457,9 @@ egrep 'OptPRACPreprc avlDel [0-9]+ 1 .* KiB$' data/prac.dat | sort -k3 -n
 echo
 egrep 'OptPRACOnln avlDel [0-9]+ 1 .* KiB$' data/prac.dat | sort -k3 -n
 echo
+echo '# Heap Sampler'
+egrep 'PRACTotl heapsampler [0-9]+ [0-9]+ .* s' data/prac.dat | sort -k3,3n -k4,4n
+echo
 echo "# End figures"
 
 fi