source: trunk/test/matrix_test.cc @ 793

Last change on this file since 793 was 793, checked in by Jari Häkkinen, 16 years ago

Fixes #193. Last test for the new functionality added to matrix.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 11.1 KB
Line 
1// $Id: matrix_test.cc 793 2007-03-11 15:40:13Z jari $
2
3/*
4  Copyright (C) The authors contributing to this file.
5
6  This file is part of the yat library, http://lev.thep.lu.se/trac/yat
7
8  The yat library is free software; you can redistribute it and/or
9  modify it under the terms of the GNU General Public License as
10  published by the Free Software Foundation; either version 2 of the
11  License, or (at your option) any later version.
12
13  The yat library is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16  General Public License for more details.
17
18  You should have received a copy of the GNU General Public License
19  along with this program; if not, write to the Free Software
20  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
21  02111-1307, USA.
22*/
23
24#include "yat/utility/matrix.h"
25
26#include <cstdio>
27#include <fstream>
28#include <iostream>
29
30class matrixwrapper
31{
32public:
33  matrixwrapper(size_t i,size_t j,double value=7)
34    : m_(i,j,value) {}
35
36  inline theplu::yat::utility::vector
37  row(const size_t& i) { return theplu::yat::utility::vector(m_,i); }
38
39  inline const theplu::yat::utility::matrix& matrix(void) const { return m_; }
40
41private:
42  theplu::yat::utility::matrix m_;
43};
44
45
46int main(const int argc,const char* argv[])
47{
48  using namespace theplu::yat;
49  std::ostream* error;
50  if (argc>1 && argv[1]==std::string("-v"))
51    error = &std::cerr;
52  else {
53    error = new std::ofstream("/dev/null");
54    if (argc>1)
55      std::cout << "matrix_test -v : for printing extra information\n";
56  }
57
58  *error << "Testing matrix class" << std::endl;
59  bool ok = true;
60  utility::matrix unit3x3(3,3);
61  for (size_t i=0; i<unit3x3.rows(); ++i)
62    unit3x3(i,i)=1;
63
64  *error << "\tcopy constructor and operator!=" << std::endl;
65  utility::matrix m(3,3,9);
66  utility::matrix m2(m);
67  if (m2!=m)
68    ok=false;
69
70  *error << "\toutput operator and istream constructor" << std::endl;
71  // Checking that the matrix output operator writes a file that the
72  // input operator can read.
73  std::ofstream my_out("data/tmp_test_matrix.txt");
74  my_out << m2;
75  my_out.close();
76  std::ifstream is("data/tmp_test_matrix.txt");
77  utility::matrix m3(is);
78  is.close();
79  if (m3!=m2)
80    ok=false;
81  std::remove("data/tmp_test_matrix.txt");
82
83  *error << "\toperator*=(double)" << std::endl;
84  utility::matrix m4(3,3,1);
85  m4 *= 9;
86  if (m4!=m) {
87    ok=false;
88    *error << "error operator*=(double)" << std::endl;
89  }
90
91  *error << "\tsub-matrix" << std::endl;
92  // Checking that sub-matrices work, i.e. mutation to the view are
93  // reflected in the viewed object.
94  is.open("data/knni_matrix.data");
95  // The stream input is a proper matrix file, with some stray empty
96  // lines and other whitespaces. The file is not expected to break
97  // things.
98  utility::matrix m5(is);
99  is.close();
100  double m5_sum=0;
101  for (size_t i=0; i<m5.rows(); ++i)
102    for (size_t j=0; j<m5.columns(); ++j)
103      m5_sum+=m5(i,j);
104  utility::matrix* m5sub=new utility::matrix(m5,3,3,3,3);
105  double m5sub_sum=0;
106  for (size_t i=0; i<m5sub->rows(); ++i)
107    for (size_t j=0; j<m5sub->columns(); ++j) {
108      m5sub_sum+=(*m5sub)(i,j);
109      (*m5sub)(i,j)=1;
110    }
111  delete m5sub;
112  double m5_sum2=0;
113  for (size_t i=0; i<m5.rows(); ++i)
114    for (size_t j=0; j<m5.columns(); ++j)
115      m5_sum2+=m5(i,j);
116  if (m5_sum2-m5_sum-9+m5sub_sum>1e-13) {
117    ok=false;
118    *error << "error sub-matrix test" << std::endl;
119  }
120
121  // Test of const view implementation (make sure no zero pointer is
122  // used). Here we use the bad style of not making the view const!
123  // If test fails with a null pointer exception it is an uncatchable
124  // core dump!
125  {
126    *error << "\tconst view implementation" << std::endl;
127    // Change the next statement to
128    // const utility::matrix m(10,10,3.0);
129    // when ticket:202 is fixes
130    utility::matrix m(10,10,3.0);
131    utility::matrix mview(m,3,3,3,3);
132    // const utility::vector vview(vv,0,5,1); // this is the proper line
133    utility::matrix m2(3,3,2.0);
134    m2.mul(mview); // should work even without const since const arg passing
135    m2.div(mview); // should work even without const since const arg passing
136  }
137
138  // checking that copy constructor creates an independent object when
139  // a non-view matrix is copied
140  {
141    *error << "\tcopy constructor" << std::endl;
142    utility::matrix m2(m5);
143    ok &= (m2.rows()==m5.rows());
144    ok &= (m2.columns()==m5.columns());
145    ok &= (m2==m5);
146    ok &= (&m2 != &m5);
147    ok &= !m2.isview();
148  }
149
150  // checking that copy constructor creates an independent object when
151  // a view matrix is copied
152  {
153    *error << "\tcopy contructor on view" << std::endl;
154    // Change the next statement to
155    // const utility::matrix m(10,10,3.0);
156    // when ticket:202 is fixes
157    utility::matrix m(10,10,3.0);
158    utility::matrix mview(m,3,3,3,3);
159    utility::matrix m3(mview);
160    ok &= (mview.rows()==m3.rows());
161    ok &= (mview.columns()==m3.columns());
162    ok &= (m3 == mview);
163    ok &= (&m3 != &mview);
164    ok &= !m3.isview();
165  }
166
167  // checking that assignment operator throws an exception if matrices
168  // differ in size
169  {
170    *error << "\tassignment operator" << std::endl;
171    // GSL will catch the error in this test there for the GSL error
172    // handling must be disabled until after the exception is
173    // catched. The GSL error handler is reinstated after the
174    // try-catch construct.
175    gsl_error_handler_t* err_handler=gsl_set_error_handler_off();
176    bool exception_happens=false;
177    try {
178      utility::matrix m(m5.rows()+1,3,0.0);
179      m=m5;
180    } catch (utility::GSL_error& err) {
181      exception_happens=true;
182    }
183    if (!exception_happens) {
184      *error << "Matrix assignment operator did not throw expected exception"
185             << std::endl;
186      ok=false;
187    }
188    gsl_set_error_handler(err_handler);
189  }
190
191  // checking that assignment operator changes the underlying object when
192  // a view is changed.
193  {
194    *error << "\tassignment operator on view" << std::endl;
195    bool this_ok(true);
196    utility::matrix mat_view(m5,3,3,3,3);
197    utility::matrix m2(3,3,12.0);
198    mat_view=m2;
199    for (u_int i=0; i<mat_view.rows(); ++i)
200      for (u_int j=0; j<mat_view.columns(); ++j)
201        if (m5(i+3,j+3)!=mat_view(i,j))
202          this_ok=false;
203    if (!this_ok) {
204      *error << "FAIL: assignemnt operator on view" << std::endl;
205      ok=false;
206    }
207  }
208
209  // checking clone functionality
210  {
211    *error << "\tclone functionality" << std::endl;
212    bool this_ok=true;
213    *error << "\t\tcloning normal matrix" << std::endl;
214    utility::matrix mat2;
215    mat2.clone(m5);
216    if (mat2.rows()!=m5.rows() || mat2.columns()!=m5.columns())
217      this_ok=false;
218    else
219      for (u_int i=0; i<m5.rows(); ++i)
220        for (u_int j=0; j<m5.columns(); ++j)
221          if (mat2(i,j)!=m5(i,j))
222            this_ok=false;
223    if (mat2.gsl_matrix_p()==m5.gsl_matrix_p())
224      this_ok=false;
225    *error << "\t\tcloning matrix view (sub matrix)" << std::endl;
226    utility::matrix* mat_view=new utility::matrix(m5,3,3,3,3);
227    utility::matrix mat_view2;
228    mat_view2.clone(*mat_view);
229    if (!mat_view2.isview())
230      this_ok=false;
231    if ( (mat_view->rows()!=mat_view2.rows()) ||
232         (mat_view->columns()!=mat_view2.columns()) )
233      this_ok=false;
234    else
235      for (u_int i=0; i<mat_view2.rows(); ++i)
236        for (u_int j=0; j<mat_view2.columns(); ++j)
237          if ((*mat_view)(i,j)!=mat_view2(i,j))
238            this_ok=false;
239    *error << "\t\tcloned matrix view independence" << std::endl;
240    delete mat_view;
241    for (u_int i=0; i<mat_view2.rows(); ++i)
242      for (u_int j=0; j<mat_view2.columns(); ++j)
243        if (m5(i+3,j+3)!=mat_view2(i,j))
244          this_ok=false;
245    if (!this_ok) {
246      *error << "FAIL: clone test" << std::endl;
247      ok=false;
248    }
249  }
250
251  *error << "\tsub-(row)vector" << std::endl;
252  // Checking that the row view works, i.e. mutation to the view are
253  // reflected in the viewed object.
254  m5_sum2=0;
255  for (size_t i=0; i<m5.rows(); ++i)
256    for (size_t j=0; j<m5.columns(); ++j)
257      m5_sum2+=m5(i,j);
258  utility::vector v5subrow(m5,3);
259  double v5subrow_sum=0;
260  for (size_t i=0; i<v5subrow.size(); ++i) {
261    v5subrow_sum+=v5subrow(i);
262    v5subrow[i]=0;
263  }
264  double m5_sum3=0;
265  for (size_t i=0; i<m5.rows(); ++i)
266    for (size_t j=0; j<m5.columns(); ++j)
267      m5_sum3+=m5(i,j);
268  if (m5_sum3-m5_sum2+v5subrow_sum>1e-13) {
269    ok=false;
270    *error << "error sub-vector test 1" << std::endl;
271  }
272
273  *error << "\tsub-(column)vector" << std::endl;
274  // Checking that the column view works, i.e. mutation to the view
275  // are reflected in the viewed object.
276  m5_sum3=0;
277  for (size_t i=0; i<m5.rows(); ++i)
278    for (size_t j=0; j<m5.columns(); ++j)
279      m5_sum3+=m5(i,j);
280  utility::vector v5subcolumn(m5,0,false);
281  double v5subcolumn_sum=0;
282  for (size_t i=0; i<v5subcolumn.size(); ++i) {
283    v5subcolumn_sum+=v5subcolumn(i);
284    v5subcolumn(i)=1;
285  }
286  double m5_sum4=0;
287  for (size_t i=0; i<m5.rows(); ++i)
288    for (size_t j=0; j<m5.columns(); ++j)
289      m5_sum4+=m5(i,j);
290  if (m5_sum4-m5_sum3-v5subcolumn.size()+v5subcolumn_sum>1e-13) {
291    ok=false;
292    *error << "error sub-vector test 2" << std::endl;
293  }
294
295  // Checking that the column view above mutates the values in the row
296  // view.
297  double v5subrow_sum2=0;
298  for (size_t i=0; i<v5subrow.size(); ++i)
299    v5subrow_sum2+=v5subrow(i);
300  if (v5subrow_sum2-v5subcolumn(3)>1e-13) {
301    ok=false;
302    *error << "error sub-vector test 3" << std::endl;
303  }
304
305  *error << "\tsub-vector and vector copying" << std::endl;
306  // Checking that a view is not inherited through the copy
307  // contructor.
308  utility::vector v6(v5subrow);
309  v6.set_all(2);
310  double v5subrow_sum3=0;
311  for (size_t i=0; i<v5subrow.size(); ++i)
312    v5subrow_sum3+=v5subrow(i);
313  if (v5subrow_sum3-v5subrow_sum2>1e-13) {
314    ok=false;
315    *error << "error sub-vector test 4" << std::endl;
316  }
317  // Checking that values in a vector is copied into a viewed matrix.
318  v5subrow.set(v6);
319  double m5_sum5=0;
320  for (size_t i=0; i<m5.rows(); ++i)
321    for (size_t j=0; j<m5.columns(); ++j)
322      m5_sum5+=m5(i,j);
323  if (m5_sum5-m5_sum4-v5subrow.size()*2+v5subrow_sum3>1e-13) {
324    ok=false;
325    *error << "error sub-vector test 5" << std::endl;
326  }
327
328  // Checking that the memberfunction matrixwrapper::row() returns a
329  // view
330  *error << "\tthat class member returns a view" << std::endl;
331  matrixwrapper mw(5,2);
332  utility::vector mwrow=mw.row(2);
333  if (mwrow.gsl_vector_p()->data != &(mw.matrix()(2,0))) {
334    ok=false;
335    *error << "error sub-vector test 7" << std::endl;
336  }
337
338  *error << "\tsub-matrix of a sub-matrix" << std::endl;
339  // Checking that a sub-matrix of a sub-matrix can be created.
340  utility::matrix sub(m5,3,3,3,3);
341  utility::matrix subsub(sub,2,1,1,2);
342  subsub(0,0)=23221121;
343  if (&sub(2,1)!=&subsub(0,0) || subsub(0,0)!=m5(5,4) || &subsub(0,0)!=&m5(5,4)){
344    ok=false;
345    *error << "error sub-matrix test 1" << std::endl;
346  }
347
348  *error << "\tmatrix::nan()" << std::endl;
349  is.open("data/sorlie_centroid_data.txt");
350  utility::matrix* m_nan = new utility::matrix(is,'\t');
351  utility::matrix m_weight;
352  utility::nan(*m_nan,m_weight);
353  is.close();
354  if (m_weight(0,0)!=1){
355    *error << "error in matrix::nan(): element(0,0) is " << m_weight(0,0)
356           << " expected 1." << std::endl;
357    ok=false;
358  }
359  if (m_weight(3,2)!=0){
360    *error << "error in matrix::nan(): element(3,2) is " << m_weight(3,2)
361           << " expected 0." << std::endl;
362    ok=false;
363  }
364  delete m_nan;
365
366  *error << "\toperator*=(matrix&)" << std::endl;
367  utility::matrix m6(unit3x3);
368  m6 *= m;
369  if (m6!=m) {
370    ok=false;
371    *error << "error operator*=(matrix) 1" << std::endl;
372  }
373  m6 *= unit3x3;
374  if (m6!=m) {
375    ok=false;
376    *error << "error operator*=(matrix) 2" << std::endl;
377  }
378
379  if (error!=&std::cerr)
380    delete error;
381
382  return (ok ? 0 : -1);
383}
Note: See TracBrowser for help on using the repository browser.