# 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

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

variables:

# 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 = nsum_x_sq[i] - sum_x[i]**2
dy = n
sum_y_sq[j] - sum_y[j]**2
corr[i][j] = ( nsum_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

On 8/21/06, Bil K. [email protected] 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

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 + e2 }
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

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?

Tomasz W. wrote:

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

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 + e2 }
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,

[email protected] 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.

Bil K. [email protected] 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?

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,

[email protected] 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,

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

Bil K. wrote:

http://fun3d.larc.nasa.gov

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

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

M. Edward (Ed) Borasky wrote:

[email protected] wrote:

[Bill 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’?”

Bil, Can I ask why you didn’t you just call your fortran library
from Ruby by wrapping it as an extension?

# I have no problem wrapping this little Fortran example:

c Demo of Ruby/Fortran

``````function hello_ruby_func(i) result(j)
integer, intent(in) :: i ! input
integer             :: j ! output
j = i + 123
end function hello_ruby_func
``````

# ============================================================= with this little Ruby/C extension which calls my Fortran function

#include “ruby.h”

VALUE Ruby2C2Fortran = Qnil;
void Init_mytest();

VALUE method_test_f(VALUE self) {
int x = 10;
int y = hello_ruby_func_(&x);
return INT2NUM(y);
}

# require ‘ruby2c2fortran’ include Ruby2C2Fortran puts test_f

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.

OTOH, there are still some good fortran libraries where it’s
hard to find equivalents even in R.

On 8/21/06, Bil K. [email protected] 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’?”

variables:

# 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 = nsum_x_sq[i] - sum_x[i]**2
dy = n
sum_y_sq[j] - sum_y[j]**2
(dxdy).zero? ? 0.0 :
( n
sumxy - sum_x[i]sum_y[j] ) / Math.sqrt(dxdy)
}

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

Bil K. wrote:

## Thanks,

Alex G. (I think) is porting RSPerl to Ruby. He has the “calling
R from Ruby piece” working, but not the “calling Ruby from R” piece,
IIRC.

My question as an R programmer would be, "what can Ruby do that R can’t
do (awkwardly, in some cases)?

Ron M [email protected] writes:

Bil, Can I ask why you didn’t you just call your fortran library
from Ruby by wrapping it as an extension?

Or even use Ruby/DL?

Ron M wrote:

OTOH, there are still some good fortran libraries where it’s
hard to find equivalents even in R.
R has the ability (on *nix, at any rate) to call both C/C++ and Fortran
libraries easily. So … you got good Fortran, call it from R or Ruby or
Python or Perl or whatever. Given a C wrapper for the Fortran and SWIG,
you can perform all sorts of miracles.

Ron M wrote:

Bil, Can I ask why you didn’t you just call your fortran library
from Ruby by wrapping it as an extension?

-ilities: simplicity, maintainability, and portability.

(Thanks for the example though!)

Regards,