1 /* movstat/qqracc.c
2 *
3 * Moving window QQR accumulator
4 *
5 * Copyright (C) 2018 Patrick Alken
6 *
7 * This program is free software; you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation; either version 3 of the License, or (at
10 * your option) any later version.
11 *
12 * This program is distributed in the hope that it will be useful, but
13 * WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 * General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with this program; if not, write to the Free Software
19 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
20 */
21
22 #include <config.h>
23 #include <gsl/gsl_math.h>
24 #include <gsl/gsl_vector.h>
25 #include <gsl/gsl_errno.h>
26 #include <gsl/gsl_movstat.h>
27 #include <gsl/gsl_sort.h>
28 #include <gsl/gsl_statistics.h>
29
30 typedef double qqracc_type_t;
31 typedef qqracc_type_t ringbuf_type_t;
32
33 #include "ringbuf.c"
34
35 typedef struct
36 {
37 qqracc_type_t *window; /* linear array for current window */
38 ringbuf *rbuf; /* ring buffer storing current window */
39 } qqracc_state_t;
40
41 static size_t
qqracc_size(const size_t n)42 qqracc_size(const size_t n)
43 {
44 size_t size = 0;
45
46 size += sizeof(qqracc_state_t);
47 size += n * sizeof(qqracc_type_t);
48 size += ringbuf_size(n);
49
50 return size;
51 }
52
53 static int
qqracc_init(const size_t n,void * vstate)54 qqracc_init(const size_t n, void * vstate)
55 {
56 qqracc_state_t * state = (qqracc_state_t *) vstate;
57
58 state->window = (qqracc_type_t *) ((unsigned char *) vstate + sizeof(qqracc_state_t));
59 state->rbuf = (ringbuf *) ((unsigned char *) state->window + n * sizeof(qqracc_type_t));
60
61 ringbuf_init(n, state->rbuf);
62
63 return GSL_SUCCESS;
64 }
65
66 static int
qqracc_insert(const qqracc_type_t x,void * vstate)67 qqracc_insert(const qqracc_type_t x, void * vstate)
68 {
69 qqracc_state_t * state = (qqracc_state_t *) vstate;
70
71 /* add new element to ring buffer */
72 ringbuf_insert(x, state->rbuf);
73
74 return GSL_SUCCESS;
75 }
76
77 static int
qqracc_delete(void * vstate)78 qqracc_delete(void * vstate)
79 {
80 qqracc_state_t * state = (qqracc_state_t *) vstate;
81
82 if (!ringbuf_is_empty(state->rbuf))
83 ringbuf_pop_back(state->rbuf);
84
85 return GSL_SUCCESS;
86 }
87
88 /* FIXME XXX: this is inefficient - could be improved by maintaining a sorted ring buffer */
89 static int
qqracc_get(void * params,qqracc_type_t * result,const void * vstate)90 qqracc_get(void * params, qqracc_type_t * result, const void * vstate)
91 {
92 const qqracc_state_t * state = (const qqracc_state_t *) vstate;
93 double q = *(double *) params;
94 size_t n = ringbuf_copy(state->window, state->rbuf);
95 double quant1, quant2;
96
97 gsl_sort(state->window, 1, n);
98
99 /* compute q-quantile and (1-q)-quantile */
100 quant1 = gsl_stats_quantile_from_sorted_data(state->window, 1, n, q);
101 quant2 = gsl_stats_quantile_from_sorted_data(state->window, 1, n, 1.0 - q);
102
103 /* compute q-quantile range */
104 *result = quant2 - quant1;
105
106 return GSL_SUCCESS;
107 }
108
109 static const gsl_movstat_accum qqr_accum_type =
110 {
111 qqracc_size,
112 qqracc_init,
113 qqracc_insert,
114 qqracc_delete,
115 qqracc_get
116 };
117
118 const gsl_movstat_accum *gsl_movstat_accum_qqr = &qqr_accum_type;
119