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}