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