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