Context Navigation

-              rfe6f20
+              re355762
+MinimiseConstrainedPotential::MinimiseConstrainedPotential(molecule::atomSet &_atoms) :
+  atoms(_atoms)
+MinimiseConstrainedPotential::MinimiseConstrainedPotential(
+    molecule::atomSet &_atoms,
+    std::map<atom*, atom *> &_PermutationMap) :
+  atoms(_atoms),
+  PermutationMap(_PermutationMap)
 {}
 …
 {}
 double MinimiseConstrainedPotential::operator()(atom **&PermutationMap, int _startstep, int _endstep, bool IsAngstroem)
+double MinimiseConstrainedPotential::operator()(int _startstep, int _endstep, bool IsAngstroem)
+{
   double Potential, OldPotential, OlderPotential;
-  PermutationMap = new atom *[atoms.size()];
-  DistanceList = new DistanceMap *[atoms.size()];
-  DistanceIterators = new DistanceMap::iterator[atoms.size()];
-  DoubleList = new int[atoms.size()];
-  StepList = new DistanceMap::iterator[atoms.size()];
   int round;
   atom *Sprinter = NULL;
 …
   // set to zero
+  for (size_t i=0;i<atoms.size();i++) {
+    PermutationMap[i] = NULL;
+    DoubleList[i] = 0;
+  }
+  PermutationMap.clear();
+  DoubleList.clear();
+  for (molecule::atomSet::const_iterator iter = atoms.begin(); iter != atoms.end(); ++iter) {
+    DistanceList[*iter].clear();
+  }
+  DistanceList.clear();
+  DistanceIterators.clear();
+  DistanceIterators.clear();
   /// Minimise the potential
 …
   DoLog(1) && (Log() << Verbose(1) << "Making the PermutationMap injective ... " << endl);
   MakeInjectivePermutation();
   delete[](DoubleList);
+  DoubleList.clear();
   // argument minimise the constrained potential in this injective PermutationMap
 …
       iter = atoms.begin();
       for (; iter != atoms.end(); ++iter) {
+        CalculateDoubleList();
         PrintPermutationMap();
         Sprinter = DistanceIterators[(*iter)->getNr()]->second;   // store initial partner
         Strider = DistanceIterators[(*iter)->getNr()];  //remember old iterator
         DistanceIterators[(*iter)->getNr()] = StepList[(*iter)->getNr()];
         if (DistanceIterators[(*iter)->getNr()] == DistanceList[(*iter)->getNr()]->end()) {// stop, before we run through the list and still on
           DistanceIterators[(*iter)->getNr()] == DistanceList[(*iter)->getNr()]->begin();
+        Sprinter = DistanceIterators[(*iter)]->second;   // store initial partner
+        Strider = DistanceIterators[(*iter)];  //remember old iterator
+        DistanceIterators[(*iter)] = StepList[(*iter)];
+        if (DistanceIterators[(*iter)] == DistanceList[(*iter)].end()) {// stop, before we run through the list and still on
+          DistanceIterators[(*iter)] == DistanceList[(*iter)].begin();
           break;
+        }
         //Log() << Verbose(2) << "Current Walker: " << *(*iter) << " with old/next candidate " << *Sprinter << "/" << *DistanceIterators[(*iter)->getNr()]->second << "." << endl;
+        //Log() << Verbose(2) << "Current Walker: " << *(*iter) << " with old/next candidate " << *Sprinter << "/" << *DistanceIterators[(*iter)]->second << "." << endl;
         // find source of the new target
         molecule::atomSet::const_iterator runner = atoms.begin();
         for (; runner != atoms.end(); ++runner) { // find the source whose toes we might be stepping on (Walker's new target should be in use by another already)
           if (PermutationMap[(*runner)->getNr()] == DistanceIterators[(*iter)->getNr()]->second) {
             //Log() << Verbose(2) << "Found the corresponding owner " << *(*runner) << " to " << *PermutationMap[(*runner)->getNr()] << "." << endl;
+          if (PermutationMap[(*runner)] == DistanceIterators[(*iter)]->second) {
+            //Log() << Verbose(2) << "Found the corresponding owner " << *(*runner) << " to " << *PermutationMap[(*runner)] << "." << endl;
             break;
+          }
 …
         if (runner != atoms.end()) { // we found the other source
           // then look in its distance list for Sprinter
           Rider = DistanceList[(*runner)->getNr()]->begin();
           for (; Rider != DistanceList[(*runner)->getNr()]->end(); Rider++)
+          Rider = DistanceList[(*runner)].begin();
+          for (; Rider != DistanceList[(*runner)].end(); Rider++)
             if (Rider->second == Sprinter)
               break;
           if (Rider != DistanceList[(*runner)->getNr()]->end()) { // if we have found one
             //Log() << Verbose(2) << "Current Other: " << *(*runner) << " with old/next candidate " << *PermutationMap[(*runner)->getNr()] << "/" << *Rider->second << "." << endl;
+          if (Rider != DistanceList[(*runner)].end()) { // if we have found one
+            //Log() << Verbose(2) << "Current Other: " << *(*runner) << " with old/next candidate " << *PermutationMap[(*runner)] << "/" << *Rider->second << "." << endl;
             // exchange both
+            PermutationMap[(*iter)->getNr()] = DistanceIterators[(*iter)->getNr()]->second; // put next farther distance into PermutationMap
+            PermutationMap[(*runner)->getNr()] = Sprinter;  // and hand the old target to its respective owner
+            PermutationMap[(*iter)] = DistanceIterators[(*iter)]->second; // put next farther distance into PermutationMap
+            PermutationMap[(*runner)] = Sprinter;  // and hand the old target to its respective owner
+            CalculateDoubleList();
             PrintPermutationMap();
             // calculate the new potential
 …
             if (Potential > OldPotential) { // we made everything worse! Undo ...
               //Log() << Verbose(3) << "Nay, made the potential worse: " << Potential << " vs. " << OldPotential << "!" << endl;
               //Log() << Verbose(3) << "Setting " << *(*runner) << "'s source to " << *DistanceIterators[(*runner)->getNr()]->second << "." << endl;
+              //Log() << Verbose(3) << "Setting " << *(*runner) << "'s source to " << *DistanceIterators[(*runner)]->second << "." << endl;
               // Undo for Runner (note, we haven't moved the iteration yet, we may use this)
               PermutationMap[(*runner)->getNr()] = DistanceIterators[(*runner)->getNr()]->second;
+              PermutationMap[(*runner)] = DistanceIterators[(*runner)]->second;
               // Undo for Walker
               DistanceIterators[(*iter)->getNr()] = Strider;  // take next farther distance target
               //Log() << Verbose(3) << "Setting " << *(*iter) << "'s source to " << *DistanceIterators[(*iter)->getNr()]->second << "." << endl;
               PermutationMap[(*iter)->getNr()] = DistanceIterators[(*iter)->getNr()]->second;
+              DistanceIterators[(*iter)] = Strider;  // take next farther distance target
+              //Log() << Verbose(3) << "Setting " << *(*iter) << "'s source to " << *DistanceIterators[(*iter)]->second << "." << endl;
+              PermutationMap[(*iter)] = DistanceIterators[(*iter)]->second;
             } else {
               DistanceIterators[(*runner)->getNr()] = Rider;  // if successful also move the pointer in the iterator list
+              DistanceIterators[(*runner)] = Rider;  // if successful also move the pointer in the iterator list
               DoLog(3) && (Log() << Verbose(3) << "Found a better permutation, new potential is " << Potential << " vs." << OldPotential << "." << endl);
               OldPotential = Potential;
 …
+          }
         } else {
           PermutationMap[(*iter)->getNr()] = DistanceIterators[(*iter)->getNr()]->second; // new target has no source!
+          PermutationMap[(*iter)] = DistanceIterators[(*iter)]->second; // new target has no source!
+        }
         StepList[(*iter)->getNr()]++; // take next farther distance target
+        StepList[(*iter)]++; // take next farther distance target
+      }
     } while (++iter != atoms.end());
 …
-  /// free memory and return with evaluated potential
-  for (int i=atoms.size(); i--;)
-    DistanceList[i]->clear();
-  delete[](DistanceList);
-  delete[](DistanceIterators);
   return ConstrainedPotential();
 };
 …
 void MinimiseConstrainedPotential::FillDistanceList()
+{
-  for (int i=atoms.size(); i--;) {
-    DistanceList[i] = new DistanceMap;    // is the distance sorted target list per atom
-    DistanceList[i]->clear();
+  }
   for (molecule::atomSet::const_iterator iter = atoms.begin(); iter != atoms.end(); ++iter) {
     for (molecule::atomSet::const_iterator runner = atoms.begin(); runner != atoms.end(); ++runner) {
       DistanceList[(*iter)->getNr()]->insert( DistancePair((*iter)->getPositionAtStep(startstep).distance((*runner)->getPositionAtStep(endstep)), (*runner)) );
+      DistanceList[(*iter)].insert( DistancePair((*iter)->getPositionAtStep(startstep).distance((*runner)->getPositionAtStep(endstep)), (*runner)) );
+    }
+  }
 …
+{
   for (molecule::atomSet::const_iterator iter = atoms.begin(); iter != atoms.end(); ++iter) {
     StepList[(*iter)->getNr()] = DistanceList[(*iter)->getNr()]->begin();    // stores the step to the next iterator that could be a possible next target
     PermutationMap[(*iter)->getNr()] = DistanceList[(*iter)->getNr()]->begin()->second;   // always pick target with the smallest distance
     DoubleList[DistanceList[(*iter)->getNr()]->begin()->second->getNr()]++;            // increase this target's source count (>1? not injective)
     DistanceIterators[(*iter)->getNr()] = DistanceList[(*iter)->getNr()]->begin();    // and remember which one we picked
     DoLog(2) && (Log() << Verbose(2) << **iter << " starts with distance " << DistanceList[(*iter)->getNr()]->begin()->first << "." << endl);
+    StepList[(*iter)] = DistanceList[(*iter)].begin();    // stores the step to the next iterator that could be a possible next target
+    PermutationMap[(*iter)] = DistanceList[(*iter)].begin()->second;   // always pick target with the smallest distance
+    DoubleList[DistanceList[(*iter)].begin()->second]++;            // increase this target's source count (>1? not injective)
+    DistanceIterators[(*iter)] = DistanceList[(*iter)].begin();    // and remember which one we picked
+    DoLog(2) && (Log() << Verbose(2) << **iter << " starts with distance " << DistanceList[(*iter)].begin()->first << "." << endl);
+  }
 };
 …
+  }
   while ((Potential) > PenaltyConstants[2]) {
+    CalculateDoubleList();
     PrintPermutationMap();
     iter++;
     if (iter == atoms.end()) // round-robin at the end
       iter = atoms.begin();
     if (DoubleList[DistanceIterators[(*iter)->getNr()]->second->getNr()] <= 1)  // no need to make those injective that aren't
+    if (DoubleList[DistanceIterators[(*iter)]->second] <= 1)  // no need to make those injective that aren't
       continue;
     // now, try finding a new one
     Potential = TryNextNearestNeighbourForInjectivePermutation((*iter), Potential);
+  }
+  for (int i=atoms.size(); i--;) // now each single entry in the DoubleList should be <=1
+    if (DoubleList[i] > 1) {
+  for (molecule::atomSet::const_iterator iter = atoms.begin(); iter != atoms.end(); ++iter) {
+    // now each single entry in the DoubleList should be <=1
+    if (DoubleList[*iter] > 1) {
       DoeLog(0) && (eLog()<< Verbose(0) << "Failed to create an injective PermutationMap!" << endl);
       performCriticalExit();
+    }
+  }
   DoLog(1) && (Log() << Verbose(1) << "done." << endl);
 };
+unsigned int MinimiseConstrainedPotential::CalculateDoubleList()
+{
+  for (molecule::atomSet::const_iterator iter = atoms.begin(); iter != atoms.end(); ++iter)
+    DoubleList[*iter] = 0;
+  unsigned int doubles = 0;
+  for (molecule::atomSet::const_iterator iter = atoms.begin(); iter != atoms.end(); ++iter)
+      DoubleList[ PermutationMap[*iter] ]++;
+  for (molecule::atomSet::const_iterator iter = atoms.begin(); iter != atoms.end(); ++iter)
+    if (DoubleList[*iter] > 1)
+      doubles++;
+  if (doubles >0)
+    DoLog(2) && (Log() << Verbose(2) << "Found " << doubles << " Doubles." << endl);
+  return doubles;
+};
 void MinimiseConstrainedPotential::PrintPermutationMap() const
+{
   stringstream zeile1, zeile2;
-  const int AtomCount = atoms.size();
-  int *DoubleList = new int[AtomCount];
-  for(int i=0;i<AtomCount;i++)
-    DoubleList[i] = 0;
   int doubles = 0;
   zeile1 << "PermutationMap: ";
   zeile2 << "                ";
+  for (int i=0;i<AtomCount;i++) {
+    DoubleList[PermutationMap[i]->getNr()]++;
+    zeile1 << i << " ";
+    zeile2 << PermutationMap[i]->getNr() << " ";
+  }
+  for (int i=0;i<AtomCount;i++)
+    if (DoubleList[i] > 1)
+    doubles++;
+  for (molecule::atomSet::const_iterator iter = atoms.begin(); iter != atoms.end(); ++iter) {
+    zeile1 << (*iter)->getName() << " ";
+    zeile2 << (PermutationMap[*iter])->getName() << " ";
+  }
+  for (molecule::atomSet::const_iterator iter = atoms.begin(); iter != atoms.end(); ++iter) {
+    std::map<atom *, unsigned int>::const_iterator value_iter = DoubleList.find(*iter);
+    if (value_iter->second > (unsigned int)1)
+      doubles++;
+  }
   if (doubles >0)
     DoLog(2) && (Log() << Verbose(2) << "Found " << doubles << " Doubles." << endl);
-  delete[](DoubleList);
 //  Log() << Verbose(2) << zeile1.str() << endl << zeile2.str() << endl;
 };
 …
   for (molecule::atomSet::const_iterator iter = atoms.begin(); iter != atoms.end(); ++iter) {
     // first term: distance to target
     Runner = PermutationMap[(*iter)->getNr()];   // find target point
+    Runner = PermutationMap[(*iter)];   // find target point
     tmp = ((*iter)->getPositionAtStep(startstep).distance(Runner->getPositionAtStep(endstep)));
     tmp *= IsAngstroem ? 1. : 1./AtomicLengthToAngstroem;
 …
+{
   double Potential = 0;
   DistanceMap::iterator NewBase = DistanceIterators[Walker->getNr()];  // store old base
+  DistanceMap::iterator NewBase = DistanceIterators[Walker];  // store old base
   do {
     NewBase++;  // take next further distance in distance to targets list that's a target of no one
   } while ((DoubleList[NewBase->second->getNr()] != 0) && (NewBase != DistanceList[Walker->getNr()]->end()));
   if (NewBase != DistanceList[Walker->getNr()]->end()) {
     PermutationMap[Walker->getNr()] = NewBase->second;
+  } while ((DoubleList[NewBase->second] != 0) && (NewBase != DistanceList[Walker].end()));
+  if (NewBase != DistanceList[Walker].end()) {
+    PermutationMap[Walker] = NewBase->second;
     Potential = fabs(ConstrainedPotential());
     if (Potential > OldPotential) { // undo
       PermutationMap[Walker->getNr()] = DistanceIterators[Walker->getNr()]->second;
+      PermutationMap[Walker] = DistanceIterators[Walker]->second;
     } else {  // do
       DoubleList[DistanceIterators[Walker->getNr()]->second->getNr()]--;  // decrease the old entry in the doubles list
       DoubleList[NewBase->second->getNr()]++;    // increase the old entry in the doubles list
       DistanceIterators[Walker->getNr()] = NewBase;
+      DoubleList[DistanceIterators[Walker]->second]--;  // decrease the old entry in the doubles list
+      DoubleList[NewBase->second]++;    // increase the old entry in the doubles list
+      DistanceIterators[Walker] = NewBase;
       OldPotential = Potential;
       DoLog(3) && (Log() << Verbose(3) << "Found a new permutation, new potential is " << OldPotential << "." << endl);
 …
   double result = 0.;
   for (molecule::atomSet::const_iterator iter = atoms.begin(); iter != atoms.end(); ++iter) {
     if ((PermutationMap[Walker->getNr()] == PermutationMap[(*iter)->getNr()]) && (Walker->getNr() < (*iter)->getNr())) {
+    if ((PermutationMap[Walker] == PermutationMap[(*iter)]) && (Walker < (*iter))) {
   //    atom *Sprinter = PermutationMap[Walker->nr];
   //        Log() << Verbose(0) << *Walker << " and " << *(*iter) << " are heading to the same target at ";
 …
       break;
     // determine normalized trajectories direction vector (n1, n2)
     Sprinter = PermutationMap[Walker->getNr()];   // find first target point
+    Sprinter = PermutationMap[Walker];   // find first target point
     trajectory1 = Sprinter->getPositionAtStep(endstep) - Walker->getPositionAtStep(startstep);
     trajectory1.Normalize();
     Norm1 = trajectory1.Norm();
     Sprinter = PermutationMap[(*iter)->getNr()];   // find second target point
+    Sprinter = PermutationMap[(*iter)];   // find second target point
     trajectory2 = Sprinter->getPositionAtStep(endstep) - (*iter)->getPositionAtStep(startstep);
     trajectory2.Normalize();
 …
       tmp = Walker->getPositionAtStep(startstep).distance((*iter)->getPositionAtStep(startstep));
     } else if (Norm1 < MYEPSILON) {
       Sprinter = PermutationMap[Walker->getNr()];   // find first target point
+      Sprinter = PermutationMap[Walker];   // find first target point
       trajectory1 = Sprinter->getPositionAtStep(endstep) - (*iter)->getPositionAtStep(startstep);
       trajectory2 *= trajectory1.ScalarProduct(trajectory2); // trajectory2 is scaled to unity, hence we don't need to divide by anything
 …
       tmp = trajectory1.Norm();  // remaining norm is distance
     } else if (Norm2 < MYEPSILON) {
       Sprinter = PermutationMap[(*iter)->getNr()];   // find second target point
+      Sprinter = PermutationMap[(*iter)];   // find second target point
       trajectory2 = Sprinter->getPositionAtStep(endstep) - Walker->getPositionAtStep(startstep);  // copy second offset
       trajectory1 *= trajectory2.ScalarProduct(trajectory1); // trajectory1 is scaled to unity, hence we don't need to divide by anything
 …
   DoLog(1) && (Log() << Verbose(1) << "Calculating forces and adding onto ForceMatrix ... " << endl);
   for(molecule::atomSet::const_iterator iter = atoms.begin(); iter != atoms.end(); ++iter) {
     atom *Sprinter = PermutationMap[(*iter)->getNr()];
+    atom *Sprinter = PermutationMap[(*iter)];
     // set forces
     for (int i=NDIM;i++;)

Note: See TracChangeset for help on using the changeset viewer.

Context Navigation

Changeset e355762 for src/Dynamics/MinimiseConstrainedPotential.cpp

Legend:

src/Dynamics/MinimiseConstrainedPotential.cpp

Download in other formats: