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