source: src/SubspaceFactorizer.cpp@ a0064e

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 a0064e was a0064e, checked in by Frederik Heber <heber@…>, 14 years ago

Moved stuff in src/Helpers and src/Patterns out to stand-alone project CodePatterns.

Makefile.am's:

  • CodePatterns added to AM_LDFLAGS and AM_CFLAGS
  • libMoleCuilderHelpers removed

Helpers/...

  • defs.hpp include prefixed with Helpers/
  • helpers.?pp removed lots of old, unused functions: bound, ask_value, ...
  • fast_functions.hpp has lots of functions removed as well.
  • all other files and unit tests moved to project CodePatterns.

Patterns/...

  • added ax_codepatterns.m4 containing AM_PATH_CODEPATTERNS to configure.ac
  • Property mode set to 100644
File size: 15.7 KB
Line 
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 * SubspaceFactorizer.cpp
10 *
11 * Created on: Nov 23, 2010
12 * Author: heber
13 */
14
15#include <cmath>
16
17#include <gsl/gsl_eigen.h>
18#include <gsl/gsl_matrix.h>
19#include <gsl/gsl_vector.h>
20#include <boost/foreach.hpp>
21#include <boost/shared_ptr.hpp>
22#include <boost/timer.hpp>
23
24#include "Helpers/Assert.hpp"
25#include "Helpers/defs.hpp"
26#include "Helpers/Log.hpp"
27#include "Helpers/toString.hpp"
28#include "Helpers/Verbose.hpp"
29#include "LinearAlgebra/Eigenspace.hpp"
30#include "LinearAlgebra/MatrixContent.hpp"
31#include "LinearAlgebra/Subspace.hpp"
32#include "LinearAlgebra/VectorContent.hpp"
33
34typedef std::set<std::set<size_t> > SetofIndexSets;
35typedef std::set<size_t> IndexSet;
36typedef std::multimap< size_t, boost::shared_ptr<Subspace> > SubspaceMap;
37typedef std::multimap< size_t, boost::shared_ptr<IndexSet> > IndexMap;
38typedef std::list< boost::shared_ptr<VectorContent> > VectorList;
39typedef std::list< std::pair<boost::shared_ptr<VectorContent>, double> > VectorValueList;
40typedef std::vector< boost::shared_ptr<VectorContent> > VectorArray;
41typedef std::vector< double > ValueArray;
42
43
44/** Iterative function to generate all power sets of indices of size \a maxelements.
45 *
46 * @param SetofSets Container for all sets
47 * @param CurrentSet pointer to current set in this container
48 * @param Indices Source set of indices
49 * @param maxelements number of elements of each set in final SetofSets
50 * @return true - generation continued, false - current set already had
51 * \a maxelements elements
52 */
53bool generatePowerSet(
54 SetofIndexSets &SetofSets,
55 SetofIndexSets::iterator &CurrentSet,
56 IndexSet &Indices,
57 const size_t maxelements)
58{
59 if (CurrentSet->size() < maxelements) {
60 // allocate the needed sets
61 const size_t size = Indices.size() - CurrentSet->size();
62 std::vector<std::set<size_t> > SetExpanded;
63 SetExpanded.reserve(size);
64
65 // copy the current set into each
66 for (size_t i=0;i<size;++i)
67 SetExpanded.push_back(*CurrentSet);
68
69 // expand each set by one index
70 size_t localindex=0;
71 BOOST_FOREACH(size_t iter, Indices) {
72 if (CurrentSet->count(iter) == 0) {
73 SetExpanded[localindex].insert(iter);
74 ++localindex;
75 }
76 }
77
78 // insert set at position of CurrentSet
79 for (size_t i=0;i<size;++i) {
80 //DoLog(1) && (Log() << Verbose(1) << "Inserting set #" << i << ": " << SetExpanded[i] << std::endl);
81 SetofSets.insert(CurrentSet, SetExpanded[i]);
82 }
83 SetExpanded.clear();
84
85 // and remove the current set
86 //SetofSets.erase(CurrentSet);
87 //CurrentSet = SetofSets.begin();
88
89 // set iterator to a valid position again
90 ++CurrentSet;
91 return true;
92 } else {
93 return false;
94 }
95}
96
97
98
99/** Prints the scalar product of each possible pair that is not orthonormal.
100 * We use class logger for printing.
101 * @param AllIndices set of all possible indices of the eigenvectors
102 * @param CurrentEigenvectors array of eigenvectors
103 * @return true - all are orthonormal to each other,
104 * false - some are not orthogonal or not normalized.
105 */
106bool checkOrthogonality(const IndexSet &AllIndices, const VectorArray &CurrentEigenvectors)
107{
108 size_t nonnormalized = 0;
109 size_t nonorthogonal = 0;
110 // check orthogonality
111 BOOST_FOREACH( size_t firstindex, AllIndices) {
112 BOOST_FOREACH( size_t secondindex, AllIndices) {
113 const double scp = (*CurrentEigenvectors[firstindex])*(*CurrentEigenvectors[secondindex]);
114 if (firstindex == secondindex) {
115 if (fabs(scp - 1.) > MYEPSILON) {
116 nonnormalized++;
117 Log() << Verbose(2) << "Vector " << firstindex << " is not normalized, off by "
118 << fabs(1.-(*CurrentEigenvectors[firstindex])*(*CurrentEigenvectors[secondindex])) << std::endl;
119 }
120 } else {
121 if (fabs(scp) > MYEPSILON) {
122 nonorthogonal++;
123 Log() << Verbose(2) << "Scalar product between " << firstindex << " and " << secondindex
124 << " is " << (*CurrentEigenvectors[firstindex])*(*CurrentEigenvectors[secondindex]) << std::endl;
125 }
126 }
127 }
128 }
129
130 if ((nonnormalized == 0) && (nonorthogonal == 0)) {
131 DoLog(1) && (DoLog(1) && (Log() << Verbose(1) << "All vectors are orthonormal to each other." << std::endl));
132 return true;
133 }
134 if ((nonnormalized == 0) && (nonorthogonal != 0))
135 DoLog(1) && (DoLog(1) && (Log() << Verbose(1) << "All vectors are normalized." << std::endl));
136 if ((nonnormalized != 0) && (nonorthogonal == 0))
137 DoLog(1) && (DoLog(1) && (Log() << Verbose(1) << "All vectors are orthogonal to each other." << std::endl));
138 return false;
139}
140
141/** Calculate the sum of the scalar product of each possible pair.
142 * @param AllIndices set of all possible indices of the eigenvectors
143 * @param CurrentEigenvectors array of eigenvectors
144 * @return sum of scalar products between all possible pairs
145 */
146double calculateOrthogonalityThreshold(const IndexSet &AllIndices, const VectorArray &CurrentEigenvectors)
147{
148 double threshold = 0.;
149 // check orthogonality
150 BOOST_FOREACH( size_t firstindex, AllIndices) {
151 BOOST_FOREACH( size_t secondindex, AllIndices) {
152 const double scp = (*CurrentEigenvectors[firstindex])*(*CurrentEigenvectors[secondindex]);
153 if (firstindex == secondindex) {
154 threshold += fabs(scp - 1.);
155 } else {
156 threshold += fabs(scp);
157 }
158 }
159 }
160 return threshold;
161}
162
163
164/** Operator for output to std::ostream operator of an IndexSet.
165 * @param ost output stream
166 * @param indexset index set to output
167 * @return ost output stream
168 */
169std::ostream & operator<<(std::ostream &ost, const IndexSet &indexset)
170{
171 ost << "{ ";
172 for (IndexSet::const_iterator iter = indexset.begin();
173 iter != indexset.end();
174 ++iter)
175 ost << *iter << " ";
176 ost << "}";
177 return ost;
178}
179
180
181int main(int argc, char **argv)
182{
183 size_t matrixdimension = 8;
184 size_t subspacelimit = 4;
185
186 if (argc < 2) {
187 std::cerr << "Usage: " << argv[0] << " <matrixdim> <subspacelimit>" << std::endl;
188 return 255;
189 } else {
190 {
191 std::stringstream s(toString(argv[1]));;
192 s >> matrixdimension;
193 }
194 {
195 std::stringstream s(toString(argv[2]));;
196 s >> subspacelimit;
197 }
198 }
199
200 MatrixContent *matrix = new MatrixContent(matrixdimension,matrixdimension);
201 matrix->setZero();
202 for (size_t i=0; i<matrixdimension ; i++) {
203 for (size_t j=0; j<= i; ++j) {
204 //const double value = 10. * rand() / (double)RAND_MAX;
205 //const double value = i==j ? 2. : 1.;
206 if (i==j)
207 matrix->set(i,i, 2.);
208 else if (j+1 == i) {
209 matrix->set(i,j, 1.);
210 matrix->set(j,i, 1.);
211 } else {
212 matrix->set(i,j, 0.);
213 matrix->set(j,i, 0.);
214 }
215 }
216 }
217
218 Eigenspace::eigenvectorset CurrentEigenvectors;
219 Eigenspace::eigenvalueset CurrentEigenvalues;
220
221 setVerbosity(3);
222
223 boost::timer Time_generatingfullspace;
224 DoLog(0) && (Log() << Verbose(0) << std::endl << std::endl);
225 // create the total index set
226 IndexSet AllIndices;
227 for (size_t i=0;i<matrixdimension;++i)
228 AllIndices.insert(i);
229 Eigenspace FullSpace(AllIndices, *matrix);
230 DoLog(1) && (Log() << Verbose(1) << "Generated full space: " << FullSpace << std::endl);
231 DoLog(0) && (Log() << Verbose(0) << "Full space generation took " << Time_generatingfullspace.elapsed() << " seconds." << std::endl);
232
233 // generate first set of eigenvectors
234 // set to first guess, i.e. the unit vectors of R^matrixdimension
235 BOOST_FOREACH( size_t index, AllIndices) {
236 boost::shared_ptr<VectorContent> EV(new VectorContent(matrixdimension));
237 EV->setZero();
238 EV->at(index) = 1.;
239 CurrentEigenvectors.push_back(EV);
240 CurrentEigenvalues.push_back(0.);
241 }
242
243 boost::timer Time_generatingsubsets;
244 DoLog(0) && (Log() << Verbose(0) << "Generating sub sets ..." << std::endl);
245 SetofIndexSets SetofSets;
246 // note that starting off empty set is unstable
247 IndexSet singleset;
248 BOOST_FOREACH(size_t iter, AllIndices) {
249 singleset.insert(iter);
250 SetofSets.insert(singleset);
251 singleset.clear();
252 }
253 SetofIndexSets::iterator CurrentSet = SetofSets.begin();
254 while (CurrentSet != SetofSets.end()) {
255 //DoLog(2) && (Log() << Verbose(2) << "Current set is " << *CurrentSet << std::endl);
256 if (!generatePowerSet(SetofSets, CurrentSet, AllIndices, subspacelimit)) {
257 // go to next set
258 ++CurrentSet;
259 }
260 }
261 DoLog(0) && (Log() << Verbose(0) << "Sub set generation took " << Time_generatingsubsets.elapsed() << " seconds." << std::endl);
262
263 // create a subspace to each set and and to respective level
264 boost::timer Time_generatingsubspaces;
265 DoLog(0) && (Log() << Verbose(0) << "Generating sub spaces ..." << std::endl);
266 SubspaceMap Dimension_to_Indexset;
267 BOOST_FOREACH(std::set<size_t> iter, SetofSets) {
268 boost::shared_ptr<Subspace> subspace(new Subspace(iter, FullSpace));
269 DoLog(1) && (Log() << Verbose(1) << "Current subspace is " << *subspace << std::endl);
270 Dimension_to_Indexset.insert( make_pair(iter.size(), boost::shared_ptr<Subspace>(subspace)) );
271 }
272
273 for (size_t dim = 1; dim<=subspacelimit;++dim) {
274 BOOST_FOREACH( SubspaceMap::value_type subspace, Dimension_to_Indexset.equal_range(dim)) {
275 if (dim != 0) { // from level 1 and onward
276 BOOST_FOREACH( SubspaceMap::value_type entry, Dimension_to_Indexset.equal_range(dim-1)) {
277 if (subspace.second->contains(*entry.second)) {
278 // if contained then add ...
279 subspace.second->addSubset(entry.second);
280 // ... and also its containees as they are all automatically contained as well
281 BOOST_FOREACH(boost::shared_ptr<Subspace> iter, entry.second->getSubIndices()) {
282 subspace.second->addSubset(iter);
283 }
284 }
285 }
286 }
287 }
288 }
289 DoLog(0) && (Log() << Verbose(0) << "Sub space generation took " << Time_generatingsubspaces.elapsed() << " seconds." << std::endl);
290
291 // create a file handle for the eigenvalues
292 std::ofstream outputvalues("eigenvalues.dat", std::ios_base::trunc);
293 ASSERT(outputvalues.good(),
294 "SubspaceFactorizerUnittest::EigenvectorTest() - failed to open eigenvalue file!");
295 outputvalues << "# iteration ";
296 BOOST_FOREACH(size_t iter, AllIndices) {
297 outputvalues << "\teigenvalue" << iter;
298 }
299 outputvalues << std::endl;
300
301 DoLog(0) && (Log() << Verbose(0) << "Solving ..." << std::endl);
302 boost::timer Time_solving;
303 size_t run=1; // counting iterations
304 double threshold = 1.; // containing threshold value
305 while ((threshold > MYEPSILON) && (run < 20)) {
306 // for every dimension
307 for (size_t dim = 1; dim <= subspacelimit;++dim) {
308 // for every index set of this dimension
309 DoLog(1) && (Log() << Verbose(1) << std::endl << std::endl);
310 DoLog(1) && (Log() << Verbose(1) << "Current dimension is " << dim << std::endl);
311 std::pair<SubspaceMap::const_iterator,SubspaceMap::const_iterator> Bounds = Dimension_to_Indexset.equal_range(dim);
312 for (SubspaceMap::const_iterator IndexsetIter = Bounds.first;
313 IndexsetIter != Bounds.second;
314 ++IndexsetIter) {
315 Subspace& subspace = *(IndexsetIter->second);
316 // show the index set
317 DoLog(2) && (Log() << Verbose(2) << std::endl);
318 DoLog(2) && (Log() << Verbose(2) << "Current subspace is " << subspace << std::endl);
319
320 // solve
321 subspace.calculateEigenSubspace();
322
323 // note that assignment to global eigenvectors all remains within subspace
324 }
325 }
326
327 // print list of similar eigenvectors
328 DoLog(2) && (Log() << Verbose(2) << std::endl);
329 BOOST_FOREACH( size_t index, AllIndices) {
330 DoLog(2) && (Log() << Verbose(2) << "Similar to " << index << "th current eigenvector " << *(CurrentEigenvectors[index]) << " are:" << std::endl);
331 BOOST_FOREACH( SubspaceMap::value_type iter, Dimension_to_Indexset) {
332 const VectorContent & CurrentEV = (iter.second)->getEigenvectorParallelToFullOne(index);
333 if (!CurrentEV.IsZero())
334 Log() << Verbose(2)
335 << "dim" << iter.first
336 << ", subspace{" << (iter.second)->getIndices()
337 << "}: "<< CurrentEV << std::endl;
338 }
339 DoLog(2) && (Log() << Verbose(2) << std::endl);
340 }
341
342 // create new CurrentEigenvectors from averaging parallel ones.
343 BOOST_FOREACH(size_t index, AllIndices) {
344 CurrentEigenvectors[index]->setZero();
345 CurrentEigenvalues[index] = 0.;
346 size_t count = 0;
347 BOOST_FOREACH( SubspaceMap::value_type iter, Dimension_to_Indexset) {
348 const VectorContent CurrentEV = (iter.second)->getEigenvectorParallelToFullOne(index);
349 *CurrentEigenvectors[index] += CurrentEV; // * (iter.second)->getEigenvalueOfEigenvectorParallelToFullOne(index);
350 CurrentEigenvalues[index] += (iter.second)->getEigenvalueOfEigenvectorParallelToFullOne(index);
351 if (!CurrentEV.IsZero())
352 count++;
353 }
354 *CurrentEigenvectors[index] *= 1./CurrentEigenvalues[index];
355 //CurrentEigenvalues[index] /= (double)count;
356 }
357
358 // check orthonormality
359 threshold = calculateOrthogonalityThreshold(AllIndices, CurrentEigenvectors);
360 bool dontOrthonormalization = checkOrthogonality(AllIndices, CurrentEigenvectors);
361
362 // orthonormalize
363 if (!dontOrthonormalization) {
364 DoLog(1) && (Log() << Verbose(1) << "Orthonormalizing ... " << std::endl);
365 for (IndexSet::const_iterator firstindex = AllIndices.begin();
366 firstindex != AllIndices.end();
367 ++firstindex) {
368 for (IndexSet::const_iterator secondindex = firstindex;
369 secondindex != AllIndices.end();
370 ++secondindex) {
371 if (*firstindex == *secondindex) {
372 (*CurrentEigenvectors[*secondindex]) *= 1./(*CurrentEigenvectors[*secondindex]).Norm();
373 } else {
374 (*CurrentEigenvectors[*secondindex]) -=
375 ((*CurrentEigenvectors[*firstindex])*(*CurrentEigenvectors[*secondindex]))
376 *(*CurrentEigenvectors[*firstindex]);
377 }
378 }
379 }
380 }
381
382// // check orthonormality again
383// checkOrthogonality(AllIndices, CurrentEigenvectors);
384
385 // put obtained eigenvectors into full space
386 FullSpace.setEigenpairs(CurrentEigenvectors, CurrentEigenvalues);
387
388 // show new ones
389 DoLog(1) && (Log() << Verbose(1) << "Resulting new eigenvectors and -values, run " << run << " are:" << std::endl);
390 outputvalues << run;
391 BOOST_FOREACH( size_t index, AllIndices) {
392 DoLog(1) && (Log() << Verbose(1) << *CurrentEigenvectors[index] << " with " << CurrentEigenvalues[index] << std::endl);
393 outputvalues << "\t" << CurrentEigenvalues[index];
394 }
395 outputvalues << std::endl;
396
397 // and next iteration
398 DoLog(0) && (Log() << Verbose(0) << "\titeration #" << run << std::endl);
399 run++;
400 }
401 DoLog(0) && (Log() << Verbose(0) << "Solving took " << Time_solving.elapsed() << " seconds." << std::endl);
402 // show final ones
403 DoLog(0) && (Log() << Verbose(0) << "Resulting new eigenvectors and -values, run " << run << " are:" << std::endl);
404 outputvalues << run;
405 BOOST_FOREACH( size_t index, AllIndices) {
406 DoLog(0) && (Log() << Verbose(0) << *CurrentEigenvectors[index] << " with " << CurrentEigenvalues[index] << std::endl);
407 outputvalues << "\t" << CurrentEigenvalues[index];
408 }
409 outputvalues << std::endl;
410 outputvalues.close();
411
412 setVerbosity(2);
413
414 DoLog(0) && (Log() << Verbose(0) << "Solving full space ..." << std::endl);
415 boost::timer Time_comparison;
416 MatrixContent tempFullspaceMatrix = FullSpace.getEigenspaceMatrix();
417 gsl_vector *eigenvalues = tempFullspaceMatrix.transformToEigenbasis();
418 tempFullspaceMatrix.sortEigenbasis(eigenvalues);
419 DoLog(0) && (Log() << Verbose(0) << "full space solution took " << Time_comparison.elapsed() << " seconds." << std::endl);
420
421 delete matrix;
422
423 return 0;
424}
Note: See TracBrowser for help on using the repository browser.