source: trunk/test/matrix_test.cc @ 1009

Last change on this file since 1009 was 1009, checked in by Peter, 14 years ago

merging branch peter-dev into trunk delta 1008:994

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