lastreams.h

basic linear algebra classes and applications (SVD,interpolation, multivariate optimization)

                                // in the stream
  const REAL& get_ref(void)
        { if( curr_el_p >= last_el_p)
           _error("Can't get() past the AREALStrideStream boundary!");
          const REAL& curr = *curr_el_p;
          return curr_el_p += stride, curr; }
  const REAL& peek_ref(void) const      // Does *not* advance the "stream ptr"
        { if( curr_el_p >= last_el_p)
           _error("Can't peek() past the AREALStrideStream boundary!");
          return *curr_el_p; }

public:
			// No public constructors!

  bool eof(void) const          { return curr_el_p >= last_el_p; }
  bool bof(void) const          { return curr_el_p == first_el_p; }

                                // Get the current elem and advance in the
                                // stream
  REAL get(void)		{ return get_ref(); }
  REAL peek(void) const         // Does *not* advance the "stream ptr"
  				{ return peek_ref(); }
          
          			// Reset the stream at the beginning
  void rewind(void)		{ curr_el_p = first_el_p; }
 
  			// Note the current position in the stream
  AREALMark tell(void) const
  	{ return AREALMark(curr_el_p - first_el_p); }
  
  			// Note the previous (that is, just gotten)
  			// position in the stream
  AREALMark tell_prev(void) const
  	{ if( curr_el_p == first_el_p )
  	    _error("The stream is at its beginng, no previous pos exists");
  	  return AREALMark(curr_el_p - first_el_p - stride); }
 
 			// Set the current position at the mark
  AREALStrideStreamIn& seek(const AREALMark mark)
  	{ assert( (bool) mark );
  	  curr_el_p = first_el_p + mark.offset;
  	  assert( curr_el_p >= first_el_p && curr_el_p < last_el_p );
  	  return *this; }
  
 			// Set the current position according to the
	  		// given offset and seek_dir.
	  		// Note: the offset can be arbitrarily large:
	  		// if it points beyond this stream, no run-time
	   		// error is generated but the EOF condition is set
	  		// (which can later be checked with eof())
  AREALStrideStreamIn& seek(const int offset, LAS::seek_dir dir = LAS::cur)
  	{ switch(dir)
	  { case LAS::beg: curr_el_p = first_el_p + offset*stride; break;
	    case LAS::cur: curr_el_p = curr_el_p + offset*stride; break;
	    case LAS::end: { const int r = (last_el_p - first_el_p) % stride;
	    		curr_el_p = (r == 0 ? last_el_p : last_el_p - r + stride)
	    				- offset*stride; } break;
	    default: assert(0 /*wrong seek_dir*/); };
	  assert( curr_el_p >= first_el_p );
  	  return *this; }
	  
  AREALStrideStreamIn& ignore(const int how_many)
  	{ return seek(how_many,LAS::cur); }
	
  	  		// Dump the current status of the stream
  ostream& dump(ostream& os) const;
};

			// This is a readable _and_ writable stream
class AREALStrideStreamOut : public AREALStrideStreamIn
{
  AREALStrideStreamOut(const AREALStrideStreamOut&);  // Not implemented and forbidden:
  void operator = (const AREALStrideStreamOut&);// no cloning/assignment allowed

protected:
				// A protected constructor
  AREALStrideStreamOut(REAL * _beg_ptr, REAL * _end_ptr,
  		       const int _stride)
        : AREALStrideStreamIn(_beg_ptr,_end_ptr,_stride)
          {}


  AREALStrideStreamOut(const AREALStrideStreamIn& proto, const IRange range,
		LAS::seek_dir dir)
  	: AREALStrideStreamIn(proto,range,dir) {}

public:
			// No public constructors!

                                // Get the current elem and advance in the
                                // stream
  REAL& get(void)		{ return const_cast<REAL&>(get_ref()); }
  REAL& peek(void) const         // Does *not* advance the "stream ptr"
  				{ return const_cast<REAL&>(peek_ref()); }
          
  	  		// Dump the current status of the stream
//  ostream& dump(ostream& os) const;
};
 


                                // LAStreams with a stride
class LAStrideStreamIn : protected Matrix::ConstReference,
		   public AREALStrideStreamIn
{
  LAStrideStreamIn(const LAStrideStreamIn&);	// Not implemented and forbidden:
  void operator = (const LAStrideStreamIn&);	// no cloning/assignment allowed

public:
  LAStrideStreamIn(const Matrix& m, const int stride)
        : Matrix::ConstReference(m),
          AREALStrideStreamIn(m.elements,m.elements+m.nelems,stride)
          {}
  LAStrideStreamIn(const LAStrideStreamIn& ls, const IRange range,
	    LAS::seek_dir dir = LAS::cur)
  	: Matrix::ConstReference(ls.ref()),
	  AREALStrideStreamIn(ls,range,dir) {}
 inline LAStrideStreamIn(const LAStrideStreamOut& ls, const IRange range,
		    LAS::seek_dir dir = LAS::cur);
 inline LAStrideStreamIn(const ConstMatrixRow& mr);
 inline LAStrideStreamIn(const ConstMatrixDiag& md);
 rowcol get_pos(const AREALMark mark) const
  	  { return mark.get_2dpos(operator const Matrix& ()); }
};

inline
LAStrideStreamIn::LAStrideStreamIn(const ConstMatrixRow& mr)
  	: Matrix::ConstReference(mr.ref()),
  	  AREALStrideStreamIn(const_cast<REAL*>(mr.row_ptr),
  	  		      const_cast<REAL*>(mr.end_ptr),mr.stride) 
  	{ }
inline
LAStrideStreamIn::LAStrideStreamIn(const ConstMatrixDiag& md)
  	: Matrix::ConstReference(md.ref()),
  	  AREALStrideStreamIn(const_cast<REAL*>(md.start_ptr),
  	  		      const_cast<REAL*>(md.end_ptr),md.stride) 
  	{ }



class LAStrideStreamOut : protected Matrix::Reference,
		    public AREALStrideStreamOut
{
  friend class LAStrideStreamIn;
  LAStrideStreamOut(const LAStrideStreamOut&);	// Not implemented and forbidden:
  void operator = (const LAStrideStreamOut&);	// no cloning/assignment allowed
public:
  LAStrideStreamOut(Matrix& m, const int stride)
        : Matrix::Reference(m),
          AREALStrideStreamOut(m.elements,m.elements+m.nelems,stride)
          {}
  LAStrideStreamOut(const LAStrideStreamOut& ls, const IRange range,
	    LAS::seek_dir dir = LAS::cur)
  	: Matrix::Reference(ls.ref()),
	  AREALStrideStreamOut(ls,range,dir) {}
  inline LAStrideStreamOut(const MatrixRow& mr);
  inline LAStrideStreamOut(const MatrixDiag& md);
  rowcol get_pos(const AREALMark mark)
  	  { return mark.get_2dpos(operator Matrix& ()); }
};

inline LAStrideStreamOut::LAStrideStreamOut(const MatrixRow& mr)
        : Matrix::Reference(const_cast<Matrix&>(mr.ref())),
          AREALStrideStreamOut(const_cast<REAL*>(mr.row_ptr),
          		      const_cast<REAL*>(mr.end_ptr),mr.stride)
          {}
inline LAStrideStreamOut::LAStrideStreamOut(const MatrixDiag& md)
        : Matrix::Reference(const_cast<Matrix&>(md.ref())),
          AREALStrideStreamOut(const_cast<REAL*>(md.start_ptr),
          		       const_cast<REAL*>(md.end_ptr),md.stride)
          {}

inline LAStrideStreamIn::LAStrideStreamIn(const LAStrideStreamOut& ls, 
					  const IRange range,
					 LAS::seek_dir dir)
  	: Matrix::ConstReference(ls.ref()),
	  AREALStrideStreamIn(ls,range,dir) {}





 			// The following pair of streams are similar
 			// to the ones above, with an exception they
 			// sequentially walk a *block* of REALs, in
 			// a column major order
class AREALBlockStreamIn
{
  const REAL * curr_el_p;
  const REAL * const first_el_p;     	// The first element pointer
  const REAL * const last_el_p;      	// The curr_el_p >= condition means EOF
  const REAL * last_col_el_p;		// The end-boundary of the curr column
  const int col_size;		     	// The size of the column in a block
  const int eoc_jump;			// The amount to forward-feed the
  				// curr_el_p to point at the beginning of the
  				// next col once we have moved past the curr col
  
  AREALBlockStreamIn(const AREALBlockStreamIn&); // Not implemented and forbidden:
  void operator = (const AREALBlockStreamIn&);// no cloning/assignment allowed

  				// Set curr_el_p to p and compute the
  				// corresponding column boundary
  void set_curr_col(const REAL * const p)
  {
    curr_el_p = p;
    last_col_el_p = curr_el_p
 	- ((curr_el_p - first_el_p) % (col_size + eoc_jump))
	+ col_size;
  }
protected:
				// A protected constructor
  AREALBlockStreamIn(const REAL * _beg_ptr, const REAL * _end_ptr,
  		      const int _col_size, const int _eoc_jump)
        : curr_el_p(_beg_ptr),
          first_el_p(_beg_ptr),
          last_el_p(_end_ptr),
	  last_col_el_p(_beg_ptr+_col_size),
          col_size(_col_size),
	  eoc_jump(_eoc_jump)
          {}

	  
				// Another protected constructor
  AREALBlockStreamIn(const REAL * orig_ptr, const DimSpecSubranged& clip)
        : first_el_p(orig_ptr + clip.q_min_offset()),
          last_el_p(orig_ptr + clip.q_max_offset()+1),
          col_size(clip.q_nrows()),
	  eoc_jump(clip.q_row_diff())
          { last_col_el_p = (curr_el_p = first_el_p)  + col_size;
	    assert(eoc_jump>=0); }


                                // Get the current elem _ref_ and advance
                                // in the stream
  const REAL& get_ref(void)
        { if( curr_el_p >= last_el_p)
           _error("Can't get() past the AREALBlockStreamIn boundary!");
          const REAL& curr = *curr_el_p++;
	  if( curr_el_p >= last_col_el_p )
	    last_col_el_p = (curr_el_p += eoc_jump) + col_size;
	  return curr; }
  const REAL& peek_ref(void) const      // Does *not* advance the "stream ptr"
        { if( curr_el_p >= last_el_p)
           _error("Can't peek() past the AREALBlockStreamIn boundary!");
          return *curr_el_p; }
public:
			// No public constructors!

  bool eof(void) const          { return curr_el_p >= last_el_p; }
  bool bof(void) const          { return curr_el_p == first_el_p; }

                                // Get the current elem and advance in the
                                // stream
  REAL get(void)		{ return get_ref(); }
  REAL peek(void) const         // Does *not* advance the "stream ptr"
  				{ return peek_ref(); }
          
          			// Reset the stream at the beginning
  void rewind(void)
  	{ last_col_el_p = (curr_el_p = first_el_p) + col_size; }
 
  			// Note the current position in the stream
  AREALMark tell(void) const
  	{ return AREALMark(curr_el_p - first_el_p); }
  
  			// Note the previous (that is, just gotten)
  			// position in the stream
  AREALMark tell_prev(void) const
  	{ if( curr_el_p == first_el_p )
  	    _error("The stream is at its beginng, no previous pos exists");
  	  return AREALMark( curr_el_p - first_el_p -
		( curr_el_p + col_size > last_col_el_p ? 1
			 : 1 + eoc_jump)); }
 
 			// Set the current position at the mark
  AREALBlockStreamIn& seek(const AREALMark mark)
  	{ assert( (bool) mark );
  	  set_curr_col(first_el_p + mark.offset);
  	  assert( curr_el_p >= first_el_p && curr_el_p < last_el_p );
  	  return *this; }
  
 			// Set the current position according to the
	  		// given offset and seek_dir.
	  		// Note: the offset can be arbitrarily large:
	  		// if it points beyond this stream, no run-time
	   		// error is generated but the EOF condition is set
	  		// (which can later be checked with eof())
  AREALBlockStreamIn& seek(const int offset, LAS::seek_dir dir = LAS::cur)
  	{ div_t res = div(offset,col_size);
	  switch(dir)
	  { case LAS::beg: curr_el_p = first_el_p + res.quot*(col_size+eoc_jump);
			   last_col_el_p = curr_el_p + col_size;
			   curr_el_p += res.rem; break;
	    case LAS::cur: curr_el_p += res.rem;
			   if( curr_el_p >= last_col_el_p )
			   	curr_el_p += eoc_jump, last_col_el_p += col_size+eoc_jump;
			   curr_el_p += res.quot*(col_size+eoc_jump);
			   last_col_el_p += res.quot*(col_size+eoc_jump);
			   break;
	    case LAS::end: last_col_el_p = last_el_p - res.quot*(col_size+eoc_jump);
			   curr_el_p = last_col_el_p - res.rem; break;
	    default: assert(0 /*wrong seek_dir*/); };
	  assert( curr_el_p >= first_el_p );
  	  return *this; }
	  
  AREALBlockStreamIn& ignore(const int how_many)
  	{ return seek(how_many,LAS::cur); }
	
  	  		// Dump the current status of the stream
  ostream& dump(ostream& os) const;
};

			// This is a readable _and_ writable stream
class AREALBlockStreamOut : public AREALBlockStreamIn
{
  AREALBlockStreamOut(const AREALBlockStreamOut&);  // Not implemented and forbidden:
  void operator = (const AREALBlockStreamOut&);// no cloning/assignment allowed

protected:
				// A protected constructor
  AREALBlockStreamOut(REAL * _beg_ptr, REAL * _end_ptr,
  		       const int _col_size, const int _eoc_jump)
        : AREALBlockStreamIn(_beg_ptr,_end_ptr,_col_size,_eoc_jump)
          {}
  AREALBlockStreamOut(REAL * orig_ptr, const DimSpecSubranged& clip)
        : AREALBlockStreamIn(orig_ptr,clip)
          {}

public:
			// No public constructors!

                                // Get the current elem and advance in the
                                // stream
  REAL& get(void)		{ return const_cast<REAL&>(get_ref()); }
  REAL& peek(void) const         // Does *not* advance the "stream ptr"
  				{ return const_cast<REAL&>(peek_ref()); }
          
};
 


                                // LAStreams over a sub-block of a matrix
class LABlockStreamIn : protected Matrix::ConstReference,
		   public AREALBlockStreamIn
{
  LABlockStreamIn(const LABlockStreamIn&);	// Not implemented and forbidden:
  void operator = (const LABlockStreamIn&);	// no cloning/assignment allowed

public:
  LABlockStreamIn(const Matrix& m, const IRange row_range,
		  const IRange col_range)
        : Matrix::ConstReference(m),
          AREALBlockStreamIn(m.elements,DimSpecSubranged(m,row_range,col_range))
          {}
  rowcol get_pos(const AREALMark mark) const
  	  { return mark.get_2dpos(operator const Matrix& ()); }
};

class LABlockStreamOut : protected Matrix::Reference,
		    public AREALBlockStreamOut
{
  LABlockStreamOut(const LABlockStreamOut&);	// Not implemented and forbidden:
  void operator = (const LABlockStreamOut&);	// no cloning/assignment allowed
public:
  LABlockStreamOut(Matrix& m, const IRange row_range,
		  const IRange col_range)
        : Matrix::Reference(m),
          AREALBlockStreamOut(m.elements,DimSpecSubranged(m,row_range,col_range))
          {}
  rowcol get_pos(const AREALMark mark)
  	  { return mark.get_2dpos(operator Matrix& ()); }
};

#endif