lib/statsample/factor/rotation.rb in statsample-0.18.0 vs lib/statsample/factor/rotation.rb in statsample-1.0.0

@@ -25,11 +25,11 @@
     attr_reader :iterations, :rotated, :component_transformation_matrix, :h2
     # Maximum number of iterations    
     attr_accessor :max_iterations
     # Maximum precision    
     attr_accessor :epsilon
-    
+    attr_accessor :use_gsl
     dirty_writer :max_iterations, :epsilon
     dirty_memoize :iterations, :rotated, :component_transformation_matrix, :h2
     
     def initialize(matrix, opts=Hash.new)
       @name=_("%s rotation") % rotation_name
@@ -39,10 +39,11 @@
       @component_transformation_matrix=nil
       @max_iterations=MAX_ITERATIONS
       @epsilon=EPSILON
       @rotated=nil
       @h2=(@matrix.collect {|c| c**2} * Matrix.column_vector([1]*@m)).column(0).to_a
+      @use_gsl=Statsample.has_gsl?
       opts.each{|k,v|
         self.send("#{k}=",v) if self.respond_to? k
       }
     end
     def report_building(g)
@@ -56,15 +57,16 @@
     def compute
       iterate
     end
     # Start iteration 
     def iterate
-      t=Matrix.identity(@m)
-      b=@matrix.dup
-      h=Matrix.diagonal(*@h2).collect {|c| Math::sqrt(c)}
+      k_matrix=@use_gsl ? GSL::Matrix : ::Matrix
+      t=k_matrix.identity(@m)
+      b=(@use_gsl ? @matrix.to_gsl : @matrix.dup)
+      h=k_matrix.diagonal(*@h2).collect {|c| Math::sqrt(c)}
       h_inverse=h.collect {|c| c!=0 ? 1/c : 0 }
-      bh=h_inverse*b
+      bh=h_inverse * b
       @not_converged=true
       @iterations=0
       while @not_converged
         break if @iterations>@max_iterations
         @iterations+=1
@@ -108,12 +110,17 @@
               t=t.to_a
               @m.times {|row_i|
                 t[row_i][i]=tx_rot[row_i]
                 t[row_i][j]=ty_rot[row_i]
               }
-              
-              bh=Matrix.rows(bh)
-              t=Matrix.rows(t)
+              #if @use_gsl
+                bh=k_matrix.[](*bh)
+                t=k_matrix.[](*t)
+              #else
+              #  bh=Matrix.rows(bh)
+              #  t=Matrix.rows(t)
+                
+              #end
             else
               num_pairs=num_pairs-1
               @not_converged=false if num_pairs==0
             end # if
           end #j