001/*-
002 *******************************************************************************
003 * Copyright (c) 2011, 2016 Diamond Light Source Ltd.
004 * All rights reserved. This program and the accompanying materials
005 * are made available under the terms of the Eclipse Public License v1.0
006 * which accompanies this distribution, and is available at
007 * http://www.eclipse.org/legal/epl-v10.html
008 *
009 * Contributors:
010 *    Peter Chang - initial API and implementation and/or initial documentation
011 *******************************************************************************/
012
013package org.eclipse.january.dataset;
014
015import java.lang.ref.SoftReference;
016import java.util.ArrayList;
017import java.util.Collections;
018import java.util.HashMap;
019import java.util.List;
020import java.util.Map;
021import java.util.TreeMap;
022
023import org.apache.commons.math3.complex.Complex;
024import org.apache.commons.math3.stat.descriptive.moment.Kurtosis;
025import org.apache.commons.math3.stat.descriptive.moment.Skewness;
026import org.eclipse.january.metadata.Dirtiable;
027import org.eclipse.january.metadata.MetadataType;
028
029
030/**
031 * Statistics class
032 * 
033 * TODO Where is mode? http://en.wikipedia.org/wiki/Mode_(statistics)
034 * 
035 */
036public class Stats {
037
038        private static class ReferencedDataset extends SoftReference<Dataset> {
039                public ReferencedDataset(Dataset d) {
040                        super(d);
041                }
042        }
043
044        private static class QStatisticsImpl<T> implements MetadataType {
045                private static final long serialVersionUID = -3296671666463190388L;
046                final static Double Q1 = 0.25;
047                final static Double Q2 = 0.5;
048                final static Double Q3 = 0.75;
049                Map<Double, T> qmap = new HashMap<Double, T>();
050                transient Map<Integer, Map<Double, ReferencedDataset>> aqmap = new HashMap<Integer, Map<Double, ReferencedDataset>>();
051                transient ReferencedDataset s; // store 0th element
052                transient Map<Integer, ReferencedDataset> smap = new HashMap<>();
053
054                @Dirtiable
055                private boolean isDirty = true;
056
057                @Override
058                public QStatisticsImpl<T> clone() {
059                        return new QStatisticsImpl<T>(this);
060                }
061
062                public QStatisticsImpl() {
063                }
064
065                private QStatisticsImpl(QStatisticsImpl<T> qstats) {
066                        if (qstats.s != null && qstats.s.get() != null) {
067                                s = new ReferencedDataset(qstats.s.get().getView(false));
068                        }
069                        qmap.putAll(qstats.qmap);
070                        for (Integer i : qstats.aqmap.keySet()) {
071                                aqmap.put(i, new HashMap<>(qstats.aqmap.get(i)));
072                        }
073                        smap.putAll(qstats.smap);
074                        isDirty = qstats.isDirty;
075                }
076
077                public void setQuantile(double q, T v) {
078                        qmap.put(q, v);
079                }
080
081                public T getQuantile(double q) {
082                        return qmap.get(q);
083                }
084
085                private Map<Double, ReferencedDataset> getMap(int axis) {
086                        Map<Double, ReferencedDataset> qm = aqmap.get(axis);
087                        if (qm == null) {
088                                qm = new HashMap<>();
089                                aqmap.put(axis, qm);
090                        }
091                        return qm;
092                }
093
094                public void setQuantile(int axis, double q, Dataset v) {
095                        Map<Double, ReferencedDataset> qm = getMap(axis);
096                        qm.put(q, new ReferencedDataset(v));
097                }
098
099                public Dataset getQuantile(int axis, double q) {
100                        Map<Double, ReferencedDataset> qm = getMap(axis);
101                        ReferencedDataset rd = qm.get(q);
102                        return rd == null ? null : rd.get();
103                }
104
105                Dataset getSortedDataset(int axis) {
106                        return smap.containsKey(axis) ? smap.get(axis).get() : null;
107                }
108
109                void setSortedDataset(int axis, Dataset v) {
110                        smap.put(axis, new ReferencedDataset(v));
111                }
112        }
113
114        // calculates statistics and returns sorted dataset (0th element if compound)
115        private static QStatisticsImpl<?> calcQuartileStats(final Dataset a) {
116                Dataset s = null;
117                final int is = a.getElementsPerItem();
118
119                if (is == 1) {
120                        s = DatasetUtils.sort(a);
121
122                        QStatisticsImpl<Double> qstats = new QStatisticsImpl<Double>();
123
124                        qstats.setQuantile(QStatisticsImpl.Q1, pQuantile(s, QStatisticsImpl.Q1));
125                        qstats.setQuantile(QStatisticsImpl.Q2, pQuantile(s, QStatisticsImpl.Q2));
126                        qstats.setQuantile(QStatisticsImpl.Q3, pQuantile(s, QStatisticsImpl.Q3));
127                        qstats.s = new ReferencedDataset(s);
128                        return qstats;
129                }
130
131                QStatisticsImpl<double[]> qstats = new QStatisticsImpl<double[]>();
132
133                Dataset w = DatasetFactory.zeros(1, a.getClass(), a.getShapeRef());
134                double[] q1 = new double[is];
135                double[] q2 = new double[is];
136                double[] q3 = new double[is];
137                qstats.setQuantile(QStatisticsImpl.Q1, q1);
138                qstats.setQuantile(QStatisticsImpl.Q2, q2);
139                qstats.setQuantile(QStatisticsImpl.Q3, q3);
140                for (int j = 0; j < is; j++) {
141                        ((CompoundDataset) a).copyElements(w, j);
142                        w.sort(null);
143
144                        q1[j] = pQuantile(w, QStatisticsImpl.Q1);
145                        q2[j] = pQuantile(w, QStatisticsImpl.Q2);
146                        q3[j] = pQuantile(w, QStatisticsImpl.Q3);
147                        if (j == 0)
148                                s = w.clone();
149                }
150                qstats.s = new ReferencedDataset(s);
151
152                return qstats;
153        }
154
155        static private QStatisticsImpl<?> getQStatistics(final Dataset a) {
156                QStatisticsImpl<?> m = a.getFirstMetadata(QStatisticsImpl.class);
157                if (m == null || m.isDirty) {
158                        m = calcQuartileStats(a);
159                        a.setMetadata(m);
160                }
161                return m;
162        }
163
164        static private QStatisticsImpl<?> getQStatistics(final Dataset a, int axis) {
165                axis = a.checkAxis(axis);
166                final int is = a.getElementsPerItem();
167                QStatisticsImpl<?> qstats = a.getFirstMetadata(QStatisticsImpl.class);
168
169                if (qstats == null || qstats.isDirty) {
170                        if (is == 1) {
171                                qstats = new QStatisticsImpl<Double>();
172                        } else {
173                                qstats = new QStatisticsImpl<double[]>();
174                        }
175                        a.setMetadata(qstats);
176                }
177
178                if (qstats.getQuantile(axis, QStatisticsImpl.Q2) == null) {
179                        if (is == 1) {
180                                Dataset s = DatasetUtils.sort(a, axis);
181
182                                qstats.setQuantile(axis, QStatisticsImpl.Q1, pQuantile(s, axis, QStatisticsImpl.Q1));
183                                qstats.setQuantile(axis, QStatisticsImpl.Q2, pQuantile(s, axis, QStatisticsImpl.Q2));
184                                qstats.setQuantile(axis, QStatisticsImpl.Q3, pQuantile(s, axis, QStatisticsImpl.Q3));
185                                qstats.setSortedDataset(axis, s);
186                        } else {
187                                Dataset w = DatasetFactory.zeros(1, a.getClass(), a.getShapeRef());
188                                CompoundDoubleDataset q1 = null, q2 = null, q3 = null;
189                                for (int j = 0; j < is; j++) {
190                                        ((CompoundDataset) a).copyElements(w, j);
191                                        w.sort(axis);
192        
193                                        final Dataset c = pQuantile(w, axis, QStatisticsImpl.Q1);
194                                        if (j == 0) {
195                                                q1 = DatasetFactory.zeros(is, CompoundDoubleDataset.class, c.getShapeRef());
196                                                q2 = DatasetFactory.zeros(is, CompoundDoubleDataset.class, c.getShapeRef());
197                                                q3 = DatasetFactory.zeros(is, CompoundDoubleDataset.class, c.getShapeRef());
198                                        }
199                                        q1.setElements(c, j);
200        
201                                        q2.setElements(pQuantile(w, axis, QStatisticsImpl.Q2), j);
202        
203                                        q3.setElements(pQuantile(w, axis, QStatisticsImpl.Q3), j);
204                                }
205                                qstats.setQuantile(axis, QStatisticsImpl.Q1, q1);
206                                qstats.setQuantile(axis, QStatisticsImpl.Q2, q2);
207                                qstats.setQuantile(axis, QStatisticsImpl.Q3, q3);
208                        }
209                }
210
211                return qstats;
212        }
213
214        // process a sorted dataset
215        private static double pQuantile(final Dataset s, final double q) {
216                double f = (s.getSize() - 1) * q; // fraction of sample number
217                if (f < 0)
218                        return Double.NaN;
219                int qpt = (int) Math.floor(f); // quantile point
220                f -= qpt;
221
222                double quantile = s.getElementDoubleAbs(qpt);
223                if (f > 0) {
224                        quantile = (1-f)*quantile + f*s.getElementDoubleAbs(qpt+1);
225                }
226                return quantile;
227        }
228
229        // process a sorted dataset and returns a double or compound double dataset
230        private static Dataset pQuantile(final Dataset s, final int axis, final double q) {
231                final int rank = s.getRank();
232                final int is = s.getElementsPerItem();
233
234                int[] oshape = s.getShape();
235
236                double f = (oshape[axis] - 1) * q; // fraction of sample number
237                int qpt = (int) Math.floor(f); // quantile point
238                f -= qpt;
239
240                oshape[axis] = 1;
241                int[] qshape = ShapeUtils.squeezeShape(oshape, false);
242                Dataset qds = DatasetFactory.zeros(is, CompoundDoubleDataset.class, qshape);
243
244                IndexIterator qiter = qds.getIterator(true);
245                int[] qpos = qiter.getPos();
246                int[] spos = oshape;
247
248                while (qiter.hasNext()) {
249                        int i = 0;
250                        for (; i < axis; i++) {
251                                spos[i] = qpos[i];
252                        }
253                        spos[i++] = qpt;
254                        for (; i < rank; i++) {
255                                spos[i] = qpos[i-1];
256                        }
257
258                        Object obj = s.getObject(spos);
259                        qds.set(obj, qpos);
260                }
261
262                if (f > 0) {
263                        qiter = qds.getIterator(true);
264                        qpos = qiter.getPos();
265                        qpt++;
266                        Dataset rds = DatasetFactory.zeros(is, CompoundDoubleDataset.class, qshape);
267                        
268                        while (qiter.hasNext()) {
269                                int i = 0;
270                                for (; i < axis; i++) {
271                                        spos[i] = qpos[i];
272                                }
273                                spos[i++] = qpt;
274                                for (; i < rank; i++) {
275                                        spos[i] = qpos[i-1];
276                                }
277
278                                Object obj = s.getObject(spos);
279                                rds.set(obj, qpos);
280                        }
281                        rds.imultiply(f);
282                        qds.imultiply(1-f);
283                        qds.iadd(rds);
284                }
285
286                return qds;
287        }
288
289        /**
290         * Calculate quantile of dataset which is defined as the inverse of the cumulative distribution function (CDF)
291         * @param a
292         * @param q
293         * @return point at which CDF has value q
294         */
295        @SuppressWarnings("unchecked")
296        public static double quantile(final Dataset a, final double q) {
297                if (q < 0 || q > 1) {
298                        throw new IllegalArgumentException("Quantile requested is outside [0,1]");
299                }
300                QStatisticsImpl<Double> qs = (QStatisticsImpl<Double>) getQStatistics(a);
301                Double qv = qs.getQuantile(q);
302                if (qv == null) {
303                        qv = pQuantile(qs.s.get(), q);
304                        qs.setQuantile(q, qv);
305                }
306                return qv;
307        }
308
309        /**
310         * Calculate quantiles of dataset which is defined as the inverse of the cumulative distribution function (CDF)
311         * @param a
312         * @param values
313         * @return points at which CDF has given values
314         */
315        @SuppressWarnings("unchecked")
316        public static double[] quantile(final Dataset a, final double... values) {
317                final double[] points  = new double[values.length];
318                QStatisticsImpl<Double> qs = (QStatisticsImpl<Double>) getQStatistics(a);
319                for (int i = 0; i < points.length; i++) {
320                        final double q = values[i];
321                        if (q < 0 || q > 1) {
322                                throw new IllegalArgumentException("Quantile requested is outside [0,1]");
323                        }
324                        Double qv = qs.getQuantile(q);
325                        if (qv == null) {
326                                qv = pQuantile(qs.s.get(), q);
327                                qs.setQuantile(q, qv);
328                        }
329                        points[i] = qv;
330                }
331
332                return points;
333        }
334
335        /**
336         * Calculate quantiles of dataset which is defined as the inverse of the cumulative distribution function (CDF)
337         * @param a
338         * @param axis
339         * @param values
340         * @return points at which CDF has given values
341         */
342        @SuppressWarnings({ "unchecked" })
343        public static Dataset[] quantile(final Dataset a, final int axis, final double... values) {
344                final Dataset[] points  = new Dataset[values.length];
345                final int is = a.getElementsPerItem();
346
347                if (is == 1) {
348                        QStatisticsImpl<Double> qs = (QStatisticsImpl<Double>) getQStatistics(a, axis);
349                        for (int i = 0; i < points.length; i++) {
350                                final double q = values[i];
351                                if (q < 0 || q > 1) {
352                                        throw new IllegalArgumentException("Quantile requested is outside [0,1]");
353                                }
354                                Dataset qv = qs.getQuantile(axis, q);
355                                if (qv == null) {
356                                        qv = pQuantile(qs.getSortedDataset(axis), axis, q);
357                                        qs.setQuantile(axis, q, qv);
358                                }
359                                points[i] = qv;
360                        }
361                } else {
362                        QStatisticsImpl<double[]> qs = (QStatisticsImpl<double[]>) getQStatistics(a);
363                        Dataset w = DatasetFactory.zeros(1, a.getClass(), a.getShapeRef());
364                        for (int j = 0; j < is; j++) {
365                                boolean copied = false;
366
367                                for (int i = 0; i < points.length; i++) {
368                                        final double q = values[i];
369                                        if (q < 0 || q > 1) {
370                                                throw new IllegalArgumentException("Quantile requested is outside [0,1]");
371                                        }
372                                        Dataset qv = qs.getQuantile(axis, q);
373                                        if (qv == null) {
374                                                if (!copied) {
375                                                        copied = true;
376                                                        ((CompoundDataset) a).copyElements(w, j);
377                                                        w.sort(axis);
378                                                }
379                                                qv = pQuantile(w, axis, q);
380                                                qs.setQuantile(axis, q, qv);
381                                                if (j == 0) {
382                                                        points[i] = DatasetFactory.zeros(is, qv.getClass(), qv.getShapeRef());
383                                                }
384                                                ((CompoundDoubleDataset) points[i]).setElements(qv, j);
385                                        }
386                                }
387                        }
388                }
389
390                return points;
391        }
392
393        /**
394         * @param a dataset
395         * @param axis
396         * @return median
397         */
398        public static Dataset median(final Dataset a, final int axis) {
399                return getQStatistics(a, axis).getQuantile(axis, QStatisticsImpl.Q2);
400        }
401
402        /**
403         * @param a dataset
404         * @return median
405         */
406        public static Object median(final Dataset a) {
407                return getQStatistics(a).getQuantile(QStatisticsImpl.Q2);
408        }
409
410        /**
411         * Interquartile range: Q3 - Q1
412         * @param a
413         * @return range
414         */
415        @SuppressWarnings("unchecked")
416        public static Object iqr(final Dataset a) {
417                final int is = a.getElementsPerItem();
418                if (is == 1) {
419                        QStatisticsImpl<Double> qs = (QStatisticsImpl<Double>) getQStatistics(a);
420                        return qs.getQuantile(QStatisticsImpl.Q3) - qs.getQuantile(QStatisticsImpl.Q1);
421                }
422
423                QStatisticsImpl<double[]> qs = (QStatisticsImpl<double[]>) getQStatistics(a);
424                double[] q1 = (double[]) qs.getQuantile(QStatisticsImpl.Q1);
425                double[] q3 = ((double[]) qs.getQuantile(QStatisticsImpl.Q3)).clone();
426                for (int j = 0; j < is; j++) {
427                        q3[j] -= q1[j];
428                }
429                return q3;
430        }
431
432        /**
433         * Interquartile range: Q3 - Q1
434         * @param a
435         * @param axis
436         * @return range
437         */
438        public static Dataset iqr(final Dataset a, final int axis) {
439                QStatisticsImpl<?> qs = getQStatistics(a, axis);
440                Dataset q3 = qs.getQuantile(axis, QStatisticsImpl.Q3);
441
442                return Maths.subtract(q3, qs.getQuantile(axis, QStatisticsImpl.Q1));
443        }
444
445        static private HigherStatisticsImpl<?> getHigherStatistic(final Dataset a, boolean[] ignoreInvalids) {
446                boolean ignoreNaNs = false;
447                boolean ignoreInfs = false;
448                if (a.hasFloatingPointElements()) {
449                        ignoreNaNs = ignoreInvalids != null && ignoreInvalids.length > 0 ? ignoreInvalids[0] : false;
450                        ignoreInfs = ignoreInvalids != null && ignoreInvalids.length > 1 ? ignoreInvalids[1] : ignoreNaNs;
451                }
452
453                HigherStatisticsImpl<?> stats = a.getFirstMetadata(HigherStatisticsImpl.class);
454                if (stats == null || stats.isDirty) {
455                        stats = calculateHigherMoments(a, ignoreNaNs, ignoreInfs);
456                        a.setMetadata(stats);
457                }
458        
459                return stats;
460        }
461
462        static private HigherStatisticsImpl<?> getHigherStatistic(final Dataset a, final boolean[] ignoreInvalids, int axis) {
463                boolean ignoreNaNs = false;
464                boolean ignoreInfs = false;
465                if (a.hasFloatingPointElements()) {
466                        ignoreNaNs = ignoreInvalids != null && ignoreInvalids.length > 0 ? ignoreInvalids[0] : false;
467                        ignoreInfs = ignoreInvalids != null && ignoreInvalids.length > 1 ? ignoreInvalids[1] : ignoreNaNs;
468                }
469
470                axis = a.checkAxis(axis);
471        
472                HigherStatisticsImpl<?> stats = a.getFirstMetadata(HigherStatisticsImpl.class);
473                if (stats == null || stats.getSkewness(axis) == null || stats.isDirty) {
474                        stats = calculateHigherMoments(a, ignoreNaNs, ignoreInfs, axis);
475                        a.setMetadata(stats);
476                }
477        
478                return stats;
479        }
480
481        private static class HigherStatisticsImpl<T> implements MetadataType {
482                private static final long serialVersionUID = -6587974784104116792L;
483                T skewness;
484                T kurtosis;
485                transient Map<Integer, ReferencedDataset> smap = new HashMap<>();
486                transient Map<Integer, ReferencedDataset> kmap = new HashMap<>();
487
488                @Dirtiable
489                private boolean isDirty = true;
490
491                @Override
492                public HigherStatisticsImpl<T> clone() {
493                        return new HigherStatisticsImpl<>(this);
494                }
495
496                public HigherStatisticsImpl() {
497                }
498
499                private HigherStatisticsImpl(HigherStatisticsImpl<T> hstats) {
500                        skewness = hstats.skewness;
501                        kurtosis = hstats.kurtosis;
502                        smap.putAll(hstats.smap);
503                        kmap.putAll(hstats.kmap);
504                        isDirty = hstats.isDirty;
505                }
506
507//              public void setSkewness(T skewness) {
508//                      this.skewness = skewness;
509//              }
510//
511//              public void setKurtosis(T kurtosis) {
512//                      this.kurtosis = kurtosis;
513//              }
514//
515//              public T getSkewness() {
516//                      return skewness;
517//              }
518//
519//              public T getKurtosis() {
520//                      return kurtosis;
521//              }
522
523                public Dataset getSkewness(int axis) {
524                        ReferencedDataset rd = smap.get(axis);
525                        return rd == null ? null : rd.get();
526                }
527
528                public Dataset getKurtosis(int axis) {
529                        ReferencedDataset rd = kmap.get(axis);
530                        return rd == null ? null : rd.get();
531                }
532
533                public void setSkewness(int axis, Dataset s) {
534                        smap.put(axis, new ReferencedDataset(s));
535                }
536
537                public void setKurtosis(int axis, Dataset k) {
538                        kmap.put(axis, new ReferencedDataset(k));
539                }
540        }
541
542        static private HigherStatisticsImpl<?> calculateHigherMoments(final Dataset a, final boolean ignoreNaNs, final boolean ignoreInfs) {
543                final int is = a.getElementsPerItem();
544                final IndexIterator iter = a.getIterator();
545
546                if (is == 1) {
547                        Skewness s = new Skewness();
548                        Kurtosis k = new Kurtosis();
549                        if (ignoreNaNs) {
550                                while (iter.hasNext()) {
551                                        final double x = a.getElementDoubleAbs(iter.index);
552                                        if (Double.isNaN(x))
553                                                continue;
554                                        s.increment(x);
555                                        k.increment(x);
556                                }
557                        } else {
558                                while (iter.hasNext()) {
559                                        final double x = a.getElementDoubleAbs(iter.index);
560                                        s.increment(x);
561                                        k.increment(x);
562                                }
563                        }
564
565                        HigherStatisticsImpl<Double> stats = new HigherStatisticsImpl<Double>();
566                        stats.skewness = s.getResult();
567                        stats.kurtosis = k.getResult();
568                        return stats;
569                }
570                final Skewness[] s = new Skewness[is];
571                final Kurtosis[] k = new Kurtosis[is];
572
573                for (int j = 0; j < is; j++) {
574                        s[j] = new Skewness();
575                        k[j] = new Kurtosis();
576                }
577                if (ignoreNaNs) {
578                        while (iter.hasNext()) {
579                                boolean skip = false;
580                                for (int j = 0; j < is; j++) {
581                                        if (Double.isNaN(a.getElementDoubleAbs(iter.index + j))) {
582                                                skip = true;
583                                                break;
584                                        }
585                                }
586                                if (!skip)
587                                        for (int j = 0; j < is; j++) {
588                                                final double val = a.getElementDoubleAbs(iter.index + j);
589                                                s[j].increment(val);
590                                                k[j].increment(val);
591                                        }
592                        }
593                } else {
594                        while (iter.hasNext()) {
595                                for (int j = 0; j < is; j++) {
596                                        final double val = a.getElementDoubleAbs(iter.index + j);
597                                        s[j].increment(val);
598                                        k[j].increment(val);
599                                }
600                        }
601                }
602                final double[] ts = new double[is];
603                final double[] tk = new double[is];
604                for (int j = 0; j < is; j++) {
605                        ts[j] = s[j].getResult();
606                        tk[j] = k[j].getResult();
607                }
608
609                HigherStatisticsImpl<double[]> stats = new HigherStatisticsImpl<double[]>();
610                stats.skewness = ts;
611                stats.kurtosis = tk;
612                return stats;
613        }
614
615        static private HigherStatisticsImpl<?> calculateHigherMoments(final Dataset a, final boolean ignoreNaNs, final boolean ignoreInfs, final int axis) {
616                final int rank = a.getRank();
617                final int is = a.getElementsPerItem();
618                final int[] oshape = a.getShape();
619                final int alen = oshape[axis];
620                oshape[axis] = 1;
621        
622                final int[] nshape = ShapeUtils.squeezeShape(oshape, false);
623                final Dataset sk;
624                final Dataset ku;
625                HigherStatisticsImpl<?> stats = null;
626        
627                if (is == 1) {
628                        if (stats == null) {
629                                stats = new HigherStatisticsImpl<Double>();
630                                a.setMetadata(stats);
631                        }
632                        sk = DatasetFactory.zeros(DoubleDataset.class, nshape);
633                        ku = DatasetFactory.zeros(DoubleDataset.class, nshape);
634                        final IndexIterator qiter = sk.getIterator(true);
635                        final int[] qpos = qiter.getPos();
636                        final int[] spos = oshape;
637        
638                        while (qiter.hasNext()) {
639                                int i = 0;
640                                for (; i < axis; i++) {
641                                        spos[i] = qpos[i];
642                                }
643                                spos[i++] = 0;
644                                for (; i < rank; i++) {
645                                        spos[i] = qpos[i - 1];
646                                }
647        
648                                Skewness s = new Skewness();
649                                Kurtosis k = new Kurtosis();
650                                if (ignoreNaNs) {
651                                        for (int j = 0; j < alen; j++) {
652                                                spos[axis] = j;
653                                                final double val = a.getDouble(spos);
654                                                if (Double.isNaN(val))
655                                                        continue;
656        
657                                                s.increment(val);
658                                                k.increment(val);
659                                        }
660                                } else {
661                                        for (int j = 0; j < alen; j++) {
662                                                spos[axis] = j;
663                                                final double val = a.getDouble(spos);
664                                                s.increment(val);
665                                                k.increment(val);
666                                        }
667                                }
668                                sk.set(s.getResult(), qpos);
669                                ku.set(k.getResult(), qpos);
670                        }
671                } else {
672                        if (stats == null) {
673                                stats = new HigherStatisticsImpl<double[]>();
674                                a.setMetadata(stats);
675                        }
676                        sk = DatasetFactory.zeros(is, CompoundDoubleDataset.class, nshape);
677                        ku = DatasetFactory.zeros(is, CompoundDoubleDataset.class, nshape);
678                        final IndexIterator qiter = sk.getIterator(true);
679                        final int[] qpos = qiter.getPos();
680                        final int[] spos = oshape;
681                        final Skewness[] s = new Skewness[is];
682                        final Kurtosis[] k = new Kurtosis[is];
683                        final double[] ts = new double[is];
684                        final double[] tk = new double[is];
685        
686                        for (int j = 0; j < is; j++) {
687                                s[j] = new Skewness();
688                                k[j] = new Kurtosis();
689                        }
690                        while (qiter.hasNext()) {
691                                int i = 0;
692                                for (; i < axis; i++) {
693                                        spos[i] = qpos[i];
694                                }
695                                spos[i++] = 0;
696                                for (; i < rank; i++) {
697                                        spos[i] = qpos[i-1];
698                                }
699        
700                                for (int j = 0; j < is; j++) {
701                                        s[j].clear();
702                                        k[j].clear();
703                                }
704                                int index = a.get1DIndex(spos);
705                                if (ignoreNaNs) {
706                                        boolean skip = false;
707                                        for (int j = 0; j < is; j++) {
708                                                if (Double.isNaN(a.getElementDoubleAbs(index + j))) {
709                                                        skip = true;
710                                                        break;
711                                                }
712                                        }
713                                        if (!skip)
714                                                for (int j = 0; j < is; j++) {
715                                                        final double val = a.getElementDoubleAbs(index + j);
716                                                        s[j].increment(val);
717                                                        k[j].increment(val);
718                                                }
719                                } else {
720                                        for (int j = 0; j < is; j++) {
721                                                final double val = a.getElementDoubleAbs(index + j);
722                                                s[j].increment(val);
723                                                k[j].increment(val);
724                                        }
725                                }
726                                for (int j = 0; j < is; j++) {
727                                        ts[j] = s[j].getResult(); 
728                                        tk[j] = k[j].getResult(); 
729                                }
730                                sk.set(ts, qpos);
731                                ku.set(tk, qpos);
732                        }
733                }
734
735                stats.setSkewness(axis, sk);
736                stats.setKurtosis(axis, ku);
737                return stats;
738        }
739
740        /**
741         * @param a dataset
742         * @param ignoreInvalids see {@link IDataset#max(boolean...)}
743         * @return skewness
744         * @since 2.0
745         */
746        public static Object skewness(final Dataset a, final boolean... ignoreInvalids) {
747                return getHigherStatistic(a, ignoreInvalids).skewness;
748        }
749
750        /**
751         * @param a dataset
752         * @param ignoreInvalids see {@link IDataset#max(boolean...)}
753         * @return kurtosis
754         * @since 2.0
755         */
756        public static Object kurtosis(final Dataset a, final boolean... ignoreInvalids) {
757                return getHigherStatistic(a, ignoreInvalids).kurtosis;
758        }
759
760        /**
761         * @param a dataset
762         * @param axis
763         * @param ignoreInvalids see {@link Dataset#max(int, boolean...)}
764         * @return skewness along axis in dataset
765         * @since 2.0
766         */
767        public static Dataset skewness(final Dataset a, final int axis, final boolean... ignoreInvalids) {
768                return getHigherStatistic(a, ignoreInvalids, axis).getSkewness(axis);
769        }
770
771        /**
772         * @param a dataset
773         * @param axis
774         * @param ignoreInvalids see {@link Dataset#max(int, boolean...)}
775         * @return kurtosis along axis in dataset
776         * @since 2.0
777         */
778        public static Dataset kurtosis(final Dataset a, final int axis, final boolean... ignoreInvalids) {
779                return getHigherStatistic(a, ignoreInvalids, axis).getKurtosis(axis);
780        }
781
782        /**
783         * @param a dataset
784         * @param ignoreInvalids see {@link IDataset#max(boolean...)}
785         * @return sum of dataset
786         * @since 2.0
787         */
788        public static Object sum(final Dataset a, final boolean... ignoreInvalids) {
789                return a.sum(ignoreInvalids);
790        }
791
792        /**
793         * @param a dataset
794         * @param dtype
795         * @param ignoreInvalids see {@link IDataset#max(boolean...)}
796         * @return typed sum of dataset
797         * @since 2.0
798         */
799        public static Object typedSum(final Dataset a, int dtype, final boolean... ignoreInvalids) {
800                if (a.isComplex()) {
801                        Complex m = (Complex) a.sum(ignoreInvalids);
802                        return m;
803                } else if (a instanceof CompoundDataset) {
804                        return  DTypeUtils.fromDoublesToBiggestPrimitives((double[]) a.sum(ignoreInvalids), dtype);
805                }
806                return DTypeUtils.fromDoubleToBiggestNumber(DTypeUtils.toReal(a.sum(ignoreInvalids)), dtype);
807        }
808
809        /**
810         * @param a dataset
811         * @param dtype
812         * @param axis
813         * @param ignoreInvalids see {@link Dataset#max(int, boolean...)}
814         * @return typed sum of items along axis in dataset
815         * @since 2.0
816         */
817        public static Dataset typedSum(final Dataset a, int dtype, int axis, final boolean... ignoreInvalids) {
818                return DatasetUtils.cast(a.sum(axis, ignoreInvalids), dtype);
819        }
820
821        /**
822         * @param a dataset
823         * @param ignoreInvalids see {@link IDataset#max(boolean...)}
824         * @return product of all items in dataset
825         * @since 2.0
826         */
827        public static Object product(final Dataset a, final boolean... ignoreInvalids) {
828                return typedProduct(a, a.getDType(), ignoreInvalids);
829        }
830
831        /**
832         * @param a dataset
833         * @param axis
834         * @param ignoreInvalids see {@link Dataset#max(int, boolean...)}
835         * @return product of items along axis in dataset
836         * @since 2.0
837         */
838        public static Dataset product(final Dataset a, final int axis, final boolean... ignoreInvalids) {
839                return typedProduct(a, a.getDType(), axis, ignoreInvalids);
840        }
841
842        /**
843         * @param a dataset
844         * @param dtype
845         * @param ignoreInvalids see {@link IDataset#max(boolean...)}
846         * @return typed product of all items in dataset
847         * @since 2.0
848         */
849        public static Object typedProduct(final Dataset a, final int dtype, final boolean... ignoreInvalids) {
850                final boolean ignoreNaNs;
851                final boolean ignoreInfs;
852                if (a.hasFloatingPointElements()) {
853                        ignoreNaNs = ignoreInvalids != null && ignoreInvalids.length > 0 ? ignoreInvalids[0] : false;
854                        ignoreInfs = ignoreInvalids != null && ignoreInvalids.length > 1 ? ignoreInvalids[1] : ignoreNaNs;
855                } else {
856                        ignoreNaNs = false;
857                        ignoreInfs = false;
858                }
859
860                if (a.isComplex()) {
861                        IndexIterator it = a.getIterator();
862                        double rv = 1, iv = 0;
863
864                        while (it.hasNext()) {
865                                final double r1 = a.getElementDoubleAbs(it.index);
866                                final double i1 = a.getElementDoubleAbs(it.index + 1);
867                                if (ignoreNaNs && (Double.isNaN(r1) || Double.isNaN(i1))) {
868                                        continue;
869                                }
870                                if (ignoreInfs  && (Double.isInfinite(r1) || Double.isInfinite(i1))) {
871                                        continue;
872                                }
873                                final double tv = r1*rv - i1*iv;
874                                iv = r1*iv + i1*rv;
875                                rv = tv;
876                                if (Double.isNaN(rv) && Double.isNaN(iv)) {
877                                        break;
878                                }
879                        }
880
881                        return new Complex(rv, iv);
882                }
883
884                IndexIterator it = a.getIterator();
885                final int is;
886                final long[] lresults;
887                final double[] dresults;
888
889                switch (dtype) {
890                case Dataset.BOOL:
891                case Dataset.INT8: case Dataset.INT16:
892                case Dataset.INT32: case Dataset.INT64:
893                        long lresult = 1;
894                        while (it.hasNext()) {
895                                lresult *= a.getElementLongAbs(it.index);
896                        }
897                        return new Long(lresult);
898                case Dataset.ARRAYINT8: case Dataset.ARRAYINT16:
899                case Dataset.ARRAYINT32: case Dataset.ARRAYINT64:
900                        lresult = 1;
901                        is = a.getElementsPerItem();
902                        lresults = new long[is];
903                        for (int j = 0; j < is; j++) {
904                                lresults[j] = 1;
905                        }
906                        while (it.hasNext()) {
907                                for (int j = 0; j < is; j++)
908                                        lresults[j] *= a.getElementLongAbs(it.index+j);
909                        }
910                        return lresults;
911                case Dataset.FLOAT32: case Dataset.FLOAT64:
912                        double dresult = 1.;
913                        while (it.hasNext()) {
914                                final double x = a.getElementDoubleAbs(it.index);
915                                if (ignoreNaNs && Double.isNaN(x)) {
916                                        continue;
917                                }
918                                if (ignoreInfs && Double.isInfinite(x)) {
919                                        continue;
920                                }
921                                dresult *= x;
922                                if (Double.isNaN(dresult)) {
923                                        break;
924                                }
925                        }
926                        return Double.valueOf(dresult);
927                case Dataset.ARRAYFLOAT32:
928                case Dataset.ARRAYFLOAT64:
929                        is = a.getElementsPerItem();
930                        double[] vals = new double[is];
931                        dresults = new double[is];
932                        for (int j = 0; j < is; j++) {
933                                dresults[j] = 1.;
934                        }
935                        while (it.hasNext()) {
936                                boolean okay = true;
937                                for (int j = 0; j < is; j++) {
938                                        final double val = a.getElementDoubleAbs(it.index + j);
939                                        if (ignoreNaNs && Double.isNaN(val)) {
940                                                okay = false;
941                                                break;
942                                        }
943                                        if (ignoreInfs && Double.isInfinite(val)) {
944                                                okay = false;
945                                                break;
946                                        }
947                                        vals[j] = val;
948                                }
949                                if (okay) {
950                                        for (int j = 0; j < is; j++) {
951                                                double val = vals[j];
952                                                dresults[j] *= val;
953                                        }
954                                }
955                        }
956                        return dresults;
957                }
958
959                return null;
960        }
961
962        /**
963         * @param a dataset
964         * @param dtype
965         * @param axis
966         * @param ignoreInvalids see {@link IDataset#max(boolean...)}
967         * @return typed product of items along axis in dataset
968         * @since 2.0
969         */
970        public static Dataset typedProduct(final Dataset a, final int dtype, int axis, final boolean... ignoreInvalids) {
971                axis = a.checkAxis(axis);
972                final int[] oshape = a.getShape();
973                final int is = a.getElementsPerItem();
974                final int alen = oshape[axis];
975                oshape[axis] = 1;
976
977                final boolean ignoreNaNs;
978                final boolean ignoreInfs;
979                if (a.hasFloatingPointElements()) {
980                        ignoreNaNs = ignoreInvalids != null && ignoreInvalids.length > 0 ? ignoreInvalids[0] : false;
981                        ignoreInfs = ignoreInvalids != null && ignoreInvalids.length > 1 ? ignoreInvalids[1] : ignoreNaNs;
982                } else {
983                        ignoreNaNs = false;
984                        ignoreInfs = false;
985                }
986                @SuppressWarnings("deprecation")
987                Dataset result = DatasetFactory.zeros(is, oshape, dtype);
988
989                IndexIterator qiter = result.getIterator(true);
990                int[] qpos = qiter.getPos();
991                int[] spos;
992
993                // TODO add getLongArray(long[], int...) to CompoundDataset
994                while (qiter.hasNext()) {
995                        spos = qpos.clone();
996
997                        if (a.isComplex()) {
998                                double rv = 1, iv = 0;
999                                switch (dtype) {
1000                                case Dataset.COMPLEX64:
1001                                        ComplexFloatDataset af = (ComplexFloatDataset) a;
1002                                        for (int j = 0; j < alen; j++) {
1003                                                spos[axis] = j;
1004                                                final float r1 = af.getReal(spos);
1005                                                final float i1 = af.getImag(spos);
1006                                                if (ignoreNaNs && (Float.isNaN(r1) || Float.isNaN(i1))) {
1007                                                        continue;
1008                                                }
1009                                                if (ignoreInfs  && (Float.isInfinite(r1) || Float.isInfinite(i1))) {
1010                                                        continue;
1011                                                }
1012                                                final double tv = r1*rv - i1*iv;
1013                                                iv = r1*iv + i1*rv;
1014                                                rv = tv;
1015                                                if (Double.isNaN(rv) && Double.isNaN(iv)) {
1016                                                        break;
1017                                                }
1018                                        }
1019                                        break;
1020                                case Dataset.COMPLEX128:
1021                                        ComplexDoubleDataset ad = (ComplexDoubleDataset) a;
1022                                        for (int j = 0; j < alen; j++) {
1023                                                spos[axis] = j;
1024                                                final double r1 = ad.getReal(spos);
1025                                                final double i1 = ad.getImag(spos);
1026                                                if (ignoreNaNs && (Double.isNaN(r1) || Double.isNaN(i1))) {
1027                                                        continue;
1028                                                }
1029                                                if (ignoreInfs  && (Double.isInfinite(r1) || Double.isInfinite(i1))) {
1030                                                        continue;
1031                                                }
1032                                                final double tv = r1*rv - i1*iv;
1033                                                iv = r1*iv + i1*rv;
1034                                                rv = tv;
1035                                                if (Double.isNaN(rv) && Double.isNaN(iv)) {
1036                                                        break;
1037                                                }
1038                                        }
1039                                        break;
1040                                }
1041
1042                                result.set(new Complex(rv, iv), qpos);
1043                        } else {
1044                                final long[] lresults;
1045                                final double[] dresults;
1046
1047                                switch (dtype) {
1048                                case Dataset.BOOL:
1049                                case Dataset.INT8: case Dataset.INT16:
1050                                case Dataset.INT32: case Dataset.INT64:
1051                                        long lresult = 1;
1052                                        for (int j = 0; j < alen; j++) {
1053                                                spos[axis] = j;
1054                                                lresult *= a.getInt(spos);
1055                                        }
1056                                        result.set(lresult, qpos);
1057                                        break;
1058                                case Dataset.ARRAYINT8:
1059                                        lresults = new long[is];
1060                                        for (int k = 0; k < is; k++) {
1061                                                lresults[k] = 1;
1062                                        }
1063                                        for (int j = 0; j < alen; j++) {
1064                                                spos[axis] = j;
1065                                                final byte[] va = (byte[]) a.getObject(spos);
1066                                                for (int k = 0; k < is; k++) {
1067                                                        lresults[k] *= va[k];
1068                                                }
1069                                        }
1070                                        result.set(lresults, qpos);
1071                                        break;
1072                                case Dataset.ARRAYINT16:
1073                                        lresults = new long[is];
1074                                        for (int k = 0; k < is; k++) {
1075                                                lresults[k] = 1;
1076                                        }
1077                                        for (int j = 0; j < alen; j++) {
1078                                                spos[axis] = j;
1079                                                final short[] va = (short[]) a.getObject(spos);
1080                                                for (int k = 0; k < is; k++) {
1081                                                        lresults[k] *= va[k];
1082                                                }
1083                                        }
1084                                        result.set(lresults, qpos);
1085                                        break;
1086                                case Dataset.ARRAYINT32:
1087                                        lresults = new long[is];
1088                                        for (int k = 0; k < is; k++) {
1089                                                lresults[k] = 1;
1090                                        }
1091                                        for (int j = 0; j < alen; j++) {
1092                                                spos[axis] = j;
1093                                                final int[] va = (int[]) a.getObject(spos);
1094                                                for (int k = 0; k < is; k++) {
1095                                                        lresults[k] *= va[k];
1096                                                }
1097                                        }
1098                                        result.set(lresults, qpos);
1099                                        break;
1100                                case Dataset.ARRAYINT64:
1101                                        lresults = new long[is];
1102                                        for (int k = 0; k < is; k++) {
1103                                                lresults[k] = 1;
1104                                        }
1105                                        for (int j = 0; j < alen; j++) {
1106                                                spos[axis] = j;
1107                                                final long[] va = (long[]) a.getObject(spos);
1108                                                for (int k = 0; k < is; k++) {
1109                                                        lresults[k] *= va[k];
1110                                                }
1111                                        }
1112                                        result.set(lresults, qpos);
1113                                        break;
1114                                case Dataset.FLOAT32: case Dataset.FLOAT64:
1115                                        double dresult = 1.;
1116                                        for (int j = 0; j < alen; j++) {
1117                                                spos[axis] = j;
1118                                                final double x = a.getDouble(spos); 
1119                                                if (ignoreNaNs && Double.isNaN(x)) {
1120                                                        continue;
1121                                                }
1122                                                if (ignoreInfs && Double.isInfinite(x)) {
1123                                                        continue;
1124                                                }
1125                                                dresult *= x;
1126                                                if (Double.isNaN(dresult)) {
1127                                                        break;
1128                                                }
1129                                        }
1130                                        result.set(dresult, qpos);
1131                                        break;
1132                                case Dataset.ARRAYFLOAT32: case Dataset.ARRAYFLOAT64:
1133                                        CompoundDataset da = (CompoundDataset) a;
1134                                        double[] dvalues = new double[is];
1135                                        dresults = new double[is];
1136                                        for (int k = 0; k < is; k++) {
1137                                                dresults[k] = 1.;
1138                                        }
1139                                        for (int j = 0; j < alen; j++) {
1140                                                spos[axis] = j;
1141                                                da.getDoubleArray(dvalues, spos);
1142                                                boolean okay = true;
1143                                                for (int k = 0; k < is; k++) {
1144                                                        final double val = dvalues[k];
1145                                                        if (ignoreNaNs && Double.isNaN(val)) {
1146                                                                okay = false;
1147                                                                break;
1148                                                        }
1149                                                        if (ignoreInfs && Double.isInfinite(val)) {
1150                                                                okay = false;
1151                                                                break;
1152                                                        }
1153                                                }
1154                                                if (okay) {
1155                                                        for (int k = 0; k < is; k++) {
1156                                                                dresults[k] *= dvalues[k];
1157                                                        }
1158                                                }
1159                                        }
1160                                        result.set(dresults, qpos);
1161                                        break;
1162                                }
1163                        }
1164                }
1165
1166                result.setShape(ShapeUtils.squeezeShape(oshape, axis));
1167                return result;
1168        }
1169
1170        /**
1171         * @param a dataset
1172         * @param ignoreInvalids see {@link IDataset#max(boolean...)}
1173         * @return cumulative product of items in flattened dataset
1174         * @since 2.0
1175         */
1176        public static Dataset cumulativeProduct(final Dataset a, final boolean... ignoreInvalids) {
1177                return cumulativeProduct(a.flatten(), 0, ignoreInvalids);
1178        }
1179
1180        /**
1181         * @param a dataset
1182         * @param axis
1183         * @param ignoreInvalids see {@link Dataset#max(int, boolean...)}
1184         * @return cumulative product of items along axis in dataset
1185         * @since 2.0
1186         */
1187        public static Dataset cumulativeProduct(final Dataset a, int axis, final boolean... ignoreInvalids) {
1188                axis = a.checkAxis(axis);
1189                int dtype = a.getDType();
1190                int[] oshape = a.getShape();
1191                int alen = oshape[axis];
1192                oshape[axis] = 1;
1193
1194                final boolean ignoreNaNs;
1195                final boolean ignoreInfs;
1196                if (a.hasFloatingPointElements()) {
1197                        ignoreNaNs = ignoreInvalids != null && ignoreInvalids.length > 0 ? ignoreInvalids[0] : false;
1198                        ignoreInfs = ignoreInvalids != null && ignoreInvalids.length > 1 ? ignoreInvalids[1] : ignoreNaNs;
1199                } else {
1200                        ignoreNaNs = false;
1201                        ignoreInfs = false;
1202                }
1203                Dataset result = DatasetFactory.zeros(a);
1204                PositionIterator pi = result.getPositionIterator(axis);
1205
1206                int[] pos = pi.getPos();
1207
1208                while (pi.hasNext()) {
1209
1210                        if (a.isComplex()) {
1211                                double rv = 1, iv = 0;
1212                                switch (dtype) {
1213                                case Dataset.COMPLEX64:
1214                                        ComplexFloatDataset af = (ComplexFloatDataset) a;
1215                                        ComplexFloatDataset rf = (ComplexFloatDataset) result;
1216                                        for (int j = 0; j < alen; j++) {
1217                                                if (!Double.isNaN(rv) || !Double.isNaN(iv)) {
1218                                                        pos[axis] = j;
1219                                                        final float r1 = af.getReal(pos);
1220                                                        final float i1 = af.getImag(pos);
1221                                                        if (ignoreNaNs && (Float.isNaN(r1) || Float.isNaN(i1))) {
1222                                                                continue;
1223                                                        }
1224                                                        if (ignoreInfs  && (Float.isInfinite(r1) || Float.isInfinite(i1))) {
1225                                                                continue;
1226                                                        }
1227                                                        final double tv = r1*rv - i1*iv;
1228                                                        iv = r1*iv + i1*rv;
1229                                                        rv = tv;
1230                                                }
1231                                                rf.set((float) rv, (float) iv, pos);
1232                                        }
1233                                        break;
1234                                case Dataset.COMPLEX128:
1235                                        ComplexDoubleDataset ad = (ComplexDoubleDataset) a;
1236                                        ComplexDoubleDataset rd = (ComplexDoubleDataset) result;
1237                                        for (int j = 0; j < alen; j++) {
1238                                                if (!Double.isNaN(rv) || !Double.isNaN(iv)) {
1239                                                        pos[axis] = j;
1240                                                        final double r1 = ad.getReal(pos);
1241                                                        final double i1 = ad.getImag(pos);
1242                                                        if (ignoreNaNs && (Double.isNaN(r1) || Double.isNaN(i1))) {
1243                                                                continue;
1244                                                        }
1245                                                        if (ignoreInfs  && (Double.isInfinite(r1) || Double.isInfinite(i1))) {
1246                                                                continue;
1247                                                        }
1248                                                        final double tv = r1*rv - i1*iv;
1249                                                        iv = r1*iv + i1*rv;
1250                                                        rv = tv;
1251                                                }
1252                                                rd.set(rv, iv, pos);
1253                                        }
1254                                        break;
1255                                }
1256                        } else {
1257                                final int is;
1258                                final long[] lresults;
1259                                final double[] dresults;
1260
1261                                switch (dtype) {
1262                                case Dataset.BOOL:
1263                                case Dataset.INT8: case Dataset.INT16:
1264                                case Dataset.INT32: case Dataset.INT64:
1265                                        long lresult = 1;
1266                                        for (int j = 0; j < alen; j++) {
1267                                                pos[axis] = j;
1268                                                lresult *= a.getInt(pos);
1269                                                result.set(lresult, pos);
1270                                        }
1271                                        break;
1272                                case Dataset.ARRAYINT8:
1273                                        is = a.getElementsPerItem();
1274                                        lresults = new long[is];
1275                                        for (int k = 0; k < is; k++) {
1276                                                lresults[k] = 1;
1277                                        }
1278                                        for (int j = 0; j < alen; j++) {
1279                                                pos[axis] = j;
1280                                                final byte[] va = (byte[]) a.getObject(pos);
1281                                                for (int k = 0; k < is; k++) {
1282                                                        lresults[k] *= va[k];
1283                                                }
1284                                                result.set(lresults, pos);
1285                                        }
1286                                        break;
1287                                case Dataset.ARRAYINT16:
1288                                        is = a.getElementsPerItem();
1289                                        lresults = new long[is];
1290                                        for (int k = 0; k < is; k++) {
1291                                                lresults[k] = 1;
1292                                        }
1293                                        for (int j = 0; j < alen; j++) {
1294                                                pos[axis] = j;
1295                                                final short[] va = (short[]) a.getObject(pos);
1296                                                for (int k = 0; k < is; k++) {
1297                                                        lresults[k] *= va[k];
1298                                                }
1299                                                result.set(lresults, pos);
1300                                        }
1301                                        break;
1302                                case Dataset.ARRAYINT32:
1303                                        is = a.getElementsPerItem();
1304                                        lresults = new long[is];
1305                                        for (int k = 0; k < is; k++) {
1306                                                lresults[k] = 1;
1307                                        }
1308                                        for (int j = 0; j < alen; j++) {
1309                                                pos[axis] = j;
1310                                                final int[] va = (int[]) a.getObject(pos);
1311                                                for (int k = 0; k < is; k++) {
1312                                                        lresults[k] *= va[k];
1313                                                }
1314                                                result.set(lresults, pos);
1315                                        }
1316                                        break;
1317                                case Dataset.ARRAYINT64:
1318                                        is = a.getElementsPerItem();
1319                                        lresults = new long[is];
1320                                        for (int k = 0; k < is; k++) {
1321                                                lresults[k] = 1;
1322                                        }
1323                                        for (int j = 0; j < alen; j++) {
1324                                                pos[axis] = j;
1325                                                final long[] va = (long[]) a.getObject(pos);
1326                                                for (int k = 0; k < is; k++) {
1327                                                        lresults[k] *= va[k];
1328                                                }
1329                                                result.set(lresults, pos);
1330                                        }
1331                                        break;
1332                                case Dataset.FLOAT32: case Dataset.FLOAT64:
1333                                        double dresult = 1.;
1334                                        for (int j = 0; j < alen; j++) {
1335                                                if (!Double.isNaN(dresult)) {
1336                                                        pos[axis] = j;
1337                                                        final double x = a.getDouble(pos);
1338                                                        if (ignoreNaNs && Double.isNaN(x)) {
1339                                                                continue;
1340                                                        }
1341                                                        if (ignoreInfs && Double.isInfinite(x)) {
1342                                                                continue;
1343                                                        }
1344                                                        dresult *= x;
1345                                                }
1346                                                result.set(dresult, pos);
1347                                        }
1348                                        break;
1349                                case Dataset.ARRAYFLOAT32: case Dataset.ARRAYFLOAT64:
1350                                        is = a.getElementsPerItem();
1351                                        CompoundDataset da = (CompoundDataset) a;
1352                                        double[] dvalues = new double[is];
1353                                        dresults = new double[is];
1354                                        for (int k = 0; k < is; k++) {
1355                                                dresults[k] = 1.;
1356                                        }
1357                                        for (int j = 0; j < alen; j++) {
1358                                                pos[axis] = j;
1359                                                da.getDoubleArray(dvalues, pos);
1360                                                boolean okay = true;
1361                                                for (int k = 0; k < is; k++) {
1362                                                        final double val = dvalues[k];
1363                                                        if (ignoreNaNs && Double.isNaN(val)) {
1364                                                                okay = false;
1365                                                                break;
1366                                                        }
1367                                                        if (ignoreInfs && Double.isInfinite(val)) {
1368                                                                okay = false;
1369                                                                break;
1370                                                        }
1371                                                }
1372                                                if (okay) {
1373                                                        for (int k = 0; k < is; k++) {
1374                                                                dresults[k] *= dvalues[k];
1375                                                        }
1376                                                }
1377                                                result.set(dresults, pos);
1378                                        }
1379                                        break;
1380                                }
1381                        }
1382                }
1383
1384                return result;
1385        }
1386
1387        /**
1388         * @param a dataset
1389         * @param ignoreInvalids see {@link IDataset#max(boolean...)}
1390         * @return cumulative sum of items in flattened dataset
1391         * @since 2.0
1392         */
1393        public static Dataset cumulativeSum(final Dataset a, final boolean... ignoreInvalids) {
1394                return cumulativeSum(a.flatten(), 0, ignoreInvalids);
1395        }
1396
1397        /**
1398         * @param a dataset
1399         * @param axis
1400         * @param ignoreInvalids see {@link Dataset#max(int, boolean...)}
1401         * @return cumulative sum of items along axis in dataset
1402         * @since 2.0
1403         */
1404        public static Dataset cumulativeSum(final Dataset a, int axis, final boolean... ignoreInvalids) {
1405                axis = a.checkAxis(axis);
1406                int dtype = a.getDType();
1407                int[] oshape = a.getShape();
1408                int alen = oshape[axis];
1409                oshape[axis] = 1;
1410
1411                final boolean ignoreNaNs;
1412                final boolean ignoreInfs;
1413                if (a.hasFloatingPointElements()) {
1414                        ignoreNaNs = ignoreInvalids != null && ignoreInvalids.length > 0 ? ignoreInvalids[0] : false;
1415                        ignoreInfs = ignoreInvalids != null && ignoreInvalids.length > 1 ? ignoreInvalids[1] : ignoreNaNs;
1416                } else {
1417                        ignoreNaNs = false;
1418                        ignoreInfs = false;
1419                }
1420                Dataset result = DatasetFactory.zeros(a);
1421                PositionIterator pi = result.getPositionIterator(axis);
1422
1423                int[] pos = pi.getPos();
1424
1425                while (pi.hasNext()) {
1426
1427                        if (a.isComplex()) {
1428                                double rv = 0, iv = 0;
1429                                switch (dtype) {
1430                                case Dataset.COMPLEX64:
1431                                        ComplexFloatDataset af = (ComplexFloatDataset) a;
1432                                        ComplexFloatDataset rf = (ComplexFloatDataset) result;
1433                                        for (int j = 0; j < alen; j++) {
1434                                                if (!Double.isNaN(rv) || !Double.isNaN(iv)) {
1435                                                        pos[axis] = j;
1436                                                        final float r1 = af.getReal(pos);
1437                                                        final float i1 = af.getImag(pos);
1438                                                        if (ignoreNaNs && (Float.isNaN(r1) || Float.isNaN(i1))) {
1439                                                                continue;
1440                                                        }
1441                                                        if (ignoreInfs  && (Float.isInfinite(r1) || Float.isInfinite(i1))) {
1442                                                                continue;
1443                                                        }
1444                                                        rv += r1;
1445                                                        iv += i1;
1446                                                }
1447                                                rf.set((float) rv, (float) iv, pos);
1448                                        }
1449                                        break;
1450                                case Dataset.COMPLEX128:
1451                                        ComplexDoubleDataset ad = (ComplexDoubleDataset) a;
1452                                        ComplexDoubleDataset rd = (ComplexDoubleDataset) result;
1453                                        for (int j = 0; j < alen; j++) {
1454                                                if (!Double.isNaN(rv) || !Double.isNaN(iv)) {
1455                                                        pos[axis] = j;
1456                                                        final double r1 = ad.getReal(pos);
1457                                                        final double i1 = ad.getImag(pos);
1458                                                        if (ignoreNaNs && (Double.isNaN(r1) || Double.isNaN(i1))) {
1459                                                                continue;
1460                                                        }
1461                                                        if (ignoreInfs  && (Double.isInfinite(r1) || Double.isInfinite(i1))) {
1462                                                                continue;
1463                                                        }
1464                                                        rv += r1;
1465                                                        iv += i1;
1466                                                }
1467                                                rd.set(rv, iv, pos);
1468                                        }
1469                                        break;
1470                                }
1471                        } else {
1472                                final int is;
1473                                final long[] lresults;
1474                                final double[] dresults;
1475
1476                                switch (dtype) {
1477                                case Dataset.BOOL:
1478                                case Dataset.INT8: case Dataset.INT16:
1479                                case Dataset.INT32: case Dataset.INT64:
1480                                        long lresult = 0;
1481                                        for (int j = 0; j < alen; j++) {
1482                                                pos[axis] = j;
1483                                                lresult += a.getInt(pos);
1484                                                result.set(lresult, pos);
1485                                        }
1486                                        break;
1487                                case Dataset.ARRAYINT8:
1488                                        is = a.getElementsPerItem();
1489                                        lresults = new long[is];
1490                                        for (int j = 0; j < alen; j++) {
1491                                                pos[axis] = j;
1492                                                final byte[] va = (byte[]) a.getObject(pos);
1493                                                for (int k = 0; k < is; k++) {
1494                                                        lresults[k] += va[k];
1495                                                }
1496                                                result.set(lresults, pos);
1497                                        }
1498                                        break;
1499                                case Dataset.ARRAYINT16:
1500                                        is = a.getElementsPerItem();
1501                                        lresults = new long[is];
1502                                        for (int j = 0; j < alen; j++) {
1503                                                pos[axis] = j;
1504                                                final short[] va = (short[]) a.getObject(pos);
1505                                                for (int k = 0; k < is; k++) {
1506                                                        lresults[k] += va[k];
1507                                                }
1508                                                result.set(lresults, pos);
1509                                        }
1510                                        break;
1511                                case Dataset.ARRAYINT32:
1512                                        is = a.getElementsPerItem();
1513                                        lresults = new long[is];
1514                                        for (int j = 0; j < alen; j++) {
1515                                                pos[axis] = j;
1516                                                final int[] va = (int[]) a.getObject(pos);
1517                                                for (int k = 0; k < is; k++) {
1518                                                        lresults[k] += va[k];
1519                                                }
1520                                                result.set(lresults, pos);
1521                                        }
1522                                        break;
1523                                case Dataset.ARRAYINT64:
1524                                        is = a.getElementsPerItem();
1525                                        lresults = new long[is];
1526                                        for (int j = 0; j < alen; j++) {
1527                                                pos[axis] = j;
1528                                                final long[] va = (long[]) a.getObject(pos);
1529                                                for (int k = 0; k < is; k++) {
1530                                                        lresults[k] += va[k];
1531                                                }
1532                                                result.set(lresults, pos);
1533                                        }
1534                                        break;
1535                                case Dataset.FLOAT32: case Dataset.FLOAT64:
1536                                        double dresult = 0.;
1537                                        for (int j = 0; j < alen; j++) {
1538                                                if (!Double.isNaN(dresult)) {
1539                                                        pos[axis] = j;
1540                                                        final double x = a.getDouble(pos);
1541                                                        if (ignoreNaNs && Double.isNaN(x)) {
1542                                                                continue;
1543                                                        }
1544                                                        if (ignoreInfs && Double.isInfinite(x)) {
1545                                                                continue;
1546                                                        }
1547                                                        dresult += x;
1548                                                }
1549                                                result.set(dresult, pos);
1550                                        }
1551                                        break;
1552                                case Dataset.ARRAYFLOAT32: case Dataset.ARRAYFLOAT64:
1553                                        is = a.getElementsPerItem();
1554                                        CompoundDataset da = (CompoundDataset) a;
1555                                        double[] dvalues = new double[is];
1556                                        dresults = new double[is];
1557                                        for (int j = 0; j < alen; j++) {
1558                                                pos[axis] = j;
1559                                                da.getDoubleArray(dvalues, pos);
1560                                                boolean okay = true;
1561                                                for (int k = 0; k < is; k++) {
1562                                                        final double val = dvalues[k];
1563                                                        if (ignoreNaNs && Double.isNaN(val)) {
1564                                                                okay = false;
1565                                                                break;
1566                                                        }
1567                                                        if (ignoreInfs && Double.isInfinite(val)) {
1568                                                                okay = false;
1569                                                                break;
1570                                                        }
1571                                                }
1572                                                if (okay) {
1573                                                        for (int k = 0; k < is; k++) {
1574                                                                dresults[k] += dvalues[k];
1575                                                        }
1576                                                }
1577                                                result.set(dresults, pos);
1578                                        }
1579                                        break;
1580                                }
1581                        }
1582                }
1583
1584                return result;
1585        }
1586
1587        /**
1588         * @param a
1589         * @return average deviation value of all items the dataset
1590         */
1591        public static Object averageDeviation(final Dataset a) {
1592                final IndexIterator it = a.getIterator();
1593                final int is = a.getElementsPerItem();
1594
1595                if (is == 1) {
1596                        double mean = (Double) a.mean();
1597                        double sum = 0.0;
1598
1599                        while (it.hasNext()) {
1600                                sum += Math.abs(a.getElementDoubleAbs(it.index) - mean);
1601                        }
1602
1603                        return sum / a.getSize();
1604                }
1605
1606                double[] means = (double[]) a.mean();
1607                double[] sums = new double[is];
1608
1609                while (it.hasNext()) {
1610                        for (int j = 0; j < is; j++)
1611                                sums[j] += Math.abs(a.getElementDoubleAbs(it.index + j) - means[j]);
1612                }
1613
1614                double n = a.getSize();
1615                for (int j = 0; j < is; j++)
1616                        sums[j] /= n;
1617
1618                return sums;
1619        }
1620
1621        /**
1622         * The residual is the sum of squared differences
1623         * @param a
1624         * @param b
1625         * @return residual value
1626         */
1627        public static double residual(final Dataset a, final Dataset b) {
1628                return a.residual(b);
1629        }
1630
1631        /**
1632         * The residual is the sum of squared differences
1633         * @param a
1634         * @param b
1635         * @param w
1636         * @return residual value
1637         */
1638        public static double weightedResidual(final Dataset a, final Dataset b, final Dataset w) {
1639                return a.residual(b, w, false);
1640        }
1641
1642        /**
1643         * Calculate approximate outlier values. These are defined as the values in the dataset
1644         * that are approximately below and above the given thresholds - in terms of percentages
1645         * of dataset size.
1646         * <p>
1647         * It approximates by limiting the number of items (given by length) used internally by
1648         * data structures - the larger this is, the more accurate will those outlier values become.
1649         * The actual thresholds used are returned in the array.
1650         * <p>
1651         * Also, the low and high values will be made distinct if possible by adjusting the thresholds
1652         * @param a
1653         * @param lo percentage threshold for lower limit
1654         * @param hi percentage threshold for higher limit
1655         * @param length maximum number of items used internally, if negative, then unlimited
1656         * @return double array with low and high values, and low and high percentage thresholds
1657         */
1658        public static double[] outlierValues(final Dataset a, double lo, double hi, final int length) {
1659                if (lo <= 0 || hi <= 0 || lo >= hi || hi >= 100  || Double.isNaN(lo)|| Double.isNaN(hi)) {
1660                        throw new IllegalArgumentException("Thresholds must be between (0,100) and in order");
1661                }
1662                final int size = a.getSize();
1663                int nl = Math.max((int) ((lo*size)/100), 1);
1664                if (length > 0 && nl > length)
1665                        nl = length;
1666                int nh = Math.max((int) (((100-hi)*size)/100), 1);
1667                if (length > 0 && nh > length)
1668                        nh = length;
1669
1670                double[] results = Math.max(nl, nh) > 640 ? outlierValuesMap(a, nl, nh) : outlierValuesList(a, nl, nh);
1671
1672                results[2] = results[2]*100./size;
1673                results[3] = 100. - results[3]*100./size;
1674                return results;
1675        }
1676
1677        static double[] outlierValuesMap(final Dataset a, int nl, int nh) {
1678                final TreeMap<Double, Integer> lMap = new TreeMap<Double, Integer>();
1679                final TreeMap<Double, Integer> hMap = new TreeMap<Double, Integer>();
1680
1681                int ml = 0;
1682                int mh = 0;
1683                IndexIterator it = a.getIterator();
1684                while (it.hasNext()) {
1685                        Double x = a.getElementDoubleAbs(it.index);
1686                        Integer i;
1687                        if (ml == nl) {
1688                                Double k = lMap.lastKey();
1689                                if (x < k) {
1690                                        i = lMap.get(k) - 1;
1691                                        if (i == 0) {
1692                                                lMap.remove(k);
1693                                        } else {
1694                                                lMap.put(k, i);
1695                                        }
1696                                        i = lMap.get(x);
1697                                        if (i == null) {
1698                                                lMap.put(x, 1);
1699                                        } else {
1700                                                lMap.put(x, i + 1);
1701                                        }
1702                                }
1703                        } else {
1704                                i = lMap.get(x);
1705                                if (i == null) {
1706                                        lMap.put(x, 1);
1707                                } else {
1708                                        lMap.put(x, i + 1);
1709                                }
1710                                ml++;
1711                        }
1712
1713                        if (mh == nh) {
1714                                Double k = hMap.firstKey();
1715                                if (x > k) {
1716                                        i = hMap.get(k) - 1;
1717                                        if (i == 0) {
1718                                                hMap.remove(k);
1719                                        } else {
1720                                                hMap.put(k, i);
1721                                        }
1722                                        i = hMap.get(x);
1723                                        if (i == null) {
1724                                                hMap.put(x, 1);
1725                                        } else {
1726                                                hMap.put(x, i+1);
1727                                        }
1728                                }
1729                        } else {
1730                                i = hMap.get(x);
1731                                if (i == null) {
1732                                        hMap.put(x, 1);
1733                                } else {
1734                                        hMap.put(x, i+1);
1735                                }
1736                                mh++;
1737                        }
1738                }
1739
1740                // Attempt to make values distinct
1741                double lx = lMap.lastKey();
1742                double hx = hMap.firstKey();
1743                if (lx >= hx) {
1744                        Double h = hMap.higherKey(lx);
1745                        if (h != null) {
1746                                hx = h;
1747                                mh--;
1748                        } else {
1749                                Double l = lMap.lowerKey(hx);
1750                                if (l != null) {
1751                                        lx = l;
1752                                        ml--;
1753                                }
1754                        }
1755                        
1756                }
1757                return new double[] {lMap.lastKey(), hMap.firstKey(), ml, mh};
1758        }
1759
1760        static double[] outlierValuesList(final Dataset a, int nl, int nh) {
1761                final List<Double> lList = new ArrayList<Double>(nl);
1762                final List<Double> hList = new ArrayList<Double>(nh);
1763//              final List<Double> lList = new LinkedList<Double>();
1764//              final List<Double> hList = new LinkedList<Double>();
1765
1766                double lx = Double.POSITIVE_INFINITY;
1767                double hx = Double.NEGATIVE_INFINITY;
1768
1769                IndexIterator it = a.getIterator();
1770                while (it.hasNext()) {
1771                        double x = a.getElementDoubleAbs(it.index);
1772                        if (x < lx) {
1773                                if (lList.size() == nl) {
1774                                        lList.remove(lx);
1775                                }
1776                                lList.add(x);
1777                                lx = Collections.max(lList);
1778                        } else if (x == lx) {
1779                                if (lList.size() < nl) {
1780                                        lList.add(x);
1781                                }
1782                        }
1783
1784                        if (x > hx) {
1785                                if (hList.size() == nh) {
1786                                        hList.remove(hx);
1787                                }
1788                                hList.add(x);
1789                                hx = Collections.min(hList);
1790                        } else if (x == hx) {
1791                                if (hList.size() < nh) {
1792                                        hList.add(x);
1793                                }
1794                        }
1795                }
1796
1797                nl = lList.size();
1798                nh = hList.size();
1799
1800                // Attempt to make values distinct
1801                if (lx >= hx) {
1802                        Collections.sort(hList);
1803                        for (double h : hList) {
1804                                if (h > hx) {
1805                                        hx = h;
1806                                        break;
1807                                }
1808                                nh--;
1809                        }
1810                        if (lx >= hx) {
1811                                Collections.sort(lList);
1812                                Collections.reverse(lList);
1813                                for (double l : lList) {
1814                                        if (l < lx) {
1815                                                lx = l;
1816                                                break;
1817                                        }
1818                                        nl--;
1819                                }
1820                        }
1821                }
1822                return new double[] {lx, hx, nl, nh};
1823        }
1824
1825        /**
1826         * See {@link #covariance(Dataset a, Dataset b, boolean rowvar, boolean bias, Integer ddof)} with b = null, rowvar = true, bias = false and ddof = null.
1827         * @param a
1828         * @return covariance array of a
1829         */
1830        public static Dataset covariance(final Dataset a) {
1831                return covariance(a, true, false, null); 
1832        }
1833
1834        /**
1835         * See {@link #covariance(Dataset a, Dataset b, boolean rowvar, boolean bias, Integer ddof)} with b = null.
1836         * @param a
1837         * @return covariance array of a
1838         * @since 2.0
1839         */
1840        public static Dataset covariance(final Dataset a, 
1841                        boolean rowvar, boolean bias, Integer ddof) {
1842                return covariance(a, null, rowvar, bias, ddof);
1843        }
1844
1845        /**
1846         * See {@link #covariance(Dataset a, Dataset b, boolean rowvar, boolean bias, Integer ddof)} with b = null, rowvar = true, bias = false and ddof = null.
1847         * @param a
1848         * @return covariance array of a concatenated with b
1849         */
1850        public static Dataset covariance(final Dataset a, final Dataset b) {
1851                return covariance(a, b, true, false, null);
1852        }
1853
1854        /**
1855         * Calculate the covariance matrix (array) of a concatenated with b. This 
1856         * method is directly based on the implementation in numpy (cov). 
1857         * @param a Array containing multiple variable and observations. Each row represents a variable, each column an observation.
1858         * @param b An extra set of variables and observations. Must be of same type as a and have a compatible shape. 
1859         * @param rowvar When true (default), each row is a variable; when false each column is a variable.
1860         * @param bias Default normalisation is (N - 1) - N is number of observations. If set true, normalisation is (N). 
1861         * @param ddof Default normalisation is (N - 1). If ddof is set, then normalisation is (N - ddof).
1862         * @return covariance array of a concatenated with b
1863         * @since 2.0
1864         */
1865        public static Dataset covariance (final Dataset a, final Dataset b, 
1866                        boolean rowvar, boolean bias, Integer ddof) {
1867                
1868                //Create a working copy of the dataset & check its rank.
1869                Dataset vars = a.clone();
1870                if (a.getRank() == 1) {
1871                        vars.setShape(1, a.getShape()[0]);
1872                }
1873                
1874                //1D of variables, so consider rows as variables
1875                if (vars.getShape()[0] == 1) {
1876                        rowvar = true;
1877                }
1878                
1879                //nr is the number of records; axis is index
1880                int nr, axis;
1881                if (rowvar) {
1882                        nr = vars.getShape()[1];
1883                        axis = 0;
1884                } else {
1885                        nr = vars.getShape()[0];
1886                        axis = 1;
1887                }
1888                
1889                //Set the reduced degrees of freedom & normalisation factor
1890                if (ddof == null) {
1891                        if (bias == false) {
1892                                ddof = 1;
1893                        } else {
1894                                ddof = 0;
1895                        }
1896                }
1897                double norm_fact = nr - ddof;
1898                if (norm_fact <= 0.) {
1899                        //TODO Some sort of warning here?
1900                        norm_fact = 0.;
1901                }
1902                
1903                //Concatenate additional set of variables with main set
1904                if (b != null) {
1905                        //Create a working copy of the dataset & check its rank.
1906                        Dataset extraVars = b.clone();
1907                        if (b.getRank() == 1) {
1908                                extraVars.setShape(1, a.getShape()[0]);
1909                        }
1910                        vars = DatasetUtils.concatenate(new Dataset[]{vars, extraVars}, axis);
1911                }
1912                
1913                //Calculate deviations & covariance matrix
1914                Dataset varsMean = vars.mean(1-axis, false);
1915                //-vars & varsMean need same shape (this is a hack!)
1916                int[] meanShape = vars.getShape();
1917                meanShape[1-axis] = 1;
1918                varsMean.setShape(meanShape);
1919                vars.isubtract(varsMean);
1920                Dataset cov;
1921                if (rowvar) {
1922                        cov = Maths.divide(LinearAlgebra.dotProduct(vars, Maths.conjugate(vars.transpose())), norm_fact).squeeze();
1923                } else {
1924                        cov = Maths.divide(LinearAlgebra.dotProduct(vars.transpose(), Maths.conjugate(vars)), norm_fact).squeeze();
1925                }
1926                return cov;
1927        }
1928}