Quantile.java
001 /*
002  * Java Genetic Algorithm Library (jenetics-2.0.2).
003  * Copyright (c) 2007-2014 Franz Wilhelmstötter
004  *
005  * Licensed under the Apache License, Version 2.0 (the "License");
006  * you may not use this file except in compliance with the License.
007  * You may obtain a copy of the License at
008  *
009  *      http://www.apache.org/licenses/LICENSE-2.0
010  *
011  * Unless required by applicable law or agreed to in writing, software
012  * distributed under the License is distributed on an "AS IS" BASIS,
013  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014  * See the License for the specific language governing permissions and
015  * limitations under the License.
016  *
017  * Author:
018  *    Franz Wilhelmstötter (franz.wilhelmstoetter@gmx.at)
019  */
020 package org.jenetics.stat;
021 
022 import static java.lang.Double.compare;
023 import static java.lang.String.format;
024 import static org.jenetics.internal.util.object.eq;
025 
026 import java.util.Arrays;
027 
028 import org.jenetics.internal.util.HashBuilder;
029 
030 import org.jenetics.util.MappedAccumulator;
031 
032 
033 /**
034  * Implementation of the quantile estimation algorithm published by
035  <p>
036  <strong>Raj JAIN and Imrich CHLAMTAC</strong>:
037  <em>
038  *     The P<sup>2</sup> Algorithm for Dynamic Calculation of Quantiles and
039  *     Histograms Without Storing Observations
040  </em>
041  <br>
042  * [<a href="http://www.cse.wustl.edu/~jain/papers/ftp/psqr.pdf">Communications
043  * of the ACM; October 1985, Volume 28, Number 10</a>]
044  <p>
045  <strong>Note that this implementation is not synchronized.</strong> If
046  * multiple threads access this object concurrently, and at least one of the
047  * threads modifies it, it must be synchronized externally.
048  *
049  @see <a href="http://en.wikipedia.org/wiki/Quantile">Wikipedia: Quantile</a>
050  *
051  @author <a href="mailto:franz.wilhelmstoetter@gmx.at">Franz Wilhelmstötter</a>
052  @since 1.0
053  @version 2.0 &mdash; <em>$Date: 2014-03-28 $</em>
054  */
055 public class Quantile<N extends Number> extends MappedAccumulator<N> {
056 
057     // The desired quantile.
058     private final double _quantile;
059 
060     // Marker heights.
061     private final double[] _q = {00000};
062 
063     // Marker positions.
064     private final double[] _n = {00000};
065 
066     // Desired marker positions.
067     private final double[] _nn = {000};
068 
069     // Desired marker position increments.
070     private final double[] _dn = {000};
071 
072     private boolean _initialized;
073 
074     /**
075      * Create a new quantile accumulator with the given value.
076      *
077      @param quantile the wished quantile value.
078      @throws IllegalArgumentException if the {@code quantile} is not in the
079      *         range {@code [0, 1]}.
080      */
081     public Quantile(final double quantile) {
082         _quantile = quantile;
083         init(quantile);
084     }
085 
086     private void init(final double quantile) {
087         if (quantile < 0.0 || quantile > 1) {
088             throw new IllegalArgumentException(format(
089                     "Quantile (%s) not in the valid range of [0, 1]", quantile
090                 ));
091         }
092 
093         Arrays.fill(_q, 0);
094         Arrays.fill(_n, 0);
095         Arrays.fill(_nn, 0);
096         Arrays.fill(_dn, 0);
097 
098         _n[0= -1.0;
099         _q[20.0;
100         _initialized = compare(quantile, 0.0== || compare(quantile, 1.0== 0;
101         _samples = 0;
102     }
103 
104     /**
105      * Reset this object to its initial state.
106      */
107     public void reset() {
108         init(_quantile);
109     }
110 
111     /**
112      * Return the computed quantile value.
113      *
114      @return the quantile value.
115      */
116     public double getValue() {
117         return _q[2];
118     }
119 
120     @Override
121     public void accumulate(final N value) {
122         if (!_initialized) {
123             initialize(value.doubleValue());
124         else {
125             update(value.doubleValue());
126         }
127 
128         ++_samples;
129     }
130 
131 
132     private void initialize(double value) {
133         if (_n[00.0) {
134             _n[00.0;
135             _q[0= value;
136         else if (_n[1== 0.0) {
137             _n[11.0;
138             _q[1= value;
139         else if (_n[2== 0.0) {
140             _n[22.0;
141             _q[2= value;
142         else if (_n[3== 0.0) {
143             _n[33.0;
144             _q[3= value;
145         else if (_n[4== 0.0) {
146             _n[44.0;
147             _q[4= value;
148         }
149 
150         if (_n[4!= 0.0) {
151             Arrays.sort(_q);
152 
153             _nn[02.0*_quantile;
154             _nn[14.0*_quantile;
155             _nn[22.0*_quantile + 2.0;
156 
157             _dn[0= _quantile/2.0;
158             _dn[1= _quantile;
159             _dn[2(1.0 + _quantile)/2.0;
160 
161             _initialized = true;
162         }
163     }
164 
165     private void update(double value) {
166         assert (_initialized);
167 
168         // If min or max, handle as special case; otherwise, ...
169         if (_quantile == 0.0) {
170             if (value < _q[2]) {
171                 _q[2= value;
172             }
173         else if (_quantile == 1.0) {
174             if (value > _q[2]) {
175                 _q[2= value;
176             }
177         else {
178             // Increment marker locations and update min and max.
179             if (value < _q[0]) {
180                 ++_n[1]; ++_n[2]; ++_n[3]; ++_n[4]; _q[0= value;
181             else if (value < _q[1]) {
182                 ++_n[1]; ++_n[2]; ++_n[3]; ++_n[4];
183             else if (value < _q[2]) {
184                 ++_n[2]; ++_n[3]; ++_n[4];
185             else if (value < _q[3]) {
186                 ++_n[3]; ++_n[4];
187             else if (value < _q[4]) {
188                 ++_n[4];
189             else {
190                 ++_n[4]; _q[4= value;
191             }
192 
193             // Increment positions of markers k + 1
194             _nn[0+= _dn[0];
195             _nn[1+= _dn[1];
196             _nn[2+= _dn[2];
197 
198             // Adjust heights of markers 0 to 2 if necessary
199             double mm = _n[11.0;
200             double mp = _n[11.0;
201             if (_nn[0>= mp && _n[2> mp) {
202                 _q[1= qPlus(mp, _n[0], _n[1], _n[2], _q[0], _q[1], _q[2]);
203                 _n[1= mp;
204             else if (_nn[0<= mm && _n[0< mm) {
205                 _q[1= qMinus(mm, _n[0], _n[1], _n[2], _q[0], _q[1], _q[2]);
206                 _n[1= mm;
207             }
208 
209             mm = _n[21.0;
210             mp = _n[21.0;
211             if (_nn[1>= mp && _n[3> mp) {
212                 _q[2= qPlus(mp, _n[1], _n[2], _n[3], _q[1], _q[2], _q[3]);
213                 _n[2= mp;
214             else if (_nn[1<= mm && _n[1< mm) {
215                 _q[2= qMinus(mm, _n[1], _n[2], _n[3], _q[1], _q[2], _q[3]);
216                 _n[2= mm;
217             }
218 
219             mm = _n[31.0;
220             mp = _n[31.0;
221             if (_nn[2>= mp && _n[4> mp) {
222                 _q[3= qPlus(mp, _n[2], _n[3], _n[4], _q[2], _q[3], _q[4]);
223                 _n[3= mp;
224             else if (_nn[2<= mm && _n[2< mm) {
225                 _q[3= qMinus(mm, _n[2], _n[3], _n[4], _q[2], _q[3], _q[4]);
226                 _n[3= mm;
227             }
228         }
229     }
230 
231     private static double qPlus(
232         final double mp,
233         final double m0,
234         final double m1,
235         final double m2,
236         final double q0,
237         final double q1,
238         final double q2
239     ) {
240         double result = q1 +
241                     ((mp - m0)*(q2 - q1)/(m2 - m1+
242                     (m2 - mp)*(q1 - q0)/(m1 - m0))/(m2 - m0);
243 
244         if (result > q2) {
245             result = q1 + (q2 - q1)/(m2 - m1);
246         }
247 
248         return result;
249     }
250 
251     private static double qMinus(
252         final double mm,
253         final double m0,
254         final double m1,
255         final double m2,
256         final double q0,
257         final double q1,
258         final double q2
259     ) {
260         double result = q1 -
261                     ((mm - m0)*(q2 - q1)/(m2 - m1+
262                     (m2 - mm)*(q1 - q0)/(m1 - m0))/(m2 - m0);
263 
264         if (q0 > result) {
265             result = q1 + (q0 - q1)/(m0 - m1);
266         }
267 
268         return result;
269     }
270 
271     @Override
272     public int hashCode() {
273         return HashBuilder.of(getClass()).
274                 and(super.hashCode()).
275                 and(_quantile).
276                 and(_dn).
277                 and(_n).
278                 and(_nn).
279                 and(_q).value();
280     }
281 
282     @Override
283     public boolean equals(final Object obj) {
284         if (obj == this) {
285             return true;
286         }
287         if (obj == null || getClass() != obj.getClass()) {
288             return false;
289         }
290 
291         final Quantile<?> quantile = (Quantile<?>)obj;
292         return super.equals(obj&&
293                 eq(_quantile, quantile._quantile&&
294                 eq(_dn, quantile._dn&&
295                 eq(_n, quantile._n&&
296                 eq(_nn, quantile._nn&&
297                 eq(_q, quantile._q);
298     }
299 
300     @Override
301     public String toString() {
302         return format(
303             "%s[samples=%d, quantile=%f]",
304             getClass().getSimpleName(), getSamples(), getValue()
305         );
306     }
307 
308     @Override
309     public Quantile<N> clone() {
310         return (Quantile<N>)super.clone();
311     }
312 
313 
314     static <N extends Number> Quantile<N> median() {
315         return new Quantile<>(0.5);
316     }
317 
318 }