lastreams.h

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

// This may look like C code, but it is really -*- C++ -*-
/*
 ************************************************************************
 *
 *			  Linear Algebra Package
 *			  
 * The following functions support a _sequential_ traversal of a matrix
 * or its parts (a row, a column, a diagonal, or an arbitrary rectangular
 * block). As the README file explains, a random access to matrix elements
 * is often an unnecessary luxury. Many Numerical Math algorithms access
 * matrices by walking them in regular sequential ways. This _serial_
 * access can be implemented far more efficiently than the arbitrary
 * access.
 * 
 * The classes below act as input/ouput "matrix" streams. As with regular
 * streams, be sure to check for eof() before attempting to get() an
 * element.
 * Note, since all the methods of the matrix streams are implemented inline,
 * the whole iteration would be inline, too.
 * Note, that if the stream is applied to the entire Matrix or its
 * rectangular block, it would be traversed in the natural, that is,
 * COLUMN major order.
 *
 * $Id: LAStreams.h,v 1.1 1998/12/06 23:39:12 oleg Exp oleg $
 *
 ************************************************************************
 */

#ifndef __GNUC__
#pragma once
#endif

#ifndef _LAStreams_h
#define _LAStreams_h 1

#if !defined(_LinAlg_h)
#include "LinAlg.h"
#endif

class LAS {			// A "namespace" to define some common constants
public:
 enum seek_dir { beg, cur, end};
 
 static inline const REAL * min(const REAL * const p1, const REAL * const p2)
  	{ return p1 > p2 ? p2 : p1; }
	
 static inline const REAL * max(const REAL * const p1, const REAL * const p2)
  	{ return p1 > p2 ? p1 : p2; }
	
};
          			// A bookmark in the stream
class AREALMark
{
  friend class AREALStreamIn;
  friend class AREALStreamOut;
  friend class AREALStrideStreamIn;
  friend class AREALStrideStreamOut;
  friend class AREALBlockStreamIn;
  friend class AREALBlockStreamOut;
  size_t offset;
  enum { invalid_offset = (size_t)(-1) };
  AREALMark(const size_t _offset) : offset(_offset) {}
public:
  AREALMark(void) : offset(invalid_offset) {}
  operator bool (void) const	{ return offset != invalid_offset; }
  bool operator ! (void) const	{ return offset == invalid_offset; }
  rowcol get_2dpos(const DimSpec dims) const
  	{ assert( (bool)(*this) && dims.q_nrows() > 0 );
  	  div_t res = div(offset,dims.q_nrows());
  	  return rowcol(res.rem+dims.q_row_lwb(),res.quot+dims.q_col_lwb()); }
  friend ostream& operator << (ostream& os, const AREALMark& mark);
};
  
 
 			// The following are the pair of "streams" over
 			// the field of REALs
 			// These classes are mixins - they provide an
 			// "abstract" stream functionality. The classes
 			// don't have public constructors, as none of
 			// mixin classes do. The classes are supposed to be
 			// "blended" with other classes to imbue them with
 			// stream properties: the derived class
 			// must then properly set the current element
 			// pointer and the last element pointer.
 			// Since the classes do not declare any virtual
 			// functions, they are lightweight (and can be
 			// inlined)

class AREALStreamIn
{
  const REAL * curr_el_p;            // curr_el_p <= last_el_p always
  const REAL * const first_el_p;     // The first element pointer
  const REAL * const last_el_p;      // The == condition means EOF
  
  AREALStreamIn(const AREALStreamIn&);	 // Not implemented and forbidden:
  void operator = (const AREALStreamIn&);// no cloning/assignment allowed

  				// Apply range:dir to yield a new first_el_p
  				// invariant: new_first_el_p >= first_el_p
  const REAL * subrange_first(const IRange range, LAS::seek_dir dir) const
  {
    if( dir == LAS::beg )
      return range.lwb <= 0 ? first_el_p : first_el_p + range.lwb;
    else if( dir == LAS::end )
      return range.lwb == - IRange::INF ? first_el_p : 
   		LAS::max(first_el_p,last_el_p - range.lwb);
    else
      return range.lwb == - IRange::INF ? curr_el_p : 
   		LAS::max(first_el_p,curr_el_p + range.lwb);
  }

  				// Apply range:dir to yield a new last_el_p
  				// invariant: new_last_el_p <= last_el_p
  const REAL * subrange_last(const IRange range, LAS::seek_dir dir) const
  {
    if( dir == LAS::beg )
      return range.upb == IRange::INF ? last_el_p : 
   		LAS::min(last_el_p,first_el_p + range.upb + 1);
    else if( dir == LAS::end )
      return last_el_p;
    else
      return range.upb == IRange::INF ? last_el_p : 
   		LAS::min(last_el_p,curr_el_p + range.upb + 1);
  }

protected:
				// A protected constructor
  AREALStreamIn(const REAL * _beg_ptr, const REAL * _end_ptr)
        : curr_el_p(_beg_ptr),
          first_el_p(_beg_ptr),
          last_el_p(_end_ptr)
          {}

	  			// Subrange the current stream by applying
	  			// the range. seek_dir determines if the
	  			// range applies to the current state
	  			// of the stream or its init position
  AREALStreamIn(const AREALStreamIn& proto, const IRange range,
		LAS::seek_dir dir)
        : first_el_p(proto.subrange_first(range,dir)),
          last_el_p(proto.subrange_last(range,dir))
          { curr_el_p = first_el_p;
	    assert(first_el_p <= last_el_p); }

                                // 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() the AREALStream boundary!");
          return *curr_el_p++; }
  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 AREALStream 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 - 1); }
 
 			// Set the current position at the mark
  AREALStreamIn& 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())
  AREALStreamIn& seek(const int offset, LAS::seek_dir dir = LAS::cur)
  	{ switch(dir)
	  { case LAS::beg: curr_el_p = first_el_p + offset; break;
	    case LAS::cur: curr_el_p = curr_el_p + offset; break;
	    case LAS::end: curr_el_p = last_el_p - offset; break;
	    default: assert(0 /*wrong seek_dir*/); };
	  if( curr_el_p > last_el_p ) curr_el_p = last_el_p;
  	  else assert( curr_el_p >= first_el_p );
  	  return *this; }
	  
  AREALStreamIn& 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 AREALStreamOut : public AREALStreamIn
{
  AREALStreamOut(const AREALStreamOut&);  // Not implemented and forbidden:
  void operator = (const AREALStreamOut&);// no cloning/assignment allowed

protected:
				// A protected constructor
  AREALStreamOut(REAL * _beg_ptr, REAL * _end_ptr)
        : AREALStreamIn(_beg_ptr,_end_ptr)
          {}

  AREALStreamOut(const AREALStreamIn& proto, const IRange range,
		LAS::seek_dir dir)
  	: AREALStreamIn(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;
};
 

                                // The following class is an STL-like iterator
                                // to access _values_ of Matrix elements
                                // one-by-one.
                                // Note this is a _read-only_ stream
class LAStreamIn : protected Matrix::ConstReference,
		   public AREALStreamIn
{
  LAStreamIn(const LAStreamIn&);	// Not implemented and forbidden:
  void operator = (const LAStreamIn&);	// no cloning/assignment allowed

public:
  LAStreamIn(const Matrix& m)
        : Matrix::ConstReference(m),
          AREALStreamIn(m.elements,m.elements+m.nelems)
          {}
  inline LAStreamIn(const ConstMatrixColumn& mc);
  LAStreamIn(const LAStreamIn& ls, const IRange range,
	    LAS::seek_dir dir = LAS::cur)
  	: Matrix::ConstReference(ls.ref()),
	  AREALStreamIn(ls,range,dir) {}
  inline LAStreamIn(const LAStreamOut& ls, const IRange range,
		    LAS::seek_dir dir = LAS::cur);
  rowcol get_pos(const AREALMark mark) const
  	  { return mark.get_2dpos(operator const Matrix& ()); }
};

inline LAStreamIn::LAStreamIn(const ConstMatrixColumn& mc)
        : Matrix::ConstReference(mc.ref()),
          AREALStreamIn(mc.col_ptr,mc.col_ptr+mc.nrows)
          {}

                                // This class provides read *and* write
                                // access to matrix elements, one-by-one
class LAStreamOut : protected Matrix::Reference,
		    public AREALStreamOut
{
  friend class LAStreamIn;
  LAStreamOut(const LAStreamOut&);	// Not implemented and forbidden:
  void operator = (const LAStreamOut&);	// no cloning/assignment allowed
public:
  LAStreamOut(Matrix& m)
        : Matrix::Reference(m),
          AREALStreamOut(m.elements,m.elements+m.nelems)
          {}
  inline LAStreamOut(const MatrixColumn& mc);
  LAStreamOut(const LAStreamOut& ls, const IRange range,
	    LAS::seek_dir dir = LAS::cur)
  	: Matrix::Reference(ls.ref()),
	  AREALStreamOut(ls,range,dir) {}
  rowcol get_pos(const AREALMark mark)
  	  { return mark.get_2dpos(operator Matrix& ()); }
};

inline LAStreamOut::LAStreamOut(const MatrixColumn& mc)
        : Matrix::Reference(const_cast<Matrix&>(mc.ref())),
          AREALStreamOut(const_cast<REAL*>(mc.col_ptr),
          		 const_cast<REAL*>(mc.col_ptr)+mc.q_nrows())
          {}

inline LAStreamIn::LAStreamIn(const LAStreamOut& ls, const IRange range,
		    LAS::seek_dir dir)
  	: Matrix::ConstReference(ls.ref()),
	  AREALStreamIn(ls,range,dir) {}

 
 			// The following pair of streams are similar
 			// to the ones above, with an exception they
 			// walk a sequence of REALs with a non-unit
 			// stride
class AREALStrideStreamIn
{
  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 int stride;

  AREALStrideStreamIn(const AREALStrideStreamIn&); // Not implemented and forbidden:
  void operator = (const AREALStrideStreamIn&);// no cloning/assignment allowed

  				// Apply range:dir to yield a new first_el_p
  				// invariant: new_first_el_p >= first_el_p
  const REAL * subrange_first(const IRange range, LAS::seek_dir dir) const
  {
    if( dir == LAS::beg )
      return range.lwb <= 0 ? first_el_p : first_el_p + stride*range.lwb;
    else if( dir == LAS::end )
     if( range.lwb == - IRange::INF )
       return first_el_p;
     else
     { const int r = (last_el_p - first_el_p) % stride;
       return LAS::max(first_el_p,
		       (r == 0 ? last_el_p : last_el_p - r + stride)
	    				- range.lwb*stride); }
    else
      return range.lwb == - IRange::INF ? curr_el_p : 
   		LAS::max(first_el_p,curr_el_p + stride*range.lwb);
  }

  				// Apply range:dir to yield a new last_el_p
  				// invariant: new_last_el_p <= last_el_p
  const REAL * subrange_last(const IRange range, LAS::seek_dir dir) const
  {
    if( dir == LAS::beg )
      return range.upb == IRange::INF ? last_el_p : 
   		LAS::min(last_el_p,first_el_p + stride*(range.upb+1));
    else if( dir == LAS::end )
      return last_el_p;
    else
      return range.upb == IRange::INF ? last_el_p : 
   		LAS::min(last_el_p,curr_el_p + stride*(range.upb+1));
  }
 
protected:
				// A protected constructor
  AREALStrideStreamIn(const REAL * _beg_ptr, const REAL * _end_ptr,
  		      const int _stride)
        : curr_el_p(_beg_ptr),
          first_el_p(_beg_ptr),
          last_el_p(_end_ptr),
          stride(_stride)
          {}

	  			// Subrange the current stream by applying
	  			// the range. seek_dir determines if the
	  			// range applies to the current state
	  			// of the stream or its init position
  AREALStrideStreamIn(const AREALStrideStreamIn& proto, const IRange range,
		LAS::seek_dir dir)
        : first_el_p(proto.subrange_first(range,dir)),
          last_el_p(proto.subrange_last(range,dir)),
	  stride(proto.stride)
          { curr_el_p = first_el_p;
	    assert(first_el_p <= last_el_p); }

                                // Get the current elem _ref_ and advance