Changeset 8450da for src/Dynamics


Ignore:
Timestamp:
Apr 10, 2018, 6:43:31 AM (7 years ago)
Author:
Frederik Heber <frederik.heber@…>
Branches:
AutomationFragmentation_failures, Candidate_v1.6.1, ChemicalSpaceEvaluator, Exclude_Hydrogens_annealWithBondGraph, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_contraction-expansion, Gui_displays_atomic_force_velocity, PythonUI_with_named_parameters, StoppableMakroAction, TremoloParser_IncreasedPrecision
Children:
6a5921, 897a01
Parents:
90ece9
git-author:
Frederik Heber <frederik.heber@…> (08/03/17 09:24:07)
git-committer:
Frederik Heber <frederik.heber@…> (04/10/18 06:43:31)
Message:

ForceAnnealing functions now have better readable time step variables.

  • _TimeStep is the parameter, Old.. and CurrentTimeStep designate the two reference time steps for calculating step width.
  • Split functions into simple step width and using BB step width.
  • TESTFIX: Removed XFAIL from all Python ForceAnnealing tests.
  • TESTFIX: 5-body annealing without bond graph has slightly changed but quality of results is the same.
File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/Dynamics/ForceAnnealing.hpp

    r90ece9 r8450da  
    9292   *
    9393   *
    94    * \param CurrentTimeStep current time step (i.e. \f$ t + \Delta t \f$ in the sense of the velocity verlet)
     94   * \param _TimeStep time step to update (i.e. \f$ t + \Delta t \f$ in the sense of the velocity verlet)
    9595   * \param offset offset in matrix file to the first force component
    9696   * \todo This is not yet checked if it is correctly working with DoConstrainedMD set >0.
    9797   */
    9898  void operator()(
    99       const int _CurrentTimeStep,
     99      const int _TimeStep,
    100100      const size_t _offset,
    101101      const bool _UseBondgraph)
    102102  {
     103    const int CurrentTimeStep = _TimeStep-1;
     104    ASSERT( CurrentTimeStep >= 0,
     105        "ForceAnnealing::operator() - a new time step (upon which we work) must already have been copied.");
     106
    103107    // make sum of forces equal zero
    104108    AtomicForceManipulator<T>::correctForceMatrixForFixedCenterOfMass(
    105109        _offset,
    106         _CurrentTimeStep-1>=0 ? _CurrentTimeStep - 1 : 0);
     110        CurrentTimeStep);
    107111
    108112    // are we in initial step? Then set static entities
    … …  
    114118
    115119      // always use atomic annealing on first step
    116       maxComponents = anneal(_CurrentTimeStep);
     120      maxComponents = anneal(_TimeStep);
    117121    } else {
    118122      ++currentStep;
    … …  
    121125      // bond graph annealing is always followed by a normal annealing
    122126      if (_UseBondgraph)
    123         maxComponents = annealWithBondGraph(_CurrentTimeStep);
     127        maxComponents = annealWithBondGraph_BarzilaiBorwein(_TimeStep);
    124128      // cannot store RemnantGradient in Atom's Force as it ruins BB stepwidth calculation
    125129      else
    126         maxComponents = anneal(_CurrentTimeStep);
     130        maxComponents = anneal_BarzilaiBorwein(_TimeStep);
    127131    }
    128132
    … …  
    163167   * We assume that forces have just been calculated.
    164168   *
    165    * \param CurrentTimeStep current time step (i.e. \f$ t + \Delta t \f$ in the sense of the velocity verlet)
     169   * \param _TimeStep time step to update (i.e. \f$ t + \Delta t \f$ in the sense of the velocity verlet)
    166170   * \return to be filled with maximum force component over all atoms
    167171   */
    168172  Vector anneal(
    169       const int CurrentTimeStep)
     173      const int _TimeStep)
    170174  {
     175    const int CurrentTimeStep = _TimeStep-1;
     176    ASSERT( CurrentTimeStep >= 0,
     177        "ForceAnnealing::anneal() - a new time step (upon which we work) must already have been copied.");
     178
     179    LOG(1, "STATUS: performing simple anneal with default stepwidth " << currentDeltat << " at step #" << currentStep);
     180
    171181    Vector maxComponents;
    172182    bool deltat_decreased = false;
    … …  
    174184        iter != AtomicForceManipulator<T>::atoms.end(); ++iter) {
    175185      // atom's force vector gives steepest descent direction
    176       const Vector oldPosition = (*iter)->getPositionAtStep(CurrentTimeStep-1 >= 0 ? CurrentTimeStep - 1 : 0);
    177186      const Vector currentPosition = (*iter)->getPositionAtStep(CurrentTimeStep);
    178       const Vector oldGradient = (*iter)->getAtomicForceAtStep(CurrentTimeStep-1 >= 0 ? CurrentTimeStep - 1 : 0);
    179187      const Vector currentGradient = (*iter)->getAtomicForceAtStep(CurrentTimeStep);
    180       LOG(4, "DEBUG: oldPosition for atom " << **iter << " is " << oldPosition);
    181       LOG(4, "DEBUG: currentPosition for atom " << **iter << " is " << currentPosition);
    182       LOG(4, "DEBUG: oldGradient for atom " << **iter << " is " << oldGradient);
    183       LOG(4, "DEBUG: currentGradient for atom " << **iter << " is " << currentGradient);
     188      LOG(4, "DEBUG: currentPosition for atom #" << (*iter)->getId() << " is " << currentPosition);
     189      LOG(4, "DEBUG: currentGradient for atom #" << (*iter)->getId() << " is " << currentGradient);
    184190//      LOG(4, "DEBUG: Force for atom " << **iter << " is " << currentGradient);
    185191
    186192      // we use Barzilai-Borwein update with position reversed to get descent
    187       const double stepwidth = getBarzilaiBorweinStepwidth(
    188           currentPosition - oldPosition, currentGradient - oldGradient);
     193      double stepwidth = currentDeltat;
    189194      Vector PositionUpdate = stepwidth * currentGradient;
    190195      LOG(3, "DEBUG: Update would be " << stepwidth << "*" << currentGradient << " = " << PositionUpdate);
    … …  
    197202      // have different sign: Check whether this is the case and one step with
    198203      // deltat to interrupt this sequence
    199       const Vector PositionDifference = currentPosition - oldPosition;
    200       if ((currentStep > 1) && (!PositionDifference.IsZero()))
     204      if (currentStep > 1) {
     205        const int OldTimeStep = CurrentTimeStep-1;
     206        ASSERT( OldTimeStep >= 0,
     207            "ForceAnnealing::anneal() - if currentStep is "+toString(currentStep)
     208            +", then there should be at least three time steps.");
     209        const Vector oldPosition = (*iter)->getPositionAtStep(OldTimeStep);
     210        const Vector PositionDifference = currentPosition - oldPosition;
     211        LOG(4, "DEBUG: oldPosition for atom #" << (*iter)->getId() << " is " << oldPosition);
     212        LOG(4, "DEBUG: PositionDifference for atom #" << (*iter)->getId() << " is " << PositionDifference);
    201213        if ((PositionUpdate.ScalarProduct(PositionDifference) < 0)
    202214            && (fabs(PositionUpdate.NormSquared()-PositionDifference.NormSquared()) < 1e-3)) {
    … …  
    210222              << ", using deltat: " << currentDeltat);
    211223          PositionUpdate = currentDeltat * currentGradient;
     224        }
    212225      }
    213226
    214227      // finally set new values
    215       (*iter)->setPosition(currentPosition + PositionUpdate);
     228      (*iter)->setPositionAtStep(_TimeStep, currentPosition + PositionUpdate);
     229    }
     230
     231    return maxComponents;
     232  }
     233
     234  /** Performs Gradient optimization on the atoms using BarzilaiBorwein step width.
     235   *
     236   * \note this can only be called when there are at least two optimization
     237   * time steps present, i.e. this must be preceeded by a simple anneal().
     238   *
     239   * We assume that forces have just been calculated.
     240   *
     241   * \param _TimeStep time step to update (i.e. \f$ t + \Delta t \f$ in the sense of the velocity verlet)
     242   * \return to be filled with maximum force component over all atoms
     243   */
     244  Vector anneal_BarzilaiBorwein(
     245      const int _TimeStep)
     246  {
     247    const int OldTimeStep = _TimeStep-2;
     248    const int CurrentTimeStep = _TimeStep-1;
     249    ASSERT( OldTimeStep >= 0,
     250        "ForceAnnealing::anneal_BarzilaiBorwein() - we need two present optimization steps to compute stepwidth.");
     251    ASSERT(currentStep > 1,
     252        "ForceAnnealing::anneal_BarzilaiBorwein() - we need two present optimization steps to compute stepwidth.");
     253
     254    LOG(1, "STATUS: performing BarzilaiBorwein anneal at step #" << currentStep);
     255
     256    Vector maxComponents;
     257    bool deltat_decreased = false;
     258    for(typename AtomSetMixin<T>::iterator iter = AtomicForceManipulator<T>::atoms.begin();
     259        iter != AtomicForceManipulator<T>::atoms.end(); ++iter) {
     260      // atom's force vector gives steepest descent direction
     261      const Vector oldPosition = (*iter)->getPositionAtStep(OldTimeStep);
     262      const Vector currentPosition = (*iter)->getPositionAtStep(CurrentTimeStep);
     263      const Vector oldGradient = (*iter)->getAtomicForceAtStep(OldTimeStep);
     264      const Vector currentGradient = (*iter)->getAtomicForceAtStep(CurrentTimeStep);
     265      LOG(4, "DEBUG: oldPosition for atom #" << (*iter)->getId() << " is " << oldPosition);
     266      LOG(4, "DEBUG: currentPosition for atom #" << (*iter)->getId() << " is " << currentPosition);
     267      LOG(4, "DEBUG: oldGradient for atom #" << (*iter)->getId() << " is " << oldGradient);
     268      LOG(4, "DEBUG: currentGradient for atom #" << (*iter)->getId() << " is " << currentGradient);
     269//      LOG(4, "DEBUG: Force for atom #" << (*iter)->getId() << " is " << currentGradient);
     270
     271      // we use Barzilai-Borwein update with position reversed to get descent
     272      const Vector PositionDifference = currentPosition - oldPosition;
     273      const Vector GradientDifference = (currentGradient - oldGradient);
     274      double stepwidth = getBarzilaiBorweinStepwidth(PositionDifference, GradientDifference);
     275      Vector PositionUpdate = stepwidth * currentGradient;
     276      LOG(3, "DEBUG: Update would be " << stepwidth << "*" << currentGradient << " = " << PositionUpdate);
     277
     278      // extract largest components for showing progress of annealing
     279      for(size_t i=0;i<NDIM;++i)
     280        maxComponents[i] = std::max(maxComponents[i], fabs(currentGradient[i]));
     281
     282//      // steps may go back and forth again (updates are of same magnitude but
     283//      // have different sign: Check whether this is the case and one step with
     284//      // deltat to interrupt this sequence
     285//      if (!PositionDifference.IsZero())
     286//        if ((PositionUpdate.ScalarProduct(PositionDifference) < 0)
     287//            && (fabs(PositionUpdate.NormSquared()-PositionDifference.NormSquared()) < 1e-3)) {
     288//          // for convergence we want a null sequence here, too
     289//          if (!deltat_decreased) {
     290//            deltat_decreased = true;
     291//            currentDeltat = .5*currentDeltat;
     292//          }
     293//          LOG(2, "DEBUG: Upgrade in other direction: " << PositionUpdate
     294//              << " > " << PositionDifference
     295//              << ", using deltat: " << currentDeltat);
     296//          PositionUpdate = currentDeltat * currentGradient;
     297//      }
     298
     299      // finally set new values
     300      (*iter)->setPositionAtStep(_TimeStep, currentPosition + PositionUpdate);
    216301    }
    217302
    … …  
    256341  }
    257342
    258   /** Performs Gradient optimization on the bonds.
     343  /** Performs Gradient optimization on the bonds with BarzilaiBorwein stepwdith.
     344   *
     345   * \note this can only be called when there are at least two optimization
     346   * time steps present, i.e. this must be preceeded by a simple anneal().
    259347   *
    260348   * We assume that forces have just been calculated. These forces are projected
    … …  
    263351   *
    264352   *
    265    * \param CurrentTimeStep current time step (i.e. \f$ t + \Delta t \f$ in the sense of the velocity verlet)
     353   * \param _TimeStep time step to update (i.e. \f$ t + \Delta t \f$ in the sense of the velocity verlet)
    266354   * \param maxComponents to be filled with maximum force component over all atoms
    267355   */
    268   Vector annealWithBondGraph(
    269       const int CurrentTimeStep)
     356  Vector annealWithBondGraph_BarzilaiBorwein(
     357      const int _TimeStep)
    270358  {
     359    const int OldTimeStep = _TimeStep-2;
     360    const int CurrentTimeStep = _TimeStep-1;
     361    ASSERT(OldTimeStep >= 0,
     362        "annealWithBondGraph_BarzilaiBorwein() - we need two present optimization steps to compute stepwidth, and the new one to update on already present.");
     363    ASSERT(currentStep > 1,
     364        "annealWithBondGraph_BarzilaiBorwein() - we need two present optimization steps to compute stepwidth.");
     365
     366    LOG(1, "STATUS: performing BarzilaiBorwein anneal on bonds at step #" << currentStep);
     367
    271368    Vector maxComponents;
    272369
    … …  
    309406        AtomicForceManipulator<T>::atoms.begin(),
    310407        AtomicForceManipulator<T>::atoms.end(),
    311         CurrentTimeStep);
     408        _TimeStep); // use time step to update here as this is the current set of bonds
    312409    const BondVectors::container_t &sorted_bonds = bv.getSorted();
    313410
    … …  
    328425    RemnantGradient_per_atom_t RemnantGradient_per_atom;
    329426    for (size_t timestep = 0; timestep <= 1; ++timestep) {
    330       const size_t CurrentStep = CurrentTimeStep-timestep-1 >= 0 ? CurrentTimeStep-timestep-1 : 0;
    331       LOG(2, "DEBUG: CurrentTimeStep is " << CurrentTimeStep
     427      const size_t ReferenceTimeStep = CurrentTimeStep-timestep;
     428      LOG(2, "DEBUG: given time step is " << _TimeStep
    332429          << ", timestep is " << timestep
    333           << ", and CurrentStep is " << CurrentStep);
     430          << ", and ReferenceTimeStep is " << ReferenceTimeStep);
    334431
    335432      for(typename AtomSetMixin<T>::const_iterator iter = AtomicForceManipulator<T>::atoms.begin();
    336433          iter != AtomicForceManipulator<T>::atoms.end(); ++iter) {
    337434        const atom &walker = *(*iter);
    338         const Vector &walkerGradient = walker.getAtomicForceAtStep(CurrentStep);
     435        const Vector &walkerGradient = walker.getAtomicForceAtStep(ReferenceTimeStep);
    339436        LOG(3, "DEBUG: Gradient of atom #" << walker.getId() << ", namely "
    340437            << walker << " is " << walkerGradient << " with magnitude of "
    … …  
    346443          // gather subset of BondVectors for the current atom
    347444          const std::vector<Vector> BondVectors =
    348               bv.getAtomsBondVectorsAtStep(walker, CurrentStep);
     445              bv.getAtomsBondVectorsAtStep(walker, ReferenceTimeStep);
    349446
    350447          // go through all its bonds and calculate what magnitude is represented
    … …  
    353450              weights_per_atom[timestep].insert(
    354451                  std::make_pair(walker.getId(),
    355                   bv.getWeightsForAtomAtStep(walker, BondVectors, CurrentStep)) );
     452                  bv.getWeightsForAtomAtStep(walker, BondVectors, ReferenceTimeStep)) );
    356453          ASSERT( inserter.second,
    357454              "ForceAnnealing::operator() - weight map for atom "+toString(walker)
    … …  
    441538        LOG(4, "DEBUG: current projected gradient for "
    442539            << (side == leftside ? "left" : "right") << " side of bond is " << currentGradient);
    443         const Vector &oldPosition = bondatom[side]->getPositionAtStep(CurrentTimeStep-2 >= 0 ? CurrentTimeStep - 2 : 0);
    444         const Vector &currentPosition = bondatom[side]->getPositionAtStep(CurrentTimeStep-1>=0 ? CurrentTimeStep - 1 : 0);
     540        const Vector &oldPosition = bondatom[side]->getPositionAtStep(OldTimeStep);
     541        const Vector &currentPosition = bondatom[side]->getPositionAtStep(CurrentTimeStep);
    445542        const Vector PositionDifference = currentPosition - oldPosition;
    446543        LOG(4, "DEBUG: old position is " << oldPosition);
    … …  
    491588      atom &walker = *(*iter);
    492589      // extract largest components for showing progress of annealing
    493       const Vector currentGradient = walker.getAtomicForceAtStep(CurrentTimeStep-1>=0 ? CurrentTimeStep-1 : 0);
     590      const Vector currentGradient = walker.getAtomicForceAtStep(CurrentTimeStep);
    494591      for(size_t i=0;i<NDIM;++i)
    495592        maxComponents[i] = std::max(maxComponents[i], fabs(currentGradient[i]));
    496 
    497       // reset force vector for next step except on final one
    498       if (currentStep != maxSteps)
    499         walker.setAtomicForce(zeroVec);
    500593    }
    501594
    … …  
    521614      LOG(3, "DEBUG: Applying update " << update << " to atom #" << atomid
    522615          << ", namely " << *walker);
    523       walker->setPosition(
    524           walker->getPositionAtStep(CurrentTimeStep-1>=0 ? CurrentTimeStep - 1 : 0)
     616      walker->setPositionAtStep(_TimeStep,
     617          walker->getPositionAtStep(CurrentTimeStep)
    525618          + update - CommonTranslation);
    526       walker->setAtomicVelocity(update);
    527619//      walker->setAtomicForce( RemnantGradient_per_atom[walker->getId()] );
    528620    }
Note: See TracChangeset for help on using the changeset viewer.