source: LinearAlgebra/src/unittests/MatrixContentUnitTest.cpp@ 6253ed

Action_Thermostats Add_AtomRandomPerturbation Add_FitFragmentPartialChargesAction Add_RotateAroundBondAction Add_SelectAtomByNameAction Added_ParseSaveFragmentResults AddingActions_SaveParseParticleParameters Adding_Graph_to_ChangeBondActions Adding_MD_integration_tests Adding_ParticleName_to_Atom Adding_StructOpt_integration_tests AtomFragments Automaking_mpqc_open AutomationFragmentation_failures Candidate_v1.5.4 Candidate_v1.6.0 Candidate_v1.6.1 ChangeBugEmailaddress ChangingTestPorts ChemicalSpaceEvaluator CombiningParticlePotentialParsing Combining_Subpackages Debian_Package_split Debian_package_split_molecuildergui_only Disabling_MemDebug Docu_Python_wait EmpiricalPotential_contain_HomologyGraph EmpiricalPotential_contain_HomologyGraph_documentation Enable_parallel_make_install Enhance_userguide Enhanced_StructuralOptimization Enhanced_StructuralOptimization_continued Example_ManyWaysToTranslateAtom Exclude_Hydrogens_annealWithBondGraph FitPartialCharges_GlobalError Fix_BoundInBox_CenterInBox_MoleculeActions Fix_ChargeSampling_PBC Fix_ChronosMutex Fix_FitPartialCharges Fix_FitPotential_needs_atomicnumbers Fix_ForceAnnealing Fix_IndependentFragmentGrids Fix_ParseParticles Fix_ParseParticles_split_forward_backward_Actions Fix_PopActions Fix_QtFragmentList_sorted_selection Fix_Restrictedkeyset_FragmentMolecule Fix_StatusMsg Fix_StepWorldTime_single_argument Fix_Verbose_Codepatterns Fix_fitting_potentials Fixes ForceAnnealing_goodresults ForceAnnealing_oldresults ForceAnnealing_tocheck ForceAnnealing_with_BondGraph ForceAnnealing_with_BondGraph_continued ForceAnnealing_with_BondGraph_continued_betteresults ForceAnnealing_with_BondGraph_contraction-expansion FragmentAction_writes_AtomFragments FragmentMolecule_checks_bonddegrees GeometryObjects Gui_Fixes Gui_displays_atomic_force_velocity ImplicitCharges IndependentFragmentGrids IndependentFragmentGrids_IndividualZeroInstances IndependentFragmentGrids_IntegrationTest IndependentFragmentGrids_Sole_NN_Calculation JobMarket_RobustOnKillsSegFaults JobMarket_StableWorkerPool JobMarket_unresolvable_hostname_fix MoreRobust_FragmentAutomation ODR_violation_mpqc_open PartialCharges_OrthogonalSummation PdbParser_setsAtomName PythonUI_with_named_parameters QtGui_reactivate_TimeChanged_changes Recreated_GuiChecks Rewrite_FitPartialCharges RotateToPrincipalAxisSystem_UndoRedo SaturateAtoms_findBestMatching SaturateAtoms_singleDegree StoppableMakroAction Subpackage_CodePatterns Subpackage_JobMarket Subpackage_LinearAlgebra Subpackage_levmar Subpackage_mpqc_open Subpackage_vmg Switchable_LogView ThirdParty_MPQC_rebuilt_buildsystem TrajectoryDependenant_MaxOrder TremoloParser_IncreasedPrecision TremoloParser_MultipleTimesteps TremoloParser_setsAtomName Ubuntu_1604_changes stable
Last change on this file since 6253ed was 94d5ac6, checked in by Frederik Heber <heber@…>, 12 years ago

FIX: As we use GSL internally, we are as of now required to use GPL v2 license.

  • GNU Scientific Library is used at every place in the code, especially the sub-package LinearAlgebra is based on it which in turn is used really everywhere in the remainder of MoleCuilder. Hence, we have to use the GPL license for the whole of MoleCuilder. In effect, GPL's COPYING was present all along and stated the terms of the GPL v2 license.
  • Hence, I added the default GPL v2 disclaimer to every source file and removed the note about a (actually missing) LICENSE file.
  • also, I added a help-redistribute action which again gives the disclaimer of the GPL v2.
  • also, I changed in the disclaimer that is printed at every program start in builder_init.cpp.
  • TEST: Added check on GPL statement present in every module to test CodeChecks project-disclaimer.
  • Property mode set to 100644
File size: 10.4 KB
Line 
1/*
2 * Project: MoleCuilder
3 * Description: creates and alters molecular systems
4 * Copyright (C) 2010-2012 University of Bonn. All rights reserved.
5 *
6 *
7 * This file is part of MoleCuilder.
8 *
9 * MoleCuilder is free software: you can redistribute it and/or modify
10 * it under the terms of the GNU General Public License as published by
11 * the Free Software Foundation, either version 2 of the License, or
12 * (at your option) any later version.
13 *
14 * MoleCuilder is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 * GNU General Public License for more details.
18 *
19 * You should have received a copy of the GNU General Public License
20 * along with MoleCuilder. If not, see <http://www.gnu.org/licenses/>.
21 */
22
23/*
24 * MatrixContentUnitTest.cpp
25 *
26 * Created on: Jan 8, 2010
27 * Author: heber
28 */
29
30// include config.h
31#ifdef HAVE_CONFIG_H
32#include <config.h>
33#endif
34
35using namespace std;
36
37#include <cppunit/CompilerOutputter.h>
38#include <cppunit/extensions/TestFactoryRegistry.h>
39#include <cppunit/ui/text/TestRunner.h>
40
41#include <cmath>
42#include <limits>
43
44// include headers that implement a archive in simple text format
45#include <boost/archive/text_oarchive.hpp>
46#include <boost/archive/text_iarchive.hpp>
47
48#include "MatrixContentUnitTest.hpp"
49
50#include "MatrixContent.hpp"
51#include "VectorContent.hpp"
52
53#ifdef HAVE_TESTRUNNER
54#include "UnitTestMain.hpp"
55#endif /*HAVE_TESTRUNNER*/
56
57/********************************************** Test classes **************************************/
58
59const std::string matrixA = " 9.962848439806617e-01 5.950413028139907e-01 3.172650036465889e-01\n\
60 2.070939109924951e-01 8.365712117773689e-01 3.980567886073226e-01\n\
61 1.021600573071414e-01 8.061359922326598e-02 8.493537275043391e-01\n\
62 6.955361514014282e-01 6.873206387346102e-02 2.649165458481005e-01";
63const std::string vectorS = "1.630290761001098e+00 7.624807744374841e-01 6.374093742147465e-01";
64
65// Registers the fixture into the 'registry'
66CPPUNIT_TEST_SUITE_REGISTRATION( MatrixContentTest );
67
68
69void MatrixContentTest::setUp()
70{
71 m = new MatrixContent(4,3);
72};
73
74void MatrixContentTest::tearDown()
75{
76 delete(m);
77};
78
79/** Unit Test for accessing matrix elements.
80 *
81 */
82void MatrixContentTest::AccessTest()
83{
84 // check whether all elements are initially zero
85 for (int i=0;i<4;i++)
86 for (int j=0;j<3;j++)
87 CPPUNIT_ASSERT_EQUAL( 0., m->at(i,j) );
88
89 // set
90 for (int i=0;i<4;i++)
91 for (int j=0;j<3;j++)
92 m->set(i,j, i*3+j );
93
94 // and check
95 double *ptr = NULL;
96 for (int i=0;i<4;i++)
97 for (int j=0;j<3;j++) {
98 CPPUNIT_ASSERT_EQUAL( (double)(i*3+j), m->at(i,j) );
99 ptr = m->Pointer(i,j);
100 CPPUNIT_ASSERT_EQUAL( (double)(i*3+j), *ptr );
101 }
102
103 // assignment
104 for (int i=0;i<4;i++)
105 for (int j=0;j<3;j++)
106 m->set(i,j, i*3+j );
107 MatrixContent *dest = new MatrixContent(4,3);
108 *dest = *m;
109 for (int i=0;i<4;i++)
110 for (int j=0;j<3;j++)
111 CPPUNIT_ASSERT_EQUAL( dest->at(i,j), m->at(i,j) );
112 delete(dest);
113
114 // out of bounds
115 //CPPUNIT_ASSERT_EQUAL(0., v->at(5,2) );
116 //CPPUNIT_ASSERT_EQUAL(0., v->at(2,17) );
117 //CPPUNIT_ASSERT_EQUAL(0., v->at(1024,140040) );
118 //CPPUNIT_ASSERT_EQUAL(0., v->at(-1,0) );
119 //CPPUNIT_ASSERT_EQUAL(0., v->at(0,-1) );
120 //CPPUNIT_ASSERT_EQUAL(0., v->at(-1,-1) );
121};
122
123/** Unit Test for initializating matrices.
124 *
125 */
126void MatrixContentTest::InitializationTest()
127{
128 // set zero
129 m->setZero();
130 for (int i=0;i<4;i++)
131 for (int j=0;j<3;j++)
132 CPPUNIT_ASSERT_EQUAL( 0., m->at(i,j) );
133
134 // set all
135 m->setValue(1.5);
136 for (int i=0;i<4;i++)
137 for (int j=0;j<3;j++)
138 CPPUNIT_ASSERT_EQUAL( 1.5, m->at(i,j) );
139
140 // set basis
141 m->setIdentity();
142 for (int i=0;i<4;i++)
143 for (int j=0;j<3;j++)
144 CPPUNIT_ASSERT_EQUAL( i == j ? 1. : 0. , m->at(i,j) );
145
146 // set from array
147 double array[] = { 1., 0., 0.,
148 0., 1., 0.,
149 0., 0., 1.,
150 0., 0., 0. };
151 m->setFromDoubleArray(array);
152 for (int i=0;i<4;i++)
153 for (int j=0;j<3;j++)
154 CPPUNIT_ASSERT_EQUAL( i == j ? 1. : 0. , m->at(i,j) );
155
156};
157
158/** Unit Test for copying matrices.
159 *
160 */
161void MatrixContentTest::CopyTest()
162{
163 // set basis
164 MatrixContent *dest = NULL;
165 for (int i=0;i<4;i++) {
166 m->setValue(i);
167 dest = new MatrixContent(m);
168 for (int j=0;j<3;j++)
169 CPPUNIT_ASSERT_EQUAL( m->at(i,j) , dest->at(i,j) );
170
171 delete(dest);
172 }
173};
174
175/** Unit Test for exchanging rows and columns.
176 *
177 */
178void MatrixContentTest::ExchangeTest()
179{
180 // set to 1,1,1,2, ...
181 for (int i=0;i<4;i++)
182 for (int j=0;j<3;j++)
183 m->set(i,j, i+1 );
184
185 // swap such that nothing happens
186 CPPUNIT_ASSERT_EQUAL( true, m->SwapColumns(1,2) );
187 for (int i=0;i<4;i++)
188 for (int j=0;j<3;j++)
189 CPPUNIT_ASSERT_EQUAL( (double)i+1., m->at(i,j) );
190
191 // swap two rows
192 CPPUNIT_ASSERT_EQUAL( true, m->SwapRows(1,2) );
193 for (int i=0;i<4;i++)
194 for (int j=0;j<3;j++)
195 switch (i) {
196 case 0:
197 CPPUNIT_ASSERT_EQUAL( 1., m->at(i,j) );
198 break;
199 case 1:
200 CPPUNIT_ASSERT_EQUAL( 3., m->at(i,j) );
201 break;
202 case 2:
203 CPPUNIT_ASSERT_EQUAL( 2., m->at(i,j) );
204 break;
205 case 3:
206 CPPUNIT_ASSERT_EQUAL( 4., m->at(i,j) );
207 break;
208 default:
209 CPPUNIT_ASSERT_EQUAL( -1., m->at(i,j) );
210 break;
211 }
212 // check that op is reversable
213 CPPUNIT_ASSERT_EQUAL( true, m->SwapRows(1,2) );
214 for (int i=0;i<4;i++)
215 for (int j=0;j<3;j++)
216 CPPUNIT_ASSERT_EQUAL( (double)i+1., m->at(i,j) );
217
218 // set to 1,2,3,1, ...
219 for (int i=0;i<4;i++)
220 for (int j=0;j<3;j++)
221 m->set(i,j, j+1. );
222
223 // swap such that nothing happens
224 CPPUNIT_ASSERT_EQUAL( true, m->SwapRows(0,2) );
225 for (int i=0;i<4;i++)
226 for (int j=0;j<3;j++)
227 CPPUNIT_ASSERT_EQUAL( (double)j+1., m->at(i,j) );
228
229 // swap two columns
230 CPPUNIT_ASSERT_EQUAL( true, m->SwapColumns(0,2) );
231 for (int i=0;i<4;i++)
232 for (int j=0;j<3;j++)
233 switch (j) {
234 case 0:
235 CPPUNIT_ASSERT_EQUAL( 3., m->at(i,j) );
236 break;
237 case 1:
238 CPPUNIT_ASSERT_EQUAL( 2., m->at(i,j) );
239 break;
240 case 2:
241 CPPUNIT_ASSERT_EQUAL( 1., m->at(i,j) );
242 break;
243 default:
244 CPPUNIT_ASSERT_EQUAL( -1., m->at(i,j) );
245 break;
246 }
247 // check that op is reversable
248 CPPUNIT_ASSERT_EQUAL( true, m->SwapColumns(0,2) );
249 for (int i=0;i<4;i++)
250 for (int j=0;j<3;j++)
251 CPPUNIT_ASSERT_EQUAL( (double)j+1., m->at(i,j) );
252
253
254 // set to 1,2,3,4, ...
255 MatrixContent *n = new MatrixContent(3,4);
256 for (int i=0;i<4;i++)
257 for (int j=0;j<3;j++) {
258 m->set(i,j, 3*i+j+1 );
259 n->set(j,i, 3*i+j+1 );
260 }
261 // transpose
262 MatrixContent res = (*m).transposed();
263 CPPUNIT_ASSERT( *n == res );
264 // second transpose
265 MatrixContent res2 = res.transposed();
266 CPPUNIT_ASSERT( *m == res2 );
267 delete n;
268};
269
270/** Unit Test for matrix properties.
271 *
272 */
273void MatrixContentTest::PropertiesTest()
274{
275 // is zero
276 m->setZero();
277 CPPUNIT_ASSERT_EQUAL( true, m->IsNull() );
278 CPPUNIT_ASSERT_EQUAL( false, m->IsPositive() );
279 CPPUNIT_ASSERT_EQUAL( false, m->IsNegative() );
280 CPPUNIT_ASSERT_EQUAL( true, m->IsNonNegative() );
281
282 // is positive
283 m->setValue(0.5);
284 CPPUNIT_ASSERT_EQUAL( false, m->IsNull() );
285 CPPUNIT_ASSERT_EQUAL( true, m->IsPositive() );
286 CPPUNIT_ASSERT_EQUAL( false, m->IsNegative() );
287 CPPUNIT_ASSERT_EQUAL( true, m->IsNonNegative() );
288
289 // is negative
290 m->setValue(-0.1);
291 CPPUNIT_ASSERT_EQUAL( false, m->IsNull() );
292 CPPUNIT_ASSERT_EQUAL( false, m->IsPositive() );
293 CPPUNIT_ASSERT_EQUAL( true, m->IsNegative() );
294 CPPUNIT_ASSERT_EQUAL( false, m->IsNonNegative() );
295
296 // is positive definite
297 CPPUNIT_ASSERT_EQUAL( false, m->IsPositiveDefinite() );
298};
299
300
301/** Unit test for reading from and writing matrix to stream
302 *
303 */
304void MatrixContentTest::ReadWriteTest()
305{
306 // set up matrix
307 for (int i=0;i<4;i++)
308 for (int j=0;j<3;j++)
309 m->set(i,j, i*3+j );
310
311 // write to stream
312 std::stringstream matrixstream;
313 m->write(matrixstream);
314
315 // parse in dimensions and check
316 std::pair<size_t, size_t> matrixdimensions = MatrixContent::preparseMatrixDimensions(matrixstream);
317 CPPUNIT_ASSERT_EQUAL( (size_t)4, matrixdimensions.first );
318 CPPUNIT_ASSERT_EQUAL( (size_t)3, matrixdimensions.second );
319 // parse in matrix and check
320 MatrixContent* n = new MatrixContent(4,3, matrixstream);
321 CPPUNIT_ASSERT_EQUAL( *m, *n );
322
323 // free matrix
324 delete n;
325}
326
327/** Unit test for MatrixContent::performSingularValueDecomposition().
328 *
329 */
330void MatrixContentTest::SVDTest()
331{
332 // setup "A", U,S,V
333 std::stringstream Astream(matrixA);
334 std::stringstream Sstream(vectorS);
335 MatrixContent A(m->getRows(), m->getColumns(), Astream);
336 VectorContent S_expected(m->getColumns(), Sstream);
337 *m = A;
338
339 // check SVD
340 VectorContent S(m->getColumns());
341 MatrixContent V(m->getColumns(), m->getColumns());
342 A.performSingularValueDecomposition(V,S);
343 MatrixContent S_diag(m->getColumns(), m->getColumns());
344 for (size_t index = 0; index < m->getColumns(); ++index)
345 S_diag.set(index, index, S[index]);
346 CPPUNIT_ASSERT_EQUAL(*m,A*S_diag*(V.transpose()));
347 for (size_t index = 0; index < m->getColumns(); ++index) {
348 //CPPUNIT_ASSERT_EQUAL(S_expected[index], S[index]);
349 CPPUNIT_ASSERT(fabs(S_expected[index] - S[index]) < numeric_limits<double>::epsilon()*10.);
350 }
351
352 // set up right-hand side
353 VectorContent b(4);
354 b[0] = 1.;
355 b[1] = 0.;
356 b[2] = 0.;
357 b[3] = 0.;
358
359 // SVD
360 VectorContent x_expected(3);
361 x_expected[0] = -0.00209169;
362 x_expected[1] = -0.325399;
363 x_expected[2] = 0.628004;
364 const VectorContent x_result(A.solveOverdeterminedLinearEquation(V,S,b));
365
366 // check result
367 for (size_t index = 0; index < m->getColumns(); ++index) {
368 CPPUNIT_ASSERT(fabs(x_expected[index] - x_result[index]) < 1.e-6);
369 }
370}
371
372/** Unit test for serialization
373 *
374 */
375void MatrixContentTest::SerializationTest()
376{
377 m->setZero();
378 int count=0;
379 for (size_t x = 0; x < 4 ; ++x)
380 for (size_t y = 0; y < 3 ; ++y)
381 m->at(x,y) = count++;
382 // write element to stream
383 std::stringstream stream;
384 boost::archive::text_oarchive oa(stream);
385 oa << m;
386
387 //std::cout << "Contents of archive is " << stream.str() << std::endl;
388
389 // create and open an archive for input
390 boost::archive::text_iarchive ia(stream);
391 // read class state from archive
392 MatrixContent *newm;
393
394 ia >> newm;
395
396 CPPUNIT_ASSERT (*m == *newm);
397}
Note: See TracBrowser for help on using the repository browser.