Quantile.java
001 /*
002  * Java Genetic Algorithm Library (jenetics-1.6.0).
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 1.3 &mdash; <em>$Date: 2014-03-01 $</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      @deprecated Use {@link #getValue()} instead.
117      */
118     @Deprecated
119     public double getQuantile() {
120         return _q[2];
121     }
122 
123     /**
124      * Return the computed quantile value.
125      *
126      @return the quantile value.
127      */
128     public double getValue() {
129         return _q[2];
130     }
131 
132     @Override
133     public void accumulate(final N value) {
134         if (!_initialized) {
135             initialize(value.doubleValue());
136         else {
137             update(value.doubleValue());
138         }
139 
140         ++_samples;
141     }
142 
143 
144     private void initialize(double value) {
145         if (_n[00.0) {
146             _n[00.0;
147             _q[0= value;
148         else if (_n[1== 0.0) {
149             _n[11.0;
150             _q[1= value;
151         else if (_n[2== 0.0) {
152             _n[22.0;
153             _q[2= value;
154         else if (_n[3== 0.0) {
155             _n[33.0;
156             _q[3= value;
157         else if (_n[4== 0.0) {
158             _n[44.0;
159             _q[4= value;
160         }
161 
162         if (_n[4!= 0.0) {
163             Arrays.sort(_q);
164 
165             _nn[02.0*_quantile;
166             _nn[14.0*_quantile;
167             _nn[22.0*_quantile + 2.0;
168 
169             _dn[0= _quantile/2.0;
170             _dn[1= _quantile;
171             _dn[2(1.0 + _quantile)/2.0;
172 
173             _initialized = true;
174         }
175     }
176 
177     private void update(double value) {
178         assert (_initialized);
179 
180         // If min or max, handle as special case; otherwise, ...
181         if (_quantile == 0.0) {
182             if (value < _q[2]) {
183                 _q[2= value;
184             }
185         else if (_quantile == 1.0) {
186             if (value > _q[2]) {
187                 _q[2= value;
188             }
189         else {
190             // Increment marker locations and update min and max.
191             if (value < _q[0]) {
192                 ++_n[1]; ++_n[2]; ++_n[3]; ++_n[4]; _q[0= value;
193             else if (value < _q[1]) {
194                 ++_n[1]; ++_n[2]; ++_n[3]; ++_n[4];
195             else if (value < _q[2]) {
196                 ++_n[2]; ++_n[3]; ++_n[4];
197             else if (value < _q[3]) {
198                 ++_n[3]; ++_n[4];
199             else if (value < _q[4]) {
200                 ++_n[4];
201             else {
202                 ++_n[4]; _q[4= value;
203             }
204 
205             // Increment positions of markers k + 1
206             _nn[0+= _dn[0];
207             _nn[1+= _dn[1];
208             _nn[2+= _dn[2];
209 
210             // Adjust heights of markers 0 to 2 if necessary
211             double mm = _n[11.0;
212             double mp = _n[11.0;
213             if (_nn[0>= mp && _n[2> mp) {
214                 _q[1= qPlus(mp, _n[0], _n[1], _n[2], _q[0], _q[1], _q[2]);
215                 _n[1= mp;
216             else if (_nn[0<= mm && _n[0< mm) {
217                 _q[1= qMinus(mm, _n[0], _n[1], _n[2], _q[0], _q[1], _q[2]);
218                 _n[1= mm;
219             }
220 
221             mm = _n[21.0;
222             mp = _n[21.0;
223             if (_nn[1>= mp && _n[3> mp) {
224                 _q[2= qPlus(mp, _n[1], _n[2], _n[3], _q[1], _q[2], _q[3]);
225                 _n[2= mp;
226             else if (_nn[1<= mm && _n[1< mm) {
227                 _q[2= qMinus(mm, _n[1], _n[2], _n[3], _q[1], _q[2], _q[3]);
228                 _n[2= mm;
229             }
230 
231             mm = _n[31.0;
232             mp = _n[31.0;
233             if (_nn[2>= mp && _n[4> mp) {
234                 _q[3= qPlus(mp, _n[2], _n[3], _n[4], _q[2], _q[3], _q[4]);
235                 _n[3= mp;
236             else if (_nn[2<= mm && _n[2< mm) {
237                 _q[3= qMinus(mm, _n[2], _n[3], _n[4], _q[2], _q[3], _q[4]);
238                 _n[3= mm;
239             }
240         }
241     }
242 
243     private static double qPlus(
244         final double mp,
245         final double m0,
246         final double m1,
247         final double m2,
248         final double q0,
249         final double q1,
250         final double q2
251     ) {
252         double result = q1 +
253                     ((mp - m0)*(q2 - q1)/(m2 - m1+
254                     (m2 - mp)*(q1 - q0)/(m1 - m0))/(m2 - m0);
255 
256         if (result > q2) {
257             result = q1 + (q2 - q1)/(m2 - m1);
258         }
259 
260         return result;
261     }
262 
263     private static double qMinus(
264         final double mm,
265         final double m0,
266         final double m1,
267         final double m2,
268         final double q0,
269         final double q1,
270         final double q2
271     ) {
272         double result = q1 -
273                     ((mm - m0)*(q2 - q1)/(m2 - m1+
274                     (m2 - mm)*(q1 - q0)/(m1 - m0))/(m2 - m0);
275 
276         if (q0 > result) {
277             result = q1 + (q0 - q1)/(m0 - m1);
278         }
279 
280         return result;
281     }
282 
283     @Override
284     public int hashCode() {
285         return HashBuilder.of(getClass()).
286                 and(super.hashCode()).
287                 and(_quantile).
288                 and(_dn).
289                 and(_n).
290                 and(_nn).
291                 and(_q).value();
292     }
293 
294     @Override
295     public boolean equals(final Object obj) {
296         if (obj == this) {
297             return true;
298         }
299         if (obj == null || getClass() != obj.getClass()) {
300             return false;
301         }
302 
303         final Quantile<?> quantile = (Quantile<?>)obj;
304         return super.equals(obj&&
305                 eq(_quantile, quantile._quantile&&
306                 eq(_dn, quantile._dn&&
307                 eq(_n, quantile._n&&
308                 eq(_nn, quantile._nn&&
309                 eq(_q, quantile._q);
310     }
311 
312     @Override
313     public String toString() {
314         return format(
315             "%s[samples=%d, quantile=%f]",
316             getClass().getSimpleName(), getSamples(), getQuantile()
317         );
318     }
319 
320     @Override
321     public Quantile<N> clone() {
322         return (Quantile<N>)super.clone();
323     }
324 
325 
326     static <N extends Number> Quantile<N> median() {
327         return new Quantile<>(0.5);
328     }
329 
330 }