    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 clazz dataset class
    796         * @param a dataset
    797         * @param ignoreInvalids see {@link IDataset#max(boolean...)}
    798         * @return typed sum of dataset
    799         * @since 2.3
    800         */
    801        public static Object typedSum(final Class<? extends Dataset> clazz, final Dataset a, 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 InterfaceUtils.fromDoublesToBiggestPrimitives(clazz, (double[]) a.sum(ignoreInvalids));
    807                }
    808                return InterfaceUtils.fromDoubleToBiggestNumber(clazz, DTypeUtils.toReal(a.sum(ignoreInvalids)));
    809        }
    810
    811        /**
    812         * @param a dataset
    813         * @param dtype
    814         * @param ignoreInvalids see {@link IDataset#max(boolean...)}
    815         * @return typed sum of dataset
    816         * @since 2.0
    817         * @deprecated Please use the class-based method {@link #typedSum(Class, Dataset, boolean...)}
    818         */
    819        @Deprecated
    820        public static Object typedSum(final Dataset a, int dtype, final boolean... ignoreInvalids) {
    821                return typedSum(DTypeUtils.getInterface(dtype), a, ignoreInvalids);
    822        }
    823
    824        /**
    825         * @param clazz dataset class
    826         * @param a dataset
    827         * @param axis
    828         * @param ignoreInvalids see {@link Dataset#max(int, boolean...)}
    829         * @return typed sum of items along axis in dataset
    830         * @since 2.3
    831         */
    832        public static <T extends Dataset> T typedSum(final Class<T> clazz, final Dataset a, int axis, final boolean... ignoreInvalids) {
    833                return a.sum(axis, ignoreInvalids).cast(clazz);
    834        }
    835
    836        /**
    837         * @param clazz dataset class
    838         * @param a dataset
    839         * @param axes
    840         * @param ignoreInvalids see {@link Dataset#max(int, boolean...)}
    841         * @return typed sum of items along axes in dataset
    842         * @since 2.3
    843         */
    844        public static <T extends Dataset> T typedSum(final Class<T> clazz, final Dataset a, int[] axes, final boolean... ignoreInvalids) {
    845                return a.sum(axes, ignoreInvalids).cast(clazz);
    846        }
    847
    848        /**
    849         * @param a dataset
    850         * @param dtype
    851         * @param axis
    852         * @param ignoreInvalids see {@link Dataset#max(int, boolean...)}
    853         * @return typed sum of items along axis in dataset
    854         * @since 2.0
    855         * @deprecated Please use the class-based method {@link #typedSum(Class, Dataset, axis, boolean...)}
    856         */
    857        @Deprecated
    858        public static Dataset typedSum(final Dataset a, int dtype, int axis, final boolean... ignoreInvalids) {
    859                return DatasetUtils.cast(a.sum(axis, ignoreInvalids), dtype);
    860        }
    861
    862        /**
    863         * @param a dataset
    864         * @param ignoreInvalids see {@link IDataset#max(boolean...)}
    865         * @return product of all items in dataset
    866         * @since 2.0
    867         */
    868        public static Object product(final Dataset a, final boolean... ignoreInvalids) {
    869                return typedProduct(a.getClass(), a, ignoreInvalids);
    870        }
    871
    872        /**
    873         * @param a dataset
    874         * @param axis
    875         * @param ignoreInvalids see {@link Dataset#max(int, boolean...)}
    876         * @return product of items along axis in dataset
    877         * @since 2.0
    878         */
    879        public static Dataset product(final Dataset a, final int axis, final boolean... ignoreInvalids) {
    880                return typedProduct(a.getClass(), a, axis, ignoreInvalids);
    881        }
    882
    883        /**
    884         * @param a dataset
    885         * @param axes
    886         * @param ignoreInvalids see {@link Dataset#max(int, boolean...)}
    887         * @return product of items along axes in dataset
    888         * @since 2.2
    889         */
    890        public static Dataset product(final Dataset a, final int[] axes, final boolean... ignoreInvalids) {
    891                return typedProduct(a.getClass(), a, axes, ignoreInvalids);
    892        }
    893
    894        /**
    895         * @param clazz dataset class
    896         * @param a dataset
    897         * @param dtype
    898         * @param ignoreInvalids see {@link IDataset#max(boolean...)}
    899         * @return typed product of all items in dataset
    900         * @since 2.3
    901         */
    902        public static Object typedProduct(final Class<? extends Dataset> clazz, final Dataset a, final boolean... ignoreInvalids) {
    903                return typedProduct(a, DTypeUtils.getDType(clazz), ignoreInvalids);
    904        }
    905
    906        /**
    907         * @param a dataset
    908         * @param dtype
    909         * @param ignoreInvalids see {@link IDataset#max(boolean...)}
    910         * @return typed product of all items in dataset
    911         * @since 2.0
    912         * @deprecated Please use the class-based method {@link #typedProduct(Class, Dataset, boolean...)}
    913         */
    914        @Deprecated
    915        public static Object typedProduct(final Dataset a, final int dtype, final boolean... ignoreInvalids) {
    916                final boolean ignoreNaNs;
    917                final boolean ignoreInfs;
    918                if (a.hasFloatingPointElements()) {
    919                        ignoreNaNs = ignoreInvalids != null && ignoreInvalids.length > 0 ? ignoreInvalids[0] : false;
    920                        ignoreInfs = ignoreInvalids != null && ignoreInvalids.length > 1 ? ignoreInvalids[1] : ignoreNaNs;
    921                } else {
    922                        ignoreNaNs = false;
    923                        ignoreInfs = false;
    924                }
    925
    926                if (a.isComplex()) {
    927                        IndexIterator it = a.getIterator();
    928                        double rv = 1, iv = 0;
    929
    930                        while (it.hasNext()) {
    931                                final double r1 = a.getElementDoubleAbs(it.index);
    932                                final double i1 = a.getElementDoubleAbs(it.index + 1);
    933                                if (ignoreNaNs && (Double.isNaN(r1) || Double.isNaN(i1))) {
    934                                        continue;
    935                                }
    936                                if (ignoreInfs  && (Double.isInfinite(r1) || Double.isInfinite(i1))) {
    937                                        continue;
    938                                }
    939                                final double tv = r1*rv - i1*iv;
    940                                iv = r1*iv + i1*rv;
    941                                rv = tv;
    942                                if (Double.isNaN(rv) && Double.isNaN(iv)) {
    943                                        break;
    944                                }
    945                        }
    946
    947                        return new Complex(rv, iv);
    948                }
    949
    950                IndexIterator it = a.getIterator();
    951                final int is;
    952                final long[] lresults;
    953                final double[] dresults;
    954
    955                switch (dtype) {
    956                case Dataset.BOOL:
    957                case Dataset.INT8: case Dataset.INT16:
    958                case Dataset.INT32: case Dataset.INT64:
    959                        long lresult = 1;
    960                        while (it.hasNext()) {
    961                                lresult *= a.getElementLongAbs(it.index);
    962                        }
    963                        return Long.valueOf(lresult);
    964                case Dataset.ARRAYINT8: case Dataset.ARRAYINT16:
    965                case Dataset.ARRAYINT32: case Dataset.ARRAYINT64:
    966                        lresult = 1;
    967                        is = a.getElementsPerItem();
    968                        lresults = new long[is];
    969                        for (int j = 0; j < is; j++) {
    970                                lresults[j] = 1;
    971                        }
    972                        while (it.hasNext()) {
    973                                for (int j = 0; j < is; j++)
    974                                        lresults[j] *= a.getElementLongAbs(it.index+j);
    975                        }
    976                        return lresults;
    977                case Dataset.FLOAT32: case Dataset.FLOAT64:
    978                        double dresult = 1.;
    979                        while (it.hasNext()) {
    980                                final double x = a.getElementDoubleAbs(it.index);
    981                                if (ignoreNaNs && Double.isNaN(x)) {
    982                                        continue;
    983                                }
    984                                if (ignoreInfs && Double.isInfinite(x)) {
    985                                        continue;
    986                                }
    987                                dresult *= x;
    988                                if (Double.isNaN(dresult)) {
    989                                        break;
    990                                }
    991                        }
    992                        return Double.valueOf(dresult);
    993                case Dataset.ARRAYFLOAT32:
    994                case Dataset.ARRAYFLOAT64:
    995                        is = a.getElementsPerItem();
    996                        double[] vals = new double[is];
    997                        dresults = new double[is];
    998                        for (int j = 0; j < is; j++) {
    999                                dresults[j] = 1.;
    1000                        }
    1001                        while (it.hasNext()) {
    1002                                boolean okay = true;
    1003                                for (int j = 0; j < is; j++) {
    1004                                        final double val = a.getElementDoubleAbs(it.index + j);
    1005                                        if (ignoreNaNs && Double.isNaN(val)) {
    1006                                                okay = false;
    1007                                                break;
    1008                                        }
    1009                                        if (ignoreInfs && Double.isInfinite(val)) {
    1010                                                okay = false;
    1011                                                break;
    1012                                        }
    1013                                        vals[j] = val;
    1014                                }
    1015                                if (okay) {
    1016                                        for (int j = 0; j < is; j++) {
    1017                                                double val = vals[j];
    1018                                                dresults[j] *= val;
    1019                                        }
    1020                                }
    1021                        }
    1022                        return dresults;
    1023                }
    1024
    1025                return null;
    1026        }
    1027
    1028        /**
    1029         * @param clazz dataset class
    1030         * @param a dataset
    1031         * @param axis
    1032         * @param ignoreInvalids see {@link IDataset#max(boolean...)}
    1033         * @return typed product of items along axis in dataset
    1034         * @since 2.3
    1035         */
    1036        public static <T extends Dataset> T typedProduct(final Class<T> clazz, final Dataset a, int axis, final boolean... ignoreInvalids) {
    1037                return typedProduct(clazz, a, new int[] {axis}, ignoreInvalids);
    1038        }
    1039
    1040        /**
    1041         * @param clazz dataset class
    1042         * @param a dataset
    1043         * @param axes
    1044         * @param ignoreInvalids see {@link IDataset#max(boolean...)}
    1045         * @return typed product of items in axes of dataset
    1046         * @since 2.3
    1047         */
    1048        @SuppressWarnings("unchecked")
    1049        public static <T extends Dataset> T typedProduct(final Class<T> clazz, final Dataset a, int[] axes, final boolean... ignoreInvalids) {
    1050                return (T) typedProduct(a, DTypeUtils.getDType(clazz), axes, ignoreInvalids);
    1051        }
    1052
    1053        /**
    1054         * @param a dataset
    1055         * @param dtype
    1056         * @param axis
    1057         * @param ignoreInvalids see {@link IDataset#max(boolean...)}
    1058         * @return typed product of items along axis in dataset
    1059         * @since 2.0
    1060         * @deprecated Please use the class-based method {@link #typedProduct(Class, Dataset, int, boolean...)}
    1061         */
    1062        @Deprecated
    1063        public static Dataset typedProduct(final Dataset a, final int dtype, int axis, final boolean... ignoreInvalids) {
    1064                return typedProduct(a, dtype, new int[] {axis}, ignoreInvalids);
    1065        }
    1066
    1067        /**
    1068         * @param a dataset
    1069         * @param dtype
    1070         * @param axes
    1071         * @param ignoreInvalids see {@link IDataset#max(boolean...)}
    1072         * @return typed product of items in axes of dataset
    1073         * @since 2.2
    1074         * @deprecated Please use the class-based method {@link #typedProduct(Class, Dataset, int[], boolean...)}
    1075         */
    1076        @Deprecated
    1077        public static Dataset typedProduct(final Dataset a, final int dtype, int[] axes, final boolean... ignoreInvalids) {
    1078                axes = ShapeUtils.checkAxes(a.getRank(), axes);
    1079                SliceNDIterator siter = new SliceNDIterator(new SliceND(a.getShapeRef()), axes);
    1080
    1081                int[] nshape = siter.getUsedSlice().getSourceShape();
    1082                final int is = a.getElementsPerItem();
    1083
    1084                final boolean ignoreNaNs;
    1085                final boolean ignoreInfs;
    1086                if (a.hasFloatingPointElements()) {
    1087                        ignoreNaNs = ignoreInvalids != null && ignoreInvalids.length > 0 ? ignoreInvalids[0] : false;
    1088                        ignoreInfs = ignoreInvalids != null && ignoreInvalids.length > 1 ? ignoreInvalids[1] : ignoreNaNs;
    1089                } else {
    1090                        ignoreNaNs = false;
    1091                        ignoreInfs = false;
    1092                }
    1093                Dataset result = DatasetFactory.zeros(is, nshape, dtype);
    1094
    1095                int[] spos = siter.getUsedPos();
    1096
    1097                // TODO add getLongArray(long[], int...) to CompoundDataset
    1098                final boolean isComplex = a.isComplex();
    1099                while (siter.hasNext()) {
    1100                        IndexIterator iter = a.getSliceIterator(siter.getCurrentSlice());
    1101                        final int[] pos = iter.getPos();
    1102
    1103                        if (isComplex) {
    1104                                double rv = 1, iv = 0;
    1105                                switch (dtype) {
    1106                                case Dataset.COMPLEX64:
    1107                                        ComplexFloatDataset af = (ComplexFloatDataset) a;
    1108                                        while (iter.hasNext()) {
    1109                                                final float r1 = af.getReal(pos);
    1110                                                final float i1 = af.getImag(pos);
    1111                                                if (ignoreNaNs && (Float.isNaN(r1) || Float.isNaN(i1))) {
    1112                                                        continue;
    1113                                                }
    1114                                                if (ignoreInfs  && (Float.isInfinite(r1) || Float.isInfinite(i1))) {
    1115                                                        continue;
    1116                                                }
    1117                                                final double tv = r1*rv - i1*iv;
    1118                                                iv = r1*iv + i1*rv;
    1119                                                rv = tv;
    1120                                                if (Double.isNaN(rv) && Double.isNaN(iv)) {
    1121                                                        break;
    1122                                                }
    1123                                        }
    1124                                        break;
    1125                                case Dataset.COMPLEX128:
    1126                                        ComplexDoubleDataset ad = (ComplexDoubleDataset) a;
    1127                                        while (iter.hasNext()) {
    1128                                                final double r1 = ad.getReal(pos);
    1129                                                final double i1 = ad.getImag(pos);
    1130                                                if (ignoreNaNs && (Double.isNaN(r1) || Double.isNaN(i1))) {
    1131                                                        continue;
    1132                                                }
    1133                                                if (ignoreInfs  && (Double.isInfinite(r1) || Double.isInfinite(i1))) {
    1134                                                        continue;
    1135                                                }
    1136                                                final double tv = r1*rv - i1*iv;
    1137                                                iv = r1*iv + i1*rv;
    1138                                                rv = tv;
    1139                                                if (Double.isNaN(rv) && Double.isNaN(iv)) {
    1140                                                        break;
    1141                                                }
    1142                                        }
    1143                                        break;
    1144                                }
    1145
    1146                                result.set(new Complex(rv, iv), spos);
    1147                        } else {
    1148                                final long[] lresults;
    1149                                final double[] dresults;
    1150
    1151                                switch (dtype) {
    1152                                case Dataset.BOOL:
    1153                                case Dataset.INT8: case Dataset.INT16:
    1154                                case Dataset.INT32: case Dataset.INT64:
    1155                                        long lresult = 1;
    1156                                        while (iter.hasNext()) {
    1157                                                lresult *= a.getInt(pos);
    1158                                        }
    1159                                        result.set(lresult, spos);
    1160                                        break;
    1161                                case Dataset.ARRAYINT8:
    1162                                        lresults = new long[is];
    1163                                        for (int k = 0; k < is; k++) {
    1164                                                lresults[k] = 1;
    1165                                        }
    1166                                        while (iter.hasNext()) {
    1167                                                final byte[] va = (byte[]) a.getObject(pos);
    1168                                                for (int k = 0; k < is; k++) {
    1169                                                        lresults[k] *= va[k];
    1170                                                }
    1171                                        }
    1172                                        result.set(lresults, spos);
    1173                                        break;
    1174                                case Dataset.ARRAYINT16:
    1175                                        lresults = new long[is];
    1176                                        for (int k = 0; k < is; k++) {
    1177                                                lresults[k] = 1;
    1178                                        }
    1179                                        while (iter.hasNext()) {
    1180                                                final short[] va = (short[]) a.getObject(pos);
    1181                                                for (int k = 0; k < is; k++) {
    1182                                                        lresults[k] *= va[k];
    1183                                                }
    1184                                        }
    1185                                        result.set(lresults, spos);
    1186                                        break;
    1187                                case Dataset.ARRAYINT32:
    1188                                        lresults = new long[is];
    1189                                        for (int k = 0; k < is; k++) {
    1190                                                lresults[k] = 1;
    1191                                        }
    1192                                        while (iter.hasNext()) {
    1193                                                final int[] va = (int[]) a.getObject(pos);
    1194                                                for (int k = 0; k < is; k++) {
    1195                                                        lresults[k] *= va[k];
    1196                                                }
    1197                                        }
    1198                                        result.set(lresults, spos);
    1199                                        break;
    1200                                case Dataset.ARRAYINT64:
    1201                                        lresults = new long[is];
    1202                                        for (int k = 0; k < is; k++) {
    1203                                                lresults[k] = 1;
    1204                                        }
    1205                                        while (iter.hasNext()) {
    1206                                                final long[] va = (long[]) a.getObject(pos);
    1207                                                for (int k = 0; k < is; k++) {
    1208                                                        lresults[k] *= va[k];
    1209                                                }
    1210                                        }
    1211                                        result.set(lresults, spos);
    1212                                        break;
    1213                                case Dataset.FLOAT32: case Dataset.FLOAT64:
    1214                                        double dresult = 1.;
    1215                                        while (iter.hasNext()) {
    1216                                                final double x = a.getElementDoubleAbs(iter.index);
    1217                                                if (ignoreNaNs && Double.isNaN(x)) {
    1218                                                        continue;
    1219                                                }
    1220                                                if (ignoreInfs && Double.isInfinite(x)) {
    1221                                                        continue;
    1222                                                }
    1223                                                dresult *= x;
    1224                                                if (Double.isNaN(dresult)) {
    1225                                                        break;
    1226                                                }
    1227                                        }
    1228                                        result.set(dresult, spos);
    1229                                        break;
    1230                                case Dataset.ARRAYFLOAT32: case Dataset.ARRAYFLOAT64:
    1231                                        CompoundDataset da = (CompoundDataset) a;
    1232                                        double[] dvalues = new double[is];
    1233                                        dresults = new double[is];
    1234                                        for (int k = 0; k < is; k++) {
    1235                                                dresults[k] = 1.;
    1236                                        }
    1237                                        while (iter.hasNext()) {
    1238                                                da.getDoubleArrayAbs(iter.index, dvalues);
    1239                                                boolean okay = true;
    1240                                                for (int k = 0; k < is; k++) {
    1241                                                        final double val = dvalues[k];
    1242                                                        if (ignoreNaNs && Double.isNaN(val)) {
    1243                                                                okay = false;
    1244                                                                break;
    1245                                                        }
    1246                                                        if (ignoreInfs && Double.isInfinite(val)) {
    1247                                                                okay = false;
    1248                                                                break;
    1249                                                        }
    1250                                                }
    1251                                                if (okay) {
    1252                                                        for (int k = 0; k < is; k++) {
    1253                                                                dresults[k] *= dvalues[k];
    1254                                                        }
    1255                                                }
    1256                                        }
    1257                                        result.set(dresults, spos);
    1258                                        break;
    1259                                }
    1260                        }
    1261                }
    1262
    1263//              result.setShape(ShapeUtils.squeezeShape(oshape, axes));
    1264                return result;
    1265        }
    1266
    1267        /**
    1268         * @param a dataset
    1269         * @param ignoreInvalids see {@link IDataset#max(boolean...)}
    1270         * @return cumulative product of items in flattened dataset
    1271         * @since 2.0
    1272         */
    1273        public static Dataset cumulativeProduct(final Dataset a, final boolean... ignoreInvalids) {
    1274                return cumulativeProduct(a.flatten(), 0, ignoreInvalids);
    1275        }
    1276
    1277        /**
    1278         * @param a dataset
    1279         * @param axis
    1280         * @param ignoreInvalids see {@link Dataset#max(int, boolean...)}
    1281         * @return cumulative product of items along axis in dataset
    1282         * @since 2.0
    1283         */
    1284        public static Dataset cumulativeProduct(final Dataset a, int axis, final boolean... ignoreInvalids) {
    1285                axis = a.checkAxis(axis);
    1286                int dtype = a.getDType();
    1287                int[] oshape = a.getShape();
    1288                int alen = oshape[axis];
    1289                oshape[axis] = 1;
    1290
    1291                final boolean ignoreNaNs;
    1292                final boolean ignoreInfs;
    1293                if (a.hasFloatingPointElements()) {
    1294                        ignoreNaNs = ignoreInvalids != null && ignoreInvalids.length > 0 ? ignoreInvalids[0] : false;
    1295                        ignoreInfs = ignoreInvalids != null && ignoreInvalids.length > 1 ? ignoreInvalids[1] : ignoreNaNs;
    1296                } else {
    1297                        ignoreNaNs = false;
    1298                        ignoreInfs = false;
    1299                }
    1300                Dataset result = DatasetFactory.zeros(a);
    1301                PositionIterator pi = result.getPositionIterator(axis);
    1302
    1303                int[] pos = pi.getPos();
    1304
    1305                while (pi.hasNext()) {
    1306
    1307                        if (a.isComplex()) {
    1308                                double rv = 1, iv = 0;
    1309                                switch (dtype) {
    1310                                case Dataset.COMPLEX64:
    1311                                        ComplexFloatDataset af = (ComplexFloatDataset) a;
    1312                                        ComplexFloatDataset rf = (ComplexFloatDataset) result;
    1313                                        for (int j = 0; j < alen; j++) {
    1314                                                if (!Double.isNaN(rv) || !Double.isNaN(iv)) {
    1315                                                        pos[axis] = j;
    1316                                                        final float r1 = af.getReal(pos);
    1317                                                        final float i1 = af.getImag(pos);
    1318                                                        if (ignoreNaNs && (Float.isNaN(r1) || Float.isNaN(i1))) {
    1319                                                                continue;
    1320                                                        }
    1321                                                        if (ignoreInfs  && (Float.isInfinite(r1) || Float.isInfinite(i1))) {
    1322                                                                continue;
    1323                                                        }
    1324                                                        final double tv = r1*rv - i1*iv;
    1325                                                        iv = r1*iv + i1*rv;
    1326                                                        rv = tv;
    1327                                                }
    1328                                                rf.set((float) rv, (float) iv, pos);
    1329                                        }
    1330                                        break;
    1331                                case Dataset.COMPLEX128:
    1332                                        ComplexDoubleDataset ad = (ComplexDoubleDataset) a;
    1333                                        ComplexDoubleDataset rd = (ComplexDoubleDataset) result;
    1334                                        for (int j = 0; j < alen; j++) {
    1335                                                if (!Double.isNaN(rv) || !Double.isNaN(iv)) {
    1336                                                        pos[axis] = j;
    1337                                                        final double r1 = ad.getReal(pos);
    1338                                                        final double i1 = ad.getImag(pos);
    1339                                                        if (ignoreNaNs && (Double.isNaN(r1) || Double.isNaN(i1))) {
    1340                                                                continue;
    1341                                                        }
    1342                                                        if (ignoreInfs  && (Double.isInfinite(r1) || Double.isInfinite(i1))) {
    1343                                                                continue;
    1344                                                        }
    1345                                                        final double tv = r1*rv - i1*iv;
    1346                                                        iv = r1*iv + i1*rv;
    1347                                                        rv = tv;
    1348                                                }
    1349                                                rd.set(rv, iv, pos);
    1350                                        }
    1351                                        break;
    1352                                }
    1353                        } else {
    1354                                final int is;
    1355                                final long[] lresults;
    1356                                final double[] dresults;
    1357
    1358                                switch (dtype) {
    1359                                case Dataset.BOOL:
    1360                                case Dataset.INT8: case Dataset.INT16:
    1361                                case Dataset.INT32: case Dataset.INT64:
    1362                                        long lresult = 1;
    1363                                        for (int j = 0; j < alen; j++) {
    1364                                                pos[axis] = j;
    1365                                                lresult *= a.getInt(pos);
    1366                                                result.set(lresult, pos);
    1367                                        }
    1368                                        break;
    1369                                case Dataset.ARRAYINT8:
    1370                                        is = a.getElementsPerItem();
    1371                                        lresults = new long[is];
    1372                                        for (int k = 0; k < is; k++) {
    1373                                                lresults[k] = 1;
    1374                                        }
    1375                                        for (int j = 0; j < alen; j++) {
    1376                                                pos[axis] = j;
    1377                                                final byte[] va = (byte[]) a.getObject(pos);
    1378                                                for (int k = 0; k < is; k++) {
    1379                                                        lresults[k] *= va[k];
    1380                                                }
    1381                                                result.set(lresults, pos);
    1382                                        }
    1383                                        break;
    1384                                case Dataset.ARRAYINT16:
    1385                                        is = a.getElementsPerItem();
    1386                                        lresults = new long[is];
    1387                                        for (int k = 0; k < is; k++) {
    1388                                                lresults[k] = 1;
    1389                                        }
    1390                                        for (int j = 0; j < alen; j++) {
    1391                                                pos[axis] = j;
    1392                                                final short[] va = (short[]) a.getObject(pos);
    1393                                                for (int k = 0; k < is; k++) {
    1394                                                        lresults[k] *= va[k];
    1395                                                }
    1396                                                result.set(lresults, pos);
    1397                                        }
    1398                                        break;
    1399                                case Dataset.ARRAYINT32:
    1400                                        is = a.getElementsPerItem();
    1401                                        lresults = new long[is];
    1402                                        for (int k = 0; k < is; k++) {
    1403                                                lresults[k] = 1;
    1404                                        }
    1405                                        for (int j = 0; j < alen; j++) {
    1406                                                pos[axis] = j;
    1407                                                final int[] va = (int[]) a.getObject(pos);
    1408                                                for (int k = 0; k < is; k++) {
    1409                                                        lresults[k] *= va[k];
    1410                                                }
    1411                                                result.set(lresults, pos);
    1412                                        }
    1413                                        break;
    1414                                case Dataset.ARRAYINT64:
    1415                                        is = a.getElementsPerItem();
    1416                                        lresults = new long[is];
    1417                                        for (int k = 0; k < is; k++) {
    1418                                                lresults[k] = 1;
    1419                                        }
    1420                                        for (int j = 0; j < alen; j++) {
    1421                                                pos[axis] = j;
    1422                                                final long[] va = (long[]) a.getObject(pos);
    1423                                                for (int k = 0; k < is; k++) {
    1424                                                        lresults[k] *= va[k];
    1425                                                }
    1426                                                result.set(lresults, pos);
    1427                                        }
    1428                                        break;
    1429                                case Dataset.FLOAT32: case Dataset.FLOAT64:
    1430                                        double dresult = 1.;
    1431                                        for (int j = 0; j < alen; j++) {
    1432                                                if (!Double.isNaN(dresult)) {
    1433                                                        pos[axis] = j;
    1434                                                        final double x = a.getDouble(pos);
    1435                                                        if (ignoreNaNs && Double.isNaN(x)) {
    1436                                                                continue;
    1437                                                        }
    1438                                                        if (ignoreInfs && Double.isInfinite(x)) {
    1439                                                                continue;
    1440                                                        }
    1441                                                        dresult *= x;
    1442                                                }
    1443                                                result.set(dresult, pos);
    1444                                        }
    1445                                        break;
    1446                                case Dataset.ARRAYFLOAT32: case Dataset.ARRAYFLOAT64:
    1447                                        is = a.getElementsPerItem();
    1448                                        CompoundDataset da = (CompoundDataset) a;
    1449                                        double[] dvalues = new double[is];
    1450                                        dresults = new double[is];
    1451                                        for (int k = 0; k < is; k++) {
    1452                                                dresults[k] = 1.;
    1453                                        }
    1454                                        for (int j = 0; j < alen; j++) {
    1455                                                pos[axis] = j;
    1456                                                da.getDoubleArray(dvalues, pos);
    1457                                                boolean okay = true;
    1458                                                for (int k = 0; k < is; k++) {
    1459                                                        final double val = dvalues[k];
    1460                                                        if (ignoreNaNs && Double.isNaN(val)) {
    1461                                                                okay = false;
    1462                                                                break;
    1463                                                        }
    1464                                                        if (ignoreInfs && Double.isInfinite(val)) {
    1465                                                                okay = false;
    1466                                                                break;
    1467                                                        }
    1468                                                }
    1469                                                if (okay) {
    1470                                                        for (int k = 0; k < is; k++) {
    1471                                                                dresults[k] *= dvalues[k];
    1472                                                        }
    1473                                                }
    1474                                                result.set(dresults, pos);
    1475                                        }
    1476                                        break;
    1477                                }
    1478                        }
    1479                }
    1480
    1481                return result;
    1482        }
    1483
    1484        /**
    1485         * @param a dataset
    1486         * @param ignoreInvalids see {@link IDataset#max(boolean...)}
    1487         * @return cumulative sum of items in flattened dataset
    1488         * @since 2.0
    1489         */
    1490        public static Dataset cumulativeSum(final Dataset a, final boolean... ignoreInvalids) {
    1491                return cumulativeSum(a.flatten(), 0, ignoreInvalids);
    1492        }
    1493
    1494        /**
    1495         * @param a dataset
    1496         * @param axis
    1497         * @param ignoreInvalids see {@link Dataset#max(int, boolean...)}
    1498         * @return cumulative sum of items along axis in dataset
    1499         * @since 2.0
    1500         */
    1501        public static Dataset cumulativeSum(final Dataset a, int axis, final boolean... ignoreInvalids) {
    1502                axis = a.checkAxis(axis);
    1503                int dtype = a.getDType();
    1504                int[] oshape = a.getShape();
    1505                int alen = oshape[axis];
    1506                oshape[axis] = 1;
    1507
    1508                final boolean ignoreNaNs;
    1509                final boolean ignoreInfs;
    1510                if (a.hasFloatingPointElements()) {
    1511                        ignoreNaNs = ignoreInvalids != null && ignoreInvalids.length > 0 ? ignoreInvalids[0] : false;
    1512                        ignoreInfs = ignoreInvalids != null && ignoreInvalids.length > 1 ? ignoreInvalids[1] : ignoreNaNs;
    1513                } else {
    1514                        ignoreNaNs = false;
    1515                        ignoreInfs = false;
    1516                }
    1517                Dataset result = DatasetFactory.zeros(a);
    1518                PositionIterator pi = result.getPositionIterator(axis);
    1519
    1520                int[] pos = pi.getPos();
    1521
    1522                while (pi.hasNext()) {
    1523
    1524                        if (a.isComplex()) {
    1525                                double rv = 0, iv = 0;
    1526                                switch (dtype) {
    1527                                case Dataset.COMPLEX64:
    1528                                        ComplexFloatDataset af = (ComplexFloatDataset) a;
    1529                                        ComplexFloatDataset rf = (ComplexFloatDataset) result;
    1530                                        for (int j = 0; j < alen; j++) {
    1531                                                if (!Double.isNaN(rv) || !Double.isNaN(iv)) {
    1532                                                        pos[axis] = j;
    1533                                                        final float r1 = af.getReal(pos);
    1534                                                        final float i1 = af.getImag(pos);
    1535                                                        if (ignoreNaNs && (Float.isNaN(r1) || Float.isNaN(i1))) {
    1536                                                                continue;
    1537                                                        }
    1538                                                        if (ignoreInfs  && (Float.isInfinite(r1) || Float.isInfinite(i1))) {
    1539                                                                continue;
    1540                                                        }
    1541                                                        rv += r1;
    1542                                                        iv += i1;
    1543                                                }
    1544                                                rf.set((float) rv, (float) iv, pos);
    1545                                        }
    1546                                        break;
    1547                                case Dataset.COMPLEX128:
    1548                                        ComplexDoubleDataset ad = (ComplexDoubleDataset) a;
    1549                                        ComplexDoubleDataset rd = (ComplexDoubleDataset) result;
    1550                                        for (int j = 0; j < alen; j++) {
    1551                                                if (!Double.isNaN(rv) || !Double.isNaN(iv)) {
    1552                                                        pos[axis] = j;
    1553                                                        final double r1 = ad.getReal(pos);
    1554                                                        final double i1 = ad.getImag(pos);
    1555                                                        if (ignoreNaNs && (Double.isNaN(r1) || Double.isNaN(i1))) {
    1556                                                                continue;
    1557                                                        }
    1558                                                        if (ignoreInfs  && (Double.isInfinite(r1) || Double.isInfinite(i1))) {
    1559                                                                continue;
    1560                                                        }
    1561                                                        rv += r1;
    1562                                                        iv += i1;
    1563                                                }
    1564                                                rd.set(rv, iv, pos);
    1565                                        }
    1566                                        break;
    1567                                }
    1568                        } else {
    1569                                final int is;
    1570                                final long[] lresults;
    1571                                final double[] dresults;
    1572
    1573                                switch (dtype) {
    1574                                case Dataset.BOOL:
    1575                                case Dataset.INT8: case Dataset.INT16:
    1576                                case Dataset.INT32: case Dataset.INT64:
    1577                                        long lresult = 0;
    1578                                        for (int j = 0; j < alen; j++) {
    1579                                                pos[axis] = j;
    1580                                                lresult += a.getInt(pos);
    1581                                                result.set(lresult, pos);
    1582                                        }
    1583                                        break;
    1584                                case Dataset.ARRAYINT8:
    1585                                        is = a.getElementsPerItem();
    1586                                        lresults = new long[is];
    1587                                        for (int j = 0; j < alen; j++) {
    1588                                                pos[axis] = j;
    1589                                                final byte[] va = (byte[]) a.getObject(pos);
    1590                                                for (int k = 0; k < is; k++) {
    1591                                                        lresults[k] += va[k];
    1592                                                }
    1593                                                result.set(lresults, pos);
    1594                                        }
    1595                                        break;
    1596                                case Dataset.ARRAYINT16:
    1597                                        is = a.getElementsPerItem();
    1598                                        lresults = new long[is];
    1599                                        for (int j = 0; j < alen; j++) {
    1600                                                pos[axis] = j;
    1601                                                final short[] va = (short[]) a.getObject(pos);
    1602                                                for (int k = 0; k < is; k++) {
    1603                                                        lresults[k] += va[k];
    1604                                                }
    1605                                                result.set(lresults, pos);
    1606                                        }
    1607                                        break;
    1608                                case Dataset.ARRAYINT32:
    1609                                        is = a.getElementsPerItem();
    1610                                        lresults = new long[is];
    1611                                        for (int j = 0; j < alen; j++) {
    1612                                                pos[axis] = j;
    1613                                                final int[] va = (int[]) a.getObject(pos);
    1614                                                for (int k = 0; k < is; k++) {
    1615                                                        lresults[k] += va[k];
    1616                                                }
    1617                                                result.set(lresults, pos);
    1618                                        }
    1619                                        break;
    1620                                case Dataset.ARRAYINT64:
    1621                                        is = a.getElementsPerItem();
    1622                                        lresults = new long[is];
    1623                                        for (int j = 0; j < alen; j++) {
    1624                                                pos[axis] = j;
    1625                                                final long[] va = (long[]) a.getObject(pos);
    1626                                                for (int k = 0; k < is; k++) {
    1627                                                        lresults[k] += va[k];
    1628                                                }
    1629                                                result.set(lresults, pos);
    1630                                        }
    1631                                        break;
    1632                                case Dataset.FLOAT32: case Dataset.FLOAT64:
    1633                                        double dresult = 0.;
    1634                                        for (int j = 0; j < alen; j++) {
    1635                                                if (!Double.isNaN(dresult)) {
    1636                                                        pos[axis] = j;
    1637                                                        final double x = a.getDouble(pos);
    1638                                                        if (ignoreNaNs && Double.isNaN(x)) {
    1639                                                                continue;
    1640                                                        }
    1641                                                        if (ignoreInfs && Double.isInfinite(x)) {
    1642                                                                continue;
    1643                                                        }
    1644                                                        dresult += x;
    1645                                                }
    1646                                                result.set(dresult, pos);
    1647                                        }
    1648                                        break;
    1649                                case Dataset.ARRAYFLOAT32: case Dataset.ARRAYFLOAT64:
    1650                                        is = a.getElementsPerItem();
    1651                                        CompoundDataset da = (CompoundDataset) a;
    1652                                        double[] dvalues = new double[is];
    1653                                        dresults = new double[is];
    1654                                        for (int j = 0; j < alen; j++) {
    1655                                                pos[axis] = j;
    1656                                                da.getDoubleArray(dvalues, pos);
    1657                                                boolean okay = true;
    1658                                                for (int k = 0; k < is; k++) {
    1659                                                        final double val = dvalues[k];
    1660                                                        if (ignoreNaNs && Double.isNaN(val)) {
    1661                                                                okay = false;
    1662                                                                break;
    1663                                                        }
    1664                                                        if (ignoreInfs && Double.isInfinite(val)) {
    1665                                                                okay = false;
    1666                                                                break;
    1667                                                        }
    1668                                                }
    1669                                                if (okay) {
    1670                                                        for (int k = 0; k < is; k++) {
    1671                                                                dresults[k] += dvalues[k];
    1672                                                        }
    1673                                                }
    1674                                                result.set(dresults, pos);
    1675                                        }
    1676                                        break;
    1677                                }
    1678                        }
    1679                }
    1680
    1681                return result;
    1682        }
    1683
    1684        /**
    1685         * @param a
    1686         * @return average deviation value of all items the dataset
    1687         */
    1688        public static Object averageDeviation(final Dataset a) {
    1689                final IndexIterator it = a.getIterator();
    1690                final int is = a.getElementsPerItem();
    1691
    1692                if (is == 1) {
    1693                        double mean = (Double) a.mean();
    1694                        double sum = 0.0;
    1695
    1696                        while (it.hasNext()) {
    1697                                sum += Math.abs(a.getElementDoubleAbs(it.index) - mean);
    1698                        }
    1699
    1700                        return sum / a.getSize();
    1701                }
    1702
    1703                double[] means = (double[]) a.mean();
    1704                double[] sums = new double[is];
    1705
    1706                while (it.hasNext()) {
    1707                        for (int j = 0; j < is; j++)
    1708                                sums[j] += Math.abs(a.getElementDoubleAbs(it.index + j) - means[j]);
    1709                }
    1710
    1711                double n = a.getSize();
    1712                for (int j = 0; j < is; j++)
    1713                        sums[j] /= n;
    1714
    1715                return sums;
    1716        }
    1717
    1718        /**
    1719         * The residual is the sum of squared differences
    1720         * @param a
    1721         * @param b
    1722         * @return residual value
    1723         */
    1724        public static double residual(final Dataset a, final Dataset b) {
    1725                return a.residual(b);
    1726        }
    1727
    1728        /**
    1729         * The residual is the sum of squared differences
    1730         * @param a
    1731         * @param b
    1732         * @param w
    1733         * @return residual value
    1734         */
    1735        public static double weightedResidual(final Dataset a, final Dataset b, final Dataset w) {
    1736                return a.residual(b, w, false);
    1737        }
    1738
    1739        /**
    1740         * Calculate approximate outlier values. These are defined as the values in the dataset
    1741         * that are approximately below and above the given thresholds - in terms of percentages
    1742         * of dataset size.
    1743         * <p>
    1744         * It approximates by limiting the number of items (given by length) used internally by
    1745         * data structures - the larger this is, the more accurate will those outlier values become.
    1746         * The actual thresholds used are returned in the array.
    1747         * <p>
    1748         * Also, the low and high values will be made distinct if possible by adjusting the thresholds
    1749         * @param a
    1750         * @param lo percentage threshold for lower limit
    1751         * @param hi percentage threshold for higher limit
    1752         * @param length maximum number of items used internally, if negative, then unlimited
    1753         * @return double array with low and high values, and low and high percentage thresholds
    1754         */
    1755        public static double[] outlierValues(final Dataset a, double lo, double hi, final int length) {
    1756                return outlierValues(a, null, true, lo, hi, length);
    1757        }
    1758
    1759        /**
    1760         * Calculate approximate outlier values. These are defined as the values in the dataset
    1761         * that are approximately below and above the given thresholds - in terms of percentages
    1762         * of dataset size.
    1763         * <p>
    1764         * It approximates by limiting the number of items (given by length) used internally by
    1765         * data structures - the larger this is, the more accurate will those outlier values become.
    1766         * The actual thresholds used are returned in the array.
    1767         * <p>
    1768         * Also, the low and high values will be made distinct if possible by adjusting the thresholds
    1769         * @param a
    1770         * @param mask can be null
    1771         * @param value value of mask to match to include for calculation
    1772         * @param lo percentage threshold for lower limit
    1773         * @param hi percentage threshold for higher limit
    1774         * @param length maximum number of items used internally, if negative, then unlimited
    1775         * @return double array with low and high values, and low and high percentage thresholds
    1776         * @since 2.1
    1777         */
    1778        public static double[] outlierValues(final Dataset a, final Dataset mask, final boolean value, double lo, double hi, final int length) {
    1779                if (lo <= 0 || hi <= 0 || lo >= hi || hi >= 100  || Double.isNaN(lo)|| Double.isNaN(hi)) {
    1780                        throw new IllegalArgumentException("Thresholds must be between (0,100) and in order");
    1781                }
    1782                final int size = a.getSize();
    1783                int nl = Math.max((int) ((lo*size)/100), 1);
    1784                if (length > 0 && nl > length)
    1785                        nl = length;
    1786                int nh = Math.max((int) (((100-hi)*size)/100), 1);
    1787                if (length > 0 && nh > length)
    1788                        nh = length;
    1789
    1790                IndexIterator it = mask == null ? a.getIterator() : a.getBooleanIterator(mask, value);
    1791                double[] results = Math.max(nl, nh) > 640 ? outlierValuesMap(a, it, nl, nh) : outlierValuesList(a, it, nl, nh);
    1792
    1793                results[2] = results[2]*100./size;
    1794                results[3] = 100. - results[3]*100./size;
    1795                return results;
    1796        }
    1797
    1798        static double[] outlierValuesMap(final Dataset a, final IndexIterator it, int nl, int nh) {
    1799                final TreeMap<Double, Integer> lMap = new TreeMap<Double, Integer>();
    1800                final TreeMap<Double, Integer> hMap = new TreeMap<Double, Integer>();
    1801
    1802                int ml = 0;
    1803                int mh = 0;
    1804                while (it.hasNext()) {
    1805                        Double x = a.getElementDoubleAbs(it.index);
    1806                        if (Double.isNaN(x)) {
    1807                                continue;
    1808                        }
    1809                        Integer i;
    1810                        if (ml == nl) {
    1811                                Double k = lMap.lastKey();
    1812                                if (x < k) {
    1813                                        i = lMap.get(k) - 1;
    1814                                        if (i == 0) {
    1815                                                lMap.remove(k);
    1816                                        } else {
    1817                                                lMap.put(k, i);
    1818                                        }
    1819                                        i = lMap.get(x);
    1820                                        if (i == null) {
    1821                                                lMap.put(x, 1);
    1822                                        } else {
    1823                                                lMap.put(x, i + 1);
    1824                                        }
    1825                                }
    1826                        } else {
    1827                                i = lMap.get(x);
    1828                                if (i == null) {
    1829                                        lMap.put(x, 1);
    1830                                } else {
    1831                                        lMap.put(x, i + 1);
    1832                                }
    1833                                ml++;
    1834                        }
    1835
    1836                        if (mh == nh) {
    1837                                Double k = hMap.firstKey();
    1838                                if (x > k) {
    1839                                        i = hMap.get(k) - 1;
    1840                                        if (i == 0) {
    1841                                                hMap.remove(k);
    1842                                        } else {
    1843                                                hMap.put(k, i);
    1844                                        }
    1845                                        i = hMap.get(x);
    1846                                        if (i == null) {
    1847                                                hMap.put(x, 1);
    1848                                        } else {
    1849                                                hMap.put(x, i+1);
    1850                                        }
    1851                                }
    1852                        } else {
    1853                                i = hMap.get(x);
    1854                                if (i == null) {
    1855                                        hMap.put(x, 1);
    1856                                } else {
    1857                                        hMap.put(x, i+1);
    1858                                }
    1859                                mh++;
    1860                        }
    1861                }
    1862
    1863                // Attempt to make values distinct
    1864                double lx = lMap.lastKey();
    1865                double hx = hMap.firstKey();
    1866                if (lx >= hx) {
    1867                        Double h = hMap.higherKey(lx);
    1868                        if (h != null) {
    1869                                hx = h;
    1870                                mh--;
    1871                        } else {
    1872                                Double l = lMap.lowerKey(hx);
    1873                                if (l != null) {
    1874                                        lx = l;
    1875                                        ml--;
    1876                                }
    1877                        }
    1878
    1879                }
    1880                return new double[] {lMap.lastKey(), hMap.firstKey(), ml, mh};
    1881        }
    1882
    1883        static double[] outlierValuesList(final Dataset a, final IndexIterator it, int nl, int nh) {
    1884                final List<Double> lList = new ArrayList<Double>(nl);
    1885                final List<Double> hList = new ArrayList<Double>(nh);
    1886
    1887                double lx = Double.POSITIVE_INFINITY;
    1888                double hx = Double.NEGATIVE_INFINITY;
    1889
    1890                while (it.hasNext()) {
    1891                        double x = a.getElementDoubleAbs(it.index);
    1892                        if (Double.isNaN(x)) {
    1893                                continue;
    1894                        }
    1895                        if (x < lx) {
    1896                                if (lList.size() == nl) {
    1897                                        lList.remove(lx);
    1898                                }
    1899                                lList.add(x);
    1900                                lx = Collections.max(lList);
    1901                        } else if (x == lx) {
    1902                                if (lList.size() < nl) {
    1903                                        lList.add(x);
    1904                                }
    1905                        }
    1906
    1907                        if (x > hx) {
    1908                                if (hList.size() == nh) {
    1909                                        hList.remove(hx);
    1910                                }
    1911                                hList.add(x);
    1912                                hx = Collections.min(hList);
    1913                        } else if (x == hx) {
    1914                                if (hList.size() < nh) {
    1915                                        hList.add(x);
    1916                                }
    1917                        }
    1918                }
    1919
    1920                nl = lList.size();
    1921                nh = hList.size();
    1922
    1923                // Attempt to make values distinct
    1924                if (lx >= hx) {
    1925                        Collections.sort(hList);
    1926                        for (double h : hList) {
    1927                                if (h > hx) {
    1928                                        hx = h;
    1929                                        break;
    1930                                }
    1931                                nh--;
    1932                        }
    1933                        if (lx >= hx) {
    1934                                Collections.sort(lList);
    1935                                Collections.reverse(lList);
    1936                                for (double l : lList) {
    1937                                        if (l < lx) {
    1938                                                lx = l;
    1939                                                break;
    1940                                        }
    1941                                        nl--;
    1942                                }
    1943                        }
    1944                }
    1945                return new double[] {lx, hx, nl, nh};
    1946        }
    1947
    1948        /**
    1949         * See {@link #covariance(Dataset a, Dataset b, boolean rowvar, boolean bias, Integer ddof)} with b = null, rowvar = true, bias = false and ddof = null.
    1950         * @param a
    1951         * @return covariance array of a
    1952         */
    1953        public static Dataset covariance(final Dataset a) {
    1954                return covariance(a, true, false, null);
    1955        }
    1956
    1957        /**
    1958         * See {@link #covariance(Dataset a, Dataset b, boolean rowvar, boolean bias, Integer ddof)} with b = null.
    1959         * @param a
    1960         * @return covariance array of a
    1961         * @since 2.0
    1962         */
    1963        public static Dataset covariance(final Dataset a,
    1964                        boolean rowvar, boolean bias, Integer ddof) {
    1965                return covariance(a, null, rowvar, bias, ddof);
    1966        }
    1967
    1968        /**
    1969         * See {@link #covariance(Dataset a, Dataset b, boolean rowvar, boolean bias, Integer ddof)} with b = null, rowvar = true, bias = false and ddof = null.
    1970         * @param a
    1971         * @return covariance array of a concatenated with b
    1972         */
    1973        public static Dataset covariance(final Dataset a, final Dataset b) {
    1974                return covariance(a, b, true, false, null);
    1975        }
    1976
    1977        /**
    1978         * Calculate the covariance matrix (array) of a concatenated with b. This
    1979         * method is directly based on the implementation in numpy (cov).
    1980         * @param a Array containing multiple variable and observations. Each row represents a variable, each column an observation.
    1981         * @param b An extra set of variables and observations. Must be of same type as a and have a compatible shape.
    1982         * @param rowvar When true (default), each row is a variable; when false each column is a variable.
    1983         * @param bias Default normalisation is (N - 1) - N is number of observations. If set true, normalisation is (N).
    1984         * @param ddof Default normalisation is (N - 1). If ddof is set, then normalisation is (N - ddof).
    1985         * @return covariance array of a concatenated with b
    1986         * @since 2.0
    1987         */
    1988        public static Dataset covariance (final Dataset a, final Dataset b,
    1989                        boolean rowvar, boolean bias, Integer ddof) {
    1990
    1991                //Create a working copy of the dataset & check its rank.
    1992                Dataset vars = a.clone();
    1993                if (a.getRank() == 1) {
    1994                        vars.setShape(1, a.getShapeRef()[0]);
    1995                }
    1996
    1997                //1D of variables, so consider rows as variables
    1998                if (vars.getShapeRef()[0] == 1) {
    1999                        rowvar = true;
    2000                }
    2001
    2002                //nr is the number of records; axis is index
    2003                int nr, axis;
    2004                if (rowvar) {
    2005                        nr = vars.getShapeRef()[1];
    2006                        axis = 0;
    2007                } else {
    2008                        nr = vars.getShapeRef()[0];
    2009                        axis = 1;
    2010                }
    2011
    2012                //Set the reduced degrees of freedom & normalisation factor
    2013                if (ddof == null) {
    2014                        if (bias == false) {
    2015                                ddof = 1;
    2016                        } else {
    2017                                ddof = 0;
    2018                        }
    2019                }
    2020                double norm_fact = nr - ddof;
    2021                if (norm_fact <= 0.) {
    2022                        //TODO Some sort of warning here?
    2023                        norm_fact = 0.;
    2024                }
    2025
    2026                //Concatenate additional set of variables with main set
    2027                if (b != null) {
    2028                        //Create a working copy of the dataset & check its rank.
    2029                        Dataset extraVars = b.clone();
    2030                        if (b.getRank() == 1) {
    2031                                extraVars.setShape(1, a.getShapeRef()[0]);
    2032                        }
    2033                        vars = DatasetUtils.concatenate(new Dataset[]{vars, extraVars}, axis);
    2034                }
    2035
    2036                //Calculate deviations & covariance matrix
    2037                Dataset varsMean = vars.mean(1-axis, false);
    2038                //-vars & varsMean need same shape (this is a hack!)
    2039                int[] meanShape = vars.getShape();
    2040                meanShape[1-axis] = 1;
    2041                varsMean.setShape(meanShape);
    2042                vars.isubtract(varsMean);
    2043                Dataset cov;
    2044                if (rowvar) {
    2045                        cov = Maths.divide(LinearAlgebra.dotProduct(vars, Maths.conjugate(vars.transpose())), norm_fact).squeeze();
    2046                } else {
    2047                        cov = Maths.divide(LinearAlgebra.dotProduct(vars.transpose(), Maths.conjugate(vars)), norm_fact).squeeze();
    2048                }
    2049                return cov;
    2050        }
    2051}