Chris@1573: /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
Chris@1573: 
Chris@1573: /*
Chris@1573:     Sonic Visualiser
Chris@1573:     An audio file viewer and annotation editor.
Chris@1573:     Centre for Digital Music, Queen Mary, University of London.
Chris@1573:     This file copyright 2007-2015 Particular Programs Ltd, 
Chris@1573:     copyright 2018 Queen Mary University of London.
Chris@1573:     
Chris@1573:     This program is free software; you can redistribute it and/or
Chris@1573:     modify it under the terms of the GNU General Public License as
Chris@1573:     published by the Free Software Foundation; either version 2 of the
Chris@1573:     License, or (at your option) any later version.  See the file
Chris@1573:     COPYING included with this distribution for more information.
Chris@1573: */
Chris@1573: 
Chris@1573: #ifndef SV_MOVING_MEDIAN_H
Chris@1573: #define SV_MOVING_MEDIAN_H
Chris@1573: 
Chris@1573: #include <bqvec/Allocators.h>
Chris@1573: #include <bqvec/VectorOps.h>
Chris@1573: 
Chris@1573: #include <algorithm>
Chris@1573: #include <iostream>
Chris@1573: #include <stdexcept>
Chris@1573: 
Chris@1573: /**
Chris@1573:  * Obtain the median (or other percentile) of a moving window across a
Chris@1573:  * time series. Construct the MovingMedian object, then push() each
Chris@1573:  * new value in the time series and get() the median of the most
Chris@1573:  * recent window. The size of the window, and the percentile
Chris@1573:  * calculated, can both be changed after construction.
Chris@1573:  *
Chris@1573:  * Note that for even-sized windows, the "median" is taken to be the
Chris@1573:  * value at the start of the second half when sorted, e.g. for size 4,
Chris@1573:  * the element at index 2 (zero-based) in the sorted window.
Chris@1573:  *
Chris@1573:  * Not thread-safe.
Chris@1573:  */
Chris@1573: template <typename T>
Chris@1573: class MovingMedian
Chris@1573: {
Chris@1573: public:
Chris@1573:     MovingMedian(int size, double percentile = 50.f) :
Chris@1573:         m_size(size),
Chris@1573:         m_percentile(percentile) {
Chris@1573:         if (size < 1) throw std::logic_error("size must be >= 1");
Chris@1573:         m_frame = breakfastquay::allocate_and_zero<T>(size);
Chris@1573: 	m_sorted = breakfastquay::allocate_and_zero<T>(size);
Chris@1573:         calculateIndex();
Chris@1573:     }
Chris@1573: 
Chris@1573:     ~MovingMedian() { 
Chris@1573:         breakfastquay::deallocate(m_frame);
Chris@1573:         breakfastquay::deallocate(m_sorted);
Chris@1573:     }
Chris@1573: 
Chris@1573:     MovingMedian(const MovingMedian &) =delete;
Chris@1573:     MovingMedian &operator=(const MovingMedian &) =delete;
Chris@1573: 
Chris@1573:     void setPercentile(double p) {
Chris@1573:         m_percentile = p;
Chris@1573:         calculateIndex();
Chris@1573:     }
Chris@1573: 
Chris@1573:     void push(T value) {
Chris@1573:         if (value != value) {
Chris@1573:             std::cerr << "WARNING: MovingMedian: NaN encountered" << std::endl;
Chris@1573:             value = T();
Chris@1573:         }
Chris@1573: 	drop(m_frame[0]);
Chris@1573:         breakfastquay::v_move(m_frame, m_frame+1, m_size-1);
Chris@1573: 	m_frame[m_size-1] = value;
Chris@1573: 	put(value);
Chris@1573:     }
Chris@1573: 
Chris@1573:     T get() const {
Chris@1573: 	return m_sorted[m_index];
Chris@1573:     }
Chris@1573: 
Chris@1573:     int size() const {
Chris@1573:         return m_size;
Chris@1573:     }
Chris@1573: 
Chris@1573:     void reset() {
Chris@1573:         breakfastquay::v_zero(m_frame, m_size);
Chris@1573:         breakfastquay::v_zero(m_sorted, m_size);
Chris@1573:     }
Chris@1573: 
Chris@1573:     void resize(int target) {
Chris@1573:         if (target == m_size) return;
Chris@1573:         int diff = std::abs(target - m_size);
Chris@1573:         if (target > m_size) { // grow
Chris@1573:             // we don't want to change the median, so fill spaces with it
Chris@1573:             T fillValue = get();
Chris@1573:             m_frame = breakfastquay::reallocate(m_frame, m_size, target);
Chris@1573:             breakfastquay::v_move(m_frame + diff, m_frame, m_size);
Chris@1573:             breakfastquay::v_set(m_frame, fillValue, diff);
Chris@1573:             m_sorted = breakfastquay::reallocate(m_sorted, m_size, target);
Chris@1573:             for (int sz = m_size + 1; sz <= target; ++sz) {
Chris@1573:                 put(m_sorted, sz, fillValue);
Chris@1573:             }
Chris@1573:         } else { // shrink
Chris@1573:             for (int i = 0; i < diff; ++i) {
Chris@1573:                 drop(m_sorted, m_size - i, m_frame[i]);
Chris@1573:             }
Chris@1573:             m_sorted = breakfastquay::reallocate(m_sorted, m_size, target);
Chris@1573:             breakfastquay::v_move(m_frame, m_frame + diff, target);
Chris@1573:             m_frame = breakfastquay::reallocate(m_frame, m_size, target);
Chris@1573:         }
Chris@1573: 
Chris@1573:         m_size = target;
Chris@1573:         calculateIndex();
Chris@1573:     }
Chris@1573: 
Chris@1573:     void checkIntegrity() const {
Chris@1573:         check();
Chris@1573:     }
Chris@1573: 
Chris@1573: private:
Chris@1573:     int m_size;
Chris@1573:     double m_percentile;
Chris@1573:     int m_index;
Chris@1573:     T *m_frame;
Chris@1573:     T *m_sorted;
Chris@1573: 
Chris@1573:     void calculateIndex() {
Chris@1573:         m_index = int((m_size * m_percentile) / 100.f);
Chris@1573:         if (m_index >= m_size) m_index = m_size-1;
Chris@1573:         if (m_index < 0) m_index = 0;
Chris@1573:     }
Chris@1573:     
Chris@1573:     void put(T value) {
Chris@1573:         put(m_sorted, m_size, value);
Chris@1573:     }
Chris@1573: 
Chris@1573:     static void put(T *const sorted, int size, T value) {
Chris@1573: 
Chris@1573:         // precondition: sorted points to size-1 sorted values,
Chris@1573:         // followed by an unused slot (i.e. only the first size-1
Chris@1573:         // values of sorted are actually sorted)
Chris@1573:         // 
Chris@1573:         // postcondition: sorted points to size sorted values
Chris@1573:         
Chris@1573: 	T *ptr = std::lower_bound(sorted, sorted + size - 1, value);
Chris@1573:         breakfastquay::v_move(ptr + 1, ptr, int(sorted + size - ptr) - 1);
Chris@1573: 	*ptr = value;
Chris@1573:     }
Chris@1573: 
Chris@1573:     void drop(T value) {
Chris@1573:         drop(m_sorted, m_size, value);
Chris@1573:     }
Chris@1573: 
Chris@1573:     static void drop(T *const sorted, int size, T value) {
Chris@1573: 
Chris@1573:         // precondition: sorted points to size sorted values, one of
Chris@1573:         // which is value
Chris@1573:         //
Chris@1573:         // postcondition: sorted points to size-1 sorted values,
Chris@1573:         // followed by a slot that has been reset to default value
Chris@1573:         // (i.e. only the first size-1 values of sorted are actually
Chris@1573:         // sorted)
Chris@1573: 
Chris@1573: 	T *ptr = std::lower_bound(sorted, sorted + size, value);
Chris@1573: 	if (*ptr != value) {
Chris@1573:             throw std::logic_error
Chris@1573:                 ("MovingMedian::drop: value being dropped is not in array");
Chris@1573:         }
Chris@1573:         breakfastquay::v_move(ptr, ptr + 1, int(sorted + size - ptr) - 1);
Chris@1573:         sorted[size-1] = T();
Chris@1573:     }
Chris@1573: 
Chris@1573:     void check() const {
Chris@1573:         bool good = true;
Chris@1573:         for (int i = 1; i < m_size; ++i) {
Chris@1573:             if (m_sorted[i] < m_sorted[i-1]) {
Chris@1573:                 std::cerr << "ERROR: MovingMedian::checkIntegrity: "
Chris@1573:                           << "mis-ordered elements in sorted array starting "
Chris@1573:                           << "at index " << i << std::endl;
Chris@1573:                 good = false;
Chris@1573:                 break;
Chris@1573:             }
Chris@1573:         }
Chris@1573:         for (int i = 0; i < m_size; ++i) {
Chris@1573:             bool found = false;
Chris@1573:             for (int j = 0; j < m_size; ++j) {
Chris@1573:                 if (m_sorted[j] == m_frame[i]) {
Chris@1573:                     found = true;
Chris@1573:                     break;
Chris@1573:                 }
Chris@1573:             }
Chris@1573:             if (!found) {
Chris@1573:                 std::cerr << "ERROR: MovingMedian::checkIntegrity: "
Chris@1573:                           << "element in frame at index " << i
Chris@1573:                           << " not found in sorted array" << std::endl;
Chris@1573:                 good = false;
Chris@1573:                 break;
Chris@1573:             }
Chris@1573:         }
Chris@1573:         for (int i = 0; i < m_size; ++i) {
Chris@1573:             bool found = false;
Chris@1573:             for (int j = 0; j < m_size; ++j) {
Chris@1573:                 if (m_sorted[i] == m_frame[j]) {
Chris@1573:                     found = true;
Chris@1573:                     break;
Chris@1573:                 }
Chris@1573:             }
Chris@1573:             if (!found) {
Chris@1573:                 std::cerr << "ERROR: MovingMedian::checkIntegrity: "
Chris@1573:                           << "element in sorted array at index " << i
Chris@1573:                           << " not found in source frame" << std::endl;
Chris@1573:                 good = false;
Chris@1573:                 break;
Chris@1573:             }
Chris@1573:         }
Chris@1573:         if (!good) {
Chris@1573:             std::cerr << "Frame contains:" << std::endl;
Chris@1573:             std::cerr << "[ ";
Chris@1573:             for (int j = 0; j < m_size; ++j) {
Chris@1573:                 std::cerr << m_frame[j] << " ";
Chris@1573:             }
Chris@1573:             std::cerr << "]" << std::endl;
Chris@1573:             std::cerr << "Sorted array contains:" << std::endl;
Chris@1573:             std::cerr << "[ ";
Chris@1573:             for (int j = 0; j < m_size; ++j) {
Chris@1573:                 std::cerr << m_sorted[j] << " ";
Chris@1573:             }
Chris@1573:             std::cerr << "]" << std::endl;
Chris@1573:             throw std::logic_error("MovingMedian failed integrity check");
Chris@1573:         }
Chris@1573:     }
Chris@1573: };
Chris@1573: 
Chris@1573: #endif
Chris@1573: