vvector.cc

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

// This may look like C code, but it is really -*- C++ -*-
/*
 ************************************************************************
 *
 *			Verify Operations on Vectors
 *
 * $Id: vvector.cc,v 4.2 1998/12/01 17:28:49 oleg Exp oleg $
 *
 ************************************************************************
 */

#include "LAStreams.h"
#include <math.h>
#include "builtin.h"
#include "iostream.h"
#include <float.h>

/*
 *------------------------------------------------------------------------
 *			Service validation functions
 */
static void verify_vector_identity(const Vector& v1, const Vector& v2)
{ verify_matrix_identity(v1,v2); }

/*
 *------------------------------------------------------------------------
 *	   Test allocation functions and compatibility check
 */

static void test_allocation(void)
{

  cout << "\n\n---> Test allocation and compatibility check\n";

  Vector v1(20);
  Vector v2(1,20);
  Vector v3(0,19);
  Vector v4(v1);
  cout << "The following vector have been allocated\n";
  v1.info(), v2.info(), v3.info(), v4.info();

  cout << "\nStatus information reported for vector v3:\n";
  cout << "  Lower bound ... " << v3.q_lwb() << "\n";
  cout << "  Upper bound ... " << v3.q_upb() << "\n";
  cout << "  No. of elements " << v3.q_no_elems() << "\n";
  cout << "  Name " << v3.q_name() << "\n";
  assert( v3.q_no_elems() == v3.q_nrows() * v3.q_ncols() );

  cout << "\nCheck vectors 1 & 2 for compatibility\n";
  assert( static_cast<const DimSpec&>(v1) == v2 );
  are_compatible(v1,v2);

  cout << "Check vectors 1 & 4 for compatibility\n";
  are_compatible(v1,v4);

  assert( !(static_cast<const DimSpec&>(v1) == v3) );
#if 0
  cout << "Check vectors 1 & 3 for compatibility\n";
  are_compatible(v1,v3);
#endif

  cout << "v2 has to be compatible with v3 after resizing to v3\n";
  v2.resize_to(v3);
  are_compatible(v2,v3);

  Vector v5(v1.q_upb()+5); v5.set_name("Vector v5");
  cout << "v1 has to be compatible with v5 after resizing to v5.upb\n";
  v5.info(); cout << endl;
  v1.resize_to(v5.q_no_elems());
  are_compatible(v1,v5);

  {
    cout << "Check that shrinking does not change remaining elements" << endl;
    Vector vb(-1,20);
    {
      LAStreamOut vbs(vb);
      for(register int i=0; !vbs.eof(); i++)
         vbs.get() = i+M_PI;
    }
    Vector v = vb;
    assert( v == vb );
    assert( v != 0 );
//    ConstMatrixDA va(v);
    v.resize_to(0,10);
    for(register int i=v.q_lwb(); i<=v.q_upb(); i++)
      assert( v(i) == vb(i-v.q_lwb()+vb.q_lwb()) );
    cout << "Check that expansion expands by zeros" << endl;
    const int old_nelems = v.q_upb() - v.q_lwb() + 1;
    v.resize_to(vb);
    assert( !(v == vb) );
    register int i;
    for(i=v.q_lwb(); i<v.q_lwb()+old_nelems; i++)
      assert( v(i) == vb(i) );
    for(; i<=v.q_upb(); i++)
      assert( v(i) == 0 );
  }
  cout << "\nDone\n";
}

/*
 *------------------------------------------------------------------------
 *		     Test uniform element operations
 */

static void test_element_op(const int vsize)
{
  const double pattern = M_PI;
  register int i;

  cout << "\n\n---> Test operations that treat each element uniformly\n";

  Vector v(-1,vsize-2);
  Vector v1(v);

  cout << "Writing zeros to v...\n";
  for(i=v.q_lwb(); i<=v.q_upb(); i++)
    v(i) = 0;
  verify_element_value(v,0);

  cout << "Clearing v1 ...\n";
  v1.clear();
  verify_element_value(v1,0);

  cout << "Comparing v1 with 0 ...\n";
  assure(v1 == 0, "v1 is not zero!");

  cout << "Writing a pattern " << pattern << " by assigning to v(i)...\n";
  for(i=v.q_lwb(); i<=v.q_upb(); i++)
    v(i) = pattern;
  verify_element_value(v,pattern);

  cout << "Writing the pattern by assigning to v1 as a whole ...\n";
  to_every(v1) = pattern;
  verify_element_value(v1,pattern);

  cout << "Comparing v and v1 ...\n";
  assure(v == v1, "v and v1 are unexpectedly different!");
  cout << "Comparing (v=0) and v1 ...\n";
  assert(!(v.clear() == v1));
  assert( of_every(v) == 0 );
  
  cout << "Asigning the pattern via a LAStreamOut..." << endl;
  for(LAStreamOut va(v); !va.eof(); )
    va.get() = pattern;
  assert( v == v1 );
  
  cout << "\nClear v and add the pattern\n";
  v.clear();
  to_every(v) += pattern;
  verify_element_value(v,pattern);
  cout << "   add the doubled pattern with the negative sign\n";
  v += -2*pattern;
  verify_element_value(v,-pattern);
  cout << "   subtract the trippled pattern with the negative sign\n";
  v -= -3*pattern;
  verify_element_value(v,2*pattern);

  cout << "\nVerify comparison operations\n";
  v = pattern;
  assert( of_every(v) == pattern && !(of_every(v) != pattern) &&
  	  of_every(v) >= pattern && of_every(v) <= pattern );
  assert( of_every(v) > 0 && of_every(v) >= 0 );
  assert( of_every(v) > -pattern && of_every(v) >= -pattern );
  assert( of_every(v) < pattern+1 && of_every(v) <= pattern+1 );
  v(v.q_upb()) += 1;
  assert( !(v==pattern) && !(v != pattern) && v != pattern-1 );
  assert( of_every(v) >= pattern && !(of_every(v) > pattern) &&
  	  !(of_every(v) >= pattern+1) );
  assert( of_every(v) <= (REAL)pattern+1 && 
  	  !(of_every(v) < pattern+1) && !(of_every(v) <= pattern) );

  cout << "\nAssign 2*pattern to v by repeating additions\n";
  v = 0; to_every(v) += pattern; v += pattern;
  cout << "Assign 2*pattern to v1 by multiplying by two \n";
  v1 = pattern; v1 *= 2;
  verify_element_value(v1,2*pattern);
  assert( v == v1 );
  cout << "Multiply v1 by one half returning it to the 1*pattern\n";
  v1 *= 1/2.;
  verify_element_value(v1,pattern);

  cout << "\nAssign -pattern to v and v1\n";
  v.clear(); to_every(v) -= pattern; v1 = -pattern;
  verify_element_value(v,-pattern);
  assert( v == v1 );
  cout << "v = sqrt(sqr(v)); v1 = abs(v1); Now v and v1 have to be the same\n";
  to_every(v).sqr();
  verify_element_value(v,pattern*pattern);
  to_every(v).sqrt();
  verify_element_value(v,pattern);
  to_every(v1).abs();
  verify_element_value(v1,pattern);
  assert( v == v1 );

  {
    cout << "\nCheck out to see that sin^2(x) + cos^2(x) = 1\n";
    for(i=v.q_lwb(); i<=v.q_upb(); i++)
      v(i) = 2*M_PI/v.q_no_elems() * i;
    class SinAction : public ElementWiseAction
    {
      void operator () (REAL& element) { element = sin(element); }
      public: SinAction(void) {}
    };
    to_every(v).apply(SinAction());
    Vector v2 = v;
    
    class CosAction : public ElementAction
    {
      double factor;
      void operation(REAL& element, const int i, const int j)
      		{ element = cos(factor*i); }
      public:
      CosAction(const int no_elems): factor(2*M_PI/no_elems) {}
    };
    v1.apply(CosAction(v.q_no_elems()));
    Vector v3 = v1;
    to_every(v).sqr();
    to_every(v1).sqr();
    to_every(v) += of_every(v1);
    verify_element_value(v,1);

    cout << "\n\tdo it again through LazyMatrix promise of a vector" << endl;
    class square_add : public LazyMatrix, public ElementAction
    {
      const Vector &v1; Vector &v2;
      void operation(REAL& element, const int i, const int j)
      		{ assert(j==1); element = v1(i)*v1(i) + v2(i)*v2(i); }
     void fill_in(Matrix& m) const { m.apply(*this); }
     public: square_add(const Vector& _v1, Vector& _v2) :
   	LazyMatrix(_v1.q_row_lwb(),_v1.q_row_upb(),1,1),
   	v1(_v1), v2(_v2) {}
    };
    Vector vres = square_add(v2,v3);
    Vector vres1 = v2; assert( !(vres1 == vres) );
    verify_element_value(vres,1);
    vres1 = square_add(v2,v3);
    verify_element_value(vres1,1);
  } 

  cout << "\nVerify constructor with initialization\n";
  Vector vi(1,5,0.0,1.0,2.0,3.0,4.0,"END");
  Vector vit(5);
  for(i=vit.q_lwb(); i<=vit.q_upb(); i++)
    vit(i) = i-1;
  verify_vector_identity(vi,vit);
  
  cout << "\nDone\n\n";
}

/*
 *------------------------------------------------------------------------
 *			Test binary vector operations
 */

static void test_binary_op(const int vsize)
{
  const double pattern = M_PI;
  register int i;

  cout << "\n---> Test Binary Vector operations\n";

  Vector v(2,vsize+1);
  Vector v1(v);
  Vector vp(v);

  for(i=v.q_lwb(); i<=v.q_upb(); i++)
    vp(i) = (i-v.q_no_elems()/2.)*pattern;
  

  cout << "\nVerify assignment of a vector to the vector\n";
  v = pattern;
  v1.clear();
  v1 = v;
  verify_element_value(v1,pattern);
  assert( v1 == v );
  assert( of_every(v) == of_every(v1) );
  
  cout << "\nAdding one vector to itself, uniform pattern " << pattern << "\n";
  v.clear(); v = pattern;
  assert( of_every(v1) == of_every(v1) );
  
  v1 = v; to_every(v1) += of_every(v1);
  verify_element_value(v1,2*pattern);
  cout << "  subtracting two vectors ...\n";
  to_every(v1) -= of_every(v);
  verify_element_value(v1,pattern);
  cout << "  subtracting the vector from itself\n";
  to_every(v1) -= of_every(v1);
  verify_element_value(v1,0);
  assert( of_every(v) > of_every(v1) );
  assert( of_every(v1) < of_every(v) );
  cout << "  adding two vectors together\n";
  to_every(v1) += of_every(v);
  verify_element_value(v1,pattern);

  cout << "\nArithmetic operations on vectors with not the same elements\n";
  cout << "   adding vp to the zero vector...\n";
  v.clear(); to_every(v) += of_every(vp);
  verify_vector_identity(v,vp);
  assert( of_every(v).max_abs(of_every(vp)) <= FLT_EPSILON );
  
//  assert( v == vp );
#if 0
  v1 = v;
  cout << "   making v = 3*vp and v1 = 3*vp, via add() and successive mult\n";
  add(v,2,vp);
  v1 += v1; v1 += vp;
  verify_vector_identity(v,v1);
  cout << "   clear both v and v1, by subtracting from itself and via add()\n";
  v1 -= v1;
  add(v,-3,vp);