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