You can write Fortran in any language

(BIl Kleb) #1

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

(Tomasz Wegrzanowski) #2

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

···

On 8/21/06, Bil Kleb <Bil.Kleb@nasa.gov> wrote:

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

(Benjohn Barnes) #3

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? :slight_smile:

(Christian Neukirchen) #4

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://chneukirchen.org

(W. James) #5

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.larc.nasa.gov

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) #6

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

···

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

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

(BIl Kleb) #7

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

Yes. Thanks.

Regards,

···

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

(M. Edward (Ed) Borasky) #8

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? :slight_smile:

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.

:slight_smile:

(BIl Kleb) #9

benjohn@fysh.org wrote:

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

Exactly. But which one? and how?

Later,

···

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

(BIl Kleb) #10

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.larc.nasa.gov

(Kroeger Simon (ext)) #11

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

(Ron M) #12

M. Edward (Ed) Borasky wrote:

···

benjohn@fysh.org 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);
}

void Init_ruby2c2fortran() {
  Ruby2C2Fortran = rb_define_module("Ruby2C2Fortran");
  rb_define_method(Ruby2C2Fortran, "test_f", method_test_f, 0);
}

and calling it as follows:

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.

(Christian Neukirchen) #13

Ron M <rm_rails@cheapcomplexdevices.com> writes:

M. Edward (Ed) Borasky 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?

Or even use Ruby/DL?

···

benjohn@fysh.org wrote:

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

(M. Edward (Ed) Borasky) #14

Bil Kleb wrote:

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.larc.nasa.gov

Alex Gutteridge (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)? :slight_smile:

(M. Edward (Ed) Borasky) #15

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.

(BIl Kleb) #16

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,

···

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