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