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