001// --- BEGIN LICENSE BLOCK --- 002/* 003 * Copyright (c) 2009, Mikio L. Braun 004 * All rights reserved. 005 * 006 * Redistribution and use in source and binary forms, with or without 007 * modification, are permitted provided that the following conditions are 008 * met: 009 * 010 * * Redistributions of source code must retain the above copyright 011 * notice, this list of conditions and the following disclaimer. 012 * 013 * * Redistributions in binary form must reproduce the above 014 * copyright notice, this list of conditions and the following 015 * disclaimer in the documentation and/or other materials provided 016 * with the distribution. 017 * 018 * * Neither the name of the Technische Universität Berlin nor the 019 * names of its contributors may be used to endorse or promote 020 * products derived from this software without specific prior 021 * written permission. 022 * 023 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 024 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 025 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR 026 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT 027 * HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, 028 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT 029 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 030 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 031 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 032 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 033 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 034 */ 035// --- END LICENSE BLOCK --- 036 037package org.jblas; 038 039/** 040 * This class provides the functions from java.lang.Math for matrices. The 041 * functions are applied to each element of the matrix. 042 * 043 * @author Mikio Braun 044 */ 045public class MatrixFunctions { 046 047 /*# 048 def mapfct(f); <<-EOS 049 for (int i = 0; i < x.length; i++) 050 x.put(i, (double) #{f}(x.get(i))); 051 return x; 052 EOS 053 end 054 055 def cmapfct(f); <<-EOS 056 for (int i = 0; i < x.length; i++) 057 x.put(i, x.get(i).#{f}()); 058 return x; 059 EOS 060 end 061 #*/ 062 063 /** 064 * Sets all elements in this matrix to their absolute values. Note 065 * that this operation is in-place. 066 * @see MatrixFunctions#abs(DoubleMatrix) 067 * @return this matrix 068 */ 069 public static DoubleMatrix absi(DoubleMatrix x) { 070 /*# mapfct('Math.abs') #*/ 071//RJPP-BEGIN------------------------------------------------------------ 072 for (int i = 0; i < x.length; i++) 073 x.put(i, (double) Math.abs(x.get(i))); 074 return x; 075//RJPP-END-------------------------------------------------------------- 076 } 077 078 public static ComplexDoubleMatrix absi(ComplexDoubleMatrix x) { 079 /*# cmapfct('abs') #*/ 080//RJPP-BEGIN------------------------------------------------------------ 081 for (int i = 0; i < x.length; i++) 082 x.put(i, x.get(i).abs()); 083 return x; 084//RJPP-END-------------------------------------------------------------- 085 } 086 087 /** 088 * Applies the trigonometric <i>arccosine</i> function element wise on this 089 * matrix. Note that this is an in-place operation. 090 * @see MatrixFunctions#acos(DoubleMatrix) 091 * @return this matrix 092 */ 093 public static DoubleMatrix acosi(DoubleMatrix x) { 094 /*# mapfct('Math.acos') #*/ 095//RJPP-BEGIN------------------------------------------------------------ 096 for (int i = 0; i < x.length; i++) 097 x.put(i, (double) Math.acos(x.get(i))); 098 return x; 099//RJPP-END-------------------------------------------------------------- 100 } 101 102 /** 103 * Applies the trigonometric <i>arcsine</i> function element wise on this 104 * matrix. Note that this is an in-place operation. 105 * @see MatrixFunctions#asin(DoubleMatrix) 106 * @return this matrix 107 */ 108 public static DoubleMatrix asini(DoubleMatrix x) { 109 /*# mapfct('Math.asin') #*/ 110//RJPP-BEGIN------------------------------------------------------------ 111 for (int i = 0; i < x.length; i++) 112 x.put(i, (double) Math.asin(x.get(i))); 113 return x; 114//RJPP-END-------------------------------------------------------------- 115 } 116 117 /** 118 * Applies the trigonometric <i>arctangend</i> function element wise on this 119 * matrix. Note that this is an in-place operation. 120 * @see MatrixFunctions#atan(DoubleMatrix) 121 * @return this matrix 122 */ 123 public static DoubleMatrix atani(DoubleMatrix x) { 124 /*# mapfct('Math.atan') #*/ 125//RJPP-BEGIN------------------------------------------------------------ 126 for (int i = 0; i < x.length; i++) 127 x.put(i, (double) Math.atan(x.get(i))); 128 return x; 129//RJPP-END-------------------------------------------------------------- 130 } 131 132 /** 133 * Applies the <i>cube root</i> function element wise on this 134 * matrix. Note that this is an in-place operation. 135 * @see MatrixFunctions#cbrt(DoubleMatrix) 136 * @return this matrix 137 */ 138 public static DoubleMatrix cbrti(DoubleMatrix x) { 139 /*# mapfct('Math.cbrt') #*/ 140//RJPP-BEGIN------------------------------------------------------------ 141 for (int i = 0; i < x.length; i++) 142 x.put(i, (double) Math.cbrt(x.get(i))); 143 return x; 144//RJPP-END-------------------------------------------------------------- 145 } 146 147 /** 148 * Element-wise round up by applying the <i>ceil</i> function on each 149 * element. Note that this is an in-place operation. 150 * @see MatrixFunctions#ceil(DoubleMatrix) 151 * @return this matrix 152 */ 153 public static DoubleMatrix ceili(DoubleMatrix x) { 154 /*# mapfct('Math.ceil') #*/ 155//RJPP-BEGIN------------------------------------------------------------ 156 for (int i = 0; i < x.length; i++) 157 x.put(i, (double) Math.ceil(x.get(i))); 158 return x; 159//RJPP-END-------------------------------------------------------------- 160 } 161 162 /** 163 * Applies the <i>cosine</i> function element-wise on this 164 * matrix. Note that this is an in-place operation. 165 * @see MatrixFunctions#cos(DoubleMatrix) 166 * @return this matrix 167 */ 168 public static DoubleMatrix cosi(DoubleMatrix x) { 169 /*# mapfct('Math.cos') #*/ 170//RJPP-BEGIN------------------------------------------------------------ 171 for (int i = 0; i < x.length; i++) 172 x.put(i, (double) Math.cos(x.get(i))); 173 return x; 174//RJPP-END-------------------------------------------------------------- 175 } 176 177 /** 178 * Applies the <i>hyperbolic cosine</i> function element-wise on this 179 * matrix. Note that this is an in-place operation. 180 * @see MatrixFunctions#cosh(DoubleMatrix) 181 * @return this matrix 182 */ 183 public static DoubleMatrix coshi(DoubleMatrix x) { 184 /*# mapfct('Math.cosh') #*/ 185//RJPP-BEGIN------------------------------------------------------------ 186 for (int i = 0; i < x.length; i++) 187 x.put(i, (double) Math.cosh(x.get(i))); 188 return x; 189//RJPP-END-------------------------------------------------------------- 190 } 191 192 /** 193 * Applies the <i>exponential</i> function element-wise on this 194 * matrix. Note that this is an in-place operation. 195 * @see MatrixFunctions#exp(DoubleMatrix) 196 * @return this matrix 197 */ 198 public static DoubleMatrix expi(DoubleMatrix x) { 199 /*# mapfct('Math.exp') #*/ 200//RJPP-BEGIN------------------------------------------------------------ 201 for (int i = 0; i < x.length; i++) 202 x.put(i, (double) Math.exp(x.get(i))); 203 return x; 204//RJPP-END-------------------------------------------------------------- 205 } 206 207 /** 208 * Element-wise round down by applying the <i>floor</i> function on each 209 * element. Note that this is an in-place operation. 210 * @see MatrixFunctions#floor(DoubleMatrix) 211 * @return this matrix 212 */ 213 public static DoubleMatrix floori(DoubleMatrix x) { 214 /*# mapfct('Math.floor') #*/ 215//RJPP-BEGIN------------------------------------------------------------ 216 for (int i = 0; i < x.length; i++) 217 x.put(i, (double) Math.floor(x.get(i))); 218 return x; 219//RJPP-END-------------------------------------------------------------- 220 } 221 222 /** 223 * Applies the <i>natural logarithm</i> function element-wise on this 224 * matrix. Note that this is an in-place operation. 225 * @see MatrixFunctions#log(DoubleMatrix) 226 * @return this matrix 227 */ 228 public static DoubleMatrix logi(DoubleMatrix x) { 229 /*# mapfct('Math.log') #*/ 230//RJPP-BEGIN------------------------------------------------------------ 231 for (int i = 0; i < x.length; i++) 232 x.put(i, (double) Math.log(x.get(i))); 233 return x; 234//RJPP-END-------------------------------------------------------------- 235 } 236 237 /** 238 * Applies the <i>logarithm with basis to 10</i> element-wise on this 239 * matrix. Note that this is an in-place operation. 240 * @see MatrixFunctions#log10(DoubleMatrix) 241 * @return this matrix 242 */ 243 public static DoubleMatrix log10i(DoubleMatrix x) { 244 /*# mapfct('Math.log10') #*/ 245//RJPP-BEGIN------------------------------------------------------------ 246 for (int i = 0; i < x.length; i++) 247 x.put(i, (double) Math.log10(x.get(i))); 248 return x; 249//RJPP-END-------------------------------------------------------------- 250 } 251 252 /** 253 * Element-wise power function. Replaces each element with its 254 * power of <tt>d</tt>.Note that this is an in-place operation. 255 * @param d the exponent 256 * @see MatrixFunctions#pow(DoubleMatrix,double) 257 * @return this matrix 258 */ 259 public static DoubleMatrix powi(DoubleMatrix x, double d) { 260 if (d == 2.0) 261 return x.muli(x); 262 else { 263 for (int i = 0; i < x.length; i++) 264 x.put(i, (double) Math.pow(x.get(i), d)); 265 return x; 266 } 267 } 268 269 public static DoubleMatrix powi(double base, DoubleMatrix x) { 270 for (int i = 0; i < x.length; i++) 271 x.put(i, (double) Math.pow(base, x.get(i))); 272 return x; 273 } 274 275 public static DoubleMatrix powi(DoubleMatrix x, DoubleMatrix e) { 276 x.checkLength(e.length); 277 for (int i = 0; i < x.length; i++) 278 x.put(i, (double) Math.pow(x.get(i), e.get(i))); 279 return x; 280 } 281 282 public static DoubleMatrix signumi(DoubleMatrix x) { 283 /*# mapfct('Math.signum') #*/ 284//RJPP-BEGIN------------------------------------------------------------ 285 for (int i = 0; i < x.length; i++) 286 x.put(i, (double) Math.signum(x.get(i))); 287 return x; 288//RJPP-END-------------------------------------------------------------- 289 } 290 291 public static DoubleMatrix sini(DoubleMatrix x) { 292 /*# mapfct('Math.sin') #*/ 293//RJPP-BEGIN------------------------------------------------------------ 294 for (int i = 0; i < x.length; i++) 295 x.put(i, (double) Math.sin(x.get(i))); 296 return x; 297//RJPP-END-------------------------------------------------------------- 298 } 299 300 public static DoubleMatrix sinhi(DoubleMatrix x) { 301 /*# mapfct('Math.sinh') #*/ 302//RJPP-BEGIN------------------------------------------------------------ 303 for (int i = 0; i < x.length; i++) 304 x.put(i, (double) Math.sinh(x.get(i))); 305 return x; 306//RJPP-END-------------------------------------------------------------- 307 } 308 public static DoubleMatrix sqrti(DoubleMatrix x) { 309 /*# mapfct('Math.sqrt') #*/ 310//RJPP-BEGIN------------------------------------------------------------ 311 for (int i = 0; i < x.length; i++) 312 x.put(i, (double) Math.sqrt(x.get(i))); 313 return x; 314//RJPP-END-------------------------------------------------------------- 315 } 316 public static DoubleMatrix tani(DoubleMatrix x) { 317 /*# mapfct('Math.tan') #*/ 318//RJPP-BEGIN------------------------------------------------------------ 319 for (int i = 0; i < x.length; i++) 320 x.put(i, (double) Math.tan(x.get(i))); 321 return x; 322//RJPP-END-------------------------------------------------------------- 323 } 324 public static DoubleMatrix tanhi(DoubleMatrix x) { 325 /*# mapfct('Math.tanh') #*/ 326//RJPP-BEGIN------------------------------------------------------------ 327 for (int i = 0; i < x.length; i++) 328 x.put(i, (double) Math.tanh(x.get(i))); 329 return x; 330//RJPP-END-------------------------------------------------------------- 331 } 332 333 /** 334 * Returns a copy of this matrix where all elements are set to their 335 * absolute values. 336 * @see MatrixFunctions#absi(DoubleMatrix) 337 * @return copy of this matrix 338 */ 339 public static DoubleMatrix abs(DoubleMatrix x) { return absi(x.dup()); } 340 341 /** 342 * Returns a copy of this matrix where the trigonometric <i>acos</i> function is applied 343 * element wise. 344 * @see MatrixFunctions#acosi(DoubleMatrix) 345 * @return copy of this matrix 346 */ 347 public static DoubleMatrix acos(DoubleMatrix x) { return acosi(x.dup()); } 348 public static DoubleMatrix asin(DoubleMatrix x) { return asini(x.dup()); } 349 public static DoubleMatrix atan(DoubleMatrix x) { return atani(x.dup()); } 350 public static DoubleMatrix cbrt(DoubleMatrix x) { return cbrti(x.dup()); } 351 public static DoubleMatrix ceil(DoubleMatrix x) { return ceili(x.dup()); } 352 public static DoubleMatrix cos(DoubleMatrix x) { return cosi(x.dup()); } 353 public static DoubleMatrix cosh(DoubleMatrix x) { return coshi(x.dup()); } 354 public static DoubleMatrix exp(DoubleMatrix x) { return expi(x.dup()); } 355 public static DoubleMatrix floor(DoubleMatrix x) { return floori(x.dup()); } 356 public static DoubleMatrix log(DoubleMatrix x) { return logi(x.dup()); } 357 public static DoubleMatrix log10(DoubleMatrix x) { return log10i(x.dup()); } 358 public static double pow(double x, double y) { return (double)Math.pow(x, y); } 359 public static DoubleMatrix pow(DoubleMatrix x, double e) { return powi(x.dup(), e); } 360 public static DoubleMatrix pow(double b, DoubleMatrix x) { return powi(b, x.dup()); } 361 public static DoubleMatrix pow(DoubleMatrix x, DoubleMatrix e) { return powi(x.dup(), e); } 362 public static DoubleMatrix signum(DoubleMatrix x) { return signumi(x.dup()); } 363 public static DoubleMatrix sin(DoubleMatrix x) { return sini(x.dup()); } 364 public static DoubleMatrix sinh(DoubleMatrix x) { return sinhi(x.dup()); } 365 public static DoubleMatrix sqrt(DoubleMatrix x) { return sqrti(x.dup()); } 366 public static DoubleMatrix tan(DoubleMatrix x) { return tani(x.dup()); } 367 public static DoubleMatrix tanh(DoubleMatrix x) { return tanhi(x.dup()); } 368 369 /*# %w{abs acos asin atan cbrt ceil cos cosh exp floor log log10 signum sin sinh sqrt tan tanh}.map do |fct| <<-EOS 370 public static double #{fct}(double x) { return (double)Math.#{fct}(x); } 371 EOS 372 end 373 #*/ 374//RJPP-BEGIN------------------------------------------------------------ 375 public static double abs(double x) { return (double)Math.abs(x); } 376 public static double acos(double x) { return (double)Math.acos(x); } 377 public static double asin(double x) { return (double)Math.asin(x); } 378 public static double atan(double x) { return (double)Math.atan(x); } 379 public static double cbrt(double x) { return (double)Math.cbrt(x); } 380 public static double ceil(double x) { return (double)Math.ceil(x); } 381 public static double cos(double x) { return (double)Math.cos(x); } 382 public static double cosh(double x) { return (double)Math.cosh(x); } 383 public static double exp(double x) { return (double)Math.exp(x); } 384 public static double floor(double x) { return (double)Math.floor(x); } 385 public static double log(double x) { return (double)Math.log(x); } 386 public static double log10(double x) { return (double)Math.log10(x); } 387 public static double signum(double x) { return (double)Math.signum(x); } 388 public static double sin(double x) { return (double)Math.sin(x); } 389 public static double sinh(double x) { return (double)Math.sinh(x); } 390 public static double sqrt(double x) { return (double)Math.sqrt(x); } 391 public static double tan(double x) { return (double)Math.tan(x); } 392 public static double tanh(double x) { return (double)Math.tanh(x); } 393//RJPP-END-------------------------------------------------------------- 394 395 /** 396 * Calculate matrix exponential of a square matrix. 397 * 398 * A scaled Pade approximation algorithm is used. 399 * The algorithm has been directly translated from Golub & Van Loan "Matrix Computations", 400 * algorithm 11.3.1. Special Horner techniques from 11.2 are also used to minimize the number 401 * of matrix multiplications. 402 * 403 * @param A square matrix 404 * @return matrix exponential of A 405 */ 406 public static DoubleMatrix expm(DoubleMatrix A) { 407 // constants for pade approximation 408 final double c0 = 1.0; 409 final double c1 = 0.5; 410 final double c2 = 0.12; 411 final double c3 = 0.01833333333333333; 412 final double c4 = 0.0019927536231884053; 413 final double c5 = 1.630434782608695E-4; 414 final double c6 = 1.0351966873706E-5; 415 final double c7 = 5.175983436853E-7; 416 final double c8 = 2.0431513566525E-8; 417 final double c9 = 6.306022705717593E-10; 418 final double c10 = 1.4837700484041396E-11; 419 final double c11 = 2.5291534915979653E-13; 420 final double c12 = 2.8101705462199615E-15; 421 final double c13 = 1.5440497506703084E-17; 422 423 int j = Math.max(0, 1 + (int) Math.floor(Math.log(A.normmax()) / Math.log(2))); 424 DoubleMatrix As = A.div((double) Math.pow(2, j)); // scaled version of A 425 int n = A.getRows(); 426 427 // calculate D and N using special Horner techniques 428 DoubleMatrix As_2 = As.mmul(As); 429 DoubleMatrix As_4 = As_2.mmul(As_2); 430 DoubleMatrix As_6 = As_4.mmul(As_2); 431 // U = c0*I + c2*A^2 + c4*A^4 + (c6*I + c8*A^2 + c10*A^4 + c12*A^6)*A^6 432 DoubleMatrix U = DoubleMatrix.eye(n).muli(c0).addi(As_2.mul(c2)).addi(As_4.mul(c4)).addi( 433 DoubleMatrix.eye(n).muli(c6).addi(As_2.mul(c8)).addi(As_4.mul(c10)).addi(As_6.mul(c12)).mmuli(As_6)); 434 // V = c1*I + c3*A^2 + c5*A^4 + (c7*I + c9*A^2 + c11*A^4 + c13*A^6)*A^6 435 DoubleMatrix V = DoubleMatrix.eye(n).muli(c1).addi(As_2.mul(c3)).addi(As_4.mul(c5)).addi( 436 DoubleMatrix.eye(n).muli(c7).addi(As_2.mul(c9)).addi(As_4.mul(c11)).addi(As_6.mul(c13)).mmuli(As_6)); 437 438 DoubleMatrix AV = As.mmuli(V); 439 DoubleMatrix N = U.add(AV); 440 DoubleMatrix D = U.subi(AV); 441 442 // solve DF = N for F 443 DoubleMatrix F = Solve.solve(D, N); 444 445 // now square j times 446 for (int k = 0; k < j; k++) { 447 F.mmuli(F); 448 } 449 450 return F; 451 } 452 453 454//STOP 455 public static DoubleMatrix floatToDouble(FloatMatrix fm) { 456 DoubleMatrix dm = new DoubleMatrix(fm.rows, fm.columns); 457 458 for (int i = 0; i < fm.length; i++) 459 dm.put(i, (double) fm.get(i)); 460 461 return dm; 462 } 463 464 public static FloatMatrix doubleToFloat(DoubleMatrix dm) { 465 FloatMatrix fm = new FloatMatrix(dm.rows, dm.columns); 466 467 for (int i = 0; i < dm.length; i++) 468 fm.put(i, (float) dm.get(i)); 469 470 return fm; 471 } 472//START 473 474//BEGIN 475 // The code below has been automatically generated. 476 // DO NOT EDIT! 477 478 /*# 479 def mapfct(f); <<-EOS 480 for (int i = 0; i < x.length; i++) 481 x.put(i, (float) #{f}(x.get(i))); 482 return x; 483 EOS 484 end 485 486 def cmapfct(f); <<-EOS 487 for (int i = 0; i < x.length; i++) 488 x.put(i, x.get(i).#{f}()); 489 return x; 490 EOS 491 end 492 #*/ 493 494 /** 495 * Sets all elements in this matrix to their absolute values. Note 496 * that this operation is in-place. 497 * @see MatrixFunctions#abs(FloatMatrix) 498 * @return this matrix 499 */ 500 public static FloatMatrix absi(FloatMatrix x) { 501 /*# mapfct('Math.abs') #*/ 502//RJPP-BEGIN------------------------------------------------------------ 503 for (int i = 0; i < x.length; i++) 504 x.put(i, (float) Math.abs(x.get(i))); 505 return x; 506//RJPP-END-------------------------------------------------------------- 507 } 508 509 public static ComplexFloatMatrix absi(ComplexFloatMatrix x) { 510 /*# cmapfct('abs') #*/ 511//RJPP-BEGIN------------------------------------------------------------ 512 for (int i = 0; i < x.length; i++) 513 x.put(i, x.get(i).abs()); 514 return x; 515//RJPP-END-------------------------------------------------------------- 516 } 517 518 /** 519 * Applies the trigonometric <i>arccosine</i> function element wise on this 520 * matrix. Note that this is an in-place operation. 521 * @see MatrixFunctions#acos(FloatMatrix) 522 * @return this matrix 523 */ 524 public static FloatMatrix acosi(FloatMatrix x) { 525 /*# mapfct('Math.acos') #*/ 526//RJPP-BEGIN------------------------------------------------------------ 527 for (int i = 0; i < x.length; i++) 528 x.put(i, (float) Math.acos(x.get(i))); 529 return x; 530//RJPP-END-------------------------------------------------------------- 531 } 532 533 /** 534 * Applies the trigonometric <i>arcsine</i> function element wise on this 535 * matrix. Note that this is an in-place operation. 536 * @see MatrixFunctions#asin(FloatMatrix) 537 * @return this matrix 538 */ 539 public static FloatMatrix asini(FloatMatrix x) { 540 /*# mapfct('Math.asin') #*/ 541//RJPP-BEGIN------------------------------------------------------------ 542 for (int i = 0; i < x.length; i++) 543 x.put(i, (float) Math.asin(x.get(i))); 544 return x; 545//RJPP-END-------------------------------------------------------------- 546 } 547 548 /** 549 * Applies the trigonometric <i>arctangend</i> function element wise on this 550 * matrix. Note that this is an in-place operation. 551 * @see MatrixFunctions#atan(FloatMatrix) 552 * @return this matrix 553 */ 554 public static FloatMatrix atani(FloatMatrix x) { 555 /*# mapfct('Math.atan') #*/ 556//RJPP-BEGIN------------------------------------------------------------ 557 for (int i = 0; i < x.length; i++) 558 x.put(i, (float) Math.atan(x.get(i))); 559 return x; 560//RJPP-END-------------------------------------------------------------- 561 } 562 563 /** 564 * Applies the <i>cube root</i> function element wise on this 565 * matrix. Note that this is an in-place operation. 566 * @see MatrixFunctions#cbrt(FloatMatrix) 567 * @return this matrix 568 */ 569 public static FloatMatrix cbrti(FloatMatrix x) { 570 /*# mapfct('Math.cbrt') #*/ 571//RJPP-BEGIN------------------------------------------------------------ 572 for (int i = 0; i < x.length; i++) 573 x.put(i, (float) Math.cbrt(x.get(i))); 574 return x; 575//RJPP-END-------------------------------------------------------------- 576 } 577 578 /** 579 * Element-wise round up by applying the <i>ceil</i> function on each 580 * element. Note that this is an in-place operation. 581 * @see MatrixFunctions#ceil(FloatMatrix) 582 * @return this matrix 583 */ 584 public static FloatMatrix ceili(FloatMatrix x) { 585 /*# mapfct('Math.ceil') #*/ 586//RJPP-BEGIN------------------------------------------------------------ 587 for (int i = 0; i < x.length; i++) 588 x.put(i, (float) Math.ceil(x.get(i))); 589 return x; 590//RJPP-END-------------------------------------------------------------- 591 } 592 593 /** 594 * Applies the <i>cosine</i> function element-wise on this 595 * matrix. Note that this is an in-place operation. 596 * @see MatrixFunctions#cos(FloatMatrix) 597 * @return this matrix 598 */ 599 public static FloatMatrix cosi(FloatMatrix x) { 600 /*# mapfct('Math.cos') #*/ 601//RJPP-BEGIN------------------------------------------------------------ 602 for (int i = 0; i < x.length; i++) 603 x.put(i, (float) Math.cos(x.get(i))); 604 return x; 605//RJPP-END-------------------------------------------------------------- 606 } 607 608 /** 609 * Applies the <i>hyperbolic cosine</i> function element-wise on this 610 * matrix. Note that this is an in-place operation. 611 * @see MatrixFunctions#cosh(FloatMatrix) 612 * @return this matrix 613 */ 614 public static FloatMatrix coshi(FloatMatrix x) { 615 /*# mapfct('Math.cosh') #*/ 616//RJPP-BEGIN------------------------------------------------------------ 617 for (int i = 0; i < x.length; i++) 618 x.put(i, (float) Math.cosh(x.get(i))); 619 return x; 620//RJPP-END-------------------------------------------------------------- 621 } 622 623 /** 624 * Applies the <i>exponential</i> function element-wise on this 625 * matrix. Note that this is an in-place operation. 626 * @see MatrixFunctions#exp(FloatMatrix) 627 * @return this matrix 628 */ 629 public static FloatMatrix expi(FloatMatrix x) { 630 /*# mapfct('Math.exp') #*/ 631//RJPP-BEGIN------------------------------------------------------------ 632 for (int i = 0; i < x.length; i++) 633 x.put(i, (float) Math.exp(x.get(i))); 634 return x; 635//RJPP-END-------------------------------------------------------------- 636 } 637 638 /** 639 * Element-wise round down by applying the <i>floor</i> function on each 640 * element. Note that this is an in-place operation. 641 * @see MatrixFunctions#floor(FloatMatrix) 642 * @return this matrix 643 */ 644 public static FloatMatrix floori(FloatMatrix x) { 645 /*# mapfct('Math.floor') #*/ 646//RJPP-BEGIN------------------------------------------------------------ 647 for (int i = 0; i < x.length; i++) 648 x.put(i, (float) Math.floor(x.get(i))); 649 return x; 650//RJPP-END-------------------------------------------------------------- 651 } 652 653 /** 654 * Applies the <i>natural logarithm</i> function element-wise on this 655 * matrix. Note that this is an in-place operation. 656 * @see MatrixFunctions#log(FloatMatrix) 657 * @return this matrix 658 */ 659 public static FloatMatrix logi(FloatMatrix x) { 660 /*# mapfct('Math.log') #*/ 661//RJPP-BEGIN------------------------------------------------------------ 662 for (int i = 0; i < x.length; i++) 663 x.put(i, (float) Math.log(x.get(i))); 664 return x; 665//RJPP-END-------------------------------------------------------------- 666 } 667 668 /** 669 * Applies the <i>logarithm with basis to 10</i> element-wise on this 670 * matrix. Note that this is an in-place operation. 671 * @see MatrixFunctions#log10(FloatMatrix) 672 * @return this matrix 673 */ 674 public static FloatMatrix log10i(FloatMatrix x) { 675 /*# mapfct('Math.log10') #*/ 676//RJPP-BEGIN------------------------------------------------------------ 677 for (int i = 0; i < x.length; i++) 678 x.put(i, (float) Math.log10(x.get(i))); 679 return x; 680//RJPP-END-------------------------------------------------------------- 681 } 682 683 /** 684 * Element-wise power function. Replaces each element with its 685 * power of <tt>d</tt>.Note that this is an in-place operation. 686 * @param d the exponent 687 * @see MatrixFunctions#pow(FloatMatrix,float) 688 * @return this matrix 689 */ 690 public static FloatMatrix powi(FloatMatrix x, float d) { 691 if (d == 2.0f) 692 return x.muli(x); 693 else { 694 for (int i = 0; i < x.length; i++) 695 x.put(i, (float) Math.pow(x.get(i), d)); 696 return x; 697 } 698 } 699 700 public static FloatMatrix powi(float base, FloatMatrix x) { 701 for (int i = 0; i < x.length; i++) 702 x.put(i, (float) Math.pow(base, x.get(i))); 703 return x; 704 } 705 706 public static FloatMatrix powi(FloatMatrix x, FloatMatrix e) { 707 x.checkLength(e.length); 708 for (int i = 0; i < x.length; i++) 709 x.put(i, (float) Math.pow(x.get(i), e.get(i))); 710 return x; 711 } 712 713 public static FloatMatrix signumi(FloatMatrix x) { 714 /*# mapfct('Math.signum') #*/ 715//RJPP-BEGIN------------------------------------------------------------ 716 for (int i = 0; i < x.length; i++) 717 x.put(i, (float) Math.signum(x.get(i))); 718 return x; 719//RJPP-END-------------------------------------------------------------- 720 } 721 722 public static FloatMatrix sini(FloatMatrix x) { 723 /*# mapfct('Math.sin') #*/ 724//RJPP-BEGIN------------------------------------------------------------ 725 for (int i = 0; i < x.length; i++) 726 x.put(i, (float) Math.sin(x.get(i))); 727 return x; 728//RJPP-END-------------------------------------------------------------- 729 } 730 731 public static FloatMatrix sinhi(FloatMatrix x) { 732 /*# mapfct('Math.sinh') #*/ 733//RJPP-BEGIN------------------------------------------------------------ 734 for (int i = 0; i < x.length; i++) 735 x.put(i, (float) Math.sinh(x.get(i))); 736 return x; 737//RJPP-END-------------------------------------------------------------- 738 } 739 public static FloatMatrix sqrti(FloatMatrix x) { 740 /*# mapfct('Math.sqrt') #*/ 741//RJPP-BEGIN------------------------------------------------------------ 742 for (int i = 0; i < x.length; i++) 743 x.put(i, (float) Math.sqrt(x.get(i))); 744 return x; 745//RJPP-END-------------------------------------------------------------- 746 } 747 public static FloatMatrix tani(FloatMatrix x) { 748 /*# mapfct('Math.tan') #*/ 749//RJPP-BEGIN------------------------------------------------------------ 750 for (int i = 0; i < x.length; i++) 751 x.put(i, (float) Math.tan(x.get(i))); 752 return x; 753//RJPP-END-------------------------------------------------------------- 754 } 755 public static FloatMatrix tanhi(FloatMatrix x) { 756 /*# mapfct('Math.tanh') #*/ 757//RJPP-BEGIN------------------------------------------------------------ 758 for (int i = 0; i < x.length; i++) 759 x.put(i, (float) Math.tanh(x.get(i))); 760 return x; 761//RJPP-END-------------------------------------------------------------- 762 } 763 764 /** 765 * Returns a copy of this matrix where all elements are set to their 766 * absolute values. 767 * @see MatrixFunctions#absi(FloatMatrix) 768 * @return copy of this matrix 769 */ 770 public static FloatMatrix abs(FloatMatrix x) { return absi(x.dup()); } 771 772 /** 773 * Returns a copy of this matrix where the trigonometric <i>acos</i> function is applied 774 * element wise. 775 * @see MatrixFunctions#acosi(FloatMatrix) 776 * @return copy of this matrix 777 */ 778 public static FloatMatrix acos(FloatMatrix x) { return acosi(x.dup()); } 779 public static FloatMatrix asin(FloatMatrix x) { return asini(x.dup()); } 780 public static FloatMatrix atan(FloatMatrix x) { return atani(x.dup()); } 781 public static FloatMatrix cbrt(FloatMatrix x) { return cbrti(x.dup()); } 782 public static FloatMatrix ceil(FloatMatrix x) { return ceili(x.dup()); } 783 public static FloatMatrix cos(FloatMatrix x) { return cosi(x.dup()); } 784 public static FloatMatrix cosh(FloatMatrix x) { return coshi(x.dup()); } 785 public static FloatMatrix exp(FloatMatrix x) { return expi(x.dup()); } 786 public static FloatMatrix floor(FloatMatrix x) { return floori(x.dup()); } 787 public static FloatMatrix log(FloatMatrix x) { return logi(x.dup()); } 788 public static FloatMatrix log10(FloatMatrix x) { return log10i(x.dup()); } 789 public static float pow(float x, float y) { return (float)Math.pow(x, y); } 790 public static FloatMatrix pow(FloatMatrix x, float e) { return powi(x.dup(), e); } 791 public static FloatMatrix pow(float b, FloatMatrix x) { return powi(b, x.dup()); } 792 public static FloatMatrix pow(FloatMatrix x, FloatMatrix e) { return powi(x.dup(), e); } 793 public static FloatMatrix signum(FloatMatrix x) { return signumi(x.dup()); } 794 public static FloatMatrix sin(FloatMatrix x) { return sini(x.dup()); } 795 public static FloatMatrix sinh(FloatMatrix x) { return sinhi(x.dup()); } 796 public static FloatMatrix sqrt(FloatMatrix x) { return sqrti(x.dup()); } 797 public static FloatMatrix tan(FloatMatrix x) { return tani(x.dup()); } 798 public static FloatMatrix tanh(FloatMatrix x) { return tanhi(x.dup()); } 799 800 /*# %w{abs acos asin atan cbrt ceil cos cosh exp floor log log10 signum sin sinh sqrt tan tanh}.map do |fct| <<-EOS 801 public static float #{fct}(float x) { return (float)Math.#{fct}(x); } 802 EOS 803 end 804 #*/ 805//RJPP-BEGIN------------------------------------------------------------ 806 public static float abs(float x) { return (float)Math.abs(x); } 807 public static float acos(float x) { return (float)Math.acos(x); } 808 public static float asin(float x) { return (float)Math.asin(x); } 809 public static float atan(float x) { return (float)Math.atan(x); } 810 public static float cbrt(float x) { return (float)Math.cbrt(x); } 811 public static float ceil(float x) { return (float)Math.ceil(x); } 812 public static float cos(float x) { return (float)Math.cos(x); } 813 public static float cosh(float x) { return (float)Math.cosh(x); } 814 public static float exp(float x) { return (float)Math.exp(x); } 815 public static float floor(float x) { return (float)Math.floor(x); } 816 public static float log(float x) { return (float)Math.log(x); } 817 public static float log10(float x) { return (float)Math.log10(x); } 818 public static float signum(float x) { return (float)Math.signum(x); } 819 public static float sin(float x) { return (float)Math.sin(x); } 820 public static float sinh(float x) { return (float)Math.sinh(x); } 821 public static float sqrt(float x) { return (float)Math.sqrt(x); } 822 public static float tan(float x) { return (float)Math.tan(x); } 823 public static float tanh(float x) { return (float)Math.tanh(x); } 824//RJPP-END-------------------------------------------------------------- 825 826 /** 827 * Calculate matrix exponential of a square matrix. 828 * 829 * A scaled Pade approximation algorithm is used. 830 * The algorithm has been directly translated from Golub & Van Loan "Matrix Computations", 831 * algorithm 11.3f.1. Special Horner techniques from 11.2f are also used to minimize the number 832 * of matrix multiplications. 833 * 834 * @param A square matrix 835 * @return matrix exponential of A 836 */ 837 public static FloatMatrix expm(FloatMatrix A) { 838 // constants for pade approximation 839 final float c0 = 1.0f; 840 final float c1 = 0.5f; 841 final float c2 = 0.12f; 842 final float c3 = 0.01833333333333333f; 843 final float c4 = 0.0019927536231884053f; 844 final float c5 = 1.630434782608695E-4f; 845 final float c6 = 1.0351966873706E-5f; 846 final float c7 = 5.175983436853E-7f; 847 final float c8 = 2.0431513566525E-8f; 848 final float c9 = 6.306022705717593E-10f; 849 final float c10 = 1.4837700484041396E-11f; 850 final float c11 = 2.5291534915979653E-13f; 851 final float c12 = 2.8101705462199615E-15f; 852 final float c13 = 1.5440497506703084E-17f; 853 854 int j = Math.max(0, 1 + (int) Math.floor(Math.log(A.normmax()) / Math.log(2))); 855 FloatMatrix As = A.div((float) Math.pow(2, j)); // scaled version of A 856 int n = A.getRows(); 857 858 // calculate D and N using special Horner techniques 859 FloatMatrix As_2 = As.mmul(As); 860 FloatMatrix As_4 = As_2.mmul(As_2); 861 FloatMatrix As_6 = As_4.mmul(As_2); 862 // U = c0*I + c2*A^2 + c4*A^4 + (c6*I + c8*A^2 + c10*A^4 + c12*A^6)*A^6 863 FloatMatrix U = FloatMatrix.eye(n).muli(c0).addi(As_2.mul(c2)).addi(As_4.mul(c4)).addi( 864 FloatMatrix.eye(n).muli(c6).addi(As_2.mul(c8)).addi(As_4.mul(c10)).addi(As_6.mul(c12)).mmuli(As_6)); 865 // V = c1*I + c3*A^2 + c5*A^4 + (c7*I + c9*A^2 + c11*A^4 + c13*A^6)*A^6 866 FloatMatrix V = FloatMatrix.eye(n).muli(c1).addi(As_2.mul(c3)).addi(As_4.mul(c5)).addi( 867 FloatMatrix.eye(n).muli(c7).addi(As_2.mul(c9)).addi(As_4.mul(c11)).addi(As_6.mul(c13)).mmuli(As_6)); 868 869 FloatMatrix AV = As.mmuli(V); 870 FloatMatrix N = U.add(AV); 871 FloatMatrix D = U.subi(AV); 872 873 // solve DF = N for F 874 FloatMatrix F = Solve.solve(D, N); 875 876 // now square j times 877 for (int k = 0; k < j; k++) { 878 F.mmuli(F); 879 } 880 881 return F; 882 } 883 884 885 886//END 887}