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