1+ use nalgebralib:: { Dyn , OMatrix } ;
12use rhai:: plugin:: * ;
23
34#[ export_module]
45pub mod matrix_functions {
5- #[ cfg( feature = "smartcore" ) ]
6- use crate :: { dense_matrix_to_vec_dynamic, if_matrix_convert_to_dense_matrix_and_do} ;
76 use crate :: {
87 if_int_convert_to_float_and_do, if_int_do_else_if_array_do, if_list_do,
98 if_matrix_convert_to_vec_array_and_do,
@@ -16,10 +15,6 @@ pub mod matrix_functions {
1615 #[ cfg( feature = "nalgebra" ) ]
1716 use nalgebralib:: DMatrix ;
1817 use rhai:: { Array , Dynamic , EvalAltResult , Map , Position , FLOAT , INT } ;
19- #[ cfg( feature = "smartcore" ) ]
20- use smartcorelib:: linalg:: basic:: arrays:: Array2 ;
21- #[ cfg( feature = "smartcore" ) ]
22- use smartcorelib:: linalg:: traits:: evd:: EVDDecomposable ;
2318 use std:: collections:: BTreeMap ;
2419
2520 /// Calculates the inverse of a matrix. Fails if the matrix if not invertible, or if the
@@ -66,58 +61,90 @@ pub mod matrix_functions {
6661 } )
6762 } )
6863 }
69- /// Calculate the eigenvalues for a matrix.
64+
65+ /// Calculate the eigenvalues and eigenvectors for a matrix. Specifically, the output is an
66+ /// object map with entries for real_eigenvalues, imaginary_eigenvalues, eigenvectors, and
67+ /// residuals.
7068 /// ```typescript
7169 /// let matrix = eye(5);
7270 /// let eig = eigs(matrix);
73- /// assert_eq( eig, #{ "eigenvectors": [[1.0, 0.0, 0.0, 0.0, 0.0],
74- /// [0.0, 1.0, 0.0, 0.0, 0.0],
75- /// [0.0, 0.0, 1.0, 0.0, 0.0],
76- /// [0.0, 0.0, 0.0, 1.0, 0 .0],
77- /// [0 .0, 0.0, 0.0, 0.0, 1.0]],
78- /// "imaginary_eigenvalues": [0.0, 0.0, 0.0, 0.0, 0.0],
79- /// "real_eigenvalues": [1.0, 1.0, 1.0, 1.0, 1.0]} );
71+ /// assert(sum( eig.residuals) < 0.000001);
72+ /// ```
73+ /// ```typescript
74+ /// let matrix = [[ 0.0, 1 .0],
75+ /// [-2 .0, -3.0]];
76+ /// let eig = eigs(matrix);
77+ /// assert(sum(eig.residuals) < 0.000001 );
8078 /// ```
81- #[ cfg( feature = "smartcore " ) ]
79+ #[ cfg( feature = "nalgebra " ) ]
8280 #[ rhai_fn( name = "eigs" , return_raw, pure) ]
83- pub fn matrix_eigs ( matrix : & mut Array ) -> Result < Map , Box < EvalAltResult > > {
84- if_matrix_convert_to_dense_matrix_and_do ( matrix, |matrix_as_dm| {
85- // Try to invert
86- let dm =
87- matrix_as_dm. evd ( matrix_as_dm. approximate_eq ( & matrix_as_dm. transpose ( ) , 0.0000001 ) ) ;
88-
89- match dm {
90- Err ( e) => {
91- Err ( EvalAltResult :: ErrorArithmetic ( format ! ( "{:?}" , e) , Position :: NONE ) . into ( ) )
81+ pub fn matrix_eigs_alt ( matrix : & mut Array ) -> Result < Map , Box < EvalAltResult > > {
82+ if_matrix_convert_to_vec_array_and_do ( matrix, |matrix_as_vec| {
83+ // Convert vec_array to omatrix
84+ let mut dm = DMatrix :: from_fn ( matrix_as_vec. len ( ) , matrix_as_vec[ 0 ] . len ( ) , |i, j| {
85+ if matrix_as_vec[ 0 ] [ 0 ] . is_float ( ) {
86+ matrix_as_vec[ i] [ j] . as_float ( ) . unwrap ( )
87+ } else {
88+ matrix_as_vec[ i] [ j] . as_int ( ) . unwrap ( ) as FLOAT
9289 }
90+ } ) ;
9391
94- Ok ( evd) => {
95- let vecs: Array = dense_matrix_to_vec_dynamic ( evd. V ) ;
96- let real_values: Array = evd
97- . d
98- . into_iter ( )
99- . map ( |x| Dynamic :: from_float ( x) )
100- . collect :: < Vec < Dynamic > > ( ) ;
101- let imaginary_values: Array = evd
102- . e
103- . into_iter ( )
104- . map ( |x| Dynamic :: from_float ( x) )
105- . collect :: < Vec < Dynamic > > ( ) ;
106-
107- let mut result = BTreeMap :: new ( ) ;
108- let mut vid = smartstring:: SmartString :: new ( ) ;
109- vid. push_str ( "eigenvectors" ) ;
110- result. insert ( vid, Dynamic :: from_array ( vecs) ) ;
111- let mut did = smartstring:: SmartString :: new ( ) ;
112- did. push_str ( "real_eigenvalues" ) ;
113- result. insert ( did, Dynamic :: from_array ( real_values) ) ;
114- let mut eid = smartstring:: SmartString :: new ( ) ;
115- eid. push_str ( "imaginary_eigenvalues" ) ;
116- result. insert ( eid, Dynamic :: from_array ( imaginary_values) ) ;
117-
118- Ok ( result)
119- }
92+ // Grab shape for later
93+ let dms = dm. shape ( ) . 1 ;
94+
95+ // Get teh eigenvalues
96+ let eigenvalues = dm. complex_eigenvalues ( ) ;
97+
98+ // Iterate through eigenvalues to get eigenvectors
99+ let mut imaginary_values = vec ! [ Dynamic :: from_float( 1.0 ) ; 0 ] ;
100+ let mut real_values = vec ! [ Dynamic :: from_float( 1.0 ) ; 0 ] ;
101+ let mut residuals = vec ! [ Dynamic :: from_float( 1.0 ) ; 0 ] ;
102+ let mut eigenvectors = DMatrix :: from_element ( dms, 0 , 0.0 ) ;
103+ for ( idx, ev) in eigenvalues. iter ( ) . enumerate ( ) {
104+ // Eigenvalue components
105+ imaginary_values. push ( Dynamic :: from_float ( ev. im ) ) ;
106+ real_values. push ( Dynamic :: from_float ( ev. re ) ) ;
107+
108+ // Get eigenvector
109+ let mut A = dm. clone ( ) - DMatrix :: from_diagonal_element ( dms, dms, ev. re ) ;
110+ A = A . insert_column ( 0 , 0.0 ) ;
111+ A = A . insert_row ( 0 , 0.0 ) ;
112+ A [ ( 0 , idx + 1 ) ] = 1.0 ;
113+ let mut b = DMatrix :: from_element ( dms + 1 , 1 , 0.0 ) ;
114+ b[ ( 0 , 0 ) ] = 1.0 ;
115+ let eigenvector = A
116+ . svd ( true , true )
117+ . solve ( & b, 1e-10 )
118+ . unwrap ( )
119+ . remove_rows ( 0 , 1 )
120+ . normalize ( ) ;
121+
122+ // Verify solution
123+ residuals. push ( Dynamic :: from_float (
124+ ( dm. clone ( ) * eigenvector. clone ( ) - ev. re * eigenvector. clone ( ) ) . amax ( ) ,
125+ ) ) ;
126+
127+ eigenvectors. extend ( eigenvector. column_iter ( ) ) ;
120128 }
129+
130+ let mut result = BTreeMap :: new ( ) ;
131+ let mut vid = smartstring:: SmartString :: new ( ) ;
132+ vid. push_str ( "eigenvectors" ) ;
133+ result. insert (
134+ vid,
135+ Dynamic :: from_array ( omatrix_to_vec_dynamic ( eigenvectors) ) ,
136+ ) ;
137+ let mut did = smartstring:: SmartString :: new ( ) ;
138+ did. push_str ( "real_eigenvalues" ) ;
139+ result. insert ( did, Dynamic :: from_array ( real_values) ) ;
140+ let mut eid = smartstring:: SmartString :: new ( ) ;
141+ eid. push_str ( "imaginary_eigenvalues" ) ;
142+ result. insert ( eid, Dynamic :: from_array ( imaginary_values) ) ;
143+ let mut rid = smartstring:: SmartString :: new ( ) ;
144+ rid. push_str ( "residuals" ) ;
145+ result. insert ( rid, Dynamic :: from_array ( residuals) ) ;
146+
147+ Ok ( result)
121148 } )
122149 }
123150
0 commit comments