    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.util.ArrayList;
    016import java.util.Arrays;
    017import java.util.Collection;
    018import java.util.Iterator;
    019import java.util.List;
    020
    021import org.eclipse.january.dataset.Comparisons.Monotonicity;
    022import org.eclipse.january.dataset.internal.GeneratedMaths;
    023
    024/**
    025 * Mathematics class
    026 */
    027public class Maths extends GeneratedMaths {
    028        /**
    029         * Unwrap result from mathematical methods if necessary
    030         *
    031         * @param o
    032         * @param a
    033         * @return a dataset if a is a dataset or an object of the same class as o
    034         */
    035        public static Object unwrap(final Dataset o, final Object a) {
    036                return a instanceof Dataset ? o : o.getObjectAbs(o.getOffset());
    037        }
    038
    039        /**
    040         * Unwrap result from mathematical methods if necessary
    041         *
    042         * @param o
    043         * @param a
    044         * @return a dataset if either a and b are datasets or an object of the same
    045         *         class as o
    046         */
    047        public static Object unwrap(final Dataset o, final Object a, final Object b) {
    048                return (a instanceof Dataset || b instanceof Dataset) ? o : o.getObjectAbs(o.getOffset());
    049        }
    050
    051        /**
    052         * Unwrap result from mathematical methods if necessary
    053         *
    054         * @param o
    055         * @param a
    056         * @return a dataset if any inputs are datasets or an object of the same
    057         *         class as o
    058         */
    059        public static Object unwrap(final Dataset o, final Object... a) {
    060                boolean isAnyDataset = false;
    061                for (Object obj : a) {
    062                        if (obj instanceof Dataset) {
    063                                isAnyDataset = true;
    064                                break;
    065                        }
    066                }
    067                return isAnyDataset ? o : o.getObjectAbs(o.getOffset());
    068        }
    069
    070        /**
    071         * @param a
    072         * @param b
    073         * @return floor divide of a and b
    074         */
    075        public static Dataset floorDivide(final Object a, final Object b) {
    076                return floorDivide(a, b, null);
    077        }
    078
    079        /**
    080         * @param a
    081         * @param b
    082         * @param o
    083         *            output can be null - in which case, a new dataset is created
    084         * @return floor divide of a and b
    085         */
    086        public static Dataset floorDivide(final Object a, final Object b, final Dataset o) {
    087                return divideTowardsFloor(a, b, o).ifloor();
    088        }
    089
    090        /**
    091         * @param a
    092         * @param b
    093         * @return floor remainder of a and b
    094         */
    095        public static Dataset floorRemainder(final Object a, final Object b) {
    096                return floorRemainder(a, b, null);
    097        }
    098
    099        /**
    100         * @param a
    101         * @param b
    102         * @param o
    103         *            output can be null - in which case, a new dataset is created
    104         * @return floor remainder of a and b
    105         */
    106        public static Dataset floorRemainder(final Object a, final Object b, final Dataset o) {
    107                Dataset da = a instanceof Dataset ? (Dataset) a : DatasetFactory.createFromObject(a);
    108                Dataset db = b instanceof Dataset ? (Dataset) b : DatasetFactory.createFromObject(b);
    109                Dataset dq = floorDivide(da, db);
    110                dq.imultiply(db);
    111                return subtract(da, dq, o);
    112        }
    113
    114        /**
    115         * Find reciprocal from dataset
    116         *
    117         * @param a
    118         * @return reciprocal dataset
    119         */
    120        public static Dataset reciprocal(final Object a) {
    121                return reciprocal(a, null);
    122        }
    123
    124        /**
    125         * Find reciprocal from dataset
    126         *
    127         * @param a
    128         * @param o
    129         *            output can be null - in which case, a new dataset is created
    130         * @return reciprocal dataset
    131         */
    132        public static Dataset reciprocal(final Object a, final Dataset o) {
    133                final Dataset da = a instanceof Dataset ? (Dataset) a : DatasetFactory.createFromObject(a);
    134                return divide(1, da, o);
    135        }
    136
    137        /**
    138         * abs - absolute value of each element
    139         *
    140         * @param a
    141         * @return dataset
    142         */
    143        public static Dataset abs(final Object a) {
    144                return abs(a, null);
    145        }
    146
    147        /**
    148         * abs - absolute value of each element
    149         *
    150         * @param a
    151         * @param o
    152         *            output can be null - in which case, a new dataset is created
    153         * @return dataset
    154         */
    155        public static Dataset abs(final Object a, final Dataset o) {
    156                final Dataset da = a instanceof Dataset ? (Dataset) a : DatasetFactory.createFromObject(a);
    157                final SingleInputBroadcastIterator it = new SingleInputBroadcastIterator(da, o, true, true, false);
    158                final Dataset result = it.getOutput();
    159                final int is = result.getElementsPerItem();
    160                final int dt = result.getDType();
    161                final int as = da.getElementsPerItem();
    162                final boolean reset = result == o && is > 1;
    163
    164                switch (dt) {
    165                case Dataset.INT8:
    166                        final byte[] oi8data = ((ByteDataset) result).data;
    167                        it.setOutputDouble(false);
    168
    169                        while (it.hasNext()) {
    170                                oi8data[it.oIndex] = (byte) Math.abs(it.aLong);
    171                        }
    172                        break;
    173                case Dataset.INT16:
    174                        final short[] oi16data = ((ShortDataset) result).data;
    175                        it.setOutputDouble(false);
    176
    177                        while (it.hasNext()) {
    178                                oi16data[it.oIndex] = (short) Math.abs(it.aLong);
    179                        }
    180                        break;
    181                case Dataset.INT32:
    182                        final int[] oi32data = ((IntegerDataset) result).data;
    183                        it.setOutputDouble(false);
    184
    185                        while (it.hasNext()) {
    186                                oi32data[it.oIndex] = (int) Math.abs(it.aLong);
    187                        }
    188                        break;
    189                case Dataset.INT64:
    190                        final long[] oi64data = ((LongDataset) result).data;
    191                        it.setOutputDouble(false);
    192
    193                        while (it.hasNext()) {
    194                                oi64data[it.oIndex] = Math.abs(it.aLong);
    195                        }
    196                        break;
    197                case Dataset.ARRAYINT8:
    198                        final byte[] oai8data = ((CompoundByteDataset) result).data;
    199                        it.setOutputDouble(false);
    200
    201                        if (is == 1) {
    202                                while (it.hasNext()) {
    203                                        oai8data[it.oIndex] = (byte) Math.abs(it.aLong);
    204                                }
    205                        } else if (as == 1) {
    206                                while (it.hasNext()) {
    207                                        byte ox = (byte) Math.abs(it.aLong);
    208                                        for (int j = 0; j < is; j++) {
    209                                                oai8data[it.oIndex + j] = ox;
    210                                        }
    211                                }
    212                        } else {
    213                                while (it.hasNext()) {
    214                                        oai8data[it.oIndex] = (byte) Math.abs(it.aLong);
    215                                        for (int j = 1; j < is; j++) {
    216                                                oai8data[it.oIndex + j] = (byte) Math.abs(da.getElementLongAbs(it.aIndex + j));
    217                                        }
    218                                }
    219                        }
    220                        break;
    221                case Dataset.ARRAYINT16:
    222                        final short[] oai16data = ((CompoundShortDataset) result).data;
    223                        it.setOutputDouble(false);
    224
    225                        if (is == 1) {
    226                                while (it.hasNext()) {
    227                                        oai16data[it.oIndex] = (short) Math.abs(it.aLong);
    228                                }
    229                        } else if (as == 1) {
    230                                while (it.hasNext()) {
    231                                        short ox = (short) Math.abs(it.aLong);
    232                                        for (int j = 0; j < is; j++) {
    233                                                oai16data[it.oIndex + j] = ox;
    234                                        }
    235                                }
    236                        } else {
    237                                while (it.hasNext()) {
    238                                        oai16data[it.oIndex] = (short) Math.abs(it.aLong);
    239                                        for (int j = 1; j < is; j++) {
    240                                                oai16data[it.oIndex + j] = (short) Math.abs(da.getElementLongAbs(it.aIndex + j));
    241                                        }
    242                                }
    243                        }
    244                        break;
    245                case Dataset.ARRAYINT32:
    246                        final int[] oai32data = ((CompoundIntegerDataset) result).data;
    247                        it.setOutputDouble(false);
    248
    249                        if (is == 1) {
    250                                while (it.hasNext()) {
    251                                        oai32data[it.oIndex] = (int) Math.abs(it.aLong);
    252                                }
    253                        } else if (as == 1) {
    254                                while (it.hasNext()) {
    255                                        int ox = (int) Math.abs(it.aLong);
    256                                        for (int j = 0; j < is; j++) {
    257                                                oai32data[it.oIndex + j] = ox;
    258                                        }
    259                                }
    260                        } else {
    261                                while (it.hasNext()) {
    262                                        oai32data[it.oIndex] = (int) Math.abs(it.aLong);
    263                                        for (int j = 1; j < is; j++) {
    264                                                oai32data[it.oIndex + j] = (int) Math.abs(da.getElementLongAbs(it.aIndex + j));
    265                                        }
    266                                }
    267                        }
    268                        break;
    269                case Dataset.ARRAYINT64:
    270                        final long[] oai64data = ((CompoundLongDataset) result).data;
    271                        it.setOutputDouble(false);
    272
    273                        if (is == 1) {
    274                                while (it.hasNext()) {
    275                                        oai64data[it.oIndex] = Math.abs(it.aLong);
    276                                }
    277                        } else if (as == 1) {
    278                                while (it.hasNext()) {
    279                                        long ox = Math.abs(it.aLong);
    280                                        for (int j = 0; j < is; j++) {
    281                                                oai64data[it.oIndex + j] = ox;
    282                                        }
    283                                }
    284                        } else {
    285                                while (it.hasNext()) {
    286                                        oai64data[it.oIndex] = Math.abs(it.aLong);
    287                                        for (int j = 1; j < is; j++) {
    288                                                oai64data[it.oIndex + j] = Math.abs(da.getElementLongAbs(it.aIndex + j));
    289                                        }
    290                                }
    291                        }
    292                        break;
    293                case Dataset.FLOAT32:
    294                        final float[] of32data = ((FloatDataset) result).data;
    295                        if (as == 1) {
    296                                while (it.hasNext()) {
    297                                        of32data[it.oIndex] = (float) (Math.abs(it.aDouble));
    298                                }
    299                        } else {
    300                                while (it.hasNext()) {
    301                                        of32data[it.oIndex] = (float) (Math.hypot(it.aDouble, da.getElementDoubleAbs(it.aIndex + 1)));
    302                                }
    303                        }
    304                        break;
    305                case Dataset.FLOAT64:
    306                        final double[] of64data = ((DoubleDataset) result).data;
    307                        if (as == 1) {
    308                                while (it.hasNext()) {
    309                                        of64data[it.oIndex] = Math.abs(it.aDouble);
    310                                }
    311                        } else {
    312                                while (it.hasNext()) {
    313                                        of64data[it.oIndex] = Math.hypot(it.aDouble, da.getElementDoubleAbs(it.aIndex + 1));
    314                                }
    315                        }
    316                        break;
    317                case Dataset.ARRAYFLOAT32:
    318                        final float[] oaf32data = ((CompoundFloatDataset) result).data;
    319                        if (is == 1) {
    320                                while (it.hasNext()) {
    321                                        oaf32data[it.oIndex] = (float) (Math.abs(it.aDouble));
    322                                }
    323                        } else if (as == 1) {
    324                                while (it.hasNext()) {
    325                                        float ox = (float) (Math.abs(it.aDouble));
    326                                        for (int j = 0; j < is; j++) {
    327                                                oaf32data[it.oIndex + j] = ox;
    328                                        }
    329                                }
    330                        } else {
    331                                while (it.hasNext()) {
    332                                        oaf32data[it.oIndex] = (float) Math.abs(it.aDouble);
    333                                        for (int j = 1; j < is; j++) {
    334                                                oaf32data[it.oIndex + j] = (float) Math.abs(da.getElementDoubleAbs(it.aIndex + j));
    335                                        }
    336                                }
    337                        }
    338                        break;
    339                case Dataset.ARRAYFLOAT64:
    340                        final double[] oaf64data = ((CompoundDoubleDataset) result).data;
    341                        if (is == 1) {
    342                                while (it.hasNext()) {
    343                                        oaf64data[it.oIndex] = Math.abs(it.aDouble);
    344                                }
    345                        } else if (as == 1) {
    346                                while (it.hasNext()) {
    347                                        final double ix = it.aDouble;
    348                                        double ox = Math.abs(ix);
    349                                        for (int j = 0; j < is; j++) {
    350                                                oaf64data[it.oIndex + j] = ox;
    351                                        }
    352                                }
    353                        } else {
    354                                while (it.hasNext()) {
    355                                        oaf64data[it.oIndex] = Math.abs(it.aDouble);
    356                                        for (int j = 1; j < is; j++) {
    357                                                oaf64data[it.oIndex + j] = Math.abs(da.getElementDoubleAbs(it.aIndex + j));
    358                                        }
    359                                }
    360                        }
    361                        break;
    362                case Dataset.COMPLEX64:
    363                        final float[] oc64data = ((ComplexFloatDataset) result).data;
    364                        if (as == 1) {
    365                                while (it.hasNext()) {
    366                                        oc64data[it.oIndex] = (float) Math.abs(it.aDouble);
    367                                        if (reset) {
    368                                                oc64data[it.oIndex + 1] = 0;
    369                                        }
    370                                }
    371                        } else {
    372                                while (it.hasNext()) {
    373                                        oc64data[it.oIndex] = (float) Math.hypot(it.aDouble, da.getElementDoubleAbs(it.aIndex + 1));
    374                                        if (reset) {
    375                                                oc64data[it.oIndex + 1] = 0;
    376                                        }
    377                                }
    378                        }
    379                        break;
    380                case Dataset.COMPLEX128:
    381                        final double[] oc128data = ((ComplexDoubleDataset) result).data;
    382                        if (as == 1) {
    383                                while (it.hasNext()) {
    384                                        oc128data[it.oIndex] = Math.abs(it.aDouble);
    385                                        if (reset) {
    386                                                oc128data[it.oIndex + 1] = 0;
    387                                        }
    388                                }
    389                        } else {
    390                                while (it.hasNext()) {
    391                                        oc128data[it.oIndex] = Math.hypot(it.aDouble, da.getElementDoubleAbs(it.aIndex + 1));
    392                                        if (reset) {
    393                                                oc128data[it.oIndex + 1] = 0;
    394                                        }
    395                                }
    396                        }
    397                        break;
    398                default:
    399                        throw new IllegalArgumentException(
    400                                        "abs supports integer, compound integer, real, compound real, complex datasets only");
    401                }
    402
    403                addFunctionName(result, "abs");
    404                return result;
    405        }
    406
    407        /**
    408         * @param a
    409         * @return a^*, complex conjugate of a
    410         */
    411        public static Dataset conjugate(final Object a) {
    412                return conjugate(a, null);
    413        }
    414
    415        /**
    416         * @param a
    417         * @param o
    418         *            output can be null - in which case, a new dataset is created
    419         * @return a^*, complex conjugate of a
    420         */
    421        public static Dataset conjugate(final Object a, final Dataset o) {
    422                final Dataset da = a instanceof Dataset ? (Dataset) a : DatasetFactory.createFromObject(a);
    423                int at = da.getDType();
    424                IndexIterator it1 = da.getIterator();
    425
    426                SingleInputBroadcastIterator it = new SingleInputBroadcastIterator(da, o, true, true, true);
    427                Dataset result = it.getOutput();
    428
    429                switch (at) {
    430                case Dataset.COMPLEX64:
    431                        float[] c64data = ((ComplexFloatDataset) result).getData();
    432
    433                        for (int i = 0; it1.hasNext();) {
    434                                c64data[i++] = (float) da.getElementDoubleAbs(it1.index);
    435                                c64data[i++] = (float) -da.getElementDoubleAbs(it1.index + 1);
    436                        }
    437                        result.setName(Operations.bracketIfNecessary(da.getName()).append("^*").toString());
    438                        break;
    439                case Dataset.COMPLEX128:
    440                        double[] c128data = ((ComplexDoubleDataset) result).getData();
    441
    442                        for (int i = 0; it1.hasNext();) {
    443                                c128data[i++] = da.getElementDoubleAbs(it1.index);
    444                                c128data[i++] = -da.getElementDoubleAbs(it1.index + 1);
    445                        }
    446                        result.setName(Operations.bracketIfNecessary(da.getName()).append("^*").toString());
    447                        break;
    448                default:
    449                        result = da;
    450                }
    451
    452                return result;
    453        }
    454
    455        /**
    456         * @param a
    457         *            side of right-angled triangle
    458         * @param b
    459         *            side of right-angled triangle
    460         * @return hypotenuse of right-angled triangle: sqrt(a^2 + a^2)
    461         */
    462        public static Dataset hypot(final Object a, final Object b) {
    463                return hypot(a, b, null);
    464        }
    465
    466        /**
    467         * @param a
    468         *            side of right-angled triangle
    469         * @param b
    470         *            side of right-angled triangle
    471         * @param o
    472         *            output can be null - in which case, a new dataset is created
    473         * @return hypotenuse of right-angled triangle: sqrt(a^2 + a^2)
    474         */
    475        public static Dataset hypot(final Object a, final Object b, final Dataset o) {
    476                final Dataset da = a instanceof Dataset ? (Dataset) a : DatasetFactory.createFromObject(a);
    477                final Dataset db = b instanceof Dataset ? (Dataset) b : DatasetFactory.createFromObject(b);
    478
    479                final BroadcastIterator it = BroadcastIterator.createIterator(da, db, o, true);
    480                it.setOutputDouble(true);
    481                final Dataset result = it.getOutput();
    482                final int is = result.getElementsPerItem();
    483                final int as = da.getElementsPerItem();
    484                final int bs = db.getElementsPerItem();
    485                final int dt = result.getDType();
    486                switch (dt) {
    487                case Dataset.BOOL:
    488                        boolean[] bdata = ((BooleanDataset) result).getData();
    489
    490                        while (it.hasNext()) {
    491                                bdata[it.oIndex] = Math.hypot(it.aDouble, it.bDouble) != 0;
    492                        }
    493                        break;
    494                case Dataset.INT8:
    495                        byte[] i8data = ((ByteDataset) result).getData();
    496
    497                        while (it.hasNext()) {
    498                                i8data[it.oIndex] = (byte) toLong(Math.hypot(it.aDouble, it.bDouble));
    499                        }
    500                        break;
    501                case Dataset.INT16:
    502                        short[] i16data = ((ShortDataset) result).getData();
    503
    504                        while (it.hasNext()) {
    505                                i16data[it.oIndex] = (short) toLong(Math.hypot(it.aDouble, it.bDouble));
    506                        }
    507                        break;
    508                case Dataset.INT32:
    509                        int[] i32data = ((IntegerDataset) result).getData();
    510
    511                        while (it.hasNext()) {
    512                                i32data[it.oIndex] = (int) toLong(Math.hypot(it.aDouble, it.bDouble));
    513                        }
    514                        break;
    515                case Dataset.INT64:
    516                        long[] i64data = ((LongDataset) result).getData();
    517
    518                        while (it.hasNext()) {
    519                                i64data[it.oIndex] = toLong(Math.hypot(it.aDouble, it.bDouble));
    520                        }
    521                        break;
    522                case Dataset.FLOAT32:
    523                        float[] f32data = ((FloatDataset) result).getData();
    524
    525                        while (it.hasNext()) {
    526                                f32data[it.oIndex] = (float) Math.hypot(it.aDouble, it.bDouble);
    527                        }
    528                        break;
    529                case Dataset.FLOAT64:
    530                        double[] f64data = ((DoubleDataset) result).getData();
    531
    532                        while (it.hasNext()) {
    533                                f64data[it.oIndex] = Math.hypot(it.aDouble, it.bDouble);
    534                        }
    535                        break;
    536                case Dataset.ARRAYINT8:
    537                        byte[] ai8data = ((CompoundByteDataset) result).getData();
    538
    539                        if (is == 1) {
    540                                while (it.hasNext()) {
    541                                        ai8data[it.oIndex] = (byte) toLong(Math.hypot(it.aDouble, it.bDouble));
    542                                }
    543                        } else if (as == 1) {
    544                                while (it.hasNext()) {
    545                                        ai8data[it.oIndex] = (byte) toLong(Math.hypot(it.aDouble, it.bDouble));
    546                                        for (int j = 1; j < is; j++) {
    547                                                ai8data[it.oIndex
    548                                                                + j] = (byte) toLong(Math.hypot(it.aDouble, db.getElementDoubleAbs(it.bIndex + j)));
    549                                        }
    550                                }
    551                        } else if (bs == 1) {
    552                                while (it.hasNext()) {
    553                                        ai8data[it.oIndex] = (byte) toLong(Math.hypot(it.aDouble, it.bDouble));
    554                                        for (int j = 1; j < is; j++) {
    555                                                ai8data[it.oIndex
    556                                                                + j] = (byte) toLong(Math.hypot(da.getElementDoubleAbs(it.aIndex + j), it.bDouble));
    557                                        }
    558                                }
    559                        } else {
    560                                while (it.hasNext()) {
    561                                        ai8data[it.oIndex] = (byte) toLong(Math.hypot(it.aDouble, it.bDouble));
    562                                        for (int j = 1; j < is; j++) {
    563                                                ai8data[it.oIndex + j] = (byte) toLong(Math.hypot(da.getElementDoubleAbs(it.aIndex + j),
    564                                                                db.getElementDoubleAbs(it.bIndex + j)));
    565                                        }
    566                                }
    567                        }
    568                        break;
    569                case Dataset.ARRAYINT16:
    570                        short[] ai16data = ((CompoundShortDataset) result).getData();
    571
    572                        if (is == 1) {
    573                                while (it.hasNext()) {
    574                                        ai16data[it.oIndex] = (short) toLong(Math.hypot(it.aDouble, it.bDouble));
    575                                }
    576                        } else if (as == 1) {
    577                                while (it.hasNext()) {
    578                                        ai16data[it.oIndex] = (short) toLong(Math.hypot(it.aDouble, it.bDouble));
    579                                        for (int j = 1; j < is; j++) {
    580                                                ai16data[it.oIndex
    581                                                                + j] = (short) toLong(Math.hypot(it.aDouble, db.getElementDoubleAbs(it.bIndex + j)));
    582                                        }
    583                                }
    584                        } else if (bs == 1) {
    585                                while (it.hasNext()) {
    586                                        ai16data[it.oIndex] = (short) toLong(Math.hypot(it.aDouble, it.bDouble));
    587                                        for (int j = 1; j < is; j++) {
    588                                                ai16data[it.oIndex
    589                                                                + j] = (short) toLong(Math.hypot(da.getElementDoubleAbs(it.aIndex + j), it.bDouble));
    590                                        }
    591                                }
    592                        } else {
    593                                while (it.hasNext()) {
    594                                        ai16data[it.oIndex] = (short) toLong(Math.hypot(it.aDouble, it.bDouble));
    595                                        for (int j = 1; j < is; j++) {
    596                                                ai16data[it.oIndex + j] = (short) toLong(Math.hypot(da.getElementDoubleAbs(it.aIndex + j),
    597                                                                db.getElementDoubleAbs(it.bIndex + j)));
    598                                        }
    599                                }
    600                        }
    601                        break;
    602                case Dataset.ARRAYINT32:
    603                        int[] ai32data = ((CompoundIntegerDataset) result).getData();
    604
    605                        if (is == 1) {
    606                                while (it.hasNext()) {
    607                                        ai32data[it.oIndex] = (int) toLong(Math.hypot(it.aDouble, it.bDouble));
    608                                }
    609                        } else if (as == 1) {
    610                                while (it.hasNext()) {
    611                                        ai32data[it.oIndex] = (int) toLong(Math.hypot(it.aDouble, it.bDouble));
    612                                        for (int j = 1; j < is; j++) {
    613                                                ai32data[it.oIndex
    614                                                                + j] = (int) toLong(Math.hypot(it.aDouble, db.getElementDoubleAbs(it.bIndex + j)));
    615                                        }
    616                                }
    617                        } else if (bs == 1) {
    618                                while (it.hasNext()) {
    619                                        ai32data[it.oIndex] = (int) toLong(Math.hypot(it.aDouble, it.bDouble));
    620                                        for (int j = 1; j < is; j++) {
    621                                                ai32data[it.oIndex
    622                                                                + j] = (int) toLong(Math.hypot(da.getElementDoubleAbs(it.aIndex + j), it.bDouble));
    623                                        }
    624                                }
    625                        } else {
    626                                while (it.hasNext()) {
    627                                        ai32data[it.oIndex] = (int) toLong(Math.hypot(it.aDouble, it.bDouble));
    628                                        for (int j = 1; j < is; j++) {
    629                                                ai32data[it.oIndex + j] = (int) toLong(Math.hypot(da.getElementDoubleAbs(it.aIndex + j),
    630                                                                db.getElementDoubleAbs(it.bIndex + j)));
    631                                        }
    632                                }
    633                        }
    634                        break;
    635                case Dataset.ARRAYINT64:
    636                        long[] ai64data = ((CompoundLongDataset) result).getData();
    637
    638                        if (is == 1) {
    639                                while (it.hasNext()) {
    640                                        ai64data[it.oIndex] = toLong(Math.hypot(it.aDouble, it.bDouble));
    641                                }
    642                        } else if (as == 1) {
    643                                while (it.hasNext()) {
    644                                        ai64data[it.oIndex] = toLong(Math.hypot(it.aDouble, it.bDouble));
    645                                        for (int j = 1; j < is; j++) {
    646                                                ai64data[it.oIndex + j] = toLong(Math.hypot(it.aDouble, db.getElementDoubleAbs(it.bIndex + j)));
    647                                        }
    648                                }
    649                        } else if (bs == 1) {
    650                                while (it.hasNext()) {
    651                                        ai64data[it.oIndex] = toLong(Math.hypot(it.aDouble, it.bDouble));
    652                                        for (int j = 1; j < is; j++) {
    653                                                ai64data[it.oIndex + j] = toLong(Math.hypot(da.getElementDoubleAbs(it.aIndex + j), it.bDouble));
    654                                        }
    655                                }
    656                        } else {
    657                                while (it.hasNext()) {
    658                                        ai64data[it.oIndex] = toLong(Math.hypot(it.aDouble, it.bDouble));
    659                                        for (int j = 1; j < is; j++) {
    660                                                ai64data[it.oIndex + j] = toLong(Math.hypot(da.getElementDoubleAbs(it.aIndex + j),
    661                                                                db.getElementDoubleAbs(it.bIndex + j)));
    662                                        }
    663                                }
    664                        }
    665                        break;
    666                case Dataset.ARRAYFLOAT32:
    667                        float[] a32data = ((CompoundFloatDataset) result).getData();
    668
    669                        if (is == 1) {
    670                                while (it.hasNext()) {
    671                                        a32data[it.oIndex] = (float) Math.hypot(it.aDouble, it.bDouble);
    672                                }
    673                        } else if (as == 1) {
    674                                while (it.hasNext()) {
    675                                        a32data[it.oIndex] = (float) Math.hypot(it.aDouble, it.bDouble);
    676                                        for (int j = 1; j < is; j++) {
    677                                                a32data[it.oIndex + j] = (float) Math.hypot(it.aDouble, db.getElementDoubleAbs(it.bIndex + j));
    678                                        }
    679                                }
    680                        } else if (bs == 1) {
    681                                while (it.hasNext()) {
    682                                        a32data[it.oIndex] = (float) Math.hypot(it.aDouble, it.bDouble);
    683                                        for (int j = 1; j < is; j++) {
    684                                                a32data[it.oIndex + j] = (float) Math.hypot(da.getElementDoubleAbs(it.aIndex + j), it.bDouble);
    685                                        }
    686                                }
    687                        } else {
    688                                while (it.hasNext()) {
    689                                        a32data[it.oIndex] = (float) Math.hypot(it.aDouble, it.bDouble);
    690                                        for (int j = 1; j < is; j++) {
    691                                                a32data[it.oIndex + j] = (float) Math.hypot(da.getElementDoubleAbs(it.aIndex + j),
    692                                                                db.getElementDoubleAbs(it.bIndex + j));
    693                                        }
    694                                }
    695                        }
    696                        break;
    697                case Dataset.ARRAYFLOAT64:
    698                        double[] a64data = ((CompoundDoubleDataset) result).getData();
    699
    700                        if (is == 1) {
    701                                while (it.hasNext()) {
    702                                        a64data[it.oIndex] = Math.hypot(it.aDouble, it.bDouble);
    703                                }
    704                        } else if (as == 1) {
    705                                while (it.hasNext()) {
    706                                        a64data[it.oIndex] = Math.hypot(it.aDouble, it.bDouble);
    707                                        for (int j = 1; j < is; j++) {
    708                                                a64data[it.oIndex + j] = Math.hypot(it.aDouble, db.getElementDoubleAbs(it.bIndex + j));
    709                                        }
    710                                }
    711                        } else if (bs == 1) {
    712                                while (it.hasNext()) {
    713                                        a64data[it.oIndex] = Math.hypot(it.aDouble, it.bDouble);
    714                                        for (int j = 1; j < is; j++) {
    715                                                a64data[it.oIndex + j] = Math.hypot(da.getElementDoubleAbs(it.aIndex + j), it.bDouble);
    716                                        }
    717                                }
    718                        } else {
    719                                while (it.hasNext()) {
    720                                        a64data[it.oIndex] = Math.hypot(it.aDouble, it.bDouble);
    721                                        for (int j = 1; j < is; j++) {
    722                                                a64data[it.oIndex + j] = Math.hypot(da.getElementDoubleAbs(it.aIndex + j),
    723                                                                db.getElementDoubleAbs(it.bIndex + j));
    724                                        }
    725                                }
    726                        }
    727                        break;
    728                default:
    729                        throw new UnsupportedOperationException("hypot does not support this dataset type");
    730                }
    731
    732                addFunctionName(da, db, result, "hypot");
    733
    734                return result;
    735        }
    736
    737        /**
    738         * @param a
    739         *            opposite side of right-angled triangle
    740         * @param b
    741         *            adjacent side of right-angled triangle
    742         * @return angle of triangle: atan(b/a)
    743         */
    744        public static Dataset arctan2(final Object a, final Object b) {
    745                return arctan2(a, b, null);
    746        }
    747
    748        /**
    749         * @param a
    750         *            opposite side of right-angled triangle
    751         * @param b
    752         *            adjacent side of right-angled triangle
    753         * @param o
    754         *            output can be null - in which case, a new dataset is created
    755         * @return angle of triangle: atan(b/a)
    756         */
    757        public static Dataset arctan2(final Object a, final Object b, final Dataset o) {
    758                final Dataset da = a instanceof Dataset ? (Dataset) a : DatasetFactory.createFromObject(a);
    759                final Dataset db = b instanceof Dataset ? (Dataset) b : DatasetFactory.createFromObject(b);
    760
    761                final BroadcastIterator it = BroadcastIterator.createIterator(da, db, o, true);
    762                it.setOutputDouble(true);
    763                final Dataset result = it.getOutput();
    764                final int is = result.getElementsPerItem();
    765                final int as = da.getElementsPerItem();
    766                final int bs = db.getElementsPerItem();
    767                final int dt = result.getDType();
    768                switch (dt) {
    769                case Dataset.BOOL:
    770                        boolean[] bdata = ((BooleanDataset) result).getData();
    771
    772                        while (it.hasNext()) {
    773                                bdata[it.oIndex] = Math.atan2(it.aDouble, it.bDouble) != 0;
    774                        }
    775                        break;
    776                case Dataset.INT8:
    777                        byte[] i8data = ((ByteDataset) result).getData();
    778
    779                        while (it.hasNext()) {
    780                                i8data[it.oIndex] = (byte) toLong(Math.atan2(it.aDouble, it.bDouble));
    781                        }
    782                        break;
    783                case Dataset.INT16:
    784                        short[] i16data = ((ShortDataset) result).getData();
    785
    786                        while (it.hasNext()) {
    787                                i16data[it.oIndex] = (short) toLong(Math.atan2(it.aDouble, it.bDouble));
    788                        }
    789                        break;
    790                case Dataset.INT32:
    791                        int[] i32data = ((IntegerDataset) result).getData();
    792
    793                        while (it.hasNext()) {
    794                                i32data[it.oIndex] = (int) toLong(Math.atan2(it.aDouble, it.bDouble));
    795                        }
    796                        break;
    797                case Dataset.INT64:
    798                        long[] i64data = ((LongDataset) result).getData();
    799
    800                        while (it.hasNext()) {
    801                                i64data[it.oIndex] = toLong(Math.atan2(it.aDouble, it.bDouble));
    802                        }
    803                        break;
    804                case Dataset.FLOAT32:
    805                        float[] f32data = ((FloatDataset) result).getData();
    806
    807                        while (it.hasNext()) {
    808                                f32data[it.oIndex] = (float) Math.atan2(it.aDouble, it.bDouble);
    809                        }
    810                        break;
    811                case Dataset.FLOAT64:
    812                        double[] f64data = ((DoubleDataset) result).getData();
    813
    814                        while (it.hasNext()) {
    815                                f64data[it.oIndex] = Math.atan2(it.aDouble, it.bDouble);
    816                        }
    817                        break;
    818                case Dataset.ARRAYINT8:
    819                        byte[] ai8data = ((CompoundByteDataset) result).getData();
    820
    821                        if (is == 1) {
    822                                while (it.hasNext()) {
    823                                        ai8data[it.oIndex] = (byte) toLong(Math.atan2(it.aDouble, it.bDouble));
    824                                }
    825                        } else if (as == 1) {
    826                                while (it.hasNext()) {
    827                                        ai8data[it.oIndex] = (byte) toLong(Math.atan2(it.aDouble, it.bDouble));
    828                                        for (int j = 1; j < is; j++) {
    829                                                ai8data[it.oIndex
    830                                                                + j] = (byte) toLong(Math.atan2(it.aDouble, db.getElementDoubleAbs(it.bIndex + j)));
    831                                        }
    832                                }
    833                        } else if (bs == 1) {
    834                                while (it.hasNext()) {
    835                                        ai8data[it.oIndex] = (byte) toLong(Math.atan2(it.aDouble, it.bDouble));
    836                                        for (int j = 1; j < is; j++) {
    837                                                ai8data[it.oIndex
    838                                                                + j] = (byte) toLong(Math.atan2(da.getElementDoubleAbs(it.aIndex + j), it.bDouble));
    839                                        }
    840                                }
    841                        } else {
    842                                while (it.hasNext()) {
    843                                        ai8data[it.oIndex] = (byte) toLong(Math.atan2(it.aDouble, it.bDouble));
    844                                        for (int j = 1; j < is; j++) {
    845                                                ai8data[it.oIndex + j] = (byte) toLong(Math.atan2(da.getElementDoubleAbs(it.aIndex + j),
    846                                                                db.getElementDoubleAbs(it.bIndex + j)));
    847                                        }
    848                                }
    849                        }
    850                        break;
    851                case Dataset.ARRAYINT16:
    852                        short[] ai16data = ((CompoundShortDataset) result).getData();
    853
    854                        if (is == 1) {
    855                                while (it.hasNext()) {
    856                                        ai16data[it.oIndex] = (short) toLong(Math.atan2(it.aDouble, it.bDouble));
    857                                }
    858                        } else if (as == 1) {
    859                                while (it.hasNext()) {
    860                                        ai16data[it.oIndex] = (short) toLong(Math.atan2(it.aDouble, it.bDouble));
    861                                        for (int j = 1; j < is; j++) {
    862                                                ai16data[it.oIndex
    863                                                                + j] = (short) toLong(Math.atan2(it.aDouble, db.getElementDoubleAbs(it.bIndex + j)));
    864                                        }
    865                                }
    866                        } else if (bs == 1) {
    867                                while (it.hasNext()) {
    868                                        ai16data[it.oIndex] = (short) toLong(Math.atan2(it.aDouble, it.bDouble));
    869                                        for (int j = 1; j < is; j++) {
    870                                                ai16data[it.oIndex
    871                                                                + j] = (short) toLong(Math.atan2(da.getElementDoubleAbs(it.aIndex + j), it.bDouble));
    872                                        }
    873                                }
    874                        } else {
    875                                while (it.hasNext()) {
    876                                        ai16data[it.oIndex] = (short) toLong(Math.atan2(it.aDouble, it.bDouble));
    877                                        for (int j = 1; j < is; j++) {
    878                                                ai16data[it.oIndex + j] = (short) toLong(Math.atan2(da.getElementDoubleAbs(it.aIndex + j),
    879                                                                db.getElementDoubleAbs(it.bIndex + j)));
    880                                        }
    881                                }
    882                        }
    883                        break;
    884                case Dataset.ARRAYINT32:
    885                        int[] ai32data = ((CompoundIntegerDataset) result).getData();
    886
    887                        if (is == 1) {
    888                                while (it.hasNext()) {
    889                                        ai32data[it.oIndex] = (int) toLong(Math.atan2(it.aDouble, it.bDouble));
    890                                }
    891                        } else if (as == 1) {
    892                                while (it.hasNext()) {
    893                                        ai32data[it.oIndex] = (int) toLong(Math.atan2(it.aDouble, it.bDouble));
    894                                        for (int j = 1; j < is; j++) {
    895                                                ai32data[it.oIndex
    896                                                                + j] = (int) toLong(Math.atan2(it.aDouble, db.getElementDoubleAbs(it.bIndex + j)));
    897                                        }
    898                                }
    899                        } else if (bs == 1) {
    900                                while (it.hasNext()) {
    901                                        ai32data[it.oIndex] = (int) toLong(Math.atan2(it.aDouble, it.bDouble));
    902                                        for (int j = 1; j < is; j++) {
    903                                                ai32data[it.oIndex
    904                                                                + j] = (int) toLong(Math.atan2(da.getElementDoubleAbs(it.aIndex + j), it.bDouble));
    905                                        }
    906                                }
    907                        } else {
    908                                while (it.hasNext()) {
    909                                        ai32data[it.oIndex] = (int) toLong(Math.atan2(it.aDouble, it.bDouble));
    910                                        for (int j = 1; j < is; j++) {
    911                                                ai32data[it.oIndex + j] = (int) toLong(Math.atan2(da.getElementDoubleAbs(it.aIndex + j),
    912                                                                db.getElementDoubleAbs(it.bIndex + j)));
    913                                        }
    914                                }
    915                        }
    916                        break;
    917                case Dataset.ARRAYINT64:
    918                        long[] ai64data = ((CompoundLongDataset) result).getData();
    919
    920                        if (is == 1) {
    921                                while (it.hasNext()) {
    922                                        ai64data[it.oIndex] = toLong(Math.atan2(it.aDouble, it.bDouble));
    923                                }
    924                        } else if (as == 1) {
    925                                while (it.hasNext()) {
    926                                        ai64data[it.oIndex] = toLong(Math.atan2(it.aDouble, it.bDouble));
    927                                        for (int j = 1; j < is; j++) {
    928                                                ai64data[it.oIndex + j] = toLong(Math.atan2(it.aDouble, db.getElementDoubleAbs(it.bIndex + j)));
    929                                        }
    930                                }
    931                        } else if (bs == 1) {
    932                                while (it.hasNext()) {
    933                                        ai64data[it.oIndex] = toLong(Math.atan2(it.aDouble, it.bDouble));
    934                                        for (int j = 1; j < is; j++) {
    935                                                ai64data[it.oIndex + j] = toLong(Math.atan2(da.getElementDoubleAbs(it.aIndex + j), it.bDouble));
    936                                        }
    937                                }
    938                        } else {
    939                                while (it.hasNext()) {
    940                                        ai64data[it.oIndex] = toLong(Math.atan2(it.aDouble, it.bDouble));
    941                                        for (int j = 1; j < is; j++) {
    942                                                ai64data[it.oIndex + j] = toLong(Math.atan2(da.getElementDoubleAbs(it.aIndex + j),
    943                                                                db.getElementDoubleAbs(it.bIndex + j)));
    944                                        }
    945                                }
    946                        }
    947                        break;
    948                case Dataset.ARRAYFLOAT32:
    949                        float[] a32data = ((CompoundFloatDataset) result).getData();
    950
    951                        if (is == 1) {
    952                                while (it.hasNext()) {
    953                                        a32data[it.oIndex] = (float) Math.atan2(it.aDouble, it.bDouble);
    954                                }
    955                        } else if (as == 1) {
    956                                while (it.hasNext()) {
    957                                        a32data[it.oIndex] = (float) Math.atan2(it.aDouble, it.bDouble);
    958                                        for (int j = 1; j < is; j++) {
    959                                                a32data[it.oIndex + j] = (float) Math.atan2(it.aDouble, db.getElementDoubleAbs(it.bIndex + j));
    960                                        }
    961                                }
    962                        } else if (bs == 1) {
    963                                while (it.hasNext()) {
    964                                        a32data[it.oIndex] = (float) Math.atan2(it.aDouble, it.bDouble);
    965                                        for (int j = 1; j < is; j++) {
    966                                                a32data[it.oIndex + j] = (float) Math.atan2(da.getElementDoubleAbs(it.aIndex + j), it.bDouble);
    967                                        }
    968                                }
    969                        } else {
    970                                while (it.hasNext()) {
    971                                        a32data[it.oIndex] = (float) Math.atan2(it.aDouble, it.bDouble);
    972                                        for (int j = 1; j < is; j++) {
    973                                                a32data[it.oIndex + j] = (float) Math.atan2(da.getElementDoubleAbs(it.aIndex + j),
    974                                                                db.getElementDoubleAbs(it.bIndex + j));
    975                                        }
    976                                }
    977                        }
    978                        break;
    979                case Dataset.ARRAYFLOAT64:
    980                        double[] a64data = ((CompoundDoubleDataset) result).getData();
    981
    982                        if (is == 1) {
    983                                while (it.hasNext()) {
    984                                        a64data[it.oIndex] = Math.atan2(it.aDouble, it.bDouble);
    985                                }
    986                        } else if (as == 1) {
    987                                while (it.hasNext()) {
    988                                        a64data[it.oIndex] = Math.atan2(it.aDouble, it.bDouble);
    989                                        for (int j = 1; j < is; j++) {
    990                                                a64data[it.oIndex + j] = Math.atan2(it.aDouble, db.getElementDoubleAbs(it.bIndex + j));
    991                                        }
    992                                }
    993                        } else if (bs == 1) {
    994                                while (it.hasNext()) {
    995                                        a64data[it.oIndex] = Math.atan2(it.aDouble, it.bDouble);
    996                                        for (int j = 1; j < is; j++) {
    997                                                a64data[it.oIndex + j] = Math.atan2(da.getElementDoubleAbs(it.aIndex + j), it.bDouble);
    998                                        }
    999                                }
    1000                        } else {
    1001                                while (it.hasNext()) {
    1002                                        a64data[it.oIndex] = Math.atan2(it.aDouble, it.bDouble);
    1003                                        for (int j = 1; j < is; j++) {
    1004                                                a64data[it.oIndex + j] = Math.atan2(da.getElementDoubleAbs(it.aIndex + j),
    1005                                                                db.getElementDoubleAbs(it.bIndex + j));
    1006                                        }
    1007                                }
    1008                        }
    1009                        break;
    1010                default:
    1011                        throw new UnsupportedOperationException("atan2 does not support multiple-element dataset");
    1012                }
    1013
    1014                addFunctionName(da, db, result, "atan2");
    1015
    1016                return result;
    1017        }
    1018
    1019        /**
    1020         * Create a dataset of the arguments from a complex dataset
    1021         *
    1022         * @param a
    1023         * @return dataset of angles in radians
    1024         */
    1025        public static Dataset angle(final Object a) {
    1026                return angle(a, false, null);
    1027        }
    1028
    1029        /**
    1030         * Create a dataset of the arguments from a complex dataset
    1031         *
    1032         * @param a
    1033         * @param inDegrees
    1034         *            if true then return angles in degrees else in radians
    1035         * @return dataset of angles
    1036         */
    1037        public static Dataset angle(final Object a, final boolean inDegrees) {
    1038                return angle(a, inDegrees, null);
    1039        }
    1040
    1041        /**
    1042         * Create a dataset of the arguments from a complex dataset
    1043         *
    1044         * @param a
    1045         * @param o
    1046         *            output can be null - in which case, a new dataset is created
    1047         * @return dataset of angles in radians
    1048         */
    1049        public static Dataset angle(final Object a, final Dataset o) {
    1050                return angle(a, false, o);
    1051        }
    1052
    1053        /**
    1054         * Create a dataset of the arguments from a complex dataset
    1055         *
    1056         * @param a
    1057         * @param inDegrees
    1058         *            if true then return angles in degrees else in radians
    1059         * @param o
    1060         *            output can be null - in which case, a new dataset is created
    1061         * @return dataset of angles
    1062         */
    1063        public static Dataset angle(final Object a, final boolean inDegrees, final Dataset o) {
    1064                final Dataset da = a instanceof Dataset ? (Dataset) a : DatasetFactory.createFromObject(a);
    1065
    1066                if (!da.isComplex()) {
    1067                        throw new UnsupportedOperationException("angle does not support this dataset type");
    1068                }
    1069
    1070                final SingleInputBroadcastIterator it = new SingleInputBroadcastIterator(da, o, true, false, false);
    1071                final Dataset result = it.getOutput();
    1072                final int is = result.getElementsPerItem();
    1073                final int dt = result.getDType();
    1074
    1075                switch (dt) {
    1076                case Dataset.INT8:
    1077                        final byte[] oi8data = ((ByteDataset) result).data;
    1078                        it.setOutputDouble(false);
    1079
    1080                        if (inDegrees) {
    1081                                while (it.hasNext()) {
    1082                                        oi8data[it.oIndex] = (byte) toLong(
    1083                                                        Math.toDegrees(Math.atan2(da.getElementDoubleAbs(it.aIndex + 1), it.aDouble)));
    1084                                }
    1085                        } else {
    1086                                while (it.hasNext()) {
    1087                                        oi8data[it.oIndex] = (byte) toLong(Math.atan2(da.getElementDoubleAbs(it.aIndex + 1), it.aDouble));
    1088                                }
    1089                        }
    1090                        break;
    1091                case Dataset.INT16:
    1092                        final short[] oi16data = ((ShortDataset) result).data;
    1093                        it.setOutputDouble(false);
    1094
    1095                        if (inDegrees) {
    1096                                while (it.hasNext()) {
    1097                                        oi16data[it.oIndex] = (short) toLong(
    1098                                                        Math.toDegrees(Math.atan2(da.getElementDoubleAbs(it.aIndex + 1), it.aDouble)));
    1099                                }
    1100                        } else {
    1101                                while (it.hasNext()) {
    1102                                        oi16data[it.oIndex] = (short) toLong(Math.atan2(da.getElementDoubleAbs(it.aIndex + 1), it.aDouble));
    1103                                }
    1104                        }
    1105                        break;
    1106                case Dataset.INT32:
    1107                        final int[] oi32data = ((IntegerDataset) result).data;
    1108                        it.setOutputDouble(false);
    1109
    1110                        if (inDegrees) {
    1111                                while (it.hasNext()) {
    1112                                        oi32data[it.oIndex] = (int) toLong(
    1113                                                        Math.toDegrees(Math.atan2(da.getElementDoubleAbs(it.aIndex + 1), it.aDouble)));
    1114                                }
    1115                        } else {
    1116                                while (it.hasNext()) {
    1117                                        oi32data[it.oIndex] = (int) toLong(Math.atan2(da.getElementDoubleAbs(it.aIndex + 1), it.aDouble));
    1118                                }
    1119                        }
    1120                        break;
    1121                case Dataset.INT64:
    1122                        final long[] oi64data = ((LongDataset) result).data;
    1123                        it.setOutputDouble(false);
    1124
    1125                        if (inDegrees) {
    1126                                while (it.hasNext()) {
    1127                                        oi64data[it.oIndex] = toLong(
    1128                                                        Math.toDegrees(Math.atan2(da.getElementDoubleAbs(it.aIndex + 1), it.aDouble)));
    1129                                }
    1130                        } else {
    1131                                while (it.hasNext()) {
    1132                                        oi64data[it.oIndex] = toLong(Math.atan2(da.getElementDoubleAbs(it.aIndex + 1), it.aDouble));
    1133                                }
    1134                        }
    1135                        break;
    1136                case Dataset.ARRAYINT8:
    1137                        final byte[] oai8data = ((CompoundByteDataset) result).data;
    1138                        it.setOutputDouble(false);
    1139
    1140                        if (inDegrees) {
    1141                                while (it.hasNext()) {
    1142                                        final byte ox = (byte) toLong(
    1143                                                        Math.toDegrees(Math.atan2(da.getElementDoubleAbs(it.aIndex + 1), it.aDouble)));
    1144                                        for (int j = 0; j < is; j++) {
    1145                                                oai8data[it.oIndex + j] = ox;
    1146                                        }
    1147                                }
    1148                        } else {
    1149                                while (it.hasNext()) {
    1150                                        final byte ox = (byte) toLong(Math.atan2(da.getElementDoubleAbs(it.aIndex + 1), it.aDouble));
    1151                                        for (int j = 0; j < is; j++) {
    1152                                                oai8data[it.oIndex + j] = ox;
    1153                                        }
    1154                                }
    1155                        }
    1156                        break;
    1157                case Dataset.ARRAYINT16:
    1158                        final short[] oai16data = ((CompoundShortDataset) result).data;
    1159                        it.setOutputDouble(false);
    1160
    1161                        if (inDegrees) {
    1162                                while (it.hasNext()) {
    1163                                        final short ox = (short) toLong(
    1164                                                        Math.toDegrees(Math.atan2(da.getElementDoubleAbs(it.aIndex + 1), it.aDouble)));
    1165                                        for (int j = 0; j < is; j++) {
    1166                                                oai16data[it.oIndex + j] = ox;
    1167                                        }
    1168                                }
    1169                        } else {
    1170                                while (it.hasNext()) {
    1171                                        final short ox = (short) toLong(Math.atan2(da.getElementDoubleAbs(it.aIndex + 1), it.aDouble));
    1172                                        for (int j = 0; j < is; j++) {
    1173                                                oai16data[it.oIndex + j] = ox;
    1174                                        }
    1175                                }
    1176                        }
    1177                        break;
    1178                case Dataset.ARRAYINT32:
    1179                        final int[] oai32data = ((CompoundIntegerDataset) result).data;
    1180                        it.setOutputDouble(false);
    1181
    1182                        if (inDegrees) {
    1183                                while (it.hasNext()) {
    1184                                        final int ox = (int) toLong(
    1185                                                        Math.toDegrees(Math.atan2(da.getElementDoubleAbs(it.aIndex + 1), it.aDouble)));
    1186                                        for (int j = 0; j < is; j++) {
    1187                                                oai32data[it.oIndex + j] = ox;
    1188                                        }
    1189                                }
    1190                        } else {
    1191                                while (it.hasNext()) {
    1192                                        final int ox = (int) toLong(Math.atan2(da.getElementDoubleAbs(it.aIndex + 1), it.aDouble));
    1193                                        for (int j = 0; j < is; j++) {
    1194                                                oai32data[it.oIndex + j] = ox;
    1195                                        }
    1196                                }
    1197                        }
    1198                        break;
    1199                case Dataset.ARRAYINT64:
    1200                        final long[] oai64data = ((CompoundLongDataset) result).data;
    1201                        it.setOutputDouble(false);
    1202
    1203                        if (inDegrees) {
    1204                                while (it.hasNext()) {
    1205                                        final long ox = toLong(
    1206                                                        Math.toDegrees(Math.atan2(da.getElementDoubleAbs(it.aIndex + 1), it.aDouble)));
    1207                                        for (int j = 0; j < is; j++) {
    1208                                                oai64data[it.oIndex + j] = ox;
    1209                                        }
    1210                                }
    1211                        } else {
    1212                                while (it.hasNext()) {
    1213                                        final long ox = toLong(Math.atan2(da.getElementDoubleAbs(it.aIndex + 1), it.aDouble));
    1214                                        for (int j = 0; j < is; j++) {
    1215                                                oai64data[it.oIndex + j] = ox;
    1216                                        }
    1217                                }
    1218                        }
    1219                        break;
    1220                case Dataset.FLOAT32:
    1221                        final float[] of32data = ((FloatDataset) result).data;
    1222
    1223                        if (inDegrees) {
    1224                                while (it.hasNext()) {
    1225                                        of32data[it.oIndex] = (float) Math
    1226                                                        .toDegrees(Math.atan2(da.getElementDoubleAbs(it.aIndex + 1), it.aDouble));
    1227                                }
    1228                        } else {
    1229                                while (it.hasNext()) {
    1230                                        of32data[it.oIndex] = (float) Math.atan2(da.getElementDoubleAbs(it.aIndex + 1), it.aDouble);
    1231                                }
    1232                        }
    1233                        break;
    1234                case Dataset.FLOAT64:
    1235                        final double[] of64data = ((DoubleDataset) result).data;
    1236
    1237                        if (inDegrees) {
    1238                                while (it.hasNext()) {
    1239                                        of64data[it.oIndex] = Math.toDegrees(Math.atan2(da.getElementDoubleAbs(it.aIndex + 1), it.aDouble));
    1240                                }
    1241                        } else {
    1242                                while (it.hasNext()) {
    1243                                        of64data[it.oIndex] = Math.atan2(da.getElementDoubleAbs(it.aIndex + 1), it.aDouble);
    1244                                }
    1245                        }
    1246                        break;
    1247                case Dataset.ARRAYFLOAT32:
    1248                        final float[] oaf32data = ((CompoundFloatDataset) result).data;
    1249
    1250                        if (inDegrees) {
    1251                                while (it.hasNext()) {
    1252                                        final float ox = (float) Math
    1253                                                        .toDegrees(Math.atan2(da.getElementDoubleAbs(it.aIndex + 1), it.aDouble));
    1254                                        for (int j = 0; j < is; j++) {
    1255                                                oaf32data[it.oIndex + j] = ox;
    1256                                        }
    1257                                }
    1258                        } else {
    1259                                while (it.hasNext()) {
    1260                                        final float ox = (float) Math.atan2(da.getElementDoubleAbs(it.aIndex + 1), it.aDouble);
    1261                                        for (int j = 0; j < is; j++) {
    1262                                                oaf32data[it.oIndex + j] = ox;
    1263                                        }
    1264                                }
    1265                        }
    1266                        break;
    1267                case Dataset.ARRAYFLOAT64:
    1268                        final double[] oaf64data = ((CompoundDoubleDataset) result).data;
    1269
    1270                        if (inDegrees) {
    1271                                while (it.hasNext()) {
    1272                                        final double ox = Math.toDegrees(Math.atan2(da.getElementDoubleAbs(it.aIndex + 1), it.aDouble));
    1273                                        for (int j = 0; j < is; j++) {
    1274                                                oaf64data[it.oIndex + j] = ox;
    1275                                        }
    1276                                }
    1277                        } else {
    1278                                while (it.hasNext()) {
    1279                                        final double ox = Math.atan2(da.getElementDoubleAbs(it.aIndex + 1), it.aDouble);
    1280                                        for (int j = 0; j < is; j++) {
    1281                                                oaf64data[it.oIndex + j] = ox;
    1282                                        }
    1283                                }
    1284                        }
    1285                        break;
    1286                case Dataset.COMPLEX64:
    1287                        final float[] oc64data = ((ComplexFloatDataset) result).data;
    1288
    1289                        if (inDegrees) {
    1290                                while (it.hasNext()) {
    1291                                        oc64data[it.oIndex] = (float) Math
    1292                                                        .toDegrees(Math.atan2(da.getElementDoubleAbs(it.aIndex + 1), it.aDouble));
    1293                                        oc64data[it.oIndex + 1] = 0;
    1294                                }
    1295                        } else {
    1296                                while (it.hasNext()) {
    1297                                        oc64data[it.oIndex] = (float) Math.atan2(da.getElementDoubleAbs(it.aIndex + 1), it.aDouble);
    1298                                        oc64data[it.oIndex + 1] = 0;
    1299                                }
    1300                        }
    1301                        break;
    1302                case Dataset.COMPLEX128:
    1303                        final double[] oc128data = ((ComplexDoubleDataset) result).data;
    1304
    1305                        if (inDegrees) {
    1306                                while (it.hasNext()) {
    1307                                        oc128data[it.oIndex] = Math
    1308                                                        .toDegrees(Math.atan2(da.getElementDoubleAbs(it.aIndex + 1), it.aDouble));
    1309                                        oc128data[it.oIndex + 1] = 0;
    1310                                }
    1311                        } else {
    1312                                while (it.hasNext()) {
    1313                                        oc128data[it.oIndex] = Math.atan2(da.getElementDoubleAbs(it.aIndex + 1), it.aDouble);
    1314                                        oc128data[it.oIndex + 1] = 0;
    1315                                }
    1316                        }
    1317                        break;
    1318                default:
    1319                        throw new IllegalArgumentException("angle does not support this dataset type");
    1320                }
    1321
    1322                addFunctionName(result, "angle");
    1323
    1324                return result;
    1325        }
    1326
    1327        /**
    1328         * Create a phase only dataset. NB it will contain NaNs if there are any
    1329         * items with zero amplitude
    1330         *
    1331         * @param a
    1332         *            dataset
    1333         * @param keepZeros
    1334         *            if true then zero items are returned as zero rather than NaNs
    1335         * @return complex dataset where items have unit amplitude
    1336         */
    1337        public static Dataset phaseAsComplexNumber(final Object a, final boolean keepZeros) {
    1338                return phaseAsComplexNumber(a, null, keepZeros);
    1339        }
    1340
    1341        /**
    1342         * Create a phase only dataset. NB it will contain NaNs if there are any
    1343         * items with zero amplitude
    1344         *
    1345         * @param a
    1346         *            dataset
    1347         * @param o
    1348         *            output can be null - in which case, a new dataset is created
    1349         * @param keepZeros
    1350         *            if true then zero items are returned as zero rather than NaNs
    1351         * @return complex dataset where items have unit amplitude
    1352         */
    1353        public static Dataset phaseAsComplexNumber(final Object a, final Dataset o, final boolean keepZeros) {
    1354                final Dataset da = a instanceof Dataset ? (Dataset) a : DatasetFactory.createFromObject(a);
    1355
    1356                if (!da.isComplex()) {
    1357                        throw new IllegalArgumentException("Input dataset is not of complex type");
    1358                }
    1359                Dataset result = o == null ? DatasetFactory.zeros(da) : o;
    1360                if (!result.isComplex()) {
    1361                        throw new IllegalArgumentException("Output dataset is not of complex type");
    1362                }
    1363                final int dt = result.getDType();
    1364                SingleInputBroadcastIterator it = new SingleInputBroadcastIterator(da, result);
    1365
    1366                switch (dt) {
    1367                case Dataset.COMPLEX64:
    1368                        float[] z64data = ((ComplexFloatDataset) result).getData();
    1369
    1370                        if (keepZeros) {
    1371                                while (it.hasNext()) {
    1372                                        double rr = it.aDouble;
    1373                                        double ri = da.getElementDoubleAbs(it.aIndex + 1);
    1374                                        double am = Math.hypot(rr, ri);
    1375                                        if (am == 0) {
    1376                                                z64data[it.oIndex] = 0;
    1377                                                z64data[it.oIndex + 1] = 0;
    1378                                        } else {
    1379                                                z64data[it.oIndex] = (float) (rr / am);
    1380                                                z64data[it.oIndex + 1] = (float) (ri / am);
    1381                                        }
    1382                                }
    1383                        } else {
    1384                                while (it.hasNext()) {
    1385                                        double rr = it.aDouble;
    1386                                        double ri = da.getElementDoubleAbs(it.aIndex + 1);
    1387                                        double am = Math.hypot(rr, ri);
    1388                                        z64data[it.oIndex] = (float) (rr / am);
    1389                                        z64data[it.oIndex + 1] = (float) (ri / am);
    1390                                }
    1391                        }
    1392                        break;
    1393                case Dataset.COMPLEX128:
    1394                        double[] z128data = ((ComplexDoubleDataset) result).getData();
    1395
    1396                        if (keepZeros) {
    1397                                while (it.hasNext()) {
    1398                                        double rr = it.aDouble;
    1399                                        double ri = da.getElementDoubleAbs(it.aIndex + 1);
    1400                                        double am = Math.hypot(rr, ri);
    1401                                        if (am == 0) {
    1402                                                z128data[it.oIndex] = 0;
    1403                                                z128data[it.oIndex + 1] = 0;
    1404                                        } else {
    1405                                                z128data[it.oIndex] = rr / am;
    1406                                                z128data[it.oIndex + 1] = ri / am;
    1407                                        }
    1408                                }
    1409                        } else {
    1410                                while (it.hasNext()) {
    1411                                        double rr = it.aDouble;
    1412                                        double ri = da.getElementDoubleAbs(it.aIndex + 1);
    1413                                        double am = Math.hypot(rr, ri);
    1414                                        z128data[it.oIndex] = rr / am;
    1415                                        z128data[it.oIndex + 1] = ri / am;
    1416                                }
    1417                        }
    1418                        break;
    1419                }
    1420
    1421                addFunctionName(result, "phase");
    1422
    1423                return result;
    1424        }
    1425
    1426        /**
    1427         * Adds all sets passed in together
    1428         *
    1429         * The first IDataset must cast to Dataset
    1430         *
    1431         * For memory efficiency sake if add(...) is called with a set of size one,
    1432         * no clone is done, the original object is returned directly. Otherwise a
    1433         * new data set is returned, the sum of those passed in.
    1434         *
    1435         * @param sets
    1436         * @param requireClone
    1437         * @return sum of collection
    1438         */
    1439        public static Dataset add(final Collection<IDataset> sets, boolean requireClone) {
    1440                if (sets.isEmpty())
    1441                        return null;
    1442                final Iterator<IDataset> it = sets.iterator();
    1443                if (sets.size() == 1)
    1444                        return DatasetUtils.convertToDataset(it.next());
    1445
    1446                Dataset sum = requireClone ? ((Dataset) it.next()).clone() : (Dataset) it.next();
    1447
    1448                while (it.hasNext()) {
    1449                        add(sum, it.next(), sum);
    1450                }
    1451
    1452                return sum;
    1453        }
    1454
    1455        /**
    1456         * Multiplies all sets passed in together
    1457         *
    1458         * The first IDataset must cast to Dataset
    1459         *
    1460         * @param sets
    1461         * @param requireClone
    1462         * @return product of collection
    1463         */
    1464        public static Dataset multiply(final Collection<IDataset> sets, boolean requireClone) {
    1465                if (sets.isEmpty())
    1466                        return null;
    1467                final Iterator<IDataset> it = sets.iterator();
    1468                if (sets.size() == 1)
    1469                        return DatasetUtils.convertToDataset(it.next());
    1470                Dataset product = requireClone ? ((Dataset) it.next()).clone() : (Dataset) it.next();
    1471
    1472                while (it.hasNext()) {
    1473                        multiply(product, it.next(), product);
    1474                }
    1475
    1476                return product;
    1477        }
    1478
    1479        /**
    1480         * Linearly interpolate values at points in a 1D dataset corresponding to
    1481         * given coordinates.
    1482         *
    1483         * @param x
    1484         *            input 1-D coordinate dataset (must be ordered)
    1485         * @param d
    1486         *            input 1-D dataset
    1487         * @param x0
    1488         *            coordinate values
    1489         * @param left
    1490         *            value to use when x0 lies left of domain. If null, then
    1491         *            interpolate to zero by using leftmost interval
    1492         * @param right
    1493         *            value to use when x0 lies right of domain. If null, then
    1494         *            interpolate to zero by using rightmost interval
    1495         * @return interpolated values
    1496         */
    1497        public static Dataset interpolate(final Dataset x, final Dataset d, final IDataset x0, Number left, Number right) {
    1498                assert x.getRank() == 1;
    1499                assert d.getRank() == 1;
    1500
    1501                DoubleDataset r = DatasetFactory.zeros(DoubleDataset.class, x0.getShape());
    1502
    1503                Monotonicity mono = Comparisons.findMonotonicity(x);
    1504                if (mono == Monotonicity.NOT_ORDERED) {
    1505                        throw new IllegalArgumentException("Dataset x must be ordered");
    1506                }
    1507                DoubleDataset dx = (DoubleDataset) DatasetUtils.cast(x, Dataset.FLOAT64);
    1508                Dataset dx0 = DatasetUtils.convertToDataset(x0);
    1509                if (x == dx) {
    1510                        dx = (DoubleDataset) x.flatten();
    1511                }
    1512                double[] xa = dx.getData();
    1513                int s = xa.length - 1;
    1514                boolean isReversed = mono == Monotonicity.STRICTLY_DECREASING || mono == Monotonicity.NONINCREASING;
    1515                if (isReversed) {
    1516                        double[] txa = xa.clone();
    1517                        for (int i = 0; i <= s; i++) { // reverse order
    1518                                txa[s - i] = xa[i];
    1519                        }
    1520                        xa = txa;
    1521                }
    1522
    1523                IndexIterator it = dx0.getIterator();
    1524                int k = -1;
    1525                while (it.hasNext()) {
    1526                        k++;
    1527                        double v = dx0.getElementDoubleAbs(it.index);
    1528                        int i = Arrays.binarySearch(xa, v);
    1529                        if (i < 0) {
    1530                                // i = -(insertion point) - 1
    1531                                if (i == -1) {
    1532                                        if (left != null) {
    1533                                                r.setAbs(k, left.doubleValue());
    1534                                                continue;
    1535                                        }
    1536                                        final double d1 = xa[0] - xa[1];
    1537                                        double t = d1 - v + xa[0];
    1538                                        if (t >= 0)
    1539                                                continue; // sets to zero
    1540                                        t /= d1;
    1541                                        r.setAbs(k, t * d.getDouble(isReversed ? s : 0));
    1542                                } else if (i == -s - 2) {
    1543                                        if (right != null) {
    1544                                                r.setAbs(k, right.doubleValue());
    1545                                                continue;
    1546                                        }
    1547                                        final double d1 = xa[s] - xa[s - 1];
    1548                                        double t = d1 - v + xa[s];
    1549                                        if (t <= 0)
    1550                                                continue; // sets to zero
    1551                                        t /= d1;
    1552                                        r.setAbs(k, t * d.getDouble(isReversed ? 0 : s));
    1553                                } else {
    1554                                        i = -i - 1;
    1555                                        double t = (xa[i] - v) / (xa[i] - xa[i - 1]);
    1556                                        if (isReversed) {
    1557                                                i = s - i;
    1558                                                r.setAbs(k, t * d.getDouble(i + 1) + (1 - t) * d.getDouble(i));
    1559                                        } else {
    1560                                                r.setAbs(k, (1 - t) * d.getDouble(i) + t * d.getDouble(i - 1));
    1561                                        }
    1562                                }
    1563                        } else {
    1564                                r.setAbs(k, d.getDouble(isReversed ? s - i : i));
    1565                        }
    1566                }
    1567                return r;
    1568        }
    1569
    1570        /**
    1571         * Linearly interpolate a value at a point in a 1D dataset. The dataset is
    1572         * considered to have zero support outside its bounds. Thus points just
    1573         * outside are interpolated from the boundary value to zero.
    1574         *
    1575         * @param d
    1576         *            input dataset
    1577         * @param x0
    1578         *            coordinate
    1579         * @return interpolated value
    1580         */
    1581        public static double interpolate(final Dataset d, final double x0) {
    1582                assert d.getRank() == 1;
    1583
    1584                final int i0 = (int) Math.floor(x0);
    1585                final int e0 = d.getSize() - 1;
    1586                if (i0 < -1 || i0 > e0)
    1587                        return 0;
    1588
    1589                final double u0 = x0 - i0;
    1590
    1591                double r = 0;
    1592                final double f1 = i0 < 0 ? 0 : d.getDouble(i0);
    1593                if (u0 > 0) {
    1594                        r = (1 - u0) * f1 + (i0 == e0 ? 0 : u0 * d.getDouble(i0 + 1));
    1595                } else {
    1596                        r = f1;
    1597                }
    1598                return r;
    1599        }
    1600
    1601        /**
    1602         * Linearly interpolate a value at a point in a 1D dataset with a mask. The
    1603         * dataset is considered to have zero support outside its bounds. Thus
    1604         * points just outside are interpolated from the boundary value to zero.
    1605         *
    1606         * @param d
    1607         *            input dataset
    1608         * @param m
    1609         *            mask dataset
    1610         * @param x0
    1611         *            coordinate
    1612         * @return interpolated value
    1613         */
    1614        public static double interpolate(final Dataset d, final Dataset m, final double x0) {
    1615                assert d.getRank() == 1;
    1616                assert m.getRank() == 1;
    1617
    1618                final int i0 = (int) Math.floor(x0);
    1619                final int e0 = d.getSize() - 1;
    1620                if (i0 < -1 || i0 > e0)
    1621                        return 0;
    1622
    1623                final double u0 = x0 - i0;
    1624
    1625                double r = 0;
    1626                final double f1 = i0 < 0 ? 0 : d.getDouble(i0) * m.getDouble(i0);
    1627                if (u0 > 0) {
    1628                        r = (1 - u0) * f1 + (i0 == e0 ? 0 : u0 * d.getDouble(i0 + 1) * m.getDouble(i0 + 1));
    1629                } else {
    1630                        r = f1;
    1631                }
    1632                return r;
    1633        }
    1634
    1635        /**
    1636         * Linearly interpolate an array of values at a point in a compound 1D
    1637         * dataset. The dataset is considered to have zero support outside its
    1638         * bounds. Thus points just outside are interpolated from the boundary value
    1639         * to zero.
    1640         *
    1641         * @param values
    1642         *            interpolated array
    1643         * @param d
    1644         *            input dataset
    1645         * @param x0
    1646         *            coordinate
    1647         */
    1648        public static void interpolate(final double[] values, final CompoundDataset d, final double x0) {
    1649                assert d.getRank() == 1;
    1650
    1651                final int is = d.getElementsPerItem();
    1652                if (is != values.length)
    1653                        throw new IllegalArgumentException("Output array length must match elements in item");
    1654                final double[] f1, f2;
    1655
    1656                final int i0 = (int) Math.floor(x0);
    1657                final int e0 = d.getSize() - 1;
    1658                if (i0 < -1 || i0 > e0) {
    1659                        Arrays.fill(values, 0);
    1660                        return;
    1661                }
    1662                final double u0 = x0 - i0;
    1663
    1664                if (u0 > 0) {
    1665                        f1 = new double[is];
    1666                        if (i0 >= 0)
    1667                                d.getDoubleArray(f1, i0);
    1668                        double t = 1 - u0;
    1669                        if (i0 == e0) {
    1670                                for (int j = 0; j < is; j++)
    1671                                        values[j] = t * f1[j];
    1672                        } else {
    1673                                f2 = new double[is];
    1674                                d.getDoubleArray(f2, i0 + 1);
    1675                                for (int j = 0; j < is; j++)
    1676                                        values[j] = t * f1[j] + u0 * f2[j];
    1677                        }
    1678                } else {
    1679                        if (i0 >= 0)
    1680                                d.getDoubleArray(values, i0);
    1681                        else
    1682                                Arrays.fill(values, 0);
    1683                }
    1684        }
    1685
    1686        /**
    1687         * Linearly interpolate a value at a point in a 2D dataset. The dataset is
    1688         * considered to have zero support outside its bounds. Thus points just
    1689         * outside are interpolated from the boundary value to zero.
    1690         *
    1691         * @param d
    1692         *            input dataset
    1693         * @param x0
    1694         *            coordinate
    1695         * @param x1
    1696         *            coordinate
    1697         * @return bilinear interpolation
    1698         */
    1699        public static double interpolate(final Dataset d, final double x0, final double x1) {
    1700                final int[] s = d.getShape();
    1701                assert s.length == 2;
    1702
    1703                final int e0 = s[0] - 1;
    1704                final int e1 = s[1] - 1;
    1705                final int i0 = (int) Math.floor(x0);
    1706                final int i1 = (int) Math.floor(x1);
    1707                final double u0 = x0 - i0;
    1708                final double u1 = x1 - i1;
    1709                if (i0 < -1 || i0 > e0 || i1 < -1 || i1 > e1)
    1710                        return 0;
    1711
    1712                // use bilinear interpolation
    1713                double r = 0;
    1714                final double f1, f2, f3, f4;
    1715                f1 = i0 < 0 || i1 < 0 ? 0 : d.getDouble(i0, i1);
    1716                if (u1 > 0) {
    1717                        if (u0 > 0) {
    1718                                if (i0 == e0) {
    1719                                        f2 = 0;
    1720                                        f4 = 0;
    1721                                        f3 = i1 == e1 ? 0 : d.getDouble(i0, i1 + 1);
    1722                                } else {
    1723                                        f2 = i1 < 0 ? 0 : d.getDouble(i0 + 1, i1);
    1724                                        if (i1 == e1) {
    1725                                                f4 = 0;
    1726                                                f3 = 0;
    1727                                        } else {
    1728                                                f4 = d.getDouble(i0 + 1, i1 + 1);
    1729                                                f3 = i0 < 0 ? 0 : d.getDouble(i0, i1 + 1);
    1730                                        }
    1731                                }
    1732                                r = (1 - u0) * (1 - u1) * f1 + u0 * (1 - u1) * f2 + (1 - u0) * u1 * f3 + u0 * u1 * f4;
    1733                        } else {
    1734                                f3 = i0 < 0 || i1 == e1 ? 0 : d.getDouble(i0, i1 + 1);
    1735                                r = (1 - u1) * f1 + u1 * f3;
    1736                        }
    1737                } else { // exactly on axis 1
    1738                        if (u0 > 0) {
    1739                                f2 = i0 == e0 || i1 < 0 ? 0 : d.getDouble(i0 + 1, i1);
    1740                                r = (1 - u0) * f1 + u0 * f2;
    1741                        } else { // exactly on axis 0
    1742                                r = f1;
    1743                        }
    1744                }
    1745                return r;
    1746        }
    1747
    1748        /**
    1749         * Linearly interpolate a value at a point in a 2D dataset with a mask. The
    1750         * dataset is considered to have zero support outside its bounds. Thus
    1751         * points just outside are interpolated from the boundary value to zero.
    1752         *
    1753         * @param d
    1754         *            input dataset
    1755         * @param m
    1756         *            mask dataset
    1757         * @param x0
    1758         *            coordinate
    1759         * @param x1
    1760         *            coordinate
    1761         * @return bilinear interpolation
    1762         */
    1763        public static double interpolate(final Dataset d, final Dataset m, final double x0, final double x1) {
    1764                if (m == null)
    1765                        return interpolate(d, x0, x1);
    1766
    1767                final int[] s = d.getShape();
    1768                assert s.length == 2;
    1769                assert m.getRank() == 2;
    1770
    1771                final int e0 = s[0] - 1;
    1772                final int e1 = s[1] - 1;
    1773                final int i0 = (int) Math.floor(x0);
    1774                final int i1 = (int) Math.floor(x1);
    1775                final double u0 = x0 - i0;
    1776                final double u1 = x1 - i1;
    1777                if (i0 < -1 || i0 > e0 || i1 < -1 || i1 > e1)
    1778                        return 0;
    1779
    1780                // use bilinear interpolation
    1781                double r = 0;
    1782                final double f1, f2, f3, f4;
    1783                f1 = i0 < 0 || i1 < 0 ? 0 : d.getDouble(i0, i1) * m.getDouble(i0, i1);
    1784                if (u1 > 0) {
    1785                        if (i0 == e0) {
    1786                                f2 = 0;
    1787                                f4 = 0;
    1788                                f3 = i1 == e1 ? 0 : d.getDouble(i0, i1 + 1) * m.getDouble(i0, i1 + 1);
    1789                        } else {
    1790                                f2 = i1 < 0 ? 0 : d.getDouble(i0 + 1, i1) * m.getDouble(i0 + 1, i1);
    1791                                if (i1 == e1) {
    1792                                        f4 = 0;
    1793                                        f3 = 0;
    1794                                } else {
    1795                                        f4 = d.getDouble(i0 + 1, i1 + 1) * m.getDouble(i0 + 1, i1 + 1);
    1796                                        f3 = i0 < 0 ? 0 : d.getDouble(i0, i1 + 1) * m.getDouble(i0, i1 + 1);
    1797                                }
    1798                        }
    1799                        r = (1 - u0) * (1 - u1) * f1 + u0 * (1 - u1) * f2 + (1 - u0) * u1 * f3 + u0 * u1 * f4;
    1800                } else { // exactly on axis 1
    1801                        if (u0 > 0) {
    1802                                f2 = i0 == e0 || i1 < 0 ? 0 : d.getDouble(i0 + 1, i1) * m.getDouble(i0 + 1, i1);
    1803                                r = (1 - u0) * f1 + u0 * f2;
    1804                        } else { // exactly on axis 0
    1805                                r = f1;
    1806                        }
    1807                }
    1808                return r;
    1809        }
    1810
    1811        /**
    1812         * Linearly interpolate an array of values at a point in a compound 2D
    1813         * dataset. The dataset is considered to have zero support outside its
    1814         * bounds. Thus points just outside are interpolated from the boundary value
    1815         * to zero.
    1816         *
    1817         * @param values
    1818         *            bilinear interpolated array
    1819         * @param d
    1820         * @param x0
    1821         * @param x1
    1822         */
    1823        public static void interpolate(final double[] values, final CompoundDataset d, final double x0, final double x1) {
    1824                final int[] s = d.getShapeRef();
    1825                assert s.length == 2;
    1826
    1827                final int is = d.getElementsPerItem();
    1828                if (is != values.length)
    1829                        throw new IllegalArgumentException("Output array length must match elements in item");
    1830
    1831                final int e0 = s[0] - 1;
    1832                final int e1 = s[1] - 1;
    1833                final int i0 = (int) Math.floor(x0);
    1834                final int i1 = (int) Math.floor(x1);
    1835                final double u0 = x0 - i0;
    1836                final double u1 = x1 - i1;
    1837                if (i0 < -1 || i0 > e0 || i1 < -1 || i1 > e1) {
    1838                        Arrays.fill(values, 0);
    1839                        return;
    1840                }
    1841                // use bilinear interpolation
    1842                double[] f1 = new double[is];
    1843                if (i0 >= 0 && i1 >= 0)
    1844                        d.getDoubleArray(f1, i0, i1);
    1845
    1846                if (u1 > 0) {
    1847                        if (u0 > 0) {
    1848                                double[] f2 = new double[is];
    1849                                double[] f3 = new double[is];
    1850                                double[] f4 = new double[is];
    1851                                if (i0 != e0) {
    1852                                        if (i1 != e1)
    1853                                                d.getDoubleArray(f3, i0 + 1, i1 + 1);
    1854                                        if (i1 >= 0)
    1855                                                d.getDoubleArray(f4, i0 + 1, i1);
    1856                                }
    1857                                if (i0 >= 0 && i1 != e1)
    1858                                        d.getDoubleArray(f2, i0, i1 + 1);
    1859                                final double t0 = 1 - u0;
    1860                                final double t1 = 1 - u1;
    1861                                final double w1 = t0 * t1;
    1862                                final double w2 = t0 * u1;
    1863                                final double w3 = u0 * u1;
    1864                                final double w4 = u0 * t1;
    1865                                for (int j = 0; j < is; j++)
    1866                                        values[j] = w1 * f1[j] + w2 * f2[j] + w3 * f3[j] + w4 * f4[j];
    1867                        } else {
    1868                                double[] f2 = new double[is];
    1869                                if (i0 >= 0 && i1 != e1)
    1870                                        d.getDoubleArray(f2, i0, i1 + 1);
    1871                                final double t1 = 1 - u1;
    1872                                for (int j = 0; j < is; j++)
    1873                                        values[j] = t1 * f1[j] + u1 * f2[j];
    1874                        }
    1875                } else { // exactly on axis 1
    1876                        if (u0 > 0) {
    1877                                double[] f4 = new double[is];
    1878                                if (i0 != e0 && i1 >= 0)
    1879                                        d.getDoubleArray(f4, i0 + 1, i1);
    1880                                final double t0 = 1 - u0;
    1881                                for (int j = 0; j < is; j++)
    1882                                        values[j] = t0 * f1[j] + u0 * f4[j];
    1883                        } else { // exactly on axis 0
    1884                                if (i0 >= 0 && i1 >= 0)
    1885                                        d.getDoubleArray(values, i0, i1);
    1886                                else
    1887                                        Arrays.fill(values, 0);
    1888                        }
    1889                }
    1890        }
    1891
    1892        /**
    1893         * Linearly interpolate a value at a point in a n-D dataset. The dataset is
    1894         * considered to have zero support outside its bounds. Thus points just
    1895         * outside are interpolated from the boundary value to zero. The number of
    1896         * coordinates must match the rank of the dataset.
    1897         *
    1898         * @param d
    1899         *            input dataset
    1900         * @param x
    1901         *            coordinates
    1902         * @return interpolated value
    1903         */
    1904        public static double interpolate(final Dataset d, final double... x) {
    1905                return interpolate(d, null, x);
    1906        }
    1907
    1908        /**
    1909         * Linearly interpolate a value at a point in a n-D dataset with a mask. The
    1910         * dataset is considered to have zero support outside its bounds. Thus
    1911         * points just outside are interpolated from the boundary value to zero. The
    1912         * number of coordinates must match the rank of the dataset.
    1913         *
    1914         * @param d
    1915         *            input dataset
    1916         * @param m
    1917         *            mask dataset (can be null)
    1918         * @param x
    1919         *            coordinates
    1920         * @return interpolated value
    1921         */
    1922        public static double interpolate(final Dataset d, final Dataset m, final double... x) {
    1923                int r = d.getRank();
    1924                if (r != x.length) {
    1925                        throw new IllegalArgumentException("Number of coordinates must be equal to rank of dataset");
    1926                }
    1927
    1928                switch (r) {
    1929                case 1:
    1930                        return m == null ? interpolate(d, x[0]) : interpolate(d, m, x[0]);
    1931                case 2:
    1932                        return m == null ? interpolate(d, x[0], x[1]) : interpolate(d, m, x[0], x[1]);
    1933                }
    1934
    1935                if (m != null && r != m.getRank()) {
    1936                        throw new IllegalArgumentException("Rank of mask dataset must be equal to rank of dataset");
    1937                }
    1938
    1939                // now do it iteratively
    1940                int[] l = new int[r]; // lower indexes
    1941                double[] f = new double[r]; // fractions
    1942                for (int i = 0; i < r; i++) {
    1943                        double xi = x[i];
    1944                        l[i] = (int) Math.floor(xi);
    1945                        f[i] = xi - l[i];
    1946                }
    1947
    1948                int[] s = d.getShape();
    1949
    1950                int n = 1 << r;
    1951                double[] results = new double[n];
    1952
    1953                // iterate over permutations {l} and {l+1}
    1954                int[] twos = new int[r];
    1955                Arrays.fill(twos, 2);
    1956                PositionIterator it = new PositionIterator(twos);
    1957                int[] ip = it.getPos();
    1958                int j = 0;
    1959                if (m == null) {
    1960                        while (it.hasNext()) {
    1961                                int[] p = l.clone();
    1962                                boolean omit = false;
    1963                                for (int i = 0; i < r; i++) {
    1964                                        int pi = p[i] + ip[i];
    1965                                        if (pi < 0 || pi >= s[i]) {
    1966                                                omit = true;
    1967                                                break;
    1968                                        }
    1969                                        p[i] = pi;
    1970                                }
    1971                                results[j++] = omit ? 0 : d.getDouble(p);
    1972                        }
    1973                } else {
    1974                        while (it.hasNext()) {
    1975                                int[] p = l.clone();
    1976                                boolean omit = false;
    1977                                for (int i = 0; i < r; i++) {
    1978                                        int pi = p[i] + ip[i];
    1979                                        if (pi < 0 || pi >= s[i]) {
    1980                                                omit = true;
    1981                                                break;
    1982                                        }
    1983                                        p[i] = pi;
    1984                                }
    1985                                results[j++] = omit ? 0 : d.getDouble(p) * m.getDouble(p);
    1986                        }
    1987                }
    1988
    1989                // reduce recursively
    1990                for (int i = r - 1; i >= 0; i--) {
    1991                        results = combine(results, f[i], 1 << i);
    1992                }
    1993                return results[0];
    1994        }
    1995
    1996        private static double[] combine(double[] values, double f, int n) {
    1997                double g = 1 - f;
    1998                double[] results = new double[n];
    1999                for (int j = 0; j < n; j++) {
    2000                        int tj = 2 * j;
    2001                        results[j] = g * values[tj] + f * values[tj + 1];
    2002                }
    2003
    2004                return results;
    2005        }
    2006
    2007        /**
    2008         * Linearly interpolate an array of values at a point in a compound n-D
    2009         * dataset. The dataset is considered to have zero support outside its
    2010         * bounds. Thus points just outside are interpolated from the boundary value
    2011         * to zero.
    2012         *
    2013         * @param values
    2014         *            linearly interpolated array
    2015         * @param d
    2016         * @param x
    2017         */
    2018        public static void interpolate(final double[] values, final CompoundDataset d, final double... x) {
    2019                int r = d.getRank();
    2020                if (r != x.length) {
    2021                        throw new IllegalArgumentException("Number of coordinates must be equal to rank of dataset");
    2022                }
    2023
    2024                switch (r) {
    2025                case 1:
    2026                        interpolate(values, d, x[0]);
    2027                        return;
    2028                case 2:
    2029                        interpolate(values, d, x[0], x[1]);
    2030                        return;
    2031                }
    2032
    2033                final int is = d.getElementsPerItem();
    2034                if (is != values.length)
    2035                        throw new IllegalArgumentException("Output array length must match elements in item");
    2036
    2037                // now do it iteratively
    2038                int[] l = new int[r]; // lower indexes
    2039                double[] f = new double[r]; // fractions
    2040                for (int i = 0; i < r; i++) {
    2041                        double xi = x[i];
    2042                        l[i] = (int) Math.floor(xi);
    2043                        f[i] = xi - l[i];
    2044                }
    2045
    2046                int[] s = d.getShape();
    2047
    2048                int n = 1 << r;
    2049                double[][] results = new double[n][is];
    2050
    2051                // iterate over permutations {l} and {l+1}
    2052                int[] twos = new int[r];
    2053                Arrays.fill(twos, 2);
    2054                PositionIterator it = new PositionIterator(twos);
    2055                int[] ip = it.getPos();
    2056                int j = 0;
    2057                while (it.hasNext()) {
    2058                        int[] p = l.clone();
    2059                        boolean omit = false;
    2060                        for (int i = 0; i < r; i++) {
    2061                                int pi = p[i] + ip[i];
    2062                                if (pi < 0 || pi >= s[i]) {
    2063                                        omit = true;
    2064                                        break;
    2065                                }
    2066                                p[i] = pi;
    2067                        }
    2068                        if (!omit) {
    2069                                d.getDoubleArray(results[j++], p);
    2070                        }
    2071                }
    2072
    2073                // reduce recursively
    2074                for (int i = r - 1; i >= 0; i--) {
    2075                        results = combineArray(is, results, f[i], 1 << i);
    2076                }
    2077                for (int k = 0; k < is; k++) {
    2078                        values[k] = results[0][k];
    2079                }
    2080        }
    2081
    2082        private static double[][] combineArray(int is, double[][] values, double f, int n) {
    2083                double g = 1 - f;
    2084                double[][] results = new double[n][is];
    2085                for (int j = 0; j < n; j++) {
    2086                        int tj = 2 * j;
    2087                        for (int k = 0; k < is; k++) {
    2088                                results[j][k] = g * values[tj][k] + f * values[tj + 1][k];
    2089                        }
    2090                }
    2091
    2092                return results;
    2093        }
    2094
    2095        /**
    2096         * Linearly interpolate a value at a point in a 1D dataset. The dataset is
    2097         * considered to have zero support outside its bounds. Thus points just
    2098         * outside are interpolated from the boundary value to zero.
    2099         *
    2100         * @param d
    2101         *            input dataset
    2102         * @param x0
    2103         *            coordinate
    2104         * @return interpolated value
    2105         * @deprecated Use {@link #interpolate(Dataset, double)}
    2106         */
    2107        @Deprecated
    2108        public static double getLinear(final IDataset d, final double x0) {
    2109                return interpolate(DatasetUtils.convertToDataset(d), x0);
    2110        }
    2111
    2112        /**
    2113         * Linearly interpolate a value at a point in a compound 1D dataset. The
    2114         * dataset is considered to have zero support outside its bounds. Thus
    2115         * points just outside are interpolated from the boundary value to zero.
    2116         *
    2117         * @param values
    2118         *            interpolated array
    2119         * @param d
    2120         *            input dataset
    2121         * @param x0
    2122         *            coordinate
    2123         * @deprecated Use {@link #interpolate(double[], CompoundDataset, double)}
    2124         */
    2125        @Deprecated
    2126        public static void getLinear(final double[] values, final CompoundDataset d, final double x0) {
    2127                interpolate(values, d, x0);
    2128        }
    2129
    2130        /**
    2131         * Interpolated a value from 2D dataset
    2132         *
    2133         * @param d
    2134         *            input dataset
    2135         * @param x0
    2136         *            coordinate
    2137         * @param x1
    2138         *            coordinate
    2139         * @return bilinear interpolation
    2140         * @deprecated Use {@link #interpolate(Dataset, double, double)}
    2141         */
    2142        @Deprecated
    2143        public static double getBilinear(final IDataset d, final double x0, final double x1) {
    2144                return interpolate(DatasetUtils.convertToDataset(d), x0, x1);
    2145        }
    2146
    2147        /**
    2148         * Interpolated a value from 2D dataset with mask
    2149         *
    2150         * @param d
    2151         *            input dataset
    2152         * @param m
    2153         *            mask dataset
    2154         * @param x0
    2155         *            coordinate
    2156         * @param x1
    2157         *            coordinate
    2158         * @return bilinear interpolation
    2159         * @deprecated Use {@link #interpolate(Dataset, Dataset, double, double)}
    2160         */
    2161        @Deprecated
    2162        public static double getBilinear(final IDataset d, final IDataset m, final double x0, final double x1) {
    2163                return interpolate(DatasetUtils.convertToDataset(d), DatasetUtils.convertToDataset(m), x0, x1);
    2164        }
    2165
    2166        /**
    2167         * Interpolated a value from 2D compound dataset
    2168         *
    2169         * @param values
    2170         *            bilinear interpolated array
    2171         * @param d
    2172         * @param x0
    2173         * @param x1
    2174         * @deprecated Use
    2175         *             {@link #interpolate(double[], CompoundDataset, double, double)}
    2176         */
    2177        @Deprecated
    2178        public static void getBilinear(final double[] values, final CompoundDataset d, final double x0, final double x1) {
    2179                interpolate(values, d, x0, x1);
    2180        }
    2181
    2182        /**
    2183         * generate binomial coefficients with negative sign:
    2184         * <p>
    2185         *
    2186         * <pre>
    2187         *  (-1)^i n! / ( i! (n-i)! )
    2188         * </pre>
    2189         *
    2190         * @param n
    2191         * @return array of coefficients
    2192         */
    2193        private static int[] bincoeff(final int n) {
    2194                final int[] b = new int[n + 1];
    2195                final int hn = n / 2;
    2196
    2197                int bc = 1;
    2198                b[0] = bc;
    2199                for (int i = 1; i <= hn; i++) {
    2200                        bc = -(bc * (n - i + 1)) / i;
    2201                        b[i] = bc;
    2202                }
    2203                if (n % 2 != 0) {
    2204                        for (int i = hn + 1; i <= n; i++) {
    2205                                b[i] = -b[n - i];
    2206                        }
    2207                } else {
    2208                        for (int i = hn + 1; i <= n; i++) {
    2209                                b[i] = b[n - i];
    2210                        }
    2211                }
    2212                return b;
    2213        }
    2214
    2215        /**
    2216         * 1st order discrete difference of dataset along flattened dataset using
    2217         * finite difference
    2218         *
    2219         * @param a
    2220         *            is 1d dataset
    2221         * @param out
    2222         *            is 1d dataset
    2223         */
    2224        private static void difference(final Dataset a, final Dataset out) {
    2225                final int isize = a.getElementsPerItem();
    2226
    2227                final IndexIterator it = a.getIterator();
    2228                if (!it.hasNext())
    2229                        return;
    2230                int oi = it.index;
    2231
    2232                switch (a.getDType()) {
    2233                case Dataset.INT8:
    2234                        final byte[] i8data = ((ByteDataset) a).data;
    2235                        final byte[] oi8data = ((ByteDataset) out).getData();
    2236                        for (int i = 0; it.hasNext();) {
    2237                                oi8data[i++] = (byte) (i8data[it.index] - i8data[oi]);
    2238                                oi = it.index;
    2239                        }
    2240                        break;
    2241                case Dataset.INT16:
    2242                        final short[] i16data = ((ShortDataset) a).data;
    2243                        final short[] oi16data = ((ShortDataset) out).getData();
    2244                        for (int i = 0; it.hasNext();) {
    2245                                oi16data[i++] = (short) (i16data[it.index] - i16data[oi]);
    2246                                oi = it.index;
    2247                        }
    2248                        break;
    2249                case Dataset.INT32:
    2250                        final int[] i32data = ((IntegerDataset) a).data;
    2251                        final int[] oi32data = ((IntegerDataset) out).getData();
    2252                        for (int i = 0; it.hasNext();) {
    2253                                oi32data[i++] = i32data[it.index] - i32data[oi];
    2254                                oi = it.index;
    2255                        }
    2256                        break;
    2257                case Dataset.INT64:
    2258                        final long[] i64data = ((LongDataset) a).data;
    2259                        final long[] oi64data = ((LongDataset) out).getData();
    2260                        for (int i = 0; it.hasNext();) {
    2261                                oi64data[i++] = i64data[it.index] - i64data[oi];
    2262                                oi = it.index;
    2263                        }
    2264                        break;
    2265                case Dataset.ARRAYINT8:
    2266                        final byte[] ai8data = ((CompoundByteDataset) a).data;
    2267                        final byte[] oai8data = ((CompoundByteDataset) out).getData();
    2268                        for (int i = 0; it.hasNext();) {
    2269                                for (int k = 0; k < isize; k++) {
    2270                                        oai8data[i++] = (byte) (ai8data[it.index + k] - ai8data[oi++]);
    2271                                }
    2272                                oi = it.index;
    2273                        }
    2274                        break;
    2275                case Dataset.ARRAYINT16:
    2276                        final short[] ai16data = ((CompoundShortDataset) a).data;
    2277                        final short[] oai16data = ((CompoundShortDataset) out).getData();
    2278                        for (int i = 0; it.hasNext();) {
    2279                                for (int k = 0; k < isize; k++) {
    2280                                        oai16data[i++] = (short) (ai16data[it.index + k] - ai16data[oi++]);
    2281                                }
    2282                                oi = it.index;
    2283                        }
    2284                        break;
    2285                case Dataset.ARRAYINT32:
    2286                        final int[] ai32data = ((CompoundIntegerDataset) a).data;
    2287                        final int[] oai32data = ((CompoundIntegerDataset) out).getData();
    2288                        for (int i = 0; it.hasNext();) {
    2289                                for (int k = 0; k < isize; k++) {
    2290                                        oai32data[i++] = ai32data[it.index + k] - ai32data[oi++];
    2291                                }
    2292                                oi = it.index;
    2293                        }
    2294                        break;
    2295                case Dataset.ARRAYINT64:
    2296                        final long[] ai64data = ((CompoundLongDataset) a).data;
    2297                        final long[] oai64data = ((CompoundLongDataset) out).getData();
    2298                        for (int i = 0; it.hasNext();) {
    2299                                for (int k = 0; k < isize; k++) {
    2300                                        oai64data[i++] = ai64data[it.index + k] - ai64data[oi++];
    2301                                }
    2302                                oi = it.index;
    2303                        }
    2304                        break;
    2305                case Dataset.FLOAT32:
    2306                        final float[] f32data = ((FloatDataset) a).data;
    2307                        final float[] of32data = ((FloatDataset) out).getData();
    2308                        for (int i = 0; it.hasNext();) {
    2309                                of32data[i++] = f32data[it.index] - f32data[oi];
    2310                                oi = it.index;
    2311                        }
    2312                        break;
    2313                case Dataset.FLOAT64:
    2314                        final double[] f64data = ((DoubleDataset) a).data;
    2315                        final double[] of64data = ((DoubleDataset) out).getData();
    2316                        for (int i = 0; it.hasNext();) {
    2317                                of64data[i++] = f64data[it.index] - f64data[oi];
    2318                                oi = it.index;
    2319                        }
    2320                        break;
    2321                case Dataset.COMPLEX64:
    2322                        final float[] c64data = ((ComplexFloatDataset) a).data;
    2323                        final float[] oc64data = ((ComplexFloatDataset) out).getData();
    2324                        for (int i = 0; it.hasNext();) {
    2325                                oc64data[i++] = c64data[it.index] - c64data[oi];
    2326                                oc64data[i++] = c64data[it.index + 1] - c64data[oi + 1];
    2327                                oi = it.index;
    2328                        }
    2329                        break;
    2330                case Dataset.COMPLEX128:
    2331                        final double[] c128data = ((ComplexDoubleDataset) a).data;
    2332                        final double[] oc128data = ((ComplexDoubleDataset) out).getData();
    2333                        for (int i = 0; it.hasNext();) {
    2334                                oc128data[i++] = c128data[it.index] - c128data[oi];
    2335                                oc128data[i++] = c128data[it.index + 1] - c128data[oi + 1];
    2336                                oi = it.index;
    2337                        }
    2338                        break;
    2339                case Dataset.ARRAYFLOAT32:
    2340                        final float[] af32data = ((CompoundFloatDataset) a).data;
    2341                        final float[] oaf32data = ((CompoundFloatDataset) out).getData();
    2342                        for (int i = 0; it.hasNext();) {
    2343                                for (int k = 0; k < isize; k++) {
    2344                                        oaf32data[i++] = af32data[it.index + k] - af32data[oi++];
    2345                                }
    2346                                oi = it.index;
    2347                        }
    2348                        break;
    2349                case Dataset.ARRAYFLOAT64:
    2350                        final double[] af64data = ((CompoundDoubleDataset) a).data;
    2351                        final double[] oaf64data = ((CompoundDoubleDataset) out).getData();
    2352                        for (int i = 0; it.hasNext();) {
    2353                                for (int k = 0; k < isize; k++) {
    2354                                        oaf64data[i++] = af64data[it.index + k] - af64data[oi++];
    2355                                }
    2356                                oi = it.index;
    2357                        }
    2358                        break;
    2359                default:
    2360                        throw new UnsupportedOperationException("difference does not support this dataset type");
    2361                }
    2362        }
    2363
    2364        /**
    2365         * Get next set of indexes
    2366         *
    2367         * @param it
    2368         * @param indexes
    2369         * @return true if there is more
    2370         */
    2371        private static boolean nextIndexes(IndexIterator it, int[] indexes) {
    2372                if (!it.hasNext())
    2373                        return false;
    2374                int m = indexes.length;
    2375                int i = 0;
    2376                for (i = 0; i < m - 1; i++) {
    2377                        indexes[i] = indexes[i + 1];
    2378                }
    2379                indexes[i] = it.index;
    2380                return true;
    2381        }
    2382
    2383        /**
    2384         * General order discrete difference of dataset along flattened dataset
    2385         * using finite difference
    2386         *
    2387         * @param a
    2388         *            is 1d dataset
    2389         * @param out
    2390         *            is 1d dataset
    2391         * @param n
    2392         *            order of difference
    2393         */
    2394        private static void difference(final Dataset a, final Dataset out, final int n) {
    2395                if (n == 1) {
    2396                        difference(a, out);
    2397                        return;
    2398                }
    2399
    2400                final int isize = a.getElementsPerItem();
    2401
    2402                final int[] coeff = bincoeff(n);
    2403                final int m = n + 1;
    2404                final int[] indexes = new int[m]; // store for index values
    2405
    2406                final IndexIterator it = a.getIterator();
    2407                for (int i = 0; i < n; i++) {
    2408                        indexes[i] = it.index;
    2409                        it.hasNext();
    2410                }
    2411                indexes[n] = it.index;
    2412
    2413                switch (a.getDType()) {
    2414                case Dataset.INT8:
    2415                        final byte[] i8data = ((ByteDataset) a).data;
    2416                        final byte[] oi8data = ((ByteDataset) out).getData();
    2417                        for (int i = 0; nextIndexes(it, indexes);) {
    2418                                int ox = 0;
    2419                                for (int j = 0; j < m; j++) {
    2420                                        ox += i8data[indexes[j]] * coeff[j];
    2421                                }
    2422                                oi8data[i++] = (byte) ox;
    2423                        }
    2424                        break;
    2425                case Dataset.INT16:
    2426                        final short[] i16data = ((ShortDataset) a).data;
    2427                        final short[] oi16data = ((ShortDataset) out).getData();
    2428                        for (int i = 0; nextIndexes(it, indexes);) {
    2429                                int ox = 0;
    2430                                for (int j = 0; j < m; j++) {
    2431                                        ox += i16data[indexes[j]] * coeff[j];
    2432                                }
    2433                                oi16data[i++] = (short) ox;
    2434                        }
    2435                        break;
    2436                case Dataset.INT32:
    2437                        final int[] i32data = ((IntegerDataset) a).data;
    2438                        final int[] oi32data = ((IntegerDataset) out).getData();
    2439                        for (int i = 0; nextIndexes(it, indexes);) {
    2440                                int ox = 0;
    2441                                for (int j = 0; j < m; j++) {
    2442                                        ox += i32data[indexes[j]] * coeff[j];
    2443                                }
    2444                                oi32data[i++] = ox;
    2445                        }
    2446                        break;
    2447                case Dataset.INT64:
    2448                        final long[] i64data = ((LongDataset) a).data;
    2449                        final long[] oi64data = ((LongDataset) out).getData();
    2450                        for (int i = 0; nextIndexes(it, indexes);) {
    2451                                long ox = 0;
    2452                                for (int j = 0; j < m; j++) {
    2453                                        ox += i64data[indexes[j]] * coeff[j];
    2454                                }
    2455                                oi64data[i++] = ox;
    2456                        }
    2457                        break;
    2458                case Dataset.ARRAYINT8:
    2459                        final byte[] ai8data = ((CompoundByteDataset) a).data;
    2460                        final byte[] oai8data = ((CompoundByteDataset) out).getData();
    2461                        int[] box = new int[isize];
    2462                        for (int i = 0; nextIndexes(it, indexes);) {
    2463                                Arrays.fill(box, 0);
    2464                                for (int j = 0; j < m; j++) {
    2465                                        double c = coeff[j];
    2466                                        int l = indexes[j];
    2467                                        for (int k = 0; k < isize; k++) {
    2468                                                box[k] += ai8data[l++] * c;
    2469                                        }
    2470                                }
    2471                                for (int k = 0; k < isize; k++) {
    2472                                        oai8data[i++] = (byte) box[k];
    2473                                }
    2474                        }
    2475                        break;
    2476                case Dataset.ARRAYINT16:
    2477                        final short[] ai16data = ((CompoundShortDataset) a).data;
    2478                        final short[] oai16data = ((CompoundShortDataset) out).getData();
    2479                        int[] sox = new int[isize];
    2480                        for (int i = 0; nextIndexes(it, indexes);) {
    2481                                Arrays.fill(sox, 0);
    2482                                for (int j = 0; j < m; j++) {
    2483                                        double c = coeff[j];
    2484                                        int l = indexes[j];
    2485                                        for (int k = 0; k < isize; k++) {
    2486                                                sox[k] += ai16data[l++] * c;
    2487                                        }
    2488                                }
    2489                                for (int k = 0; k < isize; k++) {
    2490                                        oai16data[i++] = (short) sox[k];
    2491                                }
    2492                        }
    2493                        break;
    2494                case Dataset.ARRAYINT32:
    2495                        final int[] ai32data = ((CompoundIntegerDataset) a).data;
    2496                        final int[] oai32data = ((CompoundIntegerDataset) out).getData();
    2497                        int[] iox = new int[isize];
    2498                        for (int i = 0; nextIndexes(it, indexes);) {
    2499                                Arrays.fill(iox, 0);
    2500                                for (int j = 0; j < m; j++) {
    2501                                        double c = coeff[j];
    2502                                        int l = indexes[j];
    2503                                        for (int k = 0; k < isize; k++) {
    2504                                                iox[k] += ai32data[l++] * c;
    2505                                        }
    2506                                }
    2507                                for (int k = 0; k < isize; k++) {
    2508                                        oai32data[i++] = iox[k];
    2509                                }
    2510                        }
    2511                        break;
    2512                case Dataset.ARRAYINT64:
    2513                        final long[] ai64data = ((CompoundLongDataset) a).data;
    2514                        final long[] oai64data = ((CompoundLongDataset) out).getData();
    2515                        long[] lox = new long[isize];
    2516                        for (int i = 0; nextIndexes(it, indexes);) {
    2517                                Arrays.fill(lox, 0);
    2518                                for (int j = 0; j < m; j++) {
    2519                                        double c = coeff[j];
    2520                                        int l = indexes[j];
    2521                                        for (int k = 0; k < isize; k++) {
    2522                                                lox[k] += ai64data[l++] * c;
    2523                                        }
    2524                                }
    2525                                for (int k = 0; k < isize; k++) {
    2526                                        oai64data[i++] = lox[k];
    2527                                }
    2528                        }
    2529                        break;
    2530                case Dataset.FLOAT32:
    2531                        final float[] f32data = ((FloatDataset) a).data;
    2532                        final float[] of32data = ((FloatDataset) out).getData();
    2533                        for (int i = 0; nextIndexes(it, indexes);) {
    2534                                float ox = 0;
    2535                                for (int j = 0; j < m; j++) {
    2536                                        ox += f32data[indexes[j]] * coeff[j];
    2537                                }
    2538                                of32data[i++] = ox;
    2539                        }
    2540                        break;
    2541                case Dataset.FLOAT64:
    2542                        final double[] f64data = ((DoubleDataset) a).data;
    2543                        final double[] of64data = ((DoubleDataset) out).getData();
    2544                        for (int i = 0; nextIndexes(it, indexes);) {
    2545                                double ox = 0;
    2546                                for (int j = 0; j < m; j++) {
    2547                                        ox += f64data[indexes[j]] * coeff[j];
    2548                                }
    2549                                of64data[i++] = ox;
    2550                        }
    2551                        break;
    2552                case Dataset.COMPLEX64:
    2553                        final float[] c64data = ((ComplexFloatDataset) a).data;
    2554                        final float[] oc64data = ((ComplexFloatDataset) out).getData();
    2555                        for (int i = 0; nextIndexes(it, indexes);) {
    2556                                float ox = 0;
    2557                                float oy = 0;
    2558                                for (int j = 0; j < m; j++) {
    2559                                        int l = indexes[j];
    2560                                        ox += c64data[l++] * coeff[j];
    2561                                        oy += c64data[l] * coeff[j];
    2562                                }
    2563                                oc64data[i++] = ox;
    2564                                oc64data[i++] = oy;
    2565                        }
    2566                        break;
    2567                case Dataset.COMPLEX128:
    2568                        final double[] c128data = ((ComplexDoubleDataset) a).data;
    2569                        final double[] oc128data = ((ComplexDoubleDataset) out).getData();
    2570                        for (int i = 0; nextIndexes(it, indexes);) {
    2571                                double ox = 0;
    2572                                double oy = 0;
    2573                                for (int j = 0; j < m; j++) {
    2574                                        int l = indexes[j];
    2575                                        ox += c128data[l++] * coeff[j];
    2576                                        oy += c128data[l] * coeff[j];
    2577                                }
    2578                                oc128data[i++] = ox;
    2579                                oc128data[i++] = oy;
    2580                        }
    2581                        break;
    2582                case Dataset.ARRAYFLOAT32:
    2583                        final float[] af32data = ((CompoundFloatDataset) a).data;
    2584                        final float[] oaf32data = ((CompoundFloatDataset) out).getData();
    2585                        float[] fox = new float[isize];
    2586                        for (int i = 0; nextIndexes(it, indexes);) {
    2587                                Arrays.fill(fox, 0);
    2588                                for (int j = 0; j < m; j++) {
    2589                                        double c = coeff[j];
    2590                                        int l = indexes[j];
    2591                                        for (int k = 0; k < isize; k++) {
    2592                                                fox[k] += af32data[l++] * c;
    2593                                        }
    2594                                }
    2595                                for (int k = 0; k < isize; k++) {
    2596                                        oaf32data[i++] = fox[k];
    2597                                }
    2598                        }
    2599                        break;
    2600                case Dataset.ARRAYFLOAT64:
    2601                        final double[] af64data = ((CompoundDoubleDataset) a).data;
    2602                        final double[] oaf64data = ((CompoundDoubleDataset) out).getData();
    2603                        double[] dox = new double[isize];
    2604                        for (int i = 0; nextIndexes(it, indexes);) {
    2605                                Arrays.fill(dox, 0);
    2606                                for (int j = 0; j < m; j++) {
    2607                                        double c = coeff[j];
    2608                                        int l = indexes[j];
    2609                                        for (int k = 0; k < isize; k++) {
    2610                                                dox[k] += af64data[l++] * c;
    2611                                        }
    2612                                }
    2613                                for (int k = 0; k < isize; k++) {
    2614                                        oaf64data[i++] = dox[k];
    2615                                }
    2616                        }
    2617                        break;
    2618                default:
    2619                        throw new UnsupportedOperationException("difference does not support multiple-element dataset");
    2620                }
    2621        }
    2622
    2623        /**
    2624         * Discrete difference of dataset along axis using finite difference
    2625         *
    2626         * @param a
    2627         * @param n
    2628         *            order of difference
    2629         * @param axis
    2630         * @return difference
    2631         */
    2632        public static Dataset difference(Dataset a, final int n, int axis) {
    2633                Dataset ds;
    2634                final Class<? extends Dataset> clazz = a.getClass();
    2635                final int rank = a.getRank();
    2636                final int is = a.getElementsPerItem();
    2637
    2638                if (axis < 0) {
    2639                        axis += rank;
    2640                }
    2641                if (axis < 0 || axis >= rank) {
    2642                        throw new IllegalArgumentException("Axis is out of range");
    2643                }
    2644
    2645                int[] nshape = a.getShape();
    2646                if (nshape[axis] <= n) {
    2647                        nshape[axis] = 0;
    2648                        return DatasetFactory.zeros(is, clazz, nshape);
    2649                }
    2650
    2651                nshape[axis] -= n;
    2652                ds = DatasetFactory.zeros(is, clazz, nshape);
    2653                if (rank == 1) {
    2654                        difference(DatasetUtils.convertToDataset(a), ds, n);
    2655                } else {
    2656                        final Dataset src = DatasetFactory.zeros(is, clazz, a.getShapeRef()[axis]);
    2657                        final Dataset dest = DatasetFactory.zeros(is, clazz, nshape[axis]);
    2658                        final PositionIterator pi = a.getPositionIterator(axis);
    2659                        final int[] pos = pi.getPos();
    2660                        final boolean[] hit = pi.getOmit();
    2661                        while (pi.hasNext()) {
    2662                                a.copyItemsFromAxes(pos, hit, src);
    2663                                difference(src, dest, n);
    2664                                ds.setItemsOnAxes(pos, hit, dest.getBuffer());
    2665                        }
    2666                }
    2667
    2668                return ds;
    2669        }
    2670
    2671        private static double SelectedMean(Dataset data, int Min, int Max) {
    2672
    2673                double result = 0.0;
    2674                for (int i = Min, imax = data.getSize(); i <= Max; i++) {
    2675                        // clip i appropriately, imagine that effectively the two ends
    2676                        // continue
    2677                        // straight out.
    2678                        int pos = i;
    2679                        if (pos < 0) {
    2680                                pos = 0;
    2681                        } else if (pos >= imax) {
    2682                                pos = imax - 1;
    2683                        }
    2684                        result += data.getElementDoubleAbs(pos);
    2685                }
    2686
    2687                // now the sum is complete, average the values.
    2688                result /= (Max - Min) + 1;
    2689                return result;
    2690        }
    2691
    2692        private static void SelectedMeanArray(double[] out, Dataset data, int Min, int Max) {
    2693                final int isize = out.length;
    2694                for (int j = 0; j < isize; j++)
    2695                        out[j] = 0.;
    2696
    2697                for (int i = Min, imax = data.getSize(); i <= Max; i++) {
    2698                        // clip i appropriately, imagine that effectively the two ends
    2699                        // continue
    2700                        // straight out.
    2701                        int pos = i * isize;
    2702                        if (pos < 0) {
    2703                                pos = 0;
    2704                        } else if (pos >= imax) {
    2705                                pos = imax - isize;
    2706                        }
    2707                        for (int j = 0; j < isize; j++)
    2708                                out[j] += data.getElementDoubleAbs(pos + j);
    2709                }
    2710
    2711                // now the sum is complete, average the values.
    2712                double norm = 1. / (Max - Min + 1.);
    2713                for (int j = 0; j < isize; j++)
    2714                        out[j] *= norm;
    2715
    2716        }
    2717
    2718        /**
    2719         * Calculates the derivative of a line described by two datasets (x,y) given
    2720         * a spread of n either side of the point
    2721         *
    2722         * @param x
    2723         *            The x values of the function to take the derivative of.
    2724         * @param y
    2725         *            The y values of the function to take the derivative of.
    2726         * @param n
    2727         *            The spread the derivative is calculated from, i.e. the
    2728         *            smoothing, the higher the value, the more smoothing occurs.
    2729         * @return A dataset which contains all the derivative point for point.
    2730         */
    2731        public static Dataset derivative(Dataset x, Dataset y, int n) {
    2732                if (x.getRank() != 1 || y.getRank() != 1) {
    2733                        throw new IllegalArgumentException("Only one dimensional dataset supported");
    2734                }
    2735                if (y.getSize() > x.getSize()) {
    2736                        throw new IllegalArgumentException("Length of x dataset should be greater than or equal to y's");
    2737                }
    2738                int dtype = y.getDType();
    2739                Dataset result;
    2740                switch (dtype) {
    2741                case Dataset.BOOL:
    2742                case Dataset.INT8:
    2743                case Dataset.INT16:
    2744                case Dataset.ARRAYINT8:
    2745                case Dataset.ARRAYINT16:
    2746                        result = DatasetFactory.zeros(y, FloatDataset.class);
    2747                        break;
    2748                case Dataset.INT32:
    2749                case Dataset.INT64:
    2750                case Dataset.ARRAYINT32:
    2751                case Dataset.ARRAYINT64:
    2752                        result = DatasetFactory.zeros(y, DoubleDataset.class);
    2753                        break;
    2754                case Dataset.FLOAT32:
    2755                case Dataset.FLOAT64:
    2756                case Dataset.COMPLEX64:
    2757                case Dataset.COMPLEX128:
    2758                case Dataset.ARRAYFLOAT32:
    2759                case Dataset.ARRAYFLOAT64:
    2760                        result = DatasetFactory.zeros(y);
    2761                        break;
    2762                default:
    2763                        throw new UnsupportedOperationException("derivative does not support multiple-element dataset");
    2764                }
    2765
    2766                final int isize = y.getElementsPerItem();
    2767                if (isize == 1) {
    2768                        for (int i = 0, imax = x.getSize(); i < imax; i++) {
    2769                                double LeftValue = SelectedMean(y, i - n, i - 1);
    2770                                double RightValue = SelectedMean(y, i + 1, i + n);
    2771                                double LeftPosition = SelectedMean(x, i - n, i - 1);
    2772                                double RightPosition = SelectedMean(x, i + 1, i + n);
    2773
    2774                                // now the values and positions are calculated, the derivative
    2775                                // can be
    2776                                // calculated.
    2777                                result.set(((RightValue - LeftValue) / (RightPosition - LeftPosition)), i);
    2778                        }
    2779                } else {
    2780                        double[] leftValues = new double[isize];
    2781                        double[] rightValues = new double[isize];
    2782                        for (int i = 0, imax = x.getSize(); i < imax; i++) {
    2783                                SelectedMeanArray(leftValues, y, i - n, i - 1);
    2784                                SelectedMeanArray(rightValues, y, i + 1, i + n);
    2785                                double delta = SelectedMean(x, i - n, i - 1);
    2786                                delta = 1. / (SelectedMean(x, i + 1, i + n) - delta);
    2787                                for (int j = 0; j < isize; j++) {
    2788                                        rightValues[j] -= leftValues[j];
    2789                                        rightValues[j] *= delta;
    2790                                }
    2791                                result.set(rightValues, i);
    2792                        }
    2793                }
    2794
    2795                // set the name based on the changes made
    2796                result.setName(y.getName() + "'");
    2797
    2798                return result;
    2799        }
    2800
    2801        /**
    2802         * Discrete difference of dataset along axis using finite central difference
    2803         *
    2804         * @param a
    2805         * @param axis
    2806         * @return difference
    2807         */
    2808        public static Dataset centralDifference(Dataset a, int axis) {
    2809                Dataset ds;
    2810                final Class<? extends Dataset> clazz = a.getClass();
    2811                final int rank = a.getRank();
    2812                final int is = a.getElementsPerItem();
    2813
    2814                if (axis < 0) {
    2815                        axis += rank;
    2816                }
    2817                if (axis < 0 || axis >= rank) {
    2818                        throw new IllegalArgumentException("Axis is out of range");
    2819                }
    2820
    2821                final int len = a.getShapeRef()[axis];
    2822                if (len < 2) {
    2823                        throw new IllegalArgumentException("Dataset should have a size > 1 along given axis");
    2824                }
    2825                ds = DatasetFactory.zeros(is, clazz, a.getShapeRef());
    2826                if (rank == 1) {
    2827                        centralDifference(a, ds);
    2828                } else {
    2829                        final Dataset src = DatasetFactory.zeros(is, clazz, len);
    2830                        final Dataset dest = DatasetFactory.zeros(is, clazz, len);
    2831                        final PositionIterator pi = a.getPositionIterator(axis);
    2832                        final int[] pos = pi.getPos();
    2833                        final boolean[] hit = pi.getOmit();
    2834                        while (pi.hasNext()) {
    2835                                a.copyItemsFromAxes(pos, hit, src);
    2836                                centralDifference(src, dest);
    2837                                ds.setItemsOnAxes(pos, hit, dest.getBuffer());
    2838                        }
    2839                }
    2840
    2841                return ds;
    2842        }
    2843
    2844        /**
    2845         * 1st order discrete difference of dataset along flattened dataset using
    2846         * central difference
    2847         *
    2848         * @param a
    2849         *            is 1d dataset
    2850         * @param out
    2851         *            is 1d dataset
    2852         */
    2853        private static void centralDifference(final Dataset a, final Dataset out) {
    2854                final int isize = a.getElementsPerItem();
    2855                final int dt = a.getDType();
    2856
    2857                final int nlen = (out.getShapeRef()[0] - 1) * isize;
    2858                if (nlen < 1) {
    2859                        throw new IllegalArgumentException("Dataset should have a size > 1 along given axis");
    2860                }
    2861                final IndexIterator it = a.getIterator();
    2862                if (!it.hasNext())
    2863                        return;
    2864                int oi = it.index;
    2865                if (!it.hasNext())
    2866                        return;
    2867                int pi = it.index;
    2868
    2869                switch (dt) {
    2870                case Dataset.INT8:
    2871                        final byte[] i8data = ((ByteDataset) a).data;
    2872                        final byte[] oi8data = ((ByteDataset) out).getData();
    2873                        oi8data[0] = (byte) (i8data[pi] - i8data[oi]);
    2874                        for (int i = 1; it.hasNext(); i++) {
    2875                                oi8data[i] = (byte) ((i8data[it.index] - i8data[oi]) / 2);
    2876                                oi = pi;
    2877                                pi = it.index;
    2878                        }
    2879                        oi8data[nlen] = (byte) (i8data[pi] - i8data[oi]);
    2880                        break;
    2881                case Dataset.INT16:
    2882                        final short[] i16data = ((ShortDataset) a).data;
    2883                        final short[] oi16data = ((ShortDataset) out).getData();
    2884                        oi16data[0] = (short) (i16data[pi] - i16data[oi]);
    2885                        for (int i = 1; it.hasNext(); i++) {
    2886                                oi16data[i] = (short) ((i16data[it.index] - i16data[oi]) / 2);
    2887                                oi = pi;
    2888                                pi = it.index;
    2889                        }
    2890                        oi16data[nlen] = (short) (i16data[pi] - i16data[oi]);
    2891                        break;
    2892                case Dataset.INT32:
    2893                        final int[] i32data = ((IntegerDataset) a).data;
    2894                        final int[] oi32data = ((IntegerDataset) out).getData();
    2895                        oi32data[0] = i32data[pi] - i32data[oi];
    2896                        for (int i = 1; it.hasNext(); i++) {
    2897                                oi32data[i] = (i32data[it.index] - i32data[oi]) / 2;
    2898                                oi = pi;
    2899                                pi = it.index;
    2900                        }
    2901                        oi32data[nlen] = i32data[pi] - i32data[oi];
    2902                        break;
    2903                case Dataset.INT64:
    2904                        final long[] i64data = ((LongDataset) a).data;
    2905                        final long[] oi64data = ((LongDataset) out).getData();
    2906                        oi64data[0] = i64data[pi] - i64data[oi];
    2907                        for (int i = 1; it.hasNext(); i++) {
    2908                                oi64data[i] = (i64data[it.index] - i64data[oi]) / 2;
    2909                                oi = pi;
    2910                                pi = it.index;
    2911                        }
    2912                        oi64data[nlen] = i64data[pi] - i64data[oi];
    2913                        break;
    2914                case Dataset.ARRAYINT8:
    2915                        final byte[] ai8data = ((CompoundByteDataset) a).data;
    2916                        final byte[] oai8data = ((CompoundByteDataset) out).getData();
    2917                        for (int k = 0; k < isize; k++) {
    2918                                oai8data[k] = (byte) (ai8data[pi + k] - ai8data[oi + k]);
    2919                        }
    2920                        for (int i = isize; it.hasNext();) {
    2921                                int l = it.index;
    2922                                for (int k = 0; k < isize; k++) {
    2923                                        oai8data[i++] = (byte) ((ai8data[l++] - ai8data[oi++]) / 2);
    2924                                }
    2925                                oi = pi;
    2926                                pi = it.index;
    2927                        }
    2928                        for (int k = 0; k < isize; k++) {
    2929                                oai8data[nlen + k] = (byte) (ai8data[pi + k] - ai8data[oi + k]);
    2930                        }
    2931                        break;
    2932                case Dataset.ARRAYINT16:
    2933                        final short[] ai16data = ((CompoundShortDataset) a).data;
    2934                        final short[] oai16data = ((CompoundShortDataset) out).getData();
    2935                        for (int k = 0; k < isize; k++) {
    2936                                oai16data[k] = (short) (ai16data[pi + k] - ai16data[oi + k]);
    2937                        }
    2938                        for (int i = isize; it.hasNext();) {
    2939                                int l = it.index;
    2940                                for (int k = 0; k < isize; k++) {
    2941                                        oai16data[i++] = (short) ((ai16data[l++] - ai16data[oi++]) / 2);
    2942                                }
    2943                                oi = pi;
    2944                                pi = it.index;
    2945                        }
    2946                        for (int k = 0; k < isize; k++) {
    2947                                oai16data[nlen + k] = (short) (ai16data[pi + k] - ai16data[oi + k]);
    2948                        }
    2949                        break;
    2950                case Dataset.ARRAYINT32:
    2951                        final int[] ai32data = ((CompoundIntegerDataset) a).data;
    2952                        final int[] oai32data = ((CompoundIntegerDataset) out).getData();
    2953                        for (int k = 0; k < isize; k++) {
    2954                                oai32data[k] = ai32data[pi + k] - ai32data[oi + k];
    2955                        }
    2956                        for (int i = isize; it.hasNext();) {
    2957                                int l = it.index;
    2958                                for (int k = 0; k < isize; k++) {
    2959                                        oai32data[i++] = (ai32data[l++] - ai32data[oi++]) / 2;
    2960                                }
    2961                                oi = pi;
    2962                                pi = it.index;
    2963                        }
    2964                        for (int k = 0; k < isize; k++) {
    2965                                oai32data[nlen + k] = ai32data[pi + k] - ai32data[oi + k];
    2966                        }
    2967                        break;
    2968                case Dataset.ARRAYINT64:
    2969                        final long[] ai64data = ((CompoundLongDataset) a).data;
    2970                        final long[] oai64data = ((CompoundLongDataset) out).getData();
    2971                        for (int k = 0; k < isize; k++) {
    2972                                oai64data[k] = ai64data[pi + k] - ai64data[oi + k];
    2973                        }
    2974                        for (int i = isize; it.hasNext();) {
    2975                                int l = it.index;
    2976                                for (int k = 0; k < isize; k++) {
    2977                                        oai64data[i++] = (ai64data[l++] - ai64data[oi++]) / 2;
    2978                                }
    2979                                oi = pi;
    2980                                pi = it.index;
    2981                        }
    2982                        for (int k = 0; k < isize; k++) {
    2983                                oai64data[nlen + k] = ai64data[pi + k] - ai64data[oi + k];
    2984                        }
    2985                        break;
    2986                case Dataset.FLOAT32:
    2987                        final float[] f32data = ((FloatDataset) a).data;
    2988                        final float[] of32data = ((FloatDataset) out).getData();
    2989                        of32data[0] = f32data[pi] - f32data[oi];
    2990                        for (int i = 1; it.hasNext(); i++) {
    2991                                of32data[i] = (f32data[it.index] - f32data[oi]) * 0.5f;
    2992                                oi = pi;
    2993                                pi = it.index;
    2994                        }
    2995                        of32data[nlen] = f32data[pi] - f32data[oi];
    2996                        break;
    2997                case Dataset.FLOAT64:
    2998                        final double[] f64data = ((DoubleDataset) a).data;
    2999                        final double[] of64data = ((DoubleDataset) out).getData();
    3000                        of64data[0] = f64data[pi] - f64data[oi];
    3001                        for (int i = 1; it.hasNext(); i++) {
    3002                                of64data[i] = (f64data[it.index] - f64data[oi]) * 0.5f;
    3003                                oi = pi;
    3004                                pi = it.index;
    3005                        }
    3006                        of64data[nlen] = f64data[pi] - f64data[oi];
    3007                        break;
    3008                case Dataset.COMPLEX64:
    3009                        final float[] c64data = ((ComplexFloatDataset) a).data;
    3010                        final float[] oc64data = ((ComplexFloatDataset) out).getData();
    3011                        oc64data[0] = c64data[pi] - c64data[oi];
    3012                        oc64data[1] = c64data[pi + 1] - c64data[oi + 1];
    3013                        for (int i = 2; it.hasNext();) {
    3014                                oc64data[i++] = (c64data[it.index] - c64data[oi++]) * 0.5f;
    3015                                oc64data[i++] = (c64data[it.index + 1] - c64data[oi]) * 0.5f;
    3016                                oi = pi;
    3017                                pi = it.index;
    3018                        }
    3019                        oc64data[nlen] = c64data[pi] - c64data[oi];
    3020                        oc64data[nlen + 1] = c64data[pi + 1] - c64data[oi + 1];
    3021                        break;
    3022                case Dataset.COMPLEX128:
    3023                        final double[] c128data = ((ComplexDoubleDataset) a).data;
    3024                        final double[] oc128data = ((ComplexDoubleDataset) out).getData();
    3025                        oc128data[0] = c128data[pi] - c128data[oi];
    3026                        oc128data[1] = c128data[pi + 1] - c128data[oi + 1];
    3027                        for (int i = 2; it.hasNext();) {
    3028                                oc128data[i++] = (c128data[it.index] - c128data[oi++]) * 0.5f;
    3029                                oc128data[i++] = (c128data[it.index + 1] - c128data[oi]) * 0.5f;
    3030                                oi = pi;
    3031                                pi = it.index;
    3032                        }
    3033                        oc128data[nlen] = c128data[pi] - c128data[oi];
    3034                        oc128data[nlen + 1] = c128data[pi + 1] - c128data[oi + 1];
    3035                        break;
    3036                case Dataset.ARRAYFLOAT32:
    3037                        final float[] af32data = ((CompoundFloatDataset) a).data;
    3038                        final float[] oaf32data = ((CompoundFloatDataset) out).getData();
    3039                        for (int k = 0; k < isize; k++) {
    3040                                oaf32data[k] = af32data[pi + k] - af32data[oi + k];
    3041                        }
    3042                        for (int i = isize; it.hasNext();) {
    3043                                int l = it.index;
    3044                                for (int k = 0; k < isize; k++) {
    3045                                        oaf32data[i++] = (af32data[l++] - af32data[oi++]) * 0.5f;
    3046                                }
    3047                                oi = pi;
    3048                                pi = it.index;
    3049                        }
    3050                        for (int k = 0; k < isize; k++) {
    3051                                oaf32data[nlen + k] = af32data[pi + k] - af32data[oi + k];
    3052                        }
    3053                        break;
    3054                case Dataset.ARRAYFLOAT64:
    3055                        final double[] af64data = ((CompoundDoubleDataset) a).data;
    3056                        final double[] oaf64data = ((CompoundDoubleDataset) out).getData();
    3057                        for (int k = 0; k < isize; k++) {
    3058                                oaf64data[k] = af64data[pi + k] - af64data[oi + k];
    3059                        }
    3060                        for (int i = isize; it.hasNext();) {
    3061                                int l = it.index;
    3062                                for (int k = 0; k < isize; k++) {
    3063                                        oaf64data[i++] = (af64data[l++] - af64data[oi++]) * 0.5;
    3064                                }
    3065                                oi = pi;
    3066                                pi = it.index;
    3067                        }
    3068                        for (int k = 0; k < isize; k++) {
    3069                                oaf64data[nlen + k] = af64data[pi + k] - af64data[oi + k];
    3070                        }
    3071                        break;
    3072                default:
    3073                        throw new UnsupportedOperationException("difference does not support this dataset type");
    3074                }
    3075        }
    3076
    3077        /**
    3078         * Calculate gradient (or partial derivatives) by central difference
    3079         *
    3080         * @param y
    3081         * @param x
    3082         *            one or more datasets for dependent variables
    3083         * @return a list of datasets (one for each dimension in y)
    3084         */
    3085        public static List<Dataset> gradient(Dataset y, Dataset... x) {
    3086                final int rank = y.getRank();
    3087
    3088                if (x.length > 0) {
    3089                        if (x.length != rank) {
    3090                                throw new IllegalArgumentException(
    3091                                                "Number of dependent datasets must be equal to rank of first argument");
    3092                        }
    3093                        for (int a = 0; a < rank; a++) {
    3094                                int rx = x[a].getRank();
    3095                                if (rx != rank && rx != 1) {
    3096                                        throw new IllegalArgumentException(
    3097                                                        "Dependent datasets must be 1-D or match rank of first argument");
    3098                                }
    3099                                if (rx == 1) {
    3100                                        if (y.getShapeRef()[a] != x[a].getShapeRef()[0]) {
    3101                                                throw new IllegalArgumentException("Length of dependent dataset must match axis length");
    3102                                        }
    3103                                } else {
    3104                                        y.checkCompatibility(x[a]);
    3105                                }
    3106                        }
    3107                }
    3108
    3109                List<Dataset> grad = new ArrayList<Dataset>(rank);
    3110
    3111                for (int a = 0; a < rank; a++) {
    3112                        Dataset g = centralDifference(y, a);
    3113                        grad.add(g);
    3114                }
    3115
    3116                if (x.length > 0) {
    3117                        for (int a = 0; a < rank; a++) {
    3118                                Dataset g = grad.get(a);
    3119                                Dataset dx = x[a];
    3120                                int r = dx.getRank();
    3121                                if (r == rank) {
    3122                                        g.idivide(centralDifference(dx, a));
    3123                                } else {
    3124                                        final int is = dx.getElementsPerItem();
    3125                                        final Dataset bdx = DatasetFactory.zeros(is, dx.getClass(), y.getShapeRef());
    3126                                        final PositionIterator pi = y.getPositionIterator(a);
    3127                                        final int[] pos = pi.getPos();
    3128                                        final boolean[] hit = pi.getOmit();
    3129                                        dx = centralDifference(dx, 0);
    3130
    3131                                        while (pi.hasNext()) {
    3132                                                bdx.setItemsOnAxes(pos, hit, dx.getBuffer());
    3133                                        }
    3134                                        g.idivide(bdx);
    3135                                }
    3136                        }
    3137                }
    3138                return grad;
    3139        }
    3140}