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