Nmatrix compatibility

.travis.yml

@@ -2,23 +2,15 @@ language:
   ruby
 
 rvm:
-  - 1.9
+  - 1.9.3
   - 2.0
   - 2.1
   - 2.2
-  - ruby-head
-  - rbx-2
-
-env:
-  matrix:
-    - "NARRAY=1"
-    - ""
-
-matrix:
-  allow_failures:
-    - rvm: ruby-head
-    - rvm: rbx-2
 
 before_install:
   - sudo apt-get update -qq
   - sudo apt-get install -y libgsl0-dev
+  - gem update bundler
+  - bundle install
+  - bundle exec rake compile
+  - bundle exec rake test

ChangeLog

@@ -1,5 +1,11 @@
 # markup: rd
 
+== Fri March 4 2016
+* Added compatibility for NMatrix interconversion with GSL::Vector and GSL::Matrix
+* Added NMatrix compatibility to various to make them compatible with NMatrix data types.
+* Gem is now compatible with GSL 2.1.
+* Version bump to 2.1.0
+
 == Thu Jul 03 2015
 * Ruby/GSL 1.16.0.6
   * rb-gsl and gsl are now the same gem

Gemfile

@@ -2,3 +2,4 @@ source 'https://rubygems.org'
 gemspec
 
 gem 'narray' if ENV['NARRAY']
+gem 'nmatrix' if ENV['NMATRIX']

README.md

@@ -8,6 +8,14 @@ the terms of the GNU Free Documentation License.
 Ruby/GSL is a Ruby interface to the [GNU Scientific Library](https://gnu.org/software/gsl/)
 (GSL), for numerical computing with [Ruby](http://www.ruby-lang.org/).
 
+Ruby/GSL is compatible with GSL versions upto 2.1.
+
+## Usage with GSL 2.1
+
+As of this release, GSL 2.1 has not made it's way into the Debian stable repositories. Hence, after compiling GSL 2.1 from source, you will need to set the installation location in your `LD_LIBRARY_PATH` variable. After following standard GSL 2.1 installation procedures, you should do:
+  export LD_LIBRARY_PATH=/usr/local/lib
+
+The need to do this should not arise if GSL has been installed from `apt-get`.
 
 ## Installation
 
@@ -29,6 +37,66 @@ package. Some of the example scripts in the +examples/+ directory use the
 binaries of <tt>GNU plotutils</tt> and related packages are available
 [here](http://gnuwin32.sourceforge.net/packages/plotutils.htm).
 
+## NMatrix and NArray usage
+
+Ruby/GSL works with [NMatrix](https://github.com/SciRuby/nmatrix) and [NArray](https://github.com/masa16/narray) for a variety of methods. See the docs for a detailed list.
+
+### Basic Installation
+
+In order to use rb-gsl with NMatrix you must first set the `NMATRIX` environment variable and then install rb-gsl:
+  gem install nmatrix
+  export NMATRIX=1
+  gem install rb-gsl
+
+This will compile rb-gsl with NMatrix specific functions.
+
+For using rb-gsl with NArray:
+  gem install narray
+  export NARRAY=1
+  gem install rb-gsl
+
+Note that setting both `NMATRIX` and `NARRAY` variables will lead to undefined behaviour. Only one can be used at a time.
+
+### NMatrix basic usage
+
+Convert an NMatrix 1D vector to GSL::Vector:
+``` ruby
+require 'gsl'
+
+nm = NMatrix.new([5], [1,2,3,4,5], dtype: :float64)
+#=> [1.0, 2.0, 3.0, 4.0, 5.0]
+nm.to_gslv
+# => GSL::Vector
+# [ 1.000e+00 2.000e+00 3.000e+00 4.000e+00 5.000e+00 ]
+```
+
+Convert an integer 2D NMatrix to GSL::Matrix::Int:
+``` ruby
+
+require 'gsl'
+nm = NMatrix.new([3,3], [2]*9, dtype: :int32)
+#=> 
+#[
+#  [2, 2, 2]   [2, 2, 2]   [2, 2, 2] ]
+nm.to_gslm
+#=> GSL::Matrix::Int
+#[ 2 2 2 
+#  2 2 2 
+#  2 2 2 ]
+```
+
+Convert GSL::Vector to 1D NMatrix:
+``` ruby
+g = GSL::Vector.alloc(1,2,3,4)
+# => GSL::Vector
+# [ 1.000e+00 2.000e+00 3.000e+00 4.000e+00 ]
+g.to_nm
+# => [1.0, 2.0, 3.0, 4.0]
+```
+
+`to_nm` can be used on all sorts of `GSL::Vector` and `GSL::Matrix` objects to convert them to NMatrix.
+
+For a detailed list of methods that are compatible with NMatrix, see 'nmatrix' in the [docs](https://sciruby.github.com/rb-gsl).
 
 ## Reference
 
@@ -60,6 +128,14 @@ WITHOUT ANY WARRANTY.
 Any bug reports are welcome. If you encounter bugs in Ruby/GSL, please
 report them on GitHub(https://github.com/SciRuby/rb-gsl/issues).
 
+## Testing and Contributing
+
+If you wish to make contributions, run the following commands to clone and test the gem on your local machine:
+  git clone https://github.com/SciRuby/rb-gsl.git
+  cd rb-gsl
+  bash test.sh
+
+This will run tests with and without NMatrix/NArray.
 
 ## Links
 

Rakefile

@@ -10,7 +10,13 @@ Bundler::GemHelper.install_tasks
 Rake::TestTask.new do |t|
   t.libs << 'test'
   t.libs << 'test/gsl'
-  t.test_files = FileList['test/*.rb', 'test/gsl/*.rb']
+  file_list = [ 'test/*.rb', 'test/gsl/*.rb']  
+  if ENV['NMATRIX']
+    t.libs    << 'test/gsl/nmatrix_tests' 
+    file_list << 'test/gsl/nmatrix_tests/*_test.rb'
+  end
+
+  t.test_files = FileList[*file_list]
 end
 
 spec = eval(IO.read('gsl.gemspec'))

examples/linalg/QR_solve_narray.rb

@@ -1,5 +1,6 @@
 #!/usr/bin/env ruby
 require("gsl")
+require 'narray'
 include GSL
 
 m = NMatrix[[0.18, 0.60, 0.57, 0.96], [0.41, 0.24, 0.99, 0.58],
@@ -10,4 +11,4 @@
 
 b = NVector[1.0, 2, 3, 4]
 
-p Linalg::QR.solve(qr, tau, b)
+p Linalg::QR.solve(qr, tau, b)

examples/linalg/SV_narray.rb

@@ -7,6 +7,16 @@
 
 u, v, s = Linalg::SV.decomp(m)
 
+puts "u ->"
+p u
+
+puts "v ->"
+p v
+
+puts "s ->"
+p s
+
 b = NArray[1.0, 2, 3, 4]
 
+puts "solved ->"
 p Linalg::SV.solve(u, v, s, b)

examples/linalg/chol_narray.rb

@@ -5,11 +5,14 @@
 
 m = NArray[[4.0, 2], [2, 3]]
 c = Cholesky.decomp(m)
+puts "decomp ->"
 p c
 
 b = NArray[1.0, 2]
+puts "solve ->"
 p Cholesky.solve(c, b)    # Expected [-0.125, 0.75]
 
 b = NArray[1.0, 2]
 Cholesky.svx(c, b)
+puts "svx ->"
 p b

examples/wavelet/wavelet1.rb

@@ -5,7 +5,9 @@
 nc = 20
 
 data = GSL::Vector.alloc(n)
-data.fscanf("ecg.dat")
+File.new(File.dirname(__FILE__) + '/ecg.dat').each_with_index do |l, i|
+  data[i] = l.to_f
+end
 
 w = GSL::Wavelet.alloc("daubechies", 4)
 work = GSL::Wavelet::Workspace.alloc(n)
@@ -46,5 +48,5 @@
 #data3 = GSL::Wavelet.transform_inverse(w, data2, work)
 #data3 = GSL::Wavelet.transform_inverse(w, data2)
 
-GSL::graph(nil, data, data3, "-T X -C -g 3 -x 0 #{data.size} -L 'Red: data, Green: DWT'")
+# GSL::graph(nil, data, data3, "-T X -C -g 3 -x 0 #{data.size} -L 'Red: data, Green: DWT'")
 

ext/gsl_native/array.c

@@ -68,6 +68,16 @@ double* get_vector_ptr(VALUE ary, size_t *stride, size_t *n)
     ary2 = na_change_type(ary, NA_DFLOAT);
     return NA_PTR_TYPE(ary2,double*);
 #endif
+
+#ifdef HAVE_NMATRIX_H
+  } else if (NM_IsNMatrix(ary)) {
+    *n = NM_DENSE_COUNT(ary);
+    *stride = 1;
+    if (NM_STORAGE(ary)->dtype != FLOAT64) {
+      rb_raise(rb_eTypeError, "NMatrix must be :float64");
+    }
+    return (double*)NM_DENSE_ELEMENTS(ary);
+#endif
   } else {
     rb_raise(rb_eTypeError,
              "wrong argument type %s", rb_class2name(CLASS_OF(ary)));

ext/gsl_native/bspline.c

@@ -37,12 +37,25 @@ static VALUE rb_gsl_bspline_breakpoint(VALUE obj, VALUE i)
 static VALUE rb_gsl_bspline_knots(VALUE obj, VALUE b)
 {
   gsl_bspline_workspace *w;
+  Data_Get_Struct(obj, gsl_bspline_workspace, w);
+
+#ifdef HAVE_NMATRIX_H
+  if (NM_IsNMatrix(b)) {
+    NM_DENSE_STORAGE *nm_bpts;
+    gsl_vector_view v;
+
+    nm_bpts = NM_STORAGE_DENSE(b);
+    v = gsl_vector_view_array((double*) nm_bpts->elements, NM_DENSE_COUNT(b));
+    gsl_bspline_knots(&v.vector, w);
+    return Data_Wrap_Struct(cgsl_vector_view_ro, 0, NULL, w->knots);    
+  }
+#endif
+
   gsl_vector *bpts;
   CHECK_VECTOR(b);
-  Data_Get_Struct(obj, gsl_bspline_workspace, w);
   Data_Get_Struct(b, gsl_vector, bpts);
   gsl_bspline_knots(bpts, w);
-  return Data_Wrap_Struct(cgsl_vector_view_ro, 0, NULL, w->knots);
+  return Data_Wrap_Struct(cgsl_vector_view_ro, 0, NULL, w->knots);    
 }
 static VALUE rb_gsl_bspline_knots_uniform(int argc, VALUE *argv, VALUE obj)
 {

ext/gsl_native/cheb.c

@@ -115,8 +115,29 @@ static VALUE rb_gsl_cheb_eval(VALUE obj, VALUE xx)
     return ary;
     break;
   default:
+    if (VECTOR_P(xx)) {
+      Data_Get_Struct(xx, gsl_vector, v);
+      vnew = gsl_vector_alloc(v->size);
+      for (i = 0; i < v->size; i++) {
+        gsl_vector_set(vnew, i, gsl_cheb_eval(p, gsl_vector_get(v, i)));
+      }
+      return Data_Wrap_Struct(cgsl_vector, 0, gsl_vector_free, vnew);
+    } 
+#ifdef HAVE_NMATRIX_H
+    else if (NM_IsNMatrix(xx)) {
+      NM_DENSE_STORAGE *nm;
+      double *ptr1, *ptr2;
+      nm = NM_STORAGE_DENSE(xx);
+      ptr1 = (double*) nm->elements;
+      n = NM_DENSE_COUNT(xx);
+      ary = rb_nmatrix_dense_create(FLOAT64, nm->shape, nm->dim, nm->elements, n);
+      ptr2 = (double*)NM_DENSE_ELEMENTS(ary);
+      for (i = 0; i < n; i++) ptr2[i] = gsl_cheb_eval(p, ptr1[i]);
+      return ary;
+    }
+#endif
 #ifdef HAVE_NARRAY_H
-    if (NA_IsNArray(xx)) {
+    else if (NA_IsNArray(xx)) {
       struct NARRAY *na;
       double *ptr1, *ptr2;
       GetNArray(xx, na);
@@ -128,14 +149,7 @@ static VALUE rb_gsl_cheb_eval(VALUE obj, VALUE xx)
       return ary;
     }
 #endif
-    if (VECTOR_P(xx)) {
-      Data_Get_Struct(xx, gsl_vector, v);
-      vnew = gsl_vector_alloc(v->size);
-      for (i = 0; i < v->size; i++) {
-        gsl_vector_set(vnew, i, gsl_cheb_eval(p, gsl_vector_get(v, i)));
-      }
-      return Data_Wrap_Struct(cgsl_vector, 0, gsl_vector_free, vnew);
-    } else if (MATRIX_P(xx)) {
+    else if (MATRIX_P(xx)) {
       Data_Get_Struct(xx, gsl_matrix, m);
       mnew = gsl_matrix_alloc(m->size1, m->size2);
       for (i = 0; i < m->size1; i++) {
@@ -144,7 +158,8 @@ static VALUE rb_gsl_cheb_eval(VALUE obj, VALUE xx)
         }
       }
       return Data_Wrap_Struct(cgsl_matrix, 0, gsl_matrix_free, mnew);
-    } else {
+    }
+    else {
       rb_raise(rb_eTypeError, "wrong argument type");
     }
     break;
@@ -202,6 +217,25 @@ static VALUE rb_gsl_cheb_eval_err(VALUE obj, VALUE xx)
       }
       return rb_ary_new3(2, ary, aerr);
     }
+#endif
+#ifdef HAVE_NMATRIX_H
+    if (NM_IsNMatrix(xx)) {
+      NM_DENSE_STORAGE *nm;
+      double *ptr1, *ptr2, *ptr3;
+      nm = NM_STORAGE_DENSE(xx);
+      n = NM_DENSE_COUNT(xx);
+      ptr1 = (double*) nm->elements;
+      ary = rb_nmatrix_dense_create(FLOAT64, nm->shape, nm->dim, nm->elements, n);
+      aerr = rb_nmatrix_dense_create(FLOAT64, nm->shape, nm->dim, nm->elements, n);
+      ptr2 = (double*)NM_DENSE_ELEMENTS(ary);
+      ptr3 = (double*)NM_DENSE_ELEMENTS(aerr);
+      for (i = 0; i < n; i++) {
+        gsl_cheb_eval_err(p, ptr1[i], &result, &err);
+        ptr2[i] = result;
+        ptr3[i] = err;
+      }
+      return rb_ary_new3(2, ary, aerr);
+    }
 #endif
     if (VECTOR_P(xx)) {
       Data_Get_Struct(xx, gsl_vector, v);
@@ -278,6 +312,19 @@ static VALUE rb_gsl_cheb_eval_n(VALUE obj, VALUE nn, VALUE xx)
       for (i = 0; i < n; i++) ptr2[i] = gsl_cheb_eval_n(p, order, ptr1[i]);
       return ary;
     }
+#endif
+#ifdef HAVE_NMATRIX_H
+    if (NM_IsNMatrix(xx)) {
+      NM_DENSE_STORAGE *nm;
+      double *ptr1, *ptr2;
+      nm = NM_STORAGE_DENSE(xx);
+      n = NM_DENSE_COUNT(xx);
+      ptr1 = (double*) nm->elements;
+      ary = rb_nmatrix_dense_create(FLOAT64, nm->shape, nm->dim, nm->elements, n);
+      ptr2 = (double*)NM_DENSE_ELEMENTS(ary);
+      for (i = 0; i < n; i++) ptr2[i] = gsl_cheb_eval_n(p, order, ptr1[i]);
+      return ary;
+    }
 #endif
     if (VECTOR_P(xx)) {
       Data_Get_Struct(xx, gsl_vector, v);
@@ -356,6 +403,25 @@ static VALUE rb_gsl_cheb_eval_n_err(VALUE obj, VALUE nn, VALUE xx)
       }
       return rb_ary_new3(2, ary, aerr);
     }
+#endif
+#ifdef HAVE_NMATRIX_H
+    if (NM_IsNMatrix(xx)) {
+      NM_DENSE_STORAGE *nm;
+      double *ptr1, *ptr2, *ptr3;
+      nm = NM_STORAGE_DENSE(xx);
+      n = NM_DENSE_COUNT(xx);
+      ptr1 = (double*) nm->elements;
+      ary = rb_nmatrix_dense_create(FLOAT64, nm->shape, nm->dim, nm->elements, n);
+      aerr = rb_nmatrix_dense_create(FLOAT64, nm->shape, nm->dim, nm->elements, n);
+      ptr2 = (double*)NM_DENSE_ELEMENTS(ary);
+      ptr3 = (double*)NM_DENSE_ELEMENTS(aerr);
+      for (i = 0; i < n; i++) {
+        gsl_cheb_eval_n_err(p, order, ptr1[i], &result, &err);
+        ptr2[i] = result;
+        ptr3[i] = err;
+      }
+      return rb_ary_new3(2, ary, aerr);
+    }
 #endif
     if (VECTOR_P(xx)) {
       Data_Get_Struct(xx, gsl_vector, v);