vvector.cc

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

  verify_vector_identity(v,v1);
#endif

  cout << "\nTesting element-by-element multiplications and divisions\n";
  cout << "   squaring each element with sqr() and via multiplication\n";
  v = vp; v1 = vp;
  to_every(v).sqr();
  to_every(v1) *= of_every(v1);
  verify_vector_identity(v,v1);
  cout << "   compare (v = pattern^2)/pattern with pattern\n";
  v = pattern; v1 = pattern;
  to_every(v).sqr();
  to_every(v) /= of_every(v1);
  verify_element_value(v,pattern);
  compare(v1,v,"Original vector and vector after squaring and dividing");

  cout << "\nDone\n";
}


/*
 *------------------------------------------------------------------------
 *			Verify the norm calculation
 */

static void test_norms(const int vsize)
{
  cout << "\n---> Verify norm calculations\n";

  const double pattern = 10.25;

  if( vsize % 2 == 1 )
    _error("Sorry, size of the vector to test must be even for this test\n");

  Vector v(vsize);
  Vector v1(v);

  cout << "\nAssign " << pattern << " to all the elements and check norms\n";
  v = pattern;
  cout << "  1. norm should be pattern*no_elems\n";
  assert( v.norm_1() == pattern*v.q_no_elems() );
  assert( of_every(v).sum_abs() ==  v.norm_1() );
  cout << "  Square of the 2. norm has got to be pattern^2 * no_elems\n";
  assert( v.norm_2_sqr() == sqr(pattern)*v.q_no_elems() );
  assert( of_every(v).sum_squares() ==  v.norm_2_sqr() );
  cout << "  Inf norm should be pattern itself\n";
  assert( v.norm_inf() == pattern );
  assert( of_every(v).max_abs() == pattern );
  cout << "  Scalar product of vector by itself is the sqr(2. vector norm)\n";
  assert( v.norm_2_sqr() == v * v );

  const double ap_step = 1;
  const double ap_a0   = -pattern;
  const int n = v.q_no_elems();
  cout << "\nAssign the arithm progression with 1. term " << ap_a0 <<
          "\nand the difference " << ap_step << "\n";
  { LAStreamOut va(v);
    for(double val=ap_a0; !va.eof(); val += ap_step)
     va.get() = val;
  }
  int l = min(max((int)ceil(-ap_a0/ap_step),0),n);
  double norm = (2*ap_a0 + (l+n-1)*ap_step)/2*(n-l) +
    (-2*ap_a0-(l-1)*ap_step)/2*l;
  cout << "  1. norm should be " << norm << "\n";
  assert( v.norm_1() == norm );
  assert( of_every(v).sum_abs() ==  norm );
  norm = n*( sqr(ap_a0) + ap_a0*ap_step*(n-1) + sqr(ap_step)*(n-1)*(2*n-1)/6);
  cout << "  Square of the 2. norm has got to be "
          "n*[ a0^2 + a0*q*(n-1) + q^2/6*(n-1)*(2n-1) ], or " << norm << "\n";
  assert( v.norm_2_sqr() == norm );
  assert( of_every(v).sum_squares() ==  norm );
  norm = max(abs(v(v.q_lwb())),abs(v(v.q_upb())));
  cout << "  Inf norm should be max(abs(a0),abs(a0+(n-1)*q)), ie " << norm <<
          "\n";
  assert( v.norm_inf() == norm );
  assert( of_every(v).max_abs() == norm );
  cout << "  Scalar product of vector by itself is the sqr(2. vector norm)\n";
  assert( v.norm_2_sqr() == v * v );

  v1.clear();
  compare(v,v1,"Compare the vector v with a zero vector");

  cout << "\nConstruct v1 to be orthogonal to v as v(n), -v(n-1), v(n-2)...\n";
  for(register int i=0; i<v1.q_no_elems(); i++)
    v1(i+v1.q_lwb()) = v(v.q_upb()-i) * ( i % 2 == 1 ? -1 : 1 );
  cout << "||v1|| has got to be equal ||v|| regardless of the norm def\n";
  assert( v1.norm_1() == v.norm_1() );
  assert( of_every(v).max_abs(of_every(v1)) > of_every(v).max_abs() );
  assert( v1.norm_2_sqr() == v.norm_2_sqr() );
  assert( v1.norm_inf() == v.norm_inf() );
  cout << "But the scalar product has to be zero\n";
  assert( v1 * v == 0 );

  cout << "\nDone\n";
}

/*
 *------------------------------------------------------------------------
 *	    Test operations with vectors and matrix slices
 */

static void test_matrix_slices(const int vsize)
{
  const REAL pattern = 8.625;

  cout << "\n\n---> Test operations with vectors and matrix slices\n";

  Vector vc(0,vsize);
  Vector vr(0,vsize+1);
  Matrix m(0,vsize,0,vsize+1);

  cout << "\nCheck modifying the matrix column-by-column\n";
  m = pattern;
  assert( m == pattern );
  for(register int i=m.q_col_lwb(); i <= m.q_col_upb(); i++)
  {
    to_every(MatrixColumn(m,i)) = pattern-1;
    assert( !( m == pattern ) && !( m != pattern ) );
    to_every(MatrixColumn(m,i)) *= 2;
    vc = of_every(ConstMatrixColumn(m,i));
    verify_element_value(vc,2*(pattern-1));
    assert( of_every(ConstMatrixColumn(m,i)).sum_abs() ==
    	2*(pattern-1)*m.q_nrows() );
    vc = of_every(ConstMatrixColumn(m, 
    	i+1 > m.q_col_upb() ? m.q_col_lwb() : i+1));
    verify_element_value(vc,pattern);
    to_every(MatrixColumn(m,i)) *= 0.5;
    assert( of_every(ConstMatrixColumn(m,i)) == pattern - 1 );
    assert( of_every(ConstMatrixColumn(m,i)).sum_abs(of_every(vc)) ==
    	m.q_nrows() );
    to_every(MatrixColumn(m,i)) += 1;
    assert( m == pattern );
  }

  assert( m == pattern );
  for(register int i=m.q_col_lwb(); i <= m.q_col_upb(); i++)
  {
    to_every(vc) = pattern+1;
    to_every(MatrixColumn(m,i)) = of_every(vc);
    assert( !( m == pattern ) && !( m != pattern ) );
    {
      MatrixColumn mc(m,i);
      for(register int j=m.q_row_lwb(); j<=m.q_row_upb(); j++)
        mc(j) *= 2;
      LAStreamIn m_in(mc);
      for(register int j=0; j<mc.q_nrows(); j++)
        assert( m_in.get() == 2*(pattern+1) );
      assert( m_in.eof() );
      m_in.rewind();
      for(LAStreamOut m_out(mc); !m_out.eof(); )
      {
      	assert( m_out.peek() == m_in.get() );
        m_out.get() *= 2;
      }
      assert( m_in.eof() );
    }
    vc = of_every(ConstMatrixColumn(m,i));
    verify_element_value(vc,4*(pattern+1));
    assert( of_every(ConstMatrixColumn(m,i)) == of_every(vc) );
    assert( of_every(ConstMatrixColumn(m,i)) == 4*(pattern+1) );
    assert( of_every(ConstMatrixColumn(m,i)) >= 4*(pattern+1) );
    assert( of_every(ConstMatrixColumn(m,i)) < 4*(pattern+1)+1 );
    vc = 0.25;
    to_every(MatrixColumn(m,i)) *= of_every(vc);
    to_every(MatrixColumn(m,i)) -= 1;
    vc = of_every(MatrixColumn(m,i));
    verify_element_value(vc,pattern);
    assert( m == pattern );
  }

  cout << "\nCheck modifying the matrix row-by-row\n";
  m = pattern;
  assert( m == pattern );
  for(register int i=m.q_row_lwb(); i <= m.q_row_upb(); i++)
  {
    to_every(MatrixRow(m,i)) = pattern+2;
    assert( !( m == pattern ) && !( m != pattern ) );
    vr = of_every(MatrixRow(m,i));
    verify_element_value(vr,pattern+2);
    assert( of_every(ConstMatrixRow(m,i)).sum_squares() ==
    	sqr(pattern+2)*m.q_ncols() );
    vr = of_every(ConstMatrixRow(m,i+1 > m.q_row_upb() ? m.q_row_lwb() : i+1));
    verify_element_value(vr,pattern);
    assert( of_every(ConstMatrixRow(m,i)).max_abs(of_every(vr)) == 2.0 );
    to_every(MatrixRow(m,i)) += -2;
    assert( m == pattern );
  }

  assert( m == pattern );
  cout << "\nCheck modifying the matrix row-by-row, again\n";
  for(register int i=m.q_row_lwb(); i <= m.q_row_upb(); i++)
  {
    vr = pattern-2;
    to_every(MatrixRow(m,i)) = of_every(vr);
    assert( !( m == pattern ) && !( m != pattern ) );
    {
      MatrixRow mr(m,i);
      for(register int j=m.q_col_lwb(); j<=m.q_col_upb(); j++)
        mr(j) *= 4;
      LAStrideStreamIn m_in(mr);
      for(register int j=0; j<mr.q_ncols(); j++)
        assert( m_in.get() == 4*(pattern-2) );
      assert( m_in.eof() );
      m_in.rewind();
      for(LAStrideStreamOut m_out(mr); !m_out.eof(); )
      {
      	assert( m_out.peek() == m_in.get() );
        m_out.get() *= 2;
      }
      assert( m_in.eof() );
    }
    vr = of_every(ConstMatrixRow(m,i));
    verify_element_value(vr,8*(pattern-2));
    vr = -8;
    to_every(MatrixRow(m,i)) /= of_every(vr);
    vr = 2;
    to_every(MatrixRow(m,i)) -= of_every(vr);
    assert( of_every(ConstMatrixRow(m,i)) < 0 );
    assert( of_every(ConstMatrixRow(m,i)) == -pattern );
    to_every(MatrixRow(m,i)).abs();
    vr = of_every(MatrixRow(m,i));
    verify_element_value(vr,pattern);
    assert( m == pattern );
  }

  cout << "\nCheck modifying the matrix diagonal\n";
  m = pattern;
  to_every(MatrixDiag(m)) = pattern-3;
  assert( !( m == pattern ) && !( m != pattern ) );
  vc = of_every(ConstMatrixDiag(m));
  verify_element_value(vc,pattern-3);
  assert( of_every(ConstMatrixDiag(m)) == pattern-3 );
  to_every(MatrixDiag(m)) += 3;
  assert( m == pattern );
  vc = pattern+3;
  to_every(MatrixDiag(m)) = of_every(vc);
  assert( !( m == pattern ) && !( m != pattern ) );
  {
    MatrixDiag md(m);
    for(register int j=1; j<=md.q_nrows(); j++)
      md(j) /= 3;
    LAStrideStreamIn m_in(md);
    AREALMark mark = m_in.tell();
    for(register int j=0; j<md.q_ncols(); j++)
       assert( m_in.get() == (pattern+3)/3 );
    assert( m_in.eof() );
    m_in.seek(mark);
    for(LAStrideStreamOut m_out(md); !m_out.eof(); )
    {
      assert( m_out.peek() == m_in.get() );
      m_out.get() *= 2;
    }
    assert( m_in.eof() );
  }
  vc = of_every(ConstMatrixDiag(m));
  verify_element_value(vc,(pattern+3)/1.5);
  to_every(MatrixDiag(m)) *= 1.5;
  to_every(MatrixDiag(m)) -= 3;
  vc = of_every(ConstMatrixDiag(m));
  verify_element_value(vc,pattern);
  assert( m == pattern );

  cout << "\nCheck out to see that multiplying by diagonal is a column-wise"
          "\nmatrix multiplication\n";
  Matrix mm(m);
  Matrix m1(m.q_row_lwb(),::max(m.q_row_upb(),m.q_col_upb()),
	    m.q_col_lwb(),::max(m.q_row_upb(),m.q_col_upb()));
  Vector vc1(vc), vc2(vc);
  for(register int i=m.q_row_lwb(); i<m.q_row_upb(); i++)
    to_every(MatrixRow(m,i)) = pattern+i;// Make a multiplicand
  mm = m;				// Save it

  m1 = pattern+10;
  for(register int i=vr.q_lwb(); i<=vr.q_upb(); i++)
    vr(i) = i+2;
  to_every(MatrixDiag(m1)) = of_every(vr);	// Make the other multiplicand
  assert( !(m1 == pattern+10) );

  m *= ConstMatrixDiag(m1);
  for(register int i=m.q_col_lwb(); i<=m.q_col_upb(); i++)
  {
    vc1 = of_every(MatrixColumn(mm,i));
    vc1 *= vr(i);			// Do a column-wise multiplication
    vc2 = of_every(ConstMatrixColumn(m,i));
    verify_vector_identity(vc1, vc2);
  }

  cout << "\nDone\n";
}

/*
 *------------------------------------------------------------------------
 *				Root module
 */

main()
{
  cout << "\n\n" << _Minuses << 
          "\n\t\tVerify Operations on Vectors\n";

  test_allocation();
  test_element_op(20);
  test_binary_op(20);
  test_norms(20);
  test_matrix_slices(20);

  cout << "\nAll tests passed" << endl;
}