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    }