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
037 package org.jblas;
038
039 import 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 */
047 public 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 }