1 package hep.aida.bin;
2
3 import cern.colt.list.DoubleArrayList;
4 import cern.jet.stat.Descriptive;
5 /**
6 * Static and the same as its superclass, except that it can do more: Additionally computes moments of arbitrary integer order, harmonic mean, geometric mean, etc.
7 *
8 * Constructors need to be told what functionality is required for the given use case.
9 * Only maintains aggregate measures (incrementally) - the added elements themselves are not kept.
10 *
11 * @author wolfgang.hoschek@cern.ch
12 * @version 0.9, 03-Jul-99
13 */
14 public class MightyStaticBin1D extends StaticBin1D {
15 protected boolean hasSumOfLogarithms = false;
16 protected double sumOfLogarithms = 0.0; // Sum( Log(x[i]) )
17
18 protected boolean hasSumOfInversions = false;
19 protected double sumOfInversions = 0.0; // Sum( 1/x[i] )
20
21 protected double[] sumOfPowers = null; // Sum( x[i]^3 ) .. Sum( x[i]^max_k )
22 /**
23 * Constructs and returns an empty bin with limited functionality but good performance; equivalent to <tt>MightyStaticBin1D(false,false,4)</tt>.
24 */
MightyStaticBin1D()25 public MightyStaticBin1D() {
26 this(false, false, 4);
27 }
28 /**
29 * Constructs and returns an empty bin with the given capabilities.
30 *
31 * @param hasSumOfLogarithms Tells whether {@link #sumOfLogarithms()} can return meaningful results.
32 * Set this parameter to <tt>false</tt> if measures of sum of logarithms, geometric mean and product are not required.
33 * <p>
34 * @param hasSumOfInversions Tells whether {@link #sumOfInversions()} can return meaningful results.
35 * Set this parameter to <tt>false</tt> if measures of sum of inversions, harmonic mean and sumOfPowers(-1) are not required.
36 * <p>
37 * @param maxOrderForSumOfPowers The maximum order <tt>k</tt> for which {@link #sumOfPowers(int)} can return meaningful results.
38 * Set this parameter to at least 3 if the skew is required, to at least 4 if the kurtosis is required.
39 * In general, if moments are required set this parameter at least as large as the largest required moment.
40 * This method always substitutes <tt>Math.max(2,maxOrderForSumOfPowers)</tt> for the parameter passed in.
41 * Thus, <tt>sumOfPowers(0..2)</tt> always returns meaningful results.
42 *
43 * @see #hasSumOfPowers(int)
44 * @see #moment(int,double)
45 */
MightyStaticBin1D(boolean hasSumOfLogarithms, boolean hasSumOfInversions, int maxOrderForSumOfPowers)46 public MightyStaticBin1D(boolean hasSumOfLogarithms, boolean hasSumOfInversions, int maxOrderForSumOfPowers) {
47 setMaxOrderForSumOfPowers(maxOrderForSumOfPowers);
48 this.hasSumOfLogarithms = hasSumOfLogarithms;
49 this.hasSumOfInversions = hasSumOfInversions;
50 this.clear();
51 }
52 /**
53 * Adds the part of the specified list between indexes <tt>from</tt> (inclusive) and <tt>to</tt> (inclusive) to the receiver.
54 *
55 * @param list the list of which elements shall be added.
56 * @param from the index of the first element to be added (inclusive).
57 * @param to the index of the last element to be added (inclusive).
58 * @throws IndexOutOfBoundsException if <tt>list.size()>0 && (from<0 || from>to || to>=list.size())</tt>.
59 */
addAllOfFromTo(DoubleArrayList list, int from, int to)60 public synchronized void addAllOfFromTo(DoubleArrayList list, int from, int to) {
61 super.addAllOfFromTo(list, from, to);
62
63 if (this.sumOfPowers != null) {
64 //int max_k = this.min_k + this.sumOfPowers.length-1;
65 Descriptive.incrementalUpdateSumsOfPowers(list, from, to, 3, getMaxOrderForSumOfPowers(), this.sumOfPowers);
66 }
67
68 if (this.hasSumOfInversions) {
69 this.sumOfInversions += Descriptive.sumOfInversions(list, from, to);
70 }
71
72 if (this.hasSumOfLogarithms) {
73 this.sumOfLogarithms += Descriptive.sumOfLogarithms(list, from, to);
74 }
75 }
76 /**
77 * Resets the values of all measures.
78 */
clearAllMeasures()79 protected void clearAllMeasures() {
80 super.clearAllMeasures();
81
82 this.sumOfLogarithms = 0.0;
83 this.sumOfInversions = 0.0;
84
85 if (this.sumOfPowers != null) {
86 for (int i=this.sumOfPowers.length; --i >=0; ) {
87 this.sumOfPowers[i] = 0.0;
88 }
89 }
90 }
91 /**
92 * Returns a deep copy of the receiver.
93 *
94 * @return a deep copy of the receiver.
95 */
clone()96 public synchronized Object clone() {
97 MightyStaticBin1D clone = (MightyStaticBin1D) super.clone();
98 if (this.sumOfPowers != null) clone.sumOfPowers = (double[]) clone.sumOfPowers.clone();
99 return clone;
100 }
101 /**
102 * Computes the deviations from the receiver's measures to another bin's measures.
103 * @param other the other bin to compare with
104 * @return a summary of the deviations.
105 */
compareWith(AbstractBin1D other)106 public String compareWith(AbstractBin1D other) {
107 StringBuffer buf = new StringBuffer(super.compareWith(other));
108 if (other instanceof MightyStaticBin1D) {
109 MightyStaticBin1D m = (MightyStaticBin1D) other;
110 if (hasSumOfLogarithms() && m.hasSumOfLogarithms())
111 buf.append("geometric mean: "+relError(geometricMean(),m.geometricMean()) +" %\n");
112 if (hasSumOfInversions() && m.hasSumOfInversions())
113 buf.append("harmonic mean: "+relError(harmonicMean(),m.harmonicMean()) +" %\n");
114 if (hasSumOfPowers(3) && m.hasSumOfPowers(3))
115 buf.append("skew: "+relError(skew(),m.skew()) +" %\n");
116 if (hasSumOfPowers(4) && m.hasSumOfPowers(4))
117 buf.append("kurtosis: "+relError(kurtosis(),m.kurtosis()) +" %\n");
118 buf.append("\n");
119 }
120 return buf.toString();
121 }
122 /**
123 * Returns the geometric mean, which is <tt>Product( x[i] )<sup>1.0/size()</sup></tt>.
124 *
125 * This method tries to avoid overflows at the expense of an equivalent but somewhat inefficient definition:
126 * <tt>geoMean = exp( Sum( Log(x[i]) ) / size())</tt>.
127 * Note that for a geometric mean to be meaningful, the minimum of the data sequence must not be less or equal to zero.
128 * @return the geometric mean; <tt>Double.NaN</tt> if <tt>!hasSumOfLogarithms()</tt>.
129 */
geometricMean()130 public synchronized double geometricMean() {
131 return Descriptive.geometricMean(size(), sumOfLogarithms());
132 }
133 /**
134 * Returns the maximum order <tt>k</tt> for which sums of powers are retrievable, as specified upon instance construction.
135 * @see #hasSumOfPowers(int)
136 * @see #sumOfPowers(int)
137 */
getMaxOrderForSumOfPowers()138 public synchronized int getMaxOrderForSumOfPowers() {
139 /* order 0..2 is always recorded.
140 order 0 is size()
141 order 1 is sum()
142 order 2 is sum_xx()
143 */
144 if (this.sumOfPowers == null) return 2;
145
146 return 2 + this.sumOfPowers.length;
147 }
148 /**
149 * Returns the minimum order <tt>k</tt> for which sums of powers are retrievable, as specified upon instance construction.
150 * @see #hasSumOfPowers(int)
151 * @see #sumOfPowers(int)
152 */
getMinOrderForSumOfPowers()153 public synchronized int getMinOrderForSumOfPowers() {
154 int minOrder = 0;
155 if (hasSumOfInversions()) minOrder = -1;
156 return minOrder;
157 }
158 /**
159 * Returns the harmonic mean, which is <tt>size() / Sum( 1/x[i] )</tt>.
160 * Remember: If the receiver contains at least one element of <tt>0.0</tt>, the harmonic mean is <tt>0.0</tt>.
161 * @return the harmonic mean; <tt>Double.NaN</tt> if <tt>!hasSumOfInversions()</tt>.
162 * @see #hasSumOfInversions()
163 */
harmonicMean()164 public synchronized double harmonicMean() {
165 return Descriptive.harmonicMean(size(), sumOfInversions());
166 }
167 /**
168 * Returns whether <tt>sumOfInversions()</tt> can return meaningful results.
169 * @return <tt>false</tt> if the bin was constructed with insufficient parametrization, <tt>true</tt> otherwise.
170 * See the constructors for proper parametrization.
171 */
hasSumOfInversions()172 public boolean hasSumOfInversions() {
173 return this.hasSumOfInversions;
174 }
175 /**
176 * Tells whether <tt>sumOfLogarithms()</tt> can return meaningful results.
177 * @return <tt>false</tt> if the bin was constructed with insufficient parametrization, <tt>true</tt> otherwise.
178 * See the constructors for proper parametrization.
179 */
hasSumOfLogarithms()180 public boolean hasSumOfLogarithms() {
181 return this.hasSumOfLogarithms;
182 }
183 /**
184 * Tells whether <tt>sumOfPowers(k)</tt> can return meaningful results.
185 * Defined as <tt>hasSumOfPowers(k) <==> getMinOrderForSumOfPowers() <= k && k <= getMaxOrderForSumOfPowers()</tt>.
186 * A return value of <tt>true</tt> implies that <tt>hasSumOfPowers(k-1) .. hasSumOfPowers(0)</tt> will also return <tt>true</tt>.
187 * See the constructors for proper parametrization.
188 * <p>
189 * <b>Details</b>:
190 * <tt>hasSumOfPowers(0..2)</tt> will always yield <tt>true</tt>.
191 * <tt>hasSumOfPowers(-1) <==> hasSumOfInversions()</tt>.
192 *
193 * @return <tt>false</tt> if the bin was constructed with insufficient parametrization, <tt>true</tt> otherwise.
194 * @see #getMinOrderForSumOfPowers()
195 * @see #getMaxOrderForSumOfPowers()
196 */
hasSumOfPowers(int k)197 public boolean hasSumOfPowers(int k) {
198 return getMinOrderForSumOfPowers() <= k && k <= getMaxOrderForSumOfPowers();
199 }
200 /**
201 * Returns the kurtosis (aka excess), which is <tt>-3 + moment(4,mean()) / standardDeviation()<sup>4</sup></tt>.
202 * @return the kurtosis; <tt>Double.NaN</tt> if <tt>!hasSumOfPowers(4)</tt>.
203 * @see #hasSumOfPowers(int)
204 */
kurtosis()205 public synchronized double kurtosis() {
206 return Descriptive.kurtosis( moment(4,mean()), standardDeviation() );
207 }
208 /**
209 * Returns the moment of <tt>k</tt>-th order with value <tt>c</tt>,
210 * which is <tt>Sum( (x[i]-c)<sup>k</sup> ) / size()</tt>.
211 *
212 * @param k the order; must be greater than or equal to zero.
213 * @param c any number.
214 * @throws IllegalArgumentException if <tt>k < 0</tt>.
215 * @return <tt>Double.NaN</tt> if <tt>!hasSumOfPower(k)</tt>.
216 */
moment(int k, double c)217 public synchronized double moment(int k, double c) {
218 if (k<0) throw new IllegalArgumentException("k must be >= 0");
219 //checkOrder(k);
220 if (!hasSumOfPowers(k)) return Double.NaN;
221
222 int maxOrder = Math.min(k,getMaxOrderForSumOfPowers());
223 DoubleArrayList sumOfPows = new DoubleArrayList(maxOrder+1);
224 sumOfPows.add(size());
225 sumOfPows.add(sum());
226 sumOfPows.add(sumOfSquares());
227 for (int i=3; i<=maxOrder; i++) sumOfPows.add(sumOfPowers(i));
228
229 return Descriptive.moment(k, c, size(), sumOfPows.elements());
230 }
231 /**
232 * Returns the product, which is <tt>Prod( x[i] )</tt>.
233 * In other words: <tt>x[0]*x[1]*...*x[size()-1]</tt>.
234 * @return the product; <tt>Double.NaN</tt> if <tt>!hasSumOfLogarithms()</tt>.
235 * @see #hasSumOfLogarithms()
236 */
product()237 public double product() {
238 return Descriptive.product(size(), sumOfLogarithms());
239 }
240 /**
241 * Sets the range of orders in which sums of powers are to be computed.
242 * In other words, <tt>sumOfPower(k)</tt> will return <tt>Sum( x[i]^k )</tt> if <tt>min_k <= k <= max_k || 0 <= k <= 2</tt>
243 * and throw an exception otherwise.
244 * @see #isLegalOrder(int)
245 * @see #sumOfPowers(int)
246 * @see #getRangeForSumOfPowers()
247 */
setMaxOrderForSumOfPowers(int max_k)248 protected void setMaxOrderForSumOfPowers(int max_k) {
249 //if (max_k < ) throw new IllegalArgumentException();
250
251 if (max_k <=2) {
252 this.sumOfPowers = null;
253 }
254 else {
255 this.sumOfPowers = new double[max_k - 2];
256 }
257 }
258 /**
259 * Returns the skew, which is <tt>moment(3,mean()) / standardDeviation()<sup>3</sup></tt>.
260 * @return the skew; <tt>Double.NaN</tt> if <tt>!hasSumOfPowers(3)</tt>.
261 * @see #hasSumOfPowers(int)
262 */
skew()263 public synchronized double skew() {
264 return Descriptive.skew( moment(3,mean()), standardDeviation() );
265 }
266 /**
267 * Returns the sum of inversions, which is <tt>Sum( 1 / x[i] )</tt>.
268 * @return the sum of inversions; <tt>Double.NaN</tt> if <tt>!hasSumOfInversions()</tt>.
269 * @see #hasSumOfInversions()
270 */
sumOfInversions()271 public double sumOfInversions() {
272 if (! this.hasSumOfInversions) return Double.NaN;
273 //if (! this.hasSumOfInversions) throw new IllegalOperationException("You must specify upon instance construction that the sum of inversions shall be computed.");
274 return this.sumOfInversions;
275 }
276 /**
277 * Returns the sum of logarithms, which is <tt>Sum( Log(x[i]) )</tt>.
278 * @return the sum of logarithms; <tt>Double.NaN</tt> if <tt>!hasSumOfLogarithms()</tt>.
279 * @see #hasSumOfLogarithms()
280 */
sumOfLogarithms()281 public synchronized double sumOfLogarithms() {
282 if (! this.hasSumOfLogarithms) return Double.NaN;
283 //if (! this.hasSumOfLogarithms) throw new IllegalOperationException("You must specify upon instance construction that the sum of logarithms shall be computed.");
284 return this.sumOfLogarithms;
285 }
286 /**
287 * Returns the <tt>k-th</tt> order sum of powers, which is <tt>Sum( x[i]<sup>k</sup> )</tt>.
288 * @param k the order of the powers.
289 * @return the sum of powers; <tt>Double.NaN</tt> if <tt>!hasSumOfPowers(k)</tt>.
290 * @see #hasSumOfPowers(int)
291 */
sumOfPowers(int k)292 public synchronized double sumOfPowers(int k) {
293 if (!hasSumOfPowers(k)) return Double.NaN;
294 //checkOrder(k);
295 if (k == -1) return sumOfInversions();
296 if (k == 0) return size();
297 if (k == 1) return sum();
298 if (k == 2) return sumOfSquares();
299
300 return this.sumOfPowers[k-3];
301 }
302 /**
303 * Returns a String representation of the receiver.
304 */
toString()305 public synchronized String toString() {
306 StringBuffer buf = new StringBuffer(super.toString());
307
308 if (hasSumOfLogarithms()) {
309 buf.append("Geometric mean: "+geometricMean());
310 buf.append("\nProduct: "+product()+"\n");
311 }
312
313 if (hasSumOfInversions()) {
314 buf.append("Harmonic mean: "+harmonicMean());
315 buf.append("\nSum of inversions: "+sumOfInversions()+"\n");
316 }
317
318 int maxOrder = getMaxOrderForSumOfPowers();
319 int maxPrintOrder = Math.min(6,maxOrder); // don't print tons of measures
320 if (maxOrder>2) {
321 if (maxOrder>=3) {
322 buf.append("Skew: "+skew()+"\n");
323 }
324 if (maxOrder>=4) {
325 buf.append("Kurtosis: "+kurtosis()+"\n");
326 }
327 for (int i=3; i<=maxPrintOrder; i++) {
328 buf.append("Sum of powers("+i+"): "+sumOfPowers(i)+"\n");
329 }
330 for (int k=0; k<=maxPrintOrder; k++) {
331 buf.append("Moment("+k+",0): "+moment(k,0)+"\n");
332 }
333 for (int k=0; k<=maxPrintOrder; k++) {
334 buf.append("Moment("+k+",mean()): "+moment(k,mean())+"\n");
335 }
336 }
337 return buf.toString();
338 }
339 /**
340 * @throws IllegalOperationException if <tt>! isLegalOrder(k)</tt>.
341 */
xcheckOrder(int k)342 protected void xcheckOrder(int k) {
343 //if (! isLegalOrder(k)) return Double.NaN;
344 //if (! xisLegalOrder(k)) throw new IllegalOperationException("Illegal order of sum of powers: k="+k+". Upon instance construction legal range was fixed to be "+getMinOrderForSumOfPowers()+" <= k <= "+getMaxOrderForSumOfPowers());
345 }
346 /**
347 * Returns whether two bins are equal;
348 * They are equal if the other object is of the same class or a subclass of this class and both have the same size, minimum, maximum, sum, sumOfSquares, sumOfInversions and sumOfLogarithms.
349 */
xequals(Object object)350 protected boolean xequals(Object object) {
351 if (!(object instanceof MightyStaticBin1D)) return false;
352 MightyStaticBin1D other = (MightyStaticBin1D) object;
353 return super.equals(other) && sumOfInversions()==other.sumOfInversions() && sumOfLogarithms()==other.sumOfLogarithms();
354 }
355 /**
356 * Tells whether <tt>sumOfPowers(fromK) .. sumOfPowers(toK)</tt> can return meaningful results.
357 * @return <tt>false</tt> if the bin was constructed with insufficient parametrization, <tt>true</tt> otherwise.
358 * See the constructors for proper parametrization.
359 * @throws IllegalArgumentException if <tt>fromK > toK</tt>.
360 */
xhasSumOfPowers(int fromK, int toK)361 protected boolean xhasSumOfPowers(int fromK, int toK) {
362 if (fromK > toK) throw new IllegalArgumentException("fromK must be less or equal to toK");
363 return getMinOrderForSumOfPowers() <= fromK && toK <= getMaxOrderForSumOfPowers();
364 }
365 /**
366 * Returns <tt>getMinOrderForSumOfPowers() <= k && k <= getMaxOrderForSumOfPowers()</tt>.
367 */
xisLegalOrder(int k)368 protected synchronized boolean xisLegalOrder(int k) {
369 return getMinOrderForSumOfPowers() <= k && k <= getMaxOrderForSumOfPowers();
370 }
371 }
372