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