Changeset 3243

Timestamp:

May 24, 2014, 6:26:32 PM (9 years ago)

Author:

Peter

Message:

speed up Percentiler when input range is unweighted and unsorted. Do not sort the entire range but use nth_element to find the element of interest.

File:

• trunk/yat/statistics/Percentiler.h (modified) (2 diffs)

trunk/yat/statistics/Percentiler.h

@@ -102 +102 @@
       if (first==last)
         return std::numeric_limits<double>::quiet_NaN();
-      return calculate(first, last, sorted_,
-       typename utility::weighted_iterator_traits<RandomAccessIterator>::type());
+      // just to avoid long line
+      using utility::weighted_iterator_traits;
+      typedef typename weighted_iterator_traits<RandomAccessIterator>::type tag;
+      return calculate(first, last, sorted_, tag());
     }
   private:
@@ -136 +138 @@
                          utility::unweighted_iterator_tag tag) const
   {
+    size_t n = last - first;
     // range is one value only is a special case
-    if (first+1 == last)
+    if (n == 1)
       *first;
+
+    double j = n * perc_ / 100.0;
+
+    // Some special cases
+    if (j > n-1) {
+      if (sorted)
+        return *(--last);
+      return *std::max_element(first, last);
+    }
+    if (j < 1.0) {
+      if (sorted)
+        return *first;
+      return *std::min_element(first, last);
+    }
+
+    size_t i = static_cast<size_t>(j);
+    // for border case (j integer) we take the average of adjacent bins
+    size_t k = (i==j) ? i-1 : i;
+
     if (sorted) {
-      size_t n = last - first;
-      // have to take care of this special case
-      if (perc_>= 100.0)
-        return *(--last);
-      if (perc_<= 0.0)
-        return *first;
-      double j = n * perc_ / 100.0;
-      size_t i = static_cast<size_t>(j);
-      if (i==j)
-        return (first[i]+first[i-1])/2;
-      return first[i];
-    }
-
-    std::vector<double> v_copy;
-    v_copy.reserve(std::distance(first,last));
-    std::copy(first, last, std::back_inserter(v_copy));
-    std::sort(v_copy.begin(), v_copy.end());
-    return calculate(v_copy.begin(), v_copy.end(), true, tag);
+      if (i==k)
+        return first[i];
+      return (first[i]+first[k])/2;
+    }
+
+    std::vector<double> vec(first, last);
+    // find the ith element
+    std::nth_element(vec.begin(), vec.begin()+i, vec.end());
+    // if i==k simply return vec[i]
+    if (i==k)
+      return vec[i];
+    YAT_ASSERT(k==i-1);
+    // otherwise return average of vec[i] and vec[k].
+
+    // We need to find the kth element. Since we have called
+    // nth_element above, we know that all elements in (begin+i, end)
+    // are >= vec[i] and all elements in [begin, begin+i) are <=
+    // vec[i]. Consequently, kth (k=i-1) element is the max element in
+    // range [begin, begin+i).
+    return (vec[i] + *std::max_element(vec.begin(), vec.begin()+i))/2;
   }