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}