Yale insertion refactored, including slice-setting. Now just need to

enable to spec and fix the bugs.

Author: translunar

.rspec

@@ -1,2 +1,2 @@
 --color
---format nested
+--format RSpec::Longrun::Formatter

ext/nmatrix/storage/yale/class.h

@@ -312,7 +312,7 @@ class YaleStorage {
    */
   void insert(row_iterator i, size_t j, size_t* lengths, D* const v, size_t v_size) {
     // Expensive pre-processing step: find all the information we need in order to do insertions.
-    multi_row_insertion_plan p = multi_row_insertion_plan(i, j, lengths, v, v_size);
+    multi_row_insertion_plan p = insertion_plan(i, j, lengths, v, v_size);
 
     // There are more efficient ways to do this, but this is the low hanging fruit version of the algorithm.
     // Here's the full problem: http://stackoverflow.com/questions/18753375/algorithm-for-merging-short-lists-into-a-long-vector
@@ -331,7 +331,7 @@ class YaleStorage {
       size_t v_offset = 0;
       int accum       = 0;
       for (size_t ii = 0; ii < lengths[0]; ++ii, ++i) {
-        i.insert(row_stored_nd_iterator(&i, p.pos[ii]), j, lengths[1], v, v_size, v_offset);
+        i.insert(row_stored_nd_iterator(i, p.pos[ii]), j, lengths[1], v, v_size, v_offset);
       }
     }
   }
@@ -343,18 +343,17 @@ class YaleStorage {
    * +right+ and determine what other functions to call in order to properly handle
    * it.
    */
-  row_nd_iter_pair insert(SLICE* slice, VALUE right) {
+  void insert(SLICE* slice, VALUE right) {
     if (TYPE(right) == T_DATA) {
       if (RDATA(right)->dfree == (RUBY_DATA_FUNC)nm_delete || RDATA(right)->dfree == (RUBY_DATA_FUNC)nm_delete_ref) {
         rb_raise(rb_eNotImpError, "this type of slicing not yet supported");
       } else {
         rb_raise(rb_eTypeError, "unrecognized type for slice assignment");
       }
-      row_iterator it = riend();
-      return std::make_pair(it, it.ndend());
     } else {
 
       D* v;
+
       size_t v_size = 1;
       bool v_alloc = false;
 
@@ -372,21 +371,15 @@ class YaleStorage {
 
       row_iterator i = ribegin(slice->coords[0]);
 
-
-
       if (slice->single || (slice->lengths[0] == 1 && slice->lengths[1] == 1)) { // single entry
-        row_stored_nd_iterator j = i.insert(slice->coords[1], *v);
-        if (v_alloc) xfree(v);
-        return std::make_pair(i,j);
+        i.insert(slice->coords[1], *v);
       } else if (slice->lengths[0] == 1) { // single row, multiple entries
-        row_stored_nd_iterator j = i.insert(slice->coords[1], slice->lengths[1], v, v_size);
-        if (v_alloc) xfree(v);
-        return std::make_pair(i,j);
+        i.insert(slice->coords[1], slice->lengths[1], v, v_size);
       } else { // multiple rows, unknown number of entries
-        row_nd_iter_pair ij = insert(i, slice->coords[1], slice->lengths, v, v_size);
-        if (v_alloc) xfree(v);
-        return ij;
+        insert(i, slice->coords[1], slice->lengths, v, v_size);
       }
+
+      if (v_alloc) xfree(v);
     }
   }
 
@@ -676,18 +669,15 @@ class YaleStorage {
     for (const_row_iterator it = cribegin(); it != criend(); ++it) {
       for (auto jt = it.begin(); !jt.end(); ++jt) {
         if (it.i() == jt.j()) {
-          std::cerr << "copy(ns): writing to diag pos " << it.i() << std::endl;
           if (Yield)  ns_a[it.i()] = rb_yield(~jt);
           else        ns_a[it.i()] = static_cast<E>(*jt);
         } else if (*jt != const_default_obj()) {
-          std::cerr << "copy(ns): writing to pos " << sz << std::endl;
           if (Yield)  ns_a[sz]     = rb_yield(~jt);
           else        ns_a[sz]     = static_cast<E>(*jt);
           ns.ija[sz]    = jt.j();
           ++sz;
         }
       }
-      std::cerr << "copy(ns): updating row end pointer for row " << it.i() << " to " << sz << std::endl;
       ns.ija[it.i()+1]  = sz;
     }
 
@@ -713,7 +703,7 @@ class YaleStorage {
       size_t ndnz       = count_copy_ndnz();
       size_t reserve    = shape(0) + ndnz + 1;
 
-      std::cerr << "reserve = " << reserve << std::endl;
+//      std::cerr << "reserve = " << reserve << std::endl;
 
       lhs               = YaleStorage<E>::create(xshape, reserve);
 
@@ -904,6 +894,7 @@ class YaleStorage {
   void move_right(row_stored_nd_iterator position, size_t n) {
     size_t sz = size();
     for (size_t m = 0; m < sz - position.p(); ++m) {
+      //std::cerr << "moving from " << sz-1-m << " to " << sz+n-1-m << std::endl;
       ija(sz+n-1-m) = ija(sz-1-m);
       a(sz+n-1-m)   = a(sz-1-m);
     }
@@ -1022,13 +1013,13 @@ class YaleStorage {
     }
 
     // Now update row pointers following the changed row as we copy the additional values.
-    for (size_t m = real_i + 1; m < real_shape(0); ++m) {
+    for (size_t m = real_i + 1; m <= real_shape(0); ++m) {
       new_ija[m]        = ija(m) + n;
       new_a[m]          = a(m);
     }
 
     // Copy all remaining prior to insertion/removal site
-    for (size_t m = real_shape(0); m < position.p(); ++m) {
+    for (size_t m = real_shape(0) + 1; m < position.p(); ++m) {
       new_ija[m]        = ija(m);
       new_a[m]          = a(m);
     }

ext/nmatrix/storage/yale/iterators/row.h

@@ -236,12 +236,6 @@ class row_iterator_T {
 
   inline VALUE rb_i() const { return LONG2NUM(i()); }
 
-  row_stored_nd_iterator_T<D,RefType,YaleRef> ndfind(size_t j) {
-    if (j == 0) return ndbegin();
-    size_t p = y.real_find_left_boundary_pos(p_first, p_last, j + y.offset(1));
-    return row_stored_nd_iterator_T<D,RefType,YaleRef>(*this, p);
-  }
-
   row_stored_iterator_T<D,RefType,YaleRef> begin() {  return row_stored_iterator_T<D,RefType,YaleRef>(*this, p_first);  }
   row_stored_nd_iterator_T<D,RefType,YaleRef> ndbegin() {  return row_stored_nd_iterator_T<D,RefType,YaleRef>(*this, p_first);  }
   row_stored_iterator_T<D,RefType,YaleRef> end() { return row_stored_iterator_T<D,RefType,YaleRef>(*this, p_last+1, true); }
@@ -257,6 +251,21 @@ class row_iterator_T {
     row_stored_nd_iterator_T<D,RefType,YaleRef>(*this, y.real_find_left_boundary_pos(p_first, p_last, y.offset(1)));
   }
 
+  row_stored_nd_iterator_T<D,RefType,YaleRef> ndfind(size_t j) {
+    if (j == 0) return ndbegin();
+    std::cerr << "ndfind: p_first = " << p_first << " " << p_last << std::endl;
+    size_t p = p_first > p_last ? p_first : y.real_find_left_boundary_pos(p_first, p_last, j + y.offset(1));
+    std::cerr << "ndfind(" << j << ")" << " = " << p << "\t(max_p = " << y.ija(y.real_shape(0)) << ")" << std::endl;
+    row_stored_nd_iterator iter = row_stored_nd_iterator_T<D,RefType,YaleRef>(*this, p);
+    std::cerr << "iter.end = " << std::boolalpha << iter.end() << ", p = " << iter.p() << std::endl;
+    return iter;
+  }
+
+  row_stored_iterator_T<D,RefType,YaleRef> find(size_t j) {
+    if (j == 0) return begin(); // may or may not be on the diagonal
+    else return row_stored_iterator_T<D,RefType,YaleRef>(*this, ndfind(j).p(), false); // is on the diagonal, definitely
+  }
+
   /*
    * Remove an entry from an already found non-diagonal position. Adjust this row appropriately so we can continue to
    * use it.
@@ -265,9 +274,10 @@ class row_iterator_T {
   row_stored_nd_iterator erase(row_stored_nd_iterator position) {
     size_t sz = y.size();
     if (y.capacity() / nm::yale_storage::GROWTH_CONSTANT <= sz - 1) {
-      y.update_resize_move(position, i() + offset(0), -1);
+      y.update_resize_move(position, real_i(), -1);
     } else {
       y.move_left(position, 1);
+      y.update_real_row_sizes_from(real_i(), -1);
     }
     adjust_length(-1);
     return row_stored_nd_iterator(*this, position.p()-1);
@@ -287,41 +297,20 @@ class row_iterator_T {
     }
   }
 
-  //template <typename = typename std::enable_if<!std::is_const<RefType>::value>::type>
-  template <typename T = typename std::conditional<std::is_const<RefType>::value,void,row_stored_nd_iterator>::type>
-  row_stored_nd_iterator insert(row_stored_iterator position, size_t jj, const D& val) {
-    if (position.diag()) {
-      *position = val;  // simply replace existing, regardless of whether it's 0 or not
-      ++position;
-      return row_stored_nd_iterator(*this, position.p());
-    } else {
-      row_stored_nd_iterator jt(*this, position.p());
-      return insert(jt, jj, val);
-    }
-  }
 
-  /*
-   * Insert an element in column j, using position's p() as the location to insert the new column. i and j will be the
-   * coordinates. This also does a replace if column j is already present.
-   *
-   * Returns true if a new entry was added and false if an entry was replaced.
-   *
-   * Pre-conditions:
-   *   - position.p() must be between ija(real_i) and ija(real_i+1), inclusive, where real_i = i + offset(0)
-   *   - real_i and real_j must not be equal
-   */
+
   //template <typename = typename std::enable_if<!std::is_const<RefType>::value>::type>
   row_stored_nd_iterator insert(row_stored_nd_iterator position, size_t jj, const D& val) {
     size_t sz = y.size();
+    while (!position.end() && position.j() < jj) ++position; // position is just a hint. (This loop ideally only has to happen once.)
+
     if (!position.end() && position.j() == jj) {
-      std::cerr << "insert: *position = val at " << i_ << "," << jj << "\tp=" << position.p() << std::endl;
       *position = val;      // replace existing
     } else {
+
       if (sz + 1 > y.capacity()) {
-        std::cerr << "insert: update_resize_move " << i_ << "," << jj << "\tp=" << position.p() << std::endl;
         y.update_resize_move(position, real_i(), 1);
       } else {
-        std::cerr << "insert: move_right at " << i_ << "," << jj << "\tp=" << position.p() << std::endl;
         y.move_right(position, 1);
         y.update_real_row_sizes_from(real_i(), 1);
       }
@@ -333,26 +322,44 @@ class row_iterator_T {
     return position++;
   }
 
+
+  /*
+   * This version of insert doesn't return anything. Why, when the others do?
+   *
+   * Well, mainly because j here can be a diagonal entry. Most of the inserters return the *next* element following
+   * the insertion, but to do that, we have to create a row_stored_nd_iterator, which requires at least one binary
+   * search for the location following the diagonal (and as of the writing of this, two binary searches). There's no
+   * reason to do that when we never actually *use* the return value. So instead we just have void.
+   */
   //template <typename = typename std::enable_if<!std::is_const<RefType>::value>::type>
-  row_stored_nd_iterator insert(size_t j, const D& val) {
-    return insert(ndfind(j), j, val);
+  void insert(size_t j, const D& val) {
+    if (j + y.offset(1) == real_i())  a(real_i()) = val;
+    else {
+      row_stored_nd_iterator jt = ndfind(j);
+      if (!jt.end() && jt.j() == j) {
+        if (val == y.const_default_obj()) erase(jt);          // erase
+        else                              insert(jt, j, val); // replace
+      } else { // only insert if it's not the default
+        if (val != y.const_default_obj()) insert(jt, j, val);
+      }
+    }
   }
 
 
   /*
    * Determines a plan for inserting a single row. Returns an integer giving the amount of the row change.
    */
-  int single_row_insertion_plan(row_stored_nd_iterator position, size_t jj, size_t length, D const* v, size_t v_size, const size_t& v_offset) {
+  int single_row_insertion_plan(row_stored_nd_iterator position, size_t jj, size_t length, D const* v, size_t v_size, size_t v_offset) {
     int nd_change;
-    size_t m = v_offset;
-    for (size_t jc = jj; jc < jj + length; ++jc, ++m) {
-      if (m >= v_size) m %= v_size; // reset v position.
+
+    for (size_t jc = jj; jc < jj + length; ++jc, ++v_offset) {
+      if (v_offset >= v_size) v_offset %= v_size; // reset v position.
 
       if (jc + y.offset(1) != real_i()) { // diagonal    -- no nd_change here
         if (position.j() != jc) { // not present -- do we need to add it?
-          if (v[m] != y.const_default_obj()) nd_change++;
+          if (v[v_offset] != y.const_default_obj()) nd_change++;
         } else {  // position.j() == jc
-          if (v[m] == y.const_default_obj()) nd_change--;
+          if (v[v_offset] == y.const_default_obj()) nd_change--;
           ++position; // move iterator forward.
         }
       }
@@ -378,7 +385,7 @@ class row_iterator_T {
    * Insert elements into a single row. Returns an iterator to the end of the insertion range.
    */
   row_stored_nd_iterator insert(row_stored_nd_iterator position, size_t jj, size_t length, D const* v, size_t v_size, size_t& v_offset) {
-    int nd_change = single_row_insertion_plan(position, jj, length, v, v_size);
+    int nd_change = single_row_insertion_plan(position, jj, length, v, v_size, v_offset);
 
     // First record the position, just in case our iterator becomes invalid.
     size_t pp = position.p();

ext/nmatrix/storage/yale/iterators/row_stored.h

@@ -72,6 +72,13 @@ class row_stored_iterator_T : public row_stored_nd_iterator_T<D,RefType,YaleRef,
   {
   }
 
+  /* Diagonal constructor. Puts us on the diagonal (unless end is true) */
+  /*row_stored_iterator_T(RowRef& row, bool end_, size_t j)
+  : row_stored_nd_iterator_T<D,RefType,YaleRef,RowRef>(row.ndfind(j)),
+    d_visited(false),
+    d(!end_ && j + row.offset(1) == row.real_i())
+  { }*/
+
   virtual bool diag() const {
     return d;
   }

ext/nmatrix/storage/yale/iterators/row_stored_nd.h

@@ -75,8 +75,8 @@ class row_stored_nd_iterator_T {
   }
 
   // DO NOT IMPLEMENT THESE FUNCTIONS. They prevent C++ virtual slicing
-  template <typename T> row_stored_nd_iterator_T(T const& rhs);
-  template <typename T> row_stored_nd_iterator_T<D,RefType,YaleRef,RowRef> const& operator=(T const& rhs);
+  //template <typename T> row_stored_nd_iterator_T(T const& rhs);
+  //template <typename T> row_stored_nd_iterator_T<D,RefType,YaleRef,RowRef> const& operator=(T const& rhs);
 
   // Next two functions are to ensure we can still cast between nd iterators.
   row_stored_nd_iterator_T(row_stored_nd_iterator_T<D,RefType,YaleRef,RowRef> const& rhs)
@@ -108,7 +108,7 @@ class row_stored_nd_iterator_T {
   }
 
   virtual size_t j() const {
-    if (end()) throw std::out_of_range("cannot dereference an end pointer");
+    if (end()) throw std::out_of_range("cannot dereference (get j()) for an end pointer");
     return r.ija(p_) - r.offset(1);
   }
 

ext/nmatrix/storage/yale/yale.cpp

@@ -704,38 +704,9 @@ static void set_single_cell(YALE_STORAGE* storage, size_t* coords, DType& v) {
 template <typename DType>
 void set(VALUE left, SLICE* slice, VALUE right) {
   YALE_STORAGE* storage = NM_STORAGE_YALE(left);
+  YaleStorage<DType> y(storage);
 
-  // TODO: Easily modified to accept pass dense storage elements in instead of v (below). Won't work with slices.
-  if (TYPE(right) == T_DATA) {
-    if (RDATA(right)->dfree == (RUBY_DATA_FUNC)nm_delete || RDATA(right)->dfree == (RUBY_DATA_FUNC)nm_delete_ref) {
-      rb_raise(rb_eNotImpError, "this type of slicing not yet supported");
-    } else {
-      rb_raise(rb_eTypeError, "unrecognized type for slice assignment");
-    }
-
-  } else {
-    YaleStorage<DType> s(storage);
-
-    DType* v;
-    size_t v_size = 1;
-    if (TYPE(right) == T_ARRAY) {  // Allow the user to pass in an array
-      v_size = RARRAY_LEN(right);
-      v      = ALLOC_N(DType, v_size);
-      for (size_t m = 0; m < RARRAY_LEN(right); ++m) {
-        rubyval_to_cval(rb_ary_entry(right, m), storage->dtype, &(v[m]));
-      }
-    } else {
-      v = reinterpret_cast<DType*>(rubyobj_to_cval(right, storage->dtype));
-    }
-
-    if (slice->single || (slice->lengths[0] == 1 && slice->lengths[1] == 1)) { // set a single cell
-      set_single_cell<DType>(storage, slice->coords, *v);
-    } else {
-      set_multiple_cells<DType>(storage, slice->coords, slice->lengths, v, v_size);
-    }
-
-    xfree(v);
-  }
+  y.insert(slice, right);
 }
 
 ///////////

nmatrix.gemspec

@@ -50,11 +50,11 @@ EOF
 
   gem.add_dependency 'rdoc', '>=4.0.1'
 
-  gem.add_development_dependency 'rake' #, '~>0.9'
+  gem.add_development_dependency 'rake'
   gem.add_development_dependency 'bundler'
-  gem.add_development_dependency 'rspec', '~>2.9.0'
+  gem.add_development_dependency 'rspec'
+  gem.add_development_dependency 'rspec-longrun'
   gem.add_development_dependency 'pry'
-  gem.add_development_dependency 'guard-rspec', '~>0.7.0'
   gem.add_development_dependency 'rake-compiler', '~>0.8.1'
 end
 

spec/elementwise_spec.rb

@@ -33,6 +33,7 @@
   context "yale" do
     before :each do
       @n = NMatrix.new(:yale, 3, :int64)
+      @n.extend NMatrix::YaleFunctions
       @m = NMatrix.new(:yale, 3, :int64)
       @n[0,0] = 52
       @n[0,2] = 5

spec/enum_spec.rb

@@ -73,11 +73,14 @@
           vv = []
           ii = []
           jj = []
+          require 'pry'
+          binding.pry
           @n.send :__yale_stored_nondiagonal_each_with_indices__ do |v,i,j|
             vv << v
             ii << i
             jj << j
           end
+
           vv.should == [2,3,4,5,  6,8,9,10,  11,12,14,15,  16,17,18,20]
           ii.should == [[0]*4, [1]*4, [2]*4, [4]*4].flatten
           jj.should == [1,2,3,4,  0,2,3,5,   0,1,4,5,      0,2,3,5]