[b4cf2b] | 1 | /*
|
---|
| 2 | * Project: MoleCuilder
|
---|
| 3 | * Description: creates and alters molecular systems
|
---|
| 4 | * Copyright (C) 2010 University of Bonn. All rights reserved.
|
---|
| 5 | * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
|
---|
| 6 | */
|
---|
| 7 |
|
---|
| 8 | /*
|
---|
| 9 | * SubspaceFactorizerUnittest.cpp
|
---|
| 10 | *
|
---|
| 11 | * Created on: Nov 13, 2010
|
---|
| 12 | * Author: heber
|
---|
| 13 | */
|
---|
| 14 |
|
---|
| 15 | // include config.h
|
---|
| 16 | #ifdef HAVE_CONFIG_H
|
---|
| 17 | #include <config.h>
|
---|
| 18 | #endif
|
---|
| 19 |
|
---|
| 20 | #include <cppunit/CompilerOutputter.h>
|
---|
| 21 | #include <cppunit/extensions/TestFactoryRegistry.h>
|
---|
| 22 | #include <cppunit/ui/text/TestRunner.h>
|
---|
| 23 |
|
---|
| 24 | #include <cmath>
|
---|
| 25 |
|
---|
| 26 | #include <gsl/gsl_vector.h>
|
---|
[40be55] | 27 | #include <boost/shared_ptr.hpp>
|
---|
[b4cf2b] | 28 |
|
---|
[40be55] | 29 | #include "Helpers/Assert.hpp"
|
---|
| 30 | #include "Helpers/Log.hpp"
|
---|
| 31 | #include "Helpers/toString.hpp"
|
---|
| 32 | #include "Helpers/Verbose.hpp"
|
---|
[b4cf2b] | 33 | #include "LinearAlgebra/VectorContent.hpp"
|
---|
| 34 | #include "LinearAlgebra/MatrixContent.hpp"
|
---|
| 35 |
|
---|
[40be55] | 36 | #include "SubspaceFactorizerUnittest.hpp"
|
---|
| 37 |
|
---|
[b4cf2b] | 38 | #ifdef HAVE_TESTRUNNER
|
---|
| 39 | #include "UnitTestMain.hpp"
|
---|
| 40 | #endif /*HAVE_TESTRUNNER*/
|
---|
| 41 |
|
---|
| 42 | // Registers the fixture into the 'registry'
|
---|
| 43 | CPPUNIT_TEST_SUITE_REGISTRATION( SubspaceFactorizerUnittest );
|
---|
| 44 |
|
---|
| 45 | void SubspaceFactorizerUnittest::setUp(){
|
---|
| 46 | fourbyfour = new MatrixContent(4,4);
|
---|
| 47 | fourbyfour->setZero();
|
---|
| 48 | for (int i=0; i<4 ; i++) {
|
---|
| 49 | fourbyfour->set(i,i, 2.);
|
---|
| 50 | if (i < (4-1)) {
|
---|
| 51 | fourbyfour->set(i+1,i, 1.);
|
---|
| 52 | fourbyfour->set(i,i+1, 1.);
|
---|
| 53 | }
|
---|
| 54 | }
|
---|
[40be55] | 55 | transformation = new MatrixContent**[3];
|
---|
| 56 |
|
---|
| 57 | // 1d subspace
|
---|
| 58 | transformation[0] = new MatrixContent*[4];
|
---|
| 59 | for(size_t i=0; i<4;++i) {
|
---|
| 60 | transformation[0][i] = new MatrixContent(4,4);
|
---|
| 61 | transformation[0][i]->setZero();
|
---|
| 62 | for (size_t j=0; j<1; ++j)
|
---|
| 63 | transformation[0][i]->set(i+j,i+j, 1.);
|
---|
| 64 | // std::cout << i << "th transformation matrix, " << 1 << "d subspace is "
|
---|
| 65 | // << *transformation[0][i] << std::endl;
|
---|
| 66 | }
|
---|
| 67 |
|
---|
| 68 | // 2d subspace
|
---|
| 69 | transformation[1] = new MatrixContent*[3];
|
---|
| 70 | for(size_t i=0; i<3;++i) {
|
---|
| 71 | transformation[1][i] = new MatrixContent(4,4);
|
---|
| 72 | transformation[1][i]->setZero();
|
---|
| 73 | for (size_t j=0; j<2; ++j)
|
---|
| 74 | transformation[1][i]->set(i+j,i+j, 1.);
|
---|
| 75 | // std::cout << i << "th transformation matrix, " << 2 << "d subspace is "
|
---|
| 76 | // << *transformation[1][i] << std::endl;
|
---|
| 77 | }
|
---|
| 78 |
|
---|
| 79 | // 3d subspace
|
---|
| 80 | transformation[2] = new MatrixContent*[2];
|
---|
| 81 | for(size_t i=0; i<2;++i) {
|
---|
| 82 | transformation[2][i] = new MatrixContent(4,4);
|
---|
| 83 | transformation[2][i]->setZero();
|
---|
| 84 | for (size_t j=0; j<3; ++j)
|
---|
| 85 | transformation[2][i]->set(i+j,i+j, 1.);
|
---|
| 86 | // std::cout << i << "th transformation matrix, " << 3 << "d subspace is "
|
---|
| 87 | // << *transformation[2][i] << std::endl;
|
---|
| 88 | }
|
---|
[b4cf2b] | 89 | }
|
---|
| 90 | void SubspaceFactorizerUnittest::tearDown(){
|
---|
[40be55] | 91 | // delete test matrix
|
---|
[b4cf2b] | 92 | delete fourbyfour;
|
---|
[40be55] | 93 |
|
---|
| 94 | // delete all transformations
|
---|
| 95 | for(size_t i=0; i<3;++i)
|
---|
| 96 | delete transformation[0][i];
|
---|
| 97 | delete[] transformation[0];
|
---|
| 98 | for(size_t i=0; i<3;++i)
|
---|
| 99 | delete transformation[1][i];
|
---|
| 100 | delete[] transformation[1];
|
---|
| 101 | for(size_t i=0; i<2;++i)
|
---|
| 102 | delete transformation[2][i];
|
---|
| 103 | delete[] transformation[2];
|
---|
| 104 | delete[] transformation;
|
---|
[b4cf2b] | 105 | }
|
---|
| 106 |
|
---|
| 107 | void SubspaceFactorizerUnittest::BlockTest()
|
---|
| 108 | {
|
---|
[40be55] | 109 | MatrixContent temp((*fourbyfour)&(*transformation[0][0]));
|
---|
[b4cf2b] | 110 | std::cout << "Our matrix is " << *fourbyfour << "." << std::endl;
|
---|
| 111 |
|
---|
[40be55] | 112 | std::cout << "Hadamard product of " << *fourbyfour << " with " << *transformation[0][0] << " is: " << std::endl;
|
---|
[b4cf2b] | 113 | std::cout << temp << std::endl;
|
---|
| 114 |
|
---|
| 115 | gsl_vector *eigenvalues = temp.transformToEigenbasis();
|
---|
[40be55] | 116 | VectorContent *eigenvaluesView = new VectorViewContent(gsl_vector_subvector(eigenvalues, 0, eigenvalues->size));
|
---|
[b4cf2b] | 117 | std::cout << "The resulting eigenbasis is " << temp
|
---|
[40be55] | 118 | << "\n\t with eigenvalues " << *eigenvaluesView << std::endl;
|
---|
| 119 | delete eigenvaluesView;
|
---|
[b4cf2b] | 120 | gsl_vector_free(eigenvalues);
|
---|
| 121 |
|
---|
[40be55] | 122 |
|
---|
| 123 | CPPUNIT_ASSERT_EQUAL(0,0);
|
---|
| 124 | }
|
---|
| 125 |
|
---|
| 126 | /** For given set of row and column indices, we extract the small block matrix.
|
---|
| 127 | * @param bigmatrix big matrix to extract from
|
---|
| 128 | * @param rowindexset row index set
|
---|
| 129 | * @param columnindexset column index set
|
---|
| 130 | * @return small matrix with dimension equal to the number of indices for row and column.
|
---|
| 131 | */
|
---|
| 132 | MatrixContent * getSubspaceMatrix(
|
---|
| 133 | MatrixContent &bigmatrix,
|
---|
| 134 | const IndexSet &rowindexset,
|
---|
| 135 | const IndexSet &columnindexset)
|
---|
| 136 | {
|
---|
| 137 | // check whether subsystem is big enough for both index sets
|
---|
| 138 | ASSERT(rowindexset.size() <= bigmatrix.getRows(),
|
---|
| 139 | "embedSubspaceMatrix() - bigmatrix has less rows "+toString(bigmatrix.getRows())
|
---|
| 140 | +" than needed by index set "
|
---|
| 141 | +toString(rowindexset.size())+"!");
|
---|
| 142 | ASSERT(columnindexset.size() <= bigmatrix.getColumns(),
|
---|
| 143 | "embedSubspaceMatrix() - bigmatrix has less columns "+toString(bigmatrix.getColumns())
|
---|
| 144 | +" than needed by index set "
|
---|
| 145 | +toString(columnindexset.size())+"!");
|
---|
| 146 | MatrixContent *subsystem = new MatrixContent(rowindexset.size(), columnindexset.size());
|
---|
| 147 | size_t localrow = 0; // local row indices for the subsystem
|
---|
| 148 | size_t localcolumn = 0;
|
---|
| 149 | for (IndexSet::const_iterator rowindex = rowindexset.begin();
|
---|
| 150 | rowindex != rowindexset.end();
|
---|
| 151 | ++rowindex) {
|
---|
| 152 | localcolumn = 0;
|
---|
| 153 | for (IndexSet::const_iterator columnindex = columnindexset.begin();
|
---|
| 154 | columnindex != columnindexset.end();
|
---|
| 155 | ++columnindex) {
|
---|
| 156 | ASSERT((*rowindex < bigmatrix.getRows()) && (*columnindex < bigmatrix.getColumns()),
|
---|
| 157 | "current index pair ("
|
---|
| 158 | +toString(*rowindex)+","+toString(*columnindex)
|
---|
| 159 | +") is out of bounds of bigmatrix ("
|
---|
| 160 | +toString(bigmatrix.getRows())+","+toString(bigmatrix.getColumns())
|
---|
| 161 | +")");
|
---|
| 162 | subsystem->at(localrow,localcolumn) = bigmatrix.at(*rowindex, *columnindex);
|
---|
| 163 | localcolumn++;
|
---|
| 164 | }
|
---|
| 165 | localrow++;
|
---|
| 166 | }
|
---|
| 167 | return subsystem;
|
---|
| 168 | }
|
---|
| 169 |
|
---|
| 170 | /** For a given set of row and columns indices, we embed a small block matrix into a bigger space.
|
---|
| 171 | *
|
---|
| 172 | * @param bigmatrix big matrix to place submatrix into
|
---|
| 173 | * @param rowindexset
|
---|
| 174 | * @param columnindexset
|
---|
| 175 | * @return
|
---|
| 176 | */
|
---|
| 177 | void embedSubspaceMatrix(
|
---|
| 178 | MatrixContent &bigmatrix,
|
---|
| 179 | MatrixContent &subsystem,
|
---|
| 180 | const IndexSet &rowindexset,
|
---|
| 181 | const IndexSet &columnindexset)
|
---|
| 182 | {
|
---|
| 183 | // check whether bigmatrix is at least as big as subsystem
|
---|
| 184 | ASSERT(subsystem.getRows() <= bigmatrix.getRows(),
|
---|
| 185 | "embedSubspaceMatrix() - subsystem has more rows "+toString(subsystem.getRows())
|
---|
| 186 | +" than bigmatrix "
|
---|
| 187 | +toString(bigmatrix.getRows())+"!");
|
---|
| 188 | ASSERT(subsystem.getColumns() <= bigmatrix.getColumns(),
|
---|
| 189 | "embedSubspaceMatrix() - subsystem has more columns "+toString(subsystem.getColumns())
|
---|
| 190 | +" than bigmatrix "
|
---|
| 191 | +toString(bigmatrix.getColumns())+"!");
|
---|
| 192 | // check whether subsystem is big enough for both index sets
|
---|
| 193 | ASSERT(rowindexset.size() <= subsystem.getRows(),
|
---|
| 194 | "embedSubspaceMatrix() - subsystem has less rows "+toString(subsystem.getRows())
|
---|
| 195 | +" than needed by index set "
|
---|
| 196 | +toString(rowindexset.size())+"!");
|
---|
| 197 | ASSERT(columnindexset.size() <= subsystem.getColumns(),
|
---|
| 198 | "embedSubspaceMatrix() - subsystem has less columns "+toString(subsystem.getColumns())
|
---|
| 199 | +" than needed by index set "
|
---|
| 200 | +toString(columnindexset.size())+"!");
|
---|
| 201 | size_t localrow = 0; // local row indices for the subsystem
|
---|
| 202 | size_t localcolumn = 0;
|
---|
| 203 | for (IndexSet::const_iterator rowindex = rowindexset.begin();
|
---|
| 204 | rowindex != rowindexset.end();
|
---|
| 205 | ++rowindex) {
|
---|
| 206 | localcolumn = 0;
|
---|
| 207 | for (IndexSet::const_iterator columnindex = columnindexset.begin();
|
---|
| 208 | columnindex != columnindexset.end();
|
---|
| 209 | ++columnindex) {
|
---|
| 210 | bigmatrix.at(*rowindex, *columnindex) = subsystem.at(localrow,localcolumn);
|
---|
| 211 | localcolumn++;
|
---|
| 212 | }
|
---|
| 213 | localrow++;
|
---|
| 214 | }
|
---|
| 215 | }
|
---|
| 216 |
|
---|
| 217 | /** Operator for output to std:: ostream operator of an IndexSet.
|
---|
| 218 | * @param ost output stream
|
---|
| 219 | * @param indexset index set to output
|
---|
| 220 | * @return ost output stream
|
---|
| 221 | */
|
---|
| 222 | std::ostream & operator<<(std::ostream &ost, const IndexSet &indexset)
|
---|
| 223 | {
|
---|
| 224 | ost << "{ ";
|
---|
| 225 | for (IndexSet::const_iterator iter = indexset.begin();
|
---|
| 226 | iter != indexset.end();
|
---|
| 227 | ++iter)
|
---|
| 228 | ost << *iter << " ";
|
---|
| 229 | ost << "}";
|
---|
| 230 | return ost;
|
---|
| 231 | }
|
---|
| 232 |
|
---|
| 233 | void SubspaceFactorizerUnittest::EigenvectorTest()
|
---|
| 234 | {
|
---|
| 235 | VectorArray CurrentEigenvectors;
|
---|
| 236 | VectorList *ParallelEigenvectorList = new VectorList[4];
|
---|
| 237 | IndexSet AllIndices;
|
---|
| 238 |
|
---|
| 239 | // create the total index set
|
---|
| 240 | for (size_t i=0;i<4;++i)
|
---|
| 241 | AllIndices.insert(i);
|
---|
| 242 |
|
---|
| 243 | // create all consecutive index subsets for dim 1 to 3
|
---|
| 244 | IndexMap Dimension_to_Indexset;
|
---|
| 245 | for (size_t dim = 0; dim<3;++dim) {
|
---|
| 246 | for (size_t i=0;i<4-dim;++i) {
|
---|
| 247 | IndexSet *indexset = new IndexSet;
|
---|
| 248 | for (size_t j=0;j<=dim;++j) {
|
---|
| 249 | //std::cout << "Putting " << i+j << " into " << i << "th map " << dim << std::endl;
|
---|
| 250 | CPPUNIT_ASSERT_MESSAGE("index "+toString(i+j)+" already present in "+toString(dim)+"-dim "+toString(i)+"th indexset.", indexset->count(i+j) == 0);
|
---|
| 251 | indexset->insert(i+j);
|
---|
| 252 | }
|
---|
| 253 | Dimension_to_Indexset.insert( make_pair(dim, boost::shared_ptr<IndexSet>(indexset)) );
|
---|
| 254 | // no need to free indexset, is stored in shared_ptr and
|
---|
| 255 | // will get released when Dimension_to_Indexset is destroyed
|
---|
| 256 | }
|
---|
| 257 | }
|
---|
| 258 |
|
---|
| 259 | // set to first guess, i.e. the unit vectors of R^4
|
---|
| 260 | for (IndexSet::const_iterator index = AllIndices.begin();
|
---|
| 261 | index != AllIndices.end();
|
---|
| 262 | ++index) {
|
---|
| 263 | boost::shared_ptr<VectorContent> EV(new VectorContent(4));
|
---|
| 264 | EV->setZero();
|
---|
| 265 | EV->at(*index) = 1.;
|
---|
| 266 | CurrentEigenvectors.push_back(EV);
|
---|
| 267 | }
|
---|
| 268 |
|
---|
| 269 | // for every dimension
|
---|
| 270 | for (size_t dim = 0; dim<3;++dim) {
|
---|
| 271 | // for every index set of this dimension
|
---|
| 272 | Log() << Verbose(0) << std::endl << std::endl;
|
---|
| 273 | Log() << Verbose(0) << "Current dimension is " << dim << std::endl;
|
---|
| 274 | std::pair<IndexMap::const_iterator,IndexMap::const_iterator> Bounds = Dimension_to_Indexset.equal_range(dim);
|
---|
| 275 | for (IndexMap::const_iterator IndexsetIter = Bounds.first;
|
---|
| 276 | IndexsetIter != Bounds.second;
|
---|
| 277 | ++IndexsetIter) {
|
---|
| 278 | // show the index set
|
---|
| 279 | Log() << Verbose(0) << std::endl;
|
---|
| 280 | Log() << Verbose(1) << "Current index set is " << *(IndexsetIter->second) << std::endl;
|
---|
| 281 |
|
---|
| 282 | // create transformation matrices from these
|
---|
| 283 | MatrixContent *subsystem = getSubspaceMatrix(*fourbyfour, *(IndexsetIter->second), *(IndexsetIter->second));
|
---|
| 284 | Log() << Verbose(2) << "Subsystem matrix is " << *subsystem << std::endl;
|
---|
| 285 |
|
---|
| 286 | // solve _small_ systems for eigenvalues
|
---|
| 287 | VectorContent *Eigenvalues = new VectorContent(subsystem->transformToEigenbasis());
|
---|
| 288 | Log() << Verbose(2) << "Eigenvector matrix is " << *subsystem << std::endl;
|
---|
| 289 | Log() << Verbose(2) << "Eigenvalues are " << *Eigenvalues << std::endl;
|
---|
| 290 | delete Eigenvalues;
|
---|
| 291 |
|
---|
| 292 | // blow up eigenvectors to 4dim column vector again
|
---|
| 293 | MatrixContent *Eigenvectors = new MatrixContent(4,4);
|
---|
| 294 | Eigenvectors->setZero();
|
---|
| 295 | embedSubspaceMatrix(*Eigenvectors, *subsystem, *(IndexsetIter->second), *(IndexsetIter->second));
|
---|
| 296 | Log() << Verbose(1) << "4x4 Eigenvector matrix is " << *Eigenvectors << std::endl;
|
---|
| 297 |
|
---|
| 298 | // we don't need the subsystem anymore
|
---|
| 299 | delete subsystem;
|
---|
| 300 |
|
---|
| 301 | // copy the index set, we look for one new eigenvector for each old one
|
---|
| 302 | IndexSet CorrespondenceList((*IndexsetIter->second));
|
---|
| 303 |
|
---|
| 304 | // go through all eigenvectors in this subspace
|
---|
| 305 | for (IndexSet::const_iterator iter = (*IndexsetIter->second).begin();
|
---|
| 306 | iter != (*IndexsetIter->second).end();
|
---|
| 307 | ++iter) {
|
---|
| 308 | // we rather copy the column vector, as the containing matrix is destroyed lateron
|
---|
| 309 | VectorContent *CurrentEigenvector = new VectorContent(Eigenvectors->getColumnVector(*iter));
|
---|
| 310 | Log() << Verbose(1) << "Current Eigenvector is " << *CurrentEigenvector << std::endl;
|
---|
| 311 |
|
---|
| 312 | // recognize eigenvectors parallel to existing ones
|
---|
| 313 | // (for now settle with the one we are most parallel to)
|
---|
| 314 | size_t mostparallel_index = 4;
|
---|
| 315 | double mostparallel_scalarproduct = 0.;
|
---|
| 316 | for (IndexSet::const_iterator indexiter = CorrespondenceList.begin();
|
---|
| 317 | indexiter != CorrespondenceList.end();
|
---|
| 318 | ++indexiter) {
|
---|
| 319 | Log() << Verbose(2) << "Comparing to old " << *indexiter << "th eigenvector " << *(CurrentEigenvectors[*indexiter]) << std::endl;
|
---|
| 320 | const double scalarproduct = (*(CurrentEigenvectors[*indexiter])) * (*CurrentEigenvector);
|
---|
| 321 | Log() << Verbose(2) << "SKP is " << scalarproduct << std::endl;
|
---|
| 322 | if (fabs(scalarproduct) > mostparallel_scalarproduct) {
|
---|
| 323 | mostparallel_scalarproduct = fabs(scalarproduct);
|
---|
| 324 | mostparallel_index = *indexiter;
|
---|
| 325 | }
|
---|
| 326 | }
|
---|
| 327 | if (mostparallel_index != 4) {
|
---|
| 328 | // put into std::list for later use
|
---|
| 329 | Log() << Verbose(1) << "Pushing " << *CurrentEigenvector << " into parallel list [" << mostparallel_index << "]" << std::endl;
|
---|
| 330 | ParallelEigenvectorList[mostparallel_index].push_back(boost::shared_ptr<VectorContent>(CurrentEigenvector));
|
---|
| 331 | CorrespondenceList.erase(mostparallel_index);
|
---|
| 332 | }
|
---|
| 333 | // no need to delete CurrentEigenvector as is taken care of by shared_ptr
|
---|
| 334 | }
|
---|
| 335 | // free eigenvectors
|
---|
| 336 | delete Eigenvectors;
|
---|
| 337 | }
|
---|
| 338 | }
|
---|
| 339 |
|
---|
| 340 | // print list of parallel eigenvectors
|
---|
| 341 | for (IndexSet::const_iterator index = AllIndices.begin();
|
---|
| 342 | index != AllIndices.end();
|
---|
| 343 | ++index) {
|
---|
| 344 | Log() << Verbose(0) << "Similar to " << *index << "th current eigenvector " << *(CurrentEigenvectors[*index]) << " are:" << std::endl;
|
---|
| 345 | for (VectorList::const_iterator iter = ParallelEigenvectorList[*index].begin();
|
---|
| 346 | iter != ParallelEigenvectorList[*index].end();
|
---|
| 347 | ++iter) {
|
---|
| 348 | Log() << Verbose(0) << **iter << std::endl;
|
---|
| 349 | }
|
---|
| 350 | Log() << Verbose(0) << std::endl;
|
---|
| 351 | }
|
---|
| 352 |
|
---|
| 353 | // create new CurrentEigenvectors from averaging parallel ones.
|
---|
| 354 | for (IndexSet::const_iterator index = AllIndices.begin();
|
---|
| 355 | index != AllIndices.end();
|
---|
| 356 | ++index) {
|
---|
| 357 | for (VectorList::const_iterator iter = ParallelEigenvectorList[*index].begin();
|
---|
| 358 | iter != ParallelEigenvectorList[*index].end();
|
---|
| 359 | ++iter) {
|
---|
| 360 | *CurrentEigenvectors[*index] += **iter;
|
---|
| 361 | }
|
---|
| 362 | *CurrentEigenvectors[*index] *= 1./(double)(ParallelEigenvectorList[*index].size()+1);
|
---|
| 363 | }
|
---|
| 364 |
|
---|
| 365 | // show new ones
|
---|
| 366 | Log() << Verbose(0) << "Resulting new eigenvectors are:" << std::endl;
|
---|
| 367 | for (IndexSet::const_iterator index = AllIndices.begin();
|
---|
| 368 | index != AllIndices.end();
|
---|
| 369 | ++index) {
|
---|
| 370 | Log() << Verbose(0) << *CurrentEigenvectors[*index] << std::endl;
|
---|
| 371 | }
|
---|
| 372 |
|
---|
| 373 |
|
---|
| 374 | delete[] ParallelEigenvectorList;
|
---|
| 375 |
|
---|
[b4cf2b] | 376 | CPPUNIT_ASSERT_EQUAL(0,0);
|
---|
| 377 | }
|
---|
| 378 |
|
---|