From: Bil Kleb <Bil.Kleb@...>
Date: 2006-08-21T19:35:09+09:00
Subject: You can write Fortran in any language

Here is an example used for computing correlations.
It's almost a direct translation of several hundred
lines of Fortran 77.

If you can bare to look, please help, even it is simply
to say "yo bonehead, why didn't you just require 'some_library'?"

Thanks,
--
Bil
http://fun3d.larc.nasa.gov


cat << EOF > correlation_test.rb
require 'test/unit'
require 'correlation'

class TestCorrelation < Test::Unit::TestCase
   def test_1x1_perfectly_correlated
     x = [ [0,1] ]
     y = [ [0,1] ]
     c, c2 = correlation(x,y)
     assert_equal( [ [ 1.0 ] ], c )
     assert_equal( [ [ 1.0 ] ], c2 )
   end
   def test_1x1_perfectly_anti_correlated
     x = [ [0,1] ]
     y = [ [1,0] ]
     c, c2 = correlation(x,y)
     assert_equal( [ [ -1.0 ] ], c )
     assert_equal( [ [  1.0 ] ], c2 )
   end
   def test_1x1_perfectly_uncorrelated
     x = [ [0,1] ]
     y = [ [0,0] ]
     c, c2 = correlation(x,y)
     assert_equal( [ [ 0.0 ] ], c )
     assert_equal( [ [ 0.0 ] ], c2 )
   end
   def test_2x1_perfectly_correlated
     x = [ [0,1], [0,1] ]
     y = [ [0,1] ]
     c, c2 = correlation(x,y)
     assert_equal( [ [ 1.0 ], [ 1.0 ] ], c )
     assert_equal( [ [ 0.5 ], [ 0.5 ] ], c2 )
   end
   def test_2x1_perfectly_correlated_and_uncorrelated
     x = [ [0,1], [0,0] ]
     y = [ [0,1] ]
     c, c2 = correlation(x,y)
     assert_equal( [ [ 1.0 ], [ 0.0 ] ], c )
     assert_equal( [ [ 1.0 ], [ 0.0 ] ], c2 )
   end
   def test_1x2_perfectly_anti_correlated
     x = [ [0,1] ]
     y = [ [1,0], [1,0] ]
     c, c2 = correlation(x,y)
     assert_equal( [ [ -1.0, -1.0 ] ], c )
     assert_equal( [ [  1.0,  1.0 ] ], c2 )
   end
   def test_1x2_perfectly_correlated_and_uncorrelated
     x = [ [0,1] ]
     y = [ [0,1], [0,0] ]
     c, c2 = correlation(x,y)
     assert_equal( [ [ 1.0, 0.0 ] ], c )
     assert_equal( [ [ 1.0, 0.0 ] ], c2 )
   end
end
EOF


cat << EOF > correlation.rb
# Returns an array with correlation and y-normalized correlation^2
# Assumes two nested arrays of variables and samples.  For example,
# here's a case with 3 samples of 3 input variables and 2 output variables:
#
#   x = [ [a1,a2,a3], [b1,b2,b3], [c1,c2,c3] ]
#   y = [ [r1,r2,r3], [s1,s2,s3] ]

def correlation(x,y)

   fail 'size arrays unequal' if x.first.size != y.first.size # lame

   n = x.first.size

   sum_x = Array.new(x.size){0.0}
   sum_x_sq = Array.new(x.size){0.0}
   for i in 0..x.size-1
     sum_x[i]    = x[i].inject(0) { |sum, e| sum + e }
     sum_x_sq[i] = x[i].inject(0) { |sum, e| sum + e**2 }
   end

   sum_y = Array.new(y.size){0.0}
   sum_y_sq = Array.new(y.size){0.0}
   for i in 0..y.size-1
     sum_y[i]    = y[i].inject(0) { |sum, e| sum + e }
     sum_y_sq[i] = y[i].inject(0) { |sum, e| sum + e**2 }
   end

   sum_xy = Array.new(x.size){ Array.new(y.size){0.0} }
   for k in 0..n-1
     for i in 0..x.size-1
       for j in 0..y.size-1
         sum_xy[i][j] += x[i][k]*y[j][k]
       end
     end
   end

   corr = Array.new(x.size){ Array.new(y.size){0.0} }
   for i in 0..x.size-1
     for j in 0..y.size-1
       dx = n*sum_x_sq[i] - sum_x[i]**2
       dy = n*sum_y_sq[j] - sum_y[j]**2
       corr[i][j] = ( n*sum_xy[i][j] - sum_x[i]*sum_y[j] ) / Math.sqrt(dx*dy) \
         unless dx*dy==0.0
     end
   end

   sum_corr_sq_y = Array.new(y.size){0.0}
   for j in 0..y.size-1
     for i in 0..x.size-1
       sum_corr_sq_y[j] += corr[i][j]**2
     end
   end

   corr_sq = Array.new(x.size){ Array.new(y.size){0.0} }
   for i in 0..x.size-1
     for j in 0..y.size-1
       corr_sq[i][j] = corr[i][j]**2/sum_corr_sq_y[j] \
         unless sum_corr_sq_y[j]==0.0
     end
   end

   [corr, corr_sq]

end
EOF