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