[lnkForumImage]
TotalShareware - Download Free Software
Confronta i prezzi di migliaia di prodotti.
Asp Forum
 Home | Login | Register | Search 


 

Forums >

comp.lang.ruby

You can write Fortran in any language

Bil Kleb

8/21/2006 10:22:00 AM

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.lar...


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
15 Answers

Tomasz Wegrzanowski

8/21/2006 11:04:00 AM

0

On 8/21/06, Bil Kleb <Bil.Kleb@nasa.gov> wrote:
> 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'?"
[cut]
> for i in 0..x.size-1

Well, it's a very small thing, but tripple-dot makes is somewhat more readable:

for i in 0...x.size

And I'd actually use x.size.times {|i| }, for loops looks kinda
alien in Ruby code.

Now pieces of code like:
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
can be rewritten to:
sum_x = (0...x.size).map{|i| x[i].inject(0) { |sum, e| sum + e }}
sum_x_sq = (0...x.size).map{|i| x[i].inject(0) { |sum, e| sum + e**2 }}
which looks a lot neater imho.
It would look even neater to write:
sum_x = x.map{|xi| xi.sum}
sum_x_sj = x.map{|xi| xi.map{|e|e**2}.sum}
but Ruby doesn't have Enumerable#sum by default.

Oh well, you have a few ideas of making the code nicer now :-)

--
Tomasz Wegrzanowski [ http://t-a-w.blo... ]

Bil Kleb

8/21/2006 1:08:00 PM

0

Tomasz Wegrzanowski wrote:
>
> Well, it's a very small thing, but tripple-dot makes is somewhat more
> readable:

Good point.

> 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
> can be rewritten to:
> sum_x = (0...x.size).map{|i| x[i].inject(0) { |sum, e| sum + e }}
> sum_x_sq = (0...x.size).map{|i| x[i].inject(0) { |sum, e| sum + e**2 }}
> which looks a lot neater imho.

Agreed.

How about the (i,j) beasts? E.g.,

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

Doubly embedded maps?

Is there anyway to do all this without using indices?

Maybe intercept the incoming vectors, build /scalar/
correlations for each variable combination, and then
assemble the final /matrix/ of correlations?

> It would look even neater to write:
> sum_x = x.map{|xi| xi.sum}
> sum_x_sj = x.map{|xi| xi.map{|e|e**2}.sum}
> but Ruby doesn't have Enumerable#sum by default.

Agreed. I'm happy to locally extend Enumerable...

module Enumerable
def sum
self.inject(0){ |sum,e| sum + e }
end
end

> Oh well, you have a few ideas of making the code nicer now :-)

Yes. Thanks.

Regards,
--
Bil
http://fun3d.lar...

benjohn

8/21/2006 1:18:00 PM

0

> 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'?"

I know that there are good linear maths libraries for ruby. Why didn't
you just use them? :)

M. Edward (Ed) Borasky

8/21/2006 1:29:00 PM

0

benjohn@fysh.org wrote:
>> 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'?"
>
> I know that there are good linear maths libraries for ruby. Why didn't
> you just use them? :)
>
>
>
Well ... ok ... but why even write using something as low-level as a
linear algebra library. I think this is a one-liner in R.

:)

Bil Kleb

8/21/2006 1:38:00 PM

0

M. Edward (Ed) Borasky wrote:
>
> Well ... ok ... but why even write using something as low-level as a
> linear algebra library. I think this is a one-liner in R.

I thought I should be using R, but after a brief search
for Ruby bindings, installation instructions, example usage,
etc., I gave up.

Can you point me to a "How to Use R with Ruby"?

Thanks,
--
Bil
http://fun3d.lar...

Bil Kleb

8/21/2006 1:40:00 PM

0

benjohn@fysh.org wrote:
>
> I know that there are good linear maths libraries for ruby. Why didn't
> you just use them? :)

Exactly. But which one? and how?

Later,
--
Bil
http://fun3d.lar...

Christian Neukirchen

8/21/2006 1:45:00 PM

0

Bil Kleb <Bil.Kleb@NASA.gov> writes:

> 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'?"

Wouldn't you be better off using NArray for the matrix stuff?

> Thanks,
> --
> Bil

--
Christian Neukirchen <chneukirchen@gmail.com> http://chneuk...

Kroeger, Simon (ext)

8/21/2006 2:01:00 PM

0



> From: Bil Kleb [mailto:Bil.Kleb@NASA.gov]
> Sent: Monday, August 21, 2006 3:20 PM
>
> How about the (i,j) beasts? E.g.,
>
> sum_xy = Array.new(x.size){ Array.new(y.size){0.0} }
> for k in 0..n-1
> for i in 0...x.size
> for j in 0...y.size
> sum_xy[i][j] += x[i][k]*y[j][k]
> end
> end
> end
>
> Doubly embedded maps?
>
> Is there anyway to do all this without using indices?

sum_xy = x.map do |i|
y.map do |j|
(0..n-1).inject(0) {|sum, k| sum + i[k] * j[k]}
end
end

to get rid of the last index you would need zip or
SyncEnumerator (which is slow).

> Agreed. I'm happy to locally extend Enumerable...
>
> module Enumerable
> def sum
> self.inject(0){ |sum,e| sum + e }
> end
> end

module Enumerable
def sum
inject{|sum, e| sum + e}
end
end

hmm, at least if the enumerable isn't empty ...

cheers

Simon

William James

8/21/2006 3:50:00 PM

0

Bil Kleb wrote:
> 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.lar...
>

>
> 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_x = x.map{|a| a.inject(0){|sum,e| sum + e} }
sum_x_sq = x.map{|a| a.inject(0){|sum,e| sum + e**2 } }


>
> 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

sum_xy = x.map{|xa| y.map{|ya|
xa.zip(ya).inject(0){|sum,na| sum + na.first*na.last}}}

Martin DeMello

8/21/2006 4:49:00 PM

0

On 8/21/06, Bil Kleb <Bil.Kleb@nasa.gov> wrote:
> 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'?"

# 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] ]

class Array
def sum
if block_given?
inject(0) {|s, i| s + yield(i)}
else
inject(0) {|s, i| s + i}
end
end

def sums(&blk)
map {|i| i.sum(&blk)}
end

def map2d_with_indices
a = []
each_with_index {|x,i|
a[i] = []
x.each_with_index {|y,j| a[i][j] = yield [y,i,j] }
}
a
end
end

def correlation(x,y)

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

n = x.first.size

sum_x = x.sums
sum_x_sq = x.sums {|e| e**2}

sum_y = y.sums
sum_y_sq = y.sums {|e| e**2}

sum_xy = x.map { y.map { 0.0 } }

0.upto(n-1) do |k|
x.each_index do |i|
y.each_index do |j|
sum_xy[i][j] += x[i][k]*y[j][k]
end
end
end

corr = sum_xy.map2d_with_indices {|sumxy, i, j|
dx = n*sum_x_sq[i] - sum_x[i]**2
dy = n*sum_y_sq[j] - sum_y[j]**2
(dx*dy).zero? ? 0.0 :
( n*sumxy - sum_x[i]*sum_y[j] ) / Math.sqrt(dx*dy)
}

sum_corr_sq_y = corr.transpose.sums {|e| e**2}

corr_sq = corr.map2d_with_indices {|e, i, j|
sum_corr_sq_y[j].zero? ? 0.0 : (e**2)/sum_corr_sq_y[j]
}

[corr, corr_sq]
end