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