C에서 롤링 중간 알고리즘
현재 C에서 롤링 중간 필터 (롤링 평균 필터와 유사)를 구현하는 알고리즘을 연구하고 있습니다. 문헌을 검색 한 결과 합리적으로 효율적인 두 가지 방법이 있습니다. 첫 번째는 초기 값 창을 정렬 한 다음 이진 검색을 수행하여 새 값을 삽입하고 각 반복에서 기존 값을 제거하는 것입니다.
두 번째 (Hardle and Steiger, 1995, JRSS-C, Algorithm 296)는 한쪽 끝에 maxheap, 다른쪽에 minheap, 중간에 중앙값을 가진 이중 끝 힙 구조를 만듭니다. 이것은 O (n log n)가 아닌 선형 시간 알고리즘을 생성합니다.
내 문제는 다음과 같습니다. 전자를 구현하는 것이 가능하지만 수백만 개의 시계열에서이를 실행해야하므로 효율성이 중요합니다. 후자는 구현하기가 매우 어렵습니다. R의 통계 패키지 코드의 Trunmed.c 파일에서 코드를 찾았지만 해독 할 수는 없습니다.
선형 타임 롤링 중간 알고리즘에 대해 잘 작성된 C 구현을 아는 사람이 있습니까?
편집 : Trunmed.c 코드 http://google.com/codesearch/p?hl=ko&sa=N&cd=1&ct=rc#mYw3h_Lb_e0/R-2.2.0/src/library/stats/src/Trunmed.c에 링크
src/library/stats/src/Trunmed.c독립형 C ++ 클래스 / C 서브 루틴에서 비슷한 것을 원했기 때문에 R을 몇 번 살펴 보았습니다 . 이것은 실제로 하나의 두 가지 구현 src/library/stats/man/runmed.Rd입니다. (도움말 파일의 소스)를 참조하십시오.
\details{
Apart from the end values, the result \code{y = runmed(x, k)} simply has
\code{y[j] = median(x[(j-k2):(j+k2)])} (k = 2*k2+1), computed very
efficiently.
The two algorithms are internally entirely different:
\describe{
\item{"Turlach"}{is the Härdle-Steiger
algorithm (see Ref.) as implemented by Berwin Turlach.
A tree algorithm is used, ensuring performance \eqn{O(n \log
k)}{O(n * log(k))} where \code{n <- length(x)} which is
asymptotically optimal.}
\item{"Stuetzle"}{is the (older) Stuetzle-Friedman implementation
which makes use of median \emph{updating} when one observation
enters and one leaves the smoothing window. While this performs as
\eqn{O(n \times k)}{O(n * k)} which is slower asymptotically, it is
considerably faster for small \eqn{k} or \eqn{n}.}
}
}
이것을 더 독립적 인 방식으로 재사용하는 것이 좋을 것입니다. 당신은 자원 봉사하고 있습니까? R 비트 중 일부를 도울 수 있습니다.
편집 1 : 위의 이전 버전의 Trunmed.c에 대한 링크 외에도 현재 SVN 사본은 다음과 같습니다.
편집 2 : Ryan Tibshirani에는 빠른 중앙 비닝 에 대한 C 및 Fortran 코드 가 있으며 이는 창 방식에 적합한 시작점이 될 수 있습니다.
주문 통계가있는 c ++ 데이터 구조의 현대적인 구현을 찾을 수 없으므로 MAK가 제안한 최고 코더 링크에서 두 가지 아이디어를 구현했습니다 ( Match Editorial : ScrollDown to FloatingMedian).
두 개의 멀티 세트
첫 번째 아이디어는 삽입 / 삭제 당 O (ln N)를 사용하여 데이터를 두 개의 데이터 구조 (힙, 멀티 세트 등)로 분할하므로 대량 비용없이 Quantile을 동적으로 변경할 수 없습니다. 즉, 롤링 중간 값 또는 롤링 75 %를 가질 수 있지만 동시에 둘다는 아닙니다.
세그먼트 트리
두 번째 아이디어는 삽입 / 삭제 / 쿼리에 O (ln N) 인 세그먼트 트리를 사용하지만 더 유연합니다. "N"중 가장 좋은 것은 데이터 범위의 크기입니다. 따라서 롤링 중간에 창에 백만 개의 항목이 있지만 데이터가 1..65536과 다르면 롤링 창의 이동 당 1 백만 개의 작업 만 필요합니다 !!
C ++ 코드는 Denis가 위에서 게시 한 것과 유사합니다 ( "양자화 된 데이터를위한 간단한 알고리즘").
GNU 순서 통계 트리
포기하기 직전에 stdlibc ++에 주문 통계 트리가 포함되어 있음을 알았습니다 !!!
여기에는 두 가지 중요한 작업이 있습니다.
iter = tree.find_by_order(value)
order = tree.order_of_key(value)
libstdc ++ manual policy_based_data_structures_test ( "분할 및 결합"검색)를 참조하십시오 .
c ++ 0x / c ++ 11 스타일 부분 typedef를 지원하는 컴파일러의 편의 헤더에서 사용하기 위해 트리를 래핑했습니다.
#if !defined(GNU_ORDER_STATISTIC_SET_H)
#define GNU_ORDER_STATISTIC_SET_H
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/tree_policy.hpp>
// A red-black tree table storing ints and their order
// statistics. Note that since the tree uses
// tree_order_statistics_node_update as its update policy, then it
// includes its methods by_order and order_of_key.
template <typename T>
using t_order_statistic_set = __gnu_pbds::tree<
T,
__gnu_pbds::null_type,
std::less<T>,
__gnu_pbds::rb_tree_tag,
// This policy updates nodes' metadata for order statistics.
__gnu_pbds::tree_order_statistics_node_update>;
#endif //GNU_ORDER_STATISTIC_SET_H
여기서 C 구현을 수행했습니다 . 이 질문에는 몇 가지 자세한 내용이 있습니다. C-Turlach 구현에서 롤링 중간 값 .
샘플 사용법 :
int main(int argc, char* argv[])
{
int i,v;
Mediator* m = MediatorNew(15);
for (i=0;i<30;i++)
{
v = rand()&127;
printf("Inserting %3d \n",v);
MediatorInsert(m,v);
v=MediatorMedian(m);
printf("Median = %3d.\n\n",v);
ShowTree(m);
}
}
이 증분 중앙값 추정기를 사용합니다.
median += eta * sgn(sample - median)
보다 일반적인 평균 추정기와 같은 형식입니다.
mean += eta * (sample - mean)
여기 ETA는 작은 학습율 파라미터 (예이다 0.001), 및 sgn()의 하나를 리턴 시그넘 함수이다 {-1, 0, 1}. ( eta데이터가 고정적이지 않고 시간이 지남에 따라 변경 사항을 추적하려는 경우 이와 같은 상수를 사용하십시오 . 그렇지 않으면 고정 소스의 경우 eta = 1 / n수렴 하는 것과 같은 것을 사용 n하십시오. 지금까지 샘플 수는 어디 입니까?)
또한 중간 추정량을 수정하여 임의의 Quantile에서 작동하도록했습니다. 일반적으로 Quantile 함수 는 데이터를 두 부분으로 나누는 값을 알려줍니다 : p및 1 - p. 다음은이 값을 점진적으로 추정합니다.
quantile += eta * (sgn(sample - quantile) + 2.0 * p - 1.0)
값 p은 내에 있어야합니다 [0, 1]. 이것은 본질적으로 sgn()함수의 대칭 출력 {-1, 0, 1}을 한쪽으로 기울여서 데이터 샘플을 두 개의 같지 않은 크기의 빈으로 분할합니다 (분수 p및 1 - p데이터는 각각 Quantile 추정값보다 작거나 큽니다). 의 p = 0.5경우 중앙값 추정기로 줄어 듭니다.
다음은 양자화 된 데이터에 대한 간단한 알고리즘입니다 (몇 달 후).
""" median1.py: moving median 1d for quantized, e.g. 8-bit data
Method: cache the median, so that wider windows are faster.
The code is simple -- no heaps, no trees.
Keywords: median filter, moving median, running median, numpy, scipy
See Perreault + Hebert, Median Filtering in Constant Time, 2007,
http://nomis80.org/ctmf.html: nice 6-page paper and C code,
mainly for 2d images
Example:
y = medians( x, window=window, nlevel=nlevel )
uses:
med = Median1( nlevel, window, counts=np.bincount( x[0:window] ))
med.addsub( +, - ) -- see the picture in Perreault
m = med.median() -- using cached m, summ
How it works:
picture nlevel=8, window=3 -- 3 1s in an array of 8 counters:
counts: . 1 . . 1 . 1 .
sums: 0 1 1 1 2 2 3 3
^ sums[3] < 2 <= sums[4] <=> median 4
addsub( 0, 1 ) m, summ stay the same
addsub( 5, 1 ) slide right
addsub( 5, 6 ) slide left
Updating `counts` in an `addsub` is trivial, updating `sums` is not.
But we can cache the previous median `m` and the sum to m `summ`.
The less often the median changes, the faster;
so fewer levels or *wider* windows are faster.
(Like any cache, run time varies a lot, depending on the input.)
See also:
scipy.signal.medfilt -- runtime roughly ~ window size
http://stackoverflow.com/questions/1309263/rolling-median-algorithm-in-c
"""
from __future__ import division
import numpy as np # bincount, pad0
__date__ = "2009-10-27 oct"
__author_email__ = "denis-bz-py at t-online dot de"
#...............................................................................
class Median1:
""" moving median 1d for quantized, e.g. 8-bit data """
def __init__( s, nlevel, window, counts ):
s.nlevel = nlevel # >= len(counts)
s.window = window # == sum(counts)
s.half = (window // 2) + 1 # odd or even
s.setcounts( counts )
def median( s ):
""" step up or down until sum cnt to m-1 < half <= sum to m """
if s.summ - s.cnt[s.m] < s.half <= s.summ:
return s.m
j, sumj = s.m, s.summ
if sumj <= s.half:
while j < s.nlevel - 1:
j += 1
sumj += s.cnt[j]
# print "j sumj:", j, sumj
if sumj - s.cnt[j] < s.half <= sumj: break
else:
while j > 0:
sumj -= s.cnt[j]
j -= 1
# print "j sumj:", j, sumj
if sumj - s.cnt[j] < s.half <= sumj: break
s.m, s.summ = j, sumj
return s.m
def addsub( s, add, sub ):
s.cnt[add] += 1
s.cnt[sub] -= 1
assert s.cnt[sub] >= 0, (add, sub)
if add <= s.m:
s.summ += 1
if sub <= s.m:
s.summ -= 1
def setcounts( s, counts ):
assert len(counts) <= s.nlevel, (len(counts), s.nlevel)
if len(counts) < s.nlevel:
counts = pad0__( counts, s.nlevel ) # numpy array / list
sumcounts = sum(counts)
assert sumcounts == s.window, (sumcounts, s.window)
s.cnt = counts
s.slowmedian()
def slowmedian( s ):
j, sumj = -1, 0
while sumj < s.half:
j += 1
sumj += s.cnt[j]
s.m, s.summ = j, sumj
def __str__( s ):
return ("median %d: " % s.m) + \
"".join([ (" ." if c == 0 else "%2d" % c) for c in s.cnt ])
#...............................................................................
def medianfilter( x, window, nlevel=256 ):
""" moving medians, y[j] = median( x[j:j+window] )
-> a shorter list, len(y) = len(x) - window + 1
"""
assert len(x) >= window, (len(x), window)
# np.clip( x, 0, nlevel-1, out=x )
# cf http://scipy.org/Cookbook/Rebinning
cnt = np.bincount( x[0:window] )
med = Median1( nlevel=nlevel, window=window, counts=cnt )
y = (len(x) - window + 1) * [0]
y[0] = med.median()
for j in xrange( len(x) - window ):
med.addsub( x[j+window], x[j] )
y[j+1] = med.median()
return y # list
# return np.array( y )
def pad0__( x, tolen ):
""" pad x with 0 s, numpy array or list """
n = tolen - len(x)
if n > 0:
try:
x = np.r_[ x, np.zeros( n, dtype=x[0].dtype )]
except NameError:
x += n * [0]
return x
#...............................................................................
if __name__ == "__main__":
Len = 10000
window = 3
nlevel = 256
period = 100
np.set_printoptions( 2, threshold=100, edgeitems=10 )
# print medians( np.arange(3), 3 )
sinwave = (np.sin( 2 * np.pi * np.arange(Len) / period )
+ 1) * (nlevel-1) / 2
x = np.asarray( sinwave, int )
print "x:", x
for window in ( 3, 31, 63, 127, 255 ):
if window > Len: continue
print "medianfilter: Len=%d window=%d nlevel=%d:" % (Len, window, nlevel)
y = medianfilter( x, window=window, nlevel=nlevel )
print np.array( y )
# end median1.py
두 개의 파티션을 유지함으로써 롤링 중간 값을 찾을 수 있습니다.
파티션을 유지 관리하려면 최소 힙 및 최대 힙을 사용하십시오.
최대 힙에는 중간 값보다 작은 숫자가 포함됩니다.
최소 힙에는 중간 값보다 큰 숫자가 포함됩니다.
균형 구속 조건 : 총 요소 수가 짝수이면 두 힙이 동일한 요소를 가져야합니다.
총 요소 수가 홀수이면 Max Heap은 Min Heap보다 하나 이상의 요소를 갖습니다.
중앙값 요소 : 두 파티션에 동일한 수의 요소가있는 경우 중간 값은 첫 번째 파티션의 최대 요소와 두 번째 파티션의 최소 요소의 합의 절반입니다.
그렇지 않으면 중간 값은 첫 번째 파티션의 최대 요소입니다.
연산- 1- 2 개의 힙을 가져옵니다 (최소 힙 1 개 및 최대 힙 1 개). 최대 힙에는 절반 수의 요소가 포함됩니다. 최소 힙에는 절반 수의 요소가 포함됩니다. 2- 스트림의 새로운 숫자와 Max Heap의 상단을 비교합니다. 작거나 같으면 최대 힙에 해당 숫자를 추가하십시오. 그렇지 않으면 최소 힙에 숫자를 추가하십시오. 3- 최소 힙에 최대 힙보다 많은 요소가있는 경우 그런 다음 Min Heap의 상단 요소를 제거하고 Max Heap을 추가하십시오. 최대 힙에 최소 힙보다 두 개 이상의 요소가있는 경우 그런 다음 Max Heap의 상단 요소를 제거하고 Min Heap을 추가하십시오. 4- 두 힙의 요소 수가 동일한 경우 중앙값은 최대 힙의 최대 요소와 최소 힙의 최소 요소의 합의 절반입니다. 그렇지 않으면 첫 번째 파티션에서 중간 값이 max 요소가됩니다.
public class Solution {
public static void main(String[] args) {
Scanner in = new Scanner(System.in);
RunningMedianHeaps s = new RunningMedianHeaps();
int n = in.nextInt();
for(int a_i=0; a_i < n; a_i++){
printMedian(s,in.nextInt());
}
in.close();
}
public static void printMedian(RunningMedianHeaps s, int nextNum){
s.addNumberInHeap(nextNum);
System.out.printf("%.1f\n",s.getMedian());
}
}
class RunningMedianHeaps{
PriorityQueue<Integer> minHeap = new PriorityQueue<Integer>();
PriorityQueue<Integer> maxHeap = new PriorityQueue<Integer>(Comparator.reverseOrder());
public double getMedian() {
int size = minHeap.size() + maxHeap.size();
if(size % 2 == 0)
return (maxHeap.peek()+minHeap.peek())/2.0;
return maxHeap.peek()*1.0;
}
private void balanceHeaps() {
if(maxHeap.size() < minHeap.size())
{
maxHeap.add(minHeap.poll());
}
else if(maxHeap.size() > 1+minHeap.size())
{
minHeap.add(maxHeap.poll());
}
}
public void addNumberInHeap(int num) {
if(maxHeap.size()==0 || num <= maxHeap.peek())
{
maxHeap.add(num);
}
else
{
minHeap.add(num);
}
balanceHeaps();
}
}
특정 시점의 함수로 값을 참조하는 기능이 있는 경우 부트 스트랩 을 적용하여 신뢰 구간 내에서 부트 스트랩 된 중앙값을 생성하여 대체하여 값을 샘플링 할 수 있습니다. 이를 통해 들어오는 값을 데이터 구조에 지속적으로 정렬하는 것보다 효율성이 높은 근사값을 계산할 수 있습니다.
Java에서 실행중인 중간 값이 필요한 사람들에게는 ... PriorityQueue가 당신의 친구입니다. O (log N) 삽입, O (1) 현재 중앙값 및 O (N) 제거. 데이터의 분포를 알고 있다면 이것보다 훨씬 더 잘 할 수 있습니다.
public class RunningMedian {
// Two priority queues, one of reversed order.
PriorityQueue<Integer> lower = new PriorityQueue<Integer>(10,
new Comparator<Integer>() {
public int compare(Integer arg0, Integer arg1) {
return (arg0 < arg1) ? 1 : arg0 == arg1 ? 0 : -1;
}
}), higher = new PriorityQueue<Integer>();
public void insert(Integer n) {
if (lower.isEmpty() && higher.isEmpty())
lower.add(n);
else {
if (n <= lower.peek())
lower.add(n);
else
higher.add(n);
rebalance();
}
}
void rebalance() {
if (lower.size() < higher.size() - 1)
lower.add(higher.remove());
else if (higher.size() < lower.size() - 1)
higher.add(lower.remove());
}
public Integer getMedian() {
if (lower.isEmpty() && higher.isEmpty())
return null;
else if (lower.size() == higher.size())
return (lower.peek() + higher.peek()) / 2;
else
return (lower.size() < higher.size()) ? higher.peek() : lower
.peek();
}
public void remove(Integer n) {
if (lower.remove(n) || higher.remove(n))
rebalance();
}
}
스트림의 모든 값이 (상대적으로) 작은 정의 된 범위 내의 정수인 경우 간단한 정확한 솔루션을 가진 특별한 경우가 있음을 지적 할 가치가 있습니다. 예를 들어, 모두 0과 1023 사이 여야한다고 가정합니다.이 경우 1024 개의 요소와 개수로 구성된 배열을 정의하고이 값을 모두 지우십시오. 스트림의 각 값에 대해 해당하는 구간 및 개수를 증가시킵니다. 스트림이 종료 된 후 0부터 시작하는 연속 빈을 추가하여 쉽게 카운트 / 2 최고 값을 포함하는 빈을 찾습니다. 동일한 방법을 사용하여 임의의 순위 순서 값을 찾을 수 있습니다. 빈 포화를 감지하고 실행 중에 저장소 빈의 크기를 더 큰 유형으로 "업그레이드"해야하는 경우에는 약간의 문제가 있습니다.
This special case may seem artificial, but in practice it is very common. It can also be applied as an approximation for real numbers if they lie within a range and a "good enough" level of precision is known. This would hold for pretty much any set of measurements on a group of "real world" objects. For instance, the heights or weights of a group of people. Not a big enough set? It would work just as well for the lengths or weights of all the (individual) bacteria on the planet - assuming somebody could supply the data!
원본을 잘못 읽은 것처럼 보입니다. 매우 긴 스트림의 중간 값 대신 슬라이딩 창 중간 값을 원하는 것 같습니다. 이 방법은 여전히 유효합니다. 초기 창의 첫 번째 N 스트림 값을로드 한 다음 N + 1 번째 스트림 값의 경우 해당 Bin을 증가시키면서 0 번째 스트림 값에 해당하는 Bin을 줄입니다. 이 경우 감소를 허용하기 위해 마지막 N 값을 유지해야합니다. 감소는 크기 N의 배열을 주기적으로 지정하여 효율적으로 수행 할 수 있습니다. 중앙값의 위치는 -2, -1,0,1 만 변경 될 수 있기 때문에 , 2 슬라이딩 윈도우의 각 단계에서, 모든 빈을 각 단계의 중앙값까지 합할 필요는 없으며 수정 된 측면 빈에 따라 "중앙 포인터"를 조정하면됩니다. 예를 들어 새 값과 제거되는 값이 모두 현재 중앙값 아래로 떨어지면 변경되지 않습니다 (오프셋 = 0). 메모리에 편리하게 보관하기에 N이 너무 커지면이 방법이 중단됩니다.
다음은 정확한 출력이 중요하지 않을 때 사용할 수있는 것입니다 (표시 목적 등). totalcount 및 lastmedian과 새로운 값이 필요합니다.
{
totalcount++;
newmedian=lastmedian+(newvalue>lastmedian?1:-1)*(lastmedian==0?newvalue: lastmedian/totalcount*2);
}
page_display_time과 같은 것에 대해 매우 정확한 결과를 생성합니다.
규칙 : 입력 스트림은 페이지 표시 시간 순서대로 매끄럽고 개수 (> 30 등)가 많아야하며 중앙값이 0이 아니어야합니다.
예 : 페이지로드 시간, 800 개 항목, 10ms ... 3000ms, 평균 90ms, 실제 중앙값 : 11ms
30 번의 입력 후, 중앙값 오차는 일반적으로 <= 20 % (9ms..12ms)이며 점점 줄어 듭니다. 800 개의 입력 후 오류는 + -2 %입니다.
비슷한 솔루션을 가진 또 다른 사상가가 있습니다 : Median Filter Super 효율적인 구현
다음은 자바 구현입니다.
package MedianOfIntegerStream;
import java.util.Comparator;
import java.util.HashSet;
import java.util.Iterator;
import java.util.Set;
import java.util.TreeSet;
public class MedianOfIntegerStream {
public Set<Integer> rightMinSet;
public Set<Integer> leftMaxSet;
public int numOfElements;
public MedianOfIntegerStream() {
rightMinSet = new TreeSet<Integer>();
leftMaxSet = new TreeSet<Integer>(new DescendingComparator());
numOfElements = 0;
}
public void addNumberToStream(Integer num) {
leftMaxSet.add(num);
Iterator<Integer> iterMax = leftMaxSet.iterator();
Iterator<Integer> iterMin = rightMinSet.iterator();
int maxEl = iterMax.next();
int minEl = 0;
if (iterMin.hasNext()) {
minEl = iterMin.next();
}
if (numOfElements % 2 == 0) {
if (numOfElements == 0) {
numOfElements++;
return;
} else if (maxEl > minEl) {
iterMax.remove();
if (minEl != 0) {
iterMin.remove();
}
leftMaxSet.add(minEl);
rightMinSet.add(maxEl);
}
} else {
if (maxEl != 0) {
iterMax.remove();
}
rightMinSet.add(maxEl);
}
numOfElements++;
}
public Double getMedian() {
if (numOfElements % 2 != 0)
return new Double(leftMaxSet.iterator().next());
else
return (leftMaxSet.iterator().next() + rightMinSet.iterator().next()) / 2.0;
}
private class DescendingComparator implements Comparator<Integer> {
@Override
public int compare(Integer o1, Integer o2) {
return o2 - o1;
}
}
public static void main(String[] args) {
MedianOfIntegerStream streamMedian = new MedianOfIntegerStream();
streamMedian.addNumberToStream(1);
System.out.println(streamMedian.getMedian()); // should be 1
streamMedian.addNumberToStream(5);
streamMedian.addNumberToStream(10);
streamMedian.addNumberToStream(12);
streamMedian.addNumberToStream(2);
System.out.println(streamMedian.getMedian()); // should be 5
streamMedian.addNumberToStream(3);
streamMedian.addNumberToStream(8);
streamMedian.addNumberToStream(9);
System.out.println(streamMedian.getMedian()); // should be 6.5
}
}
평활 평균을 요구하는 경우 빠르고 쉬운 방법은 최신 값에 x를 곱하고 평균 값에 (1-x)를 곱하는 것입니다. 이것은 새로운 평균이됩니다.
편집 : 사용자가 요청한 것이 아니며 통계적으로 유효하지는 않지만 많은 용도로 충분합니다.
검색을 위해 여기에 (공감대에도 불구하고) 남겨 두겠습니다!
참고 URL : https://stackoverflow.com/questions/1309263/rolling-median-algorithm-in-c
'programing tip' 카테고리의 다른 글
| 각도기를 사용하여 요소가 보이는지 확인하는 방법은 무엇입니까? (0) | 2020.08.07 |
|---|---|
| Entity Framework에 개체가 있는지 확인하는 가장 좋은 방법은 무엇입니까? (0) | 2020.08.07 |
| 속기 배경 속성의 배경 크기 (CSS3) (0) | 2020.08.06 |
| Visual Studio에서 Interop 유형 포함 참과 거짓의 차이점은 무엇입니까? (0) | 2020.08.06 |
| 바이트에 대한 적절한 최대 절전 모드 주석 (0) | 2020.08.06 |