Electroneum
rolling_median.h
Go to the documentation of this file.
1 // Copyright (c) 2019, The Monero Project
2 //
3 // All rights reserved.
4 //
5 // Redistribution and use in source and binary forms, with or without modification, are
6 // permitted provided that the following conditions are met:
7 //
8 // 1. Redistributions of source code must retain the above copyright notice, this list of
9 // conditions and the following disclaimer.
10 //
11 // 2. Redistributions in binary form must reproduce the above copyright notice, this list
12 // of conditions and the following disclaimer in the documentation and/or other
13 // materials provided with the distribution.
14 //
15 // 3. Neither the name of the copyright holder nor the names of its contributors may be
16 // used to endorse or promote products derived from this software without specific
17 // prior written permission.
18 //
19 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
20 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
21 // MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL
22 // THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
23 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
24 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
25 // INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
26 // STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF
27 // THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
28 //
29 // Adapted from source by AShelly:
30 // Copyright (c) 2011 ashelly.myopenid.com, licenced under the MIT licence
31 // https://stackoverflow.com/questions/5527437/rolling-median-in-c-turlach-implementation
32 // https://stackoverflow.com/questions/1309263/rolling-median-algorithm-in-c
33 // https://ideone.com/XPbl6
34 
35 #pragma once
36 
37 #include <stdlib.h>
38 #include <stdint.h>
39 
40 namespace epee
41 {
42 namespace misc_utils
43 {
44 
45 template<typename Item>
47 {
48 private:
49  Item* data; //circular queue of values
50  int* pos; //index into `heap` for each value
51  int* heap; //max/median/min heap holding indexes into `data`.
52  int N; //allocated size.
53  int idx; //position in circular queue
54  int minCt; //count of items in min heap
55  int maxCt; //count of items in max heap
56  int sz; //count of items in heap
57 
58 private:
59 
60  //returns true if heap[i] < heap[j]
61  bool mmless(int i, int j) const
62  {
63  return data[heap[i]] < data[heap[j]];
64  }
65 
66  //swaps items i&j in heap, maintains indexes
67  bool mmexchange(int i, int j)
68  {
69  const int t = heap[i];
70  heap[i] = heap[j];
71  heap[j] = t;
72  pos[heap[i]] = i;
73  pos[heap[j]] = j;
74  return 1;
75  }
76 
77  //swaps items i&j if i<j; returns true if swapped
78  bool mmCmpExch(int i, int j)
79  {
80  return mmless(i, j) && mmexchange(i, j);
81  }
82 
83  //maintains minheap property for all items below i.
84  void minSortDown(int i)
85  {
86  for (i *= 2; i <= minCt; i *= 2)
87  {
88  if (i < minCt && mmless(i + 1, i))
89  ++i;
90  if (!mmCmpExch(i, i / 2))
91  break;
92  }
93  }
94 
95  //maintains maxheap property for all items below i. (negative indexes)
96  void maxSortDown(int i)
97  {
98  for (i *= 2; i >= -maxCt; i *= 2)
99  {
100  if (i > -maxCt && mmless(i, i - 1))
101  --i;
102  if (!mmCmpExch(i / 2, i))
103  break;
104  }
105  }
106 
107  //maintains minheap property for all items above i, including median
108  //returns true if median changed
109  bool minSortUp(int i)
110  {
111  while (i > 0 && mmCmpExch(i, i / 2))
112  i /= 2;
113  return i == 0;
114  }
115 
116  //maintains maxheap property for all items above i, including median
117  //returns true if median changed
118  bool maxSortUp(int i)
119  {
120  while (i < 0 && mmCmpExch(i / 2, i))
121  i /= 2;
122  return i == 0;
123  }
124 
125 protected:
126  rolling_median_t &operator=(const rolling_median_t&) = delete;
127  rolling_median_t(const rolling_median_t&) = delete;
128 
129 public:
130  //creates new rolling_median_t: to calculate `nItems` running median.
131  rolling_median_t(size_t N): N(N)
132  {
133  int size = N * (sizeof(Item) + sizeof(int) * 2);
134  data = (Item*)malloc(size);
135  pos = (int*) (data + N);
136  heap = pos + N + (N / 2); //points to middle of storage.
137  clear();
138  }
139 
141  {
142  free(data);
143  memcpy(this, &m, sizeof(rolling_median_t));
144  m.data = NULL;
145  }
147  {
148  free(data);
149  memcpy(this, &m, sizeof(rolling_median_t));
150  m.data = NULL;
151  return *this;
152  }
153 
155  {
156  free(data);
157  }
158 
159  void clear()
160  {
161  idx = 0;
162  minCt = 0;
163  maxCt = 0;
164  sz = 0;
165  int nItems = N;
166  while (nItems--) //set up initial heap fill pattern: median,max,min,max,...
167  {
168  pos[nItems] = ((nItems + 1) / 2) * ((nItems & 1) ? -1 : 1);
169  heap[pos[nItems]] = nItems;
170  }
171  }
172 
173  int size() const
174  {
175  return sz;
176  }
177 
178  //Inserts item, maintains median in O(lg nItems)
179  void insert(Item v)
180  {
181  int p = pos[idx];
182  Item old = data[idx];
183  data[idx] = v;
184  idx = (idx + 1) % N;
185  sz = std::min<int>(sz + 1, N);
186  if (p > 0) //new item is in minHeap
187  {
188  if (minCt < (N - 1) / 2)
189  {
190  ++minCt;
191  }
192  else if (v > old)
193  {
194  minSortDown(p);
195  return;
196  }
197  if (minSortUp(p) && mmCmpExch(0, -1))
198  maxSortDown(-1);
199  }
200  else if (p < 0) //new item is in maxheap
201  {
202  if (maxCt < N / 2)
203  {
204  ++maxCt;
205  }
206  else if (v < old)
207  {
208  maxSortDown(p);
209  return;
210  }
211  if (maxSortUp(p) && minCt && mmCmpExch(1, 0))
212  minSortDown(1);
213  }
214  else //new item is at median
215  {
216  if (maxCt && maxSortUp(-1))
217  maxSortDown(-1);
218  if (minCt && minSortUp(1))
219  minSortDown(1);
220  }
221  }
222 
223  //returns median item (or average of 2 when item count is even)
224  Item median() const
225  {
226  Item v = data[heap[0]];
227  if (minCt < maxCt)
228  {
229  v = (v + data[heap[-1]]) / 2;
230  }
231  return v;
232  }
233 };
234 
235 }
236 }
rolling_median_t & operator=(const rolling_median_t &)=delete
void * memcpy(void *a, const void *b, size_t c)
rolling_median_t(rolling_median_t &&m)
rolling_median_t(const rolling_median_t &)=delete
rolling_median_t & operator=(rolling_median_t &&m)