Changeset 40be55

Timestamp:

Dec 4, 2010, 11:56:27 PM (14 years ago)

Author:

Frederik Heber <heber@…>

Branches:

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

Children:

3dd9c7

Parents:

0fd3f2

git-author:

Frederik Heber <heber@…> (11/16/10 15:43:58)

git-committer:

Frederik Heber <heber@…> (12/04/10 23:56:27)

Message:

Current stable version of SubspaceFactorizer.

can:

• functions embedSubspaceMatrix() and getSubspaceMatrix() to embed and to obtain the small block matrix from/in a big one
• we construct the index set
• we construct the block matrices
• we solve the eigensystems
• we transfer the eigenvalues back
• we associate to original eigenvectors by parallelness, uniquely within a subset of indices.
• we construct the new set of eigenvectors by averaging.

can't:

• resulting ev are not orthonormalized or checked to be
• no iteration yet

Location:

src/unittests

Files:

• SubspaceFactorizerUnittest.cpp (modified) (2 diffs)
• SubspaceFactorizerUnittest.hpp (modified) (3 diffs)

src/unittests/SubspaceFactorizerUnittest.cpp

@@ -25 +25 @@
 
 #include <gsl/gsl_vector.h>
-
-#include "SubspaceFactorizerUnittest.hpp"
+#include <boost/shared_ptr.hpp>
+
+#include "Helpers/Assert.hpp"
+#include "Helpers/Log.hpp"
+#include "Helpers/toString.hpp"
+#include "Helpers/Verbose.hpp"
 #include "LinearAlgebra/VectorContent.hpp"
 #include "LinearAlgebra/MatrixContent.hpp"
+
+#include "SubspaceFactorizerUnittest.hpp"
 
 #ifdef HAVE_TESTRUNNER
@@ -47 +53 @@
     }
   }
-  permU1_1 = new MatrixContent(4,1);
-  permU1_1->setZero();
-  permU1_1->set(0,0, 1.);
-  permU1_2 = new MatrixContent(4,1);
-  permU1_2->setZero();
-  permU1_2->set(1,0, 1.);
-  permU1_3 = new MatrixContent(4,1);
-  permU1_3->setZero();
-  permU1_3->set(2,0, 1.);
-  permU1_4 = new MatrixContent(4,1);
-  permU1_4->setZero();
-  permU1_4->set(3,0, 1.);
-
-  permU2_1 = new MatrixContent(4,2);
-  permU2_1->setZero();
-  permU2_1->set(0,0, 1.);
-  permU2_1->set(1,1, 1.);
-  permU2_2 = new MatrixContent(4,2);
-  permU2_2->setZero();
-  permU2_2->set(1,0, 1.);
-  permU2_2->set(2,1, 1.);
-  permU2_3 = new MatrixContent(4,2);
-  permU2_3->setZero();
-  permU2_3->set(2,0, 1.);
-  permU2_3->set(3,1, 1.);
-
-  permU3_1 = new MatrixContent(4,3);
-  permU3_1->setZero();
-  permU3_1->set(0,0, 1.);
-  permU3_1->set(1,1, 1.);
-  permU3_1->set(2,2, 1.);
-  permU3_2 = new MatrixContent(4,3);
-  permU3_2->setZero();
-  permU3_2->set(1,0, 1.);
-  permU3_2->set(2,1, 1.);
-  permU3_2->set(3,2, 1.);
+  transformation = new MatrixContent**[3];
+
+  // 1d subspace
+  transformation[0] = new MatrixContent*[4];
+  for(size_t i=0; i<4;++i) {
+    transformation[0][i] = new MatrixContent(4,4);
+    transformation[0][i]->setZero();
+    for (size_t j=0; j<1; ++j)
+      transformation[0][i]->set(i+j,i+j, 1.);
+//    std::cout << i << "th transformation matrix, " << 1 << "d subspace is "
+//        << *transformation[0][i] << std::endl;
+  }
+
+  // 2d subspace
+  transformation[1] = new MatrixContent*[3];
+  for(size_t i=0; i<3;++i) {
+    transformation[1][i] = new MatrixContent(4,4);
+    transformation[1][i]->setZero();
+    for (size_t j=0; j<2; ++j)
+      transformation[1][i]->set(i+j,i+j, 1.);
+//    std::cout << i << "th transformation matrix, " << 2 << "d subspace is "
+//        << *transformation[1][i] << std::endl;
+  }
+
+  // 3d subspace
+  transformation[2] = new MatrixContent*[2];
+  for(size_t i=0; i<2;++i) {
+    transformation[2][i] = new MatrixContent(4,4);
+    transformation[2][i]->setZero();
+    for (size_t j=0; j<3; ++j)
+      transformation[2][i]->set(i+j,i+j, 1.);
+//    std::cout << i << "th transformation matrix, " << 3 << "d subspace is "
+//        << *transformation[2][i] << std::endl;
+  }
 }
 void SubspaceFactorizerUnittest::tearDown(){
+  // delete test matrix
   delete fourbyfour;
-  delete permU1_1;
-  delete permU1_2;
-  delete permU1_3;
-  delete permU1_4;
-  delete permU2_1;
-  delete permU2_2;
-  delete permU2_3;
-  delete permU3_1;
-  delete permU3_2;
+
+  // delete all transformations
+  for(size_t i=0; i<3;++i)
+    delete transformation[0][i];
+  delete[] transformation[0];
+  for(size_t i=0; i<3;++i)
+    delete transformation[1][i];
+  delete[] transformation[1];
+  for(size_t i=0; i<2;++i)
+    delete transformation[2][i];
+  delete[] transformation[2];
+  delete[] transformation;
 }
 
 void SubspaceFactorizerUnittest::BlockTest()
 {
-  MatrixContent temp((*fourbyfour)*(*permU1_1));
+  MatrixContent temp((*fourbyfour)&(*transformation[0][0]));
   std::cout << "Our matrix is " << *fourbyfour << "." << std::endl;
 
-  std::cout << "Multiplying " << *fourbyfour << " by " << *permU1_1 << " is: " << std::endl;
+  std::cout << "Hadamard product of " << *fourbyfour << " with " << *transformation[0][0] << " is: " << std::endl;
   std::cout << temp << std::endl;
 
   gsl_vector *eigenvalues = temp.transformToEigenbasis();
+  VectorContent *eigenvaluesView = new VectorViewContent(gsl_vector_subvector(eigenvalues, 0, eigenvalues->size));
   std::cout << "The resulting eigenbasis is " << temp
-      << " with eigenvalues " << gsl_vector_get(eigenvalues, 0)
-      << std::endl;
+      << "\n\t with eigenvalues " << *eigenvaluesView << std::endl;
+  delete eigenvaluesView;
   gsl_vector_free(eigenvalues);
 
+
   CPPUNIT_ASSERT_EQUAL(0,0);
 }
 
+/** For given set of row and column indices, we extract the small block matrix.
+ * @param bigmatrix big matrix to extract from
+ * @param rowindexset row index set
+ * @param columnindexset column index set
+ * @return small matrix with dimension equal to the number of indices for row and column.
+ */
+MatrixContent * getSubspaceMatrix(
+    MatrixContent &bigmatrix,
+    const IndexSet &rowindexset,
+    const IndexSet &columnindexset)
+{
+  // check whether subsystem is big enough for both index sets
+  ASSERT(rowindexset.size() <= bigmatrix.getRows(),
+      "embedSubspaceMatrix() - bigmatrix has less rows "+toString(bigmatrix.getRows())
+      +" than needed by index set "
+      +toString(rowindexset.size())+"!");
+  ASSERT(columnindexset.size() <= bigmatrix.getColumns(),
+      "embedSubspaceMatrix() - bigmatrix has less columns "+toString(bigmatrix.getColumns())
+      +" than needed by index set "
+      +toString(columnindexset.size())+"!");
+  MatrixContent *subsystem =  new MatrixContent(rowindexset.size(), columnindexset.size());
+  size_t localrow = 0; // local row indices for the subsystem
+  size_t localcolumn = 0;
+  for (IndexSet::const_iterator rowindex = rowindexset.begin();
+      rowindex != rowindexset.end();
+      ++rowindex) {
+    localcolumn = 0;
+    for (IndexSet::const_iterator columnindex = columnindexset.begin();
+        columnindex != columnindexset.end();
+        ++columnindex) {
+      ASSERT((*rowindex < bigmatrix.getRows()) && (*columnindex < bigmatrix.getColumns()),
+          "current index pair ("
+          +toString(*rowindex)+","+toString(*columnindex)
+          +") is out of bounds of bigmatrix ("
+          +toString(bigmatrix.getRows())+","+toString(bigmatrix.getColumns())
+          +")");
+      subsystem->at(localrow,localcolumn) = bigmatrix.at(*rowindex, *columnindex);
+      localcolumn++;
+    }
+    localrow++;
+  }
+  return subsystem;
+}
+
+/** For a given set of row and columns indices, we embed a small block matrix into a bigger space.
+ *
+ * @param bigmatrix big matrix to place submatrix into
+ * @param rowindexset
+ * @param columnindexset
+ * @return
+ */
+void embedSubspaceMatrix(
+    MatrixContent &bigmatrix,
+    MatrixContent &subsystem,
+    const IndexSet &rowindexset,
+    const IndexSet &columnindexset)
+{
+  // check whether bigmatrix is at least as big as subsystem
+  ASSERT(subsystem.getRows() <= bigmatrix.getRows(),
+      "embedSubspaceMatrix() - subsystem has more rows "+toString(subsystem.getRows())
+      +" than bigmatrix "
+      +toString(bigmatrix.getRows())+"!");
+  ASSERT(subsystem.getColumns() <= bigmatrix.getColumns(),
+      "embedSubspaceMatrix() - subsystem has more columns "+toString(subsystem.getColumns())
+      +" than bigmatrix "
+      +toString(bigmatrix.getColumns())+"!");
+  // check whether subsystem is big enough for both index sets
+  ASSERT(rowindexset.size() <= subsystem.getRows(),
+      "embedSubspaceMatrix() - subsystem has less rows "+toString(subsystem.getRows())
+      +" than needed by index set "
+      +toString(rowindexset.size())+"!");
+  ASSERT(columnindexset.size() <= subsystem.getColumns(),
+      "embedSubspaceMatrix() - subsystem has less columns "+toString(subsystem.getColumns())
+      +" than needed by index set "
+      +toString(columnindexset.size())+"!");
+  size_t localrow = 0; // local row indices for the subsystem
+  size_t localcolumn = 0;
+  for (IndexSet::const_iterator rowindex = rowindexset.begin();
+      rowindex != rowindexset.end();
+      ++rowindex) {
+    localcolumn = 0;
+    for (IndexSet::const_iterator columnindex = columnindexset.begin();
+        columnindex != columnindexset.end();
+        ++columnindex) {
+      bigmatrix.at(*rowindex, *columnindex) = subsystem.at(localrow,localcolumn);
+      localcolumn++;
+    }
+    localrow++;
+  }
+}
+
+/** Operator for output to std:: ostream operator of an IndexSet.
+ * @param ost output stream
+ * @param indexset index set to output
+ * @return ost output stream
+ */
+std::ostream & operator<<(std::ostream &ost, const IndexSet &indexset)
+{
+  ost << "{ ";
+  for (IndexSet::const_iterator iter = indexset.begin();
+      iter != indexset.end();
+      ++iter)
+    ost << *iter << " ";
+  ost << "}";
+  return ost;
+}
+
+void SubspaceFactorizerUnittest::EigenvectorTest()
+{
+  VectorArray CurrentEigenvectors;
+  VectorList *ParallelEigenvectorList = new VectorList[4];
+  IndexSet AllIndices;
+
+  // create the total index set
+  for (size_t i=0;i<4;++i)
+    AllIndices.insert(i);
+
+  // create all consecutive index subsets for dim 1 to 3
+  IndexMap Dimension_to_Indexset;
+  for (size_t dim = 0; dim<3;++dim) {
+    for (size_t i=0;i<4-dim;++i) {
+      IndexSet *indexset = new IndexSet;
+      for (size_t j=0;j<=dim;++j) {
+        //std::cout << "Putting " << i+j << " into " << i << "th map " << dim << std::endl;
+        CPPUNIT_ASSERT_MESSAGE("index "+toString(i+j)+" already present in "+toString(dim)+"-dim "+toString(i)+"th indexset.", indexset->count(i+j) == 0);
+        indexset->insert(i+j);
+      }
+      Dimension_to_Indexset.insert( make_pair(dim, boost::shared_ptr<IndexSet>(indexset)) );
+      // no need to free indexset, is stored in shared_ptr and
+      // will get released when Dimension_to_Indexset is destroyed
+    }
+  }
+
+  // set to first guess, i.e. the unit vectors of R^4
+  for (IndexSet::const_iterator index = AllIndices.begin();
+      index != AllIndices.end();
+      ++index) {
+    boost::shared_ptr<VectorContent> EV(new VectorContent(4));
+    EV->setZero();
+    EV->at(*index) = 1.;
+    CurrentEigenvectors.push_back(EV);
+  }
+
+  // for every dimension
+  for (size_t dim = 0; dim<3;++dim) {
+    // for every index set of this dimension
+    Log() << Verbose(0) << std::endl << std::endl;
+    Log() << Verbose(0) << "Current dimension is " << dim << std::endl;
+    std::pair<IndexMap::const_iterator,IndexMap::const_iterator> Bounds = Dimension_to_Indexset.equal_range(dim);
+    for (IndexMap::const_iterator IndexsetIter = Bounds.first;
+        IndexsetIter != Bounds.second;
+        ++IndexsetIter) {
+      // show the index set
+      Log() << Verbose(0) << std::endl;
+      Log() << Verbose(1) << "Current index set is " << *(IndexsetIter->second) << std::endl;
+
+      // create transformation matrices from these
+      MatrixContent *subsystem = getSubspaceMatrix(*fourbyfour, *(IndexsetIter->second), *(IndexsetIter->second));
+      Log() << Verbose(2) << "Subsystem matrix is " << *subsystem << std::endl;
+
+      // solve _small_ systems for eigenvalues
+      VectorContent *Eigenvalues = new VectorContent(subsystem->transformToEigenbasis());
+      Log() << Verbose(2) << "Eigenvector matrix is " << *subsystem << std::endl;
+      Log() << Verbose(2) << "Eigenvalues are " << *Eigenvalues << std::endl;
+      delete Eigenvalues;
+
+      // blow up eigenvectors to 4dim column vector again
+      MatrixContent *Eigenvectors = new MatrixContent(4,4);
+      Eigenvectors->setZero();
+      embedSubspaceMatrix(*Eigenvectors, *subsystem, *(IndexsetIter->second), *(IndexsetIter->second));
+      Log() << Verbose(1) << "4x4 Eigenvector matrix is " << *Eigenvectors << std::endl;
+
+      // we don't need the subsystem anymore
+      delete subsystem;
+
+      // copy the index set, we look for one new eigenvector for each old one
+      IndexSet CorrespondenceList((*IndexsetIter->second));
+
+      // go through all eigenvectors in this subspace
+      for (IndexSet::const_iterator iter = (*IndexsetIter->second).begin();
+          iter != (*IndexsetIter->second).end();
+          ++iter) {
+        // we rather copy the column vector, as the containing matrix is destroyed lateron
+        VectorContent *CurrentEigenvector = new VectorContent(Eigenvectors->getColumnVector(*iter));
+        Log() << Verbose(1) << "Current Eigenvector is " << *CurrentEigenvector << std::endl;
+
+        // recognize eigenvectors parallel to existing ones
+        // (for now settle with the one we are most parallel to)
+        size_t mostparallel_index = 4;
+        double mostparallel_scalarproduct = 0.;
+        for (IndexSet::const_iterator indexiter = CorrespondenceList.begin();
+            indexiter != CorrespondenceList.end();
+            ++indexiter) {
+          Log() << Verbose(2) << "Comparing to old " << *indexiter << "th eigenvector " << *(CurrentEigenvectors[*indexiter]) << std::endl;
+          const double scalarproduct = (*(CurrentEigenvectors[*indexiter])) * (*CurrentEigenvector);
+          Log() << Verbose(2) << "SKP is " << scalarproduct << std::endl;
+          if (fabs(scalarproduct) > mostparallel_scalarproduct) {
+            mostparallel_scalarproduct = fabs(scalarproduct);
+            mostparallel_index = *indexiter;
+          }
+        }
+        if (mostparallel_index != 4) {
+          // put into std::list for later use
+          Log() << Verbose(1) << "Pushing " << *CurrentEigenvector << " into parallel list [" << mostparallel_index << "]" << std::endl;
+          ParallelEigenvectorList[mostparallel_index].push_back(boost::shared_ptr<VectorContent>(CurrentEigenvector));
+          CorrespondenceList.erase(mostparallel_index);
+        }
+        // no need to delete CurrentEigenvector as is taken care of by shared_ptr
+      }
+      // free eigenvectors
+      delete Eigenvectors;
+    }
+  }
+
+  // print list of parallel eigenvectors
+  for (IndexSet::const_iterator index = AllIndices.begin();
+      index != AllIndices.end();
+      ++index) {
+    Log() << Verbose(0) << "Similar to " << *index << "th current eigenvector " << *(CurrentEigenvectors[*index]) << " are:" << std::endl;
+    for (VectorList::const_iterator iter = ParallelEigenvectorList[*index].begin();
+        iter != ParallelEigenvectorList[*index].end();
+        ++iter) {
+      Log() << Verbose(0) << **iter << std::endl;
+    }
+    Log() << Verbose(0) << std::endl;
+  }
+
+  // create new CurrentEigenvectors from averaging parallel ones.
+  for (IndexSet::const_iterator index = AllIndices.begin();
+      index != AllIndices.end();
+      ++index) {
+    for (VectorList::const_iterator iter = ParallelEigenvectorList[*index].begin();
+        iter != ParallelEigenvectorList[*index].end();
+        ++iter) {
+      *CurrentEigenvectors[*index] += **iter;
+    }
+    *CurrentEigenvectors[*index] *= 1./(double)(ParallelEigenvectorList[*index].size()+1);
+  }
+
+  // show new ones
+  Log() << Verbose(0) << "Resulting new eigenvectors are:" << std::endl;
+  for (IndexSet::const_iterator index = AllIndices.begin();
+      index != AllIndices.end();
+      ++index) {
+    Log() << Verbose(0) << *CurrentEigenvectors[*index] << std::endl;
+  }
+
+
+  delete[] ParallelEigenvectorList;
+
+  CPPUNIT_ASSERT_EQUAL(0,0);
+}
+

src/unittests/SubspaceFactorizerUnittest.hpp

@@ -11 +11 @@
 #include <cppunit/extensions/HelperMacros.h>
 
+typedef std::set<size_t> IndexSet;
+typedef std::multimap< size_t, boost::shared_ptr<IndexSet> > IndexMap;
+typedef std::list< boost::shared_ptr<VectorContent> > VectorList;
+typedef std::vector< boost::shared_ptr<VectorContent> > VectorArray;
+
 class MatrixContent;
 
@@ -17 +22 @@
   CPPUNIT_TEST_SUITE( SubspaceFactorizerUnittest) ;
   CPPUNIT_TEST ( BlockTest );
+  CPPUNIT_TEST ( EigenvectorTest );
   CPPUNIT_TEST_SUITE_END();
 
@@ -24 +30 @@
 
   void BlockTest();
+  void EigenvectorTest();
 
   MatrixContent *fourbyfour;
 
-  MatrixContent *permU1_1;
-  MatrixContent *permU1_2;
-  MatrixContent *permU1_3;
-  MatrixContent *permU1_4;
-  MatrixContent *permU2_1;
-  MatrixContent *permU2_2;
-  MatrixContent *permU2_3;
-  MatrixContent *permU3_1;
-  MatrixContent *permU3_2;
+  MatrixContent ***transformation;
 };
 
+MatrixContent * getSubspaceMatrix(MatrixContent &bigmatrix, const IndexSet &rowindexset, const IndexSet &columnindexset);
+void embedSubspaceMatrix(MatrixContent &bigmatrix, MatrixContent &subsystem, const IndexSet &rowindexset, const IndexSet &columnindexset);
+std::ostream & operator<<(std::ostream &ost, const IndexSet &indexset);
+
 #endif /* SUBSPACEFACTORIZERUNITTEST_HPP_ */