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