Loading...
Searching...
No Matches
eonc::ConjugateGradients Class Reference

Decleration of the Conjugate Gradients optimizer. More...

#include <ConjugateGradients.h>

Inheritance diagram for eonc::ConjugateGradients:

Public Member Functions

 ConjugateGradients (std::shared_ptr< ObjectiveFunction > a_objf, const Parameters &a_params)
 Conjugate Gradients optimizer constructor.
 ~ConjugateGradients ()=default
 Conjugant Gradient deconstructor.
int step (double a_maxMove) override
 Calls the next step in the algorithm.
int run (size_t a_maxIterations, double a_maxMove) override
 Runs the conjugate gradient.
Eigen::VectorXd getStep ()
 Gets the direction of the next step.
Public Member Functions inherited from eonc::Optimizer
 Optimizer (std::shared_ptr< ObjectiveFunction > a_objf, const OptimizerConfig &a_config)
 Optimizer (std::shared_ptr< ObjectiveFunction > a_objf, OptType a_optype, const OptimizerConfig &a_config)
 Optimizer (std::shared_ptr< ObjectiveFunction > a_objf, const Parameters &a_params)
 Optimizer (std::shared_ptr< ObjectiveFunction > a_objf, OptType a_optype, const Parameters &a_params)
virtual ~Optimizer ()

Private Member Functions

int single_step (double a_maxMove)
 Steps the conjugate gradient.
int line_search (double a_maxMove)
 Steps the conjugate gradient.

Private Attributes

Eigen::VectorXd m_direction
 Current step direction of the conjugate gradient.
Eigen::VectorXd m_directionOld
 Algorithms previous step direction.
Eigen::VectorXd m_directionNorm
 Normalised version of the current direction vector.
Eigen::VectorXd m_force
 Current force vector.
Eigen::VectorXd m_forceOld
 Previous force vector.
eonc::log::FileScoped m_log {"cg", "_cg.log"}
size_t m_cg_i
 Counts the number of discrete steps until algorithm convergence.

Additional Inherited Members

Protected Attributes inherited from eonc::Optimizer
const OptimizerConfig m_optConfig
std::shared_ptr< ObjectiveFunctionm_objf

Detailed Description

Decleration of the Conjugate Gradients optimizer.

Definition at line 38 of file ConjugateGradients.h.

Constructor & Destructor Documentation

◆ ConjugateGradients()

eonc::ConjugateGradients::ConjugateGradients ( std::shared_ptr< ObjectiveFunction > a_objf,
const Parameters & a_params )
inline

Conjugate Gradients optimizer constructor.

Parameters
std::shared_ptr<ObjectiveFunction>m_objf that tells the optimizer how to run
constParameters& m_params defined by the config.init file

Definition at line 46 of file ConjugateGradients.h.

48 : Optimizer(a_objf, OptType::CG, a_params),
49 m_directionOld{(a_objf->getPositions()).setZero()},
50 m_forceOld{(a_objf->getPositions()).setZero()}, // use setZero instead
51 m_cg_i{0} {}
Eigen::VectorXd m_directionOld
Algorithms previous step direction.
size_t m_cg_i
Counts the number of discrete steps until algorithm convergence.
Eigen::VectorXd m_forceOld
Previous force vector.
Optimizer(std::shared_ptr< ObjectiveFunction > a_objf, const OptimizerConfig &a_config)
Definition Optimizer.h:71

◆ ~ConjugateGradients()

eonc::ConjugateGradients::~ConjugateGradients ( )
default

Conjugant Gradient deconstructor.

Member Function Documentation

◆ getStep()

Eigen::VectorXd ConjugateGradients::getStep ( )

Gets the direction of the next step.

Definition at line 17 of file ConjugateGradients.cpp.

17 {
18 double a = std::fabs(m_force.dot(m_forceOld));
19 double b = m_forceOld.squaredNorm();
20 double gamma = 0.0;
21 if (a < 0.5 * b) {
22 // Polak-Ribiere way to determine how much to mix in of old direction
23 gamma = eonc::safemath::safe_div(m_force.dot(m_force - m_forceOld), b, 0.0);
24 } else {
25 gamma = 0;
26 }
29 eonc::safemath::safe_normalize_inplace(m_directionNorm);
32
33 // Only if value for max nr of iteration before reset
34 if (m_optConfig.opts.cg.max_iter_before_reset > 0 &&
35 m_optConfig.opts.cg.max_iter_before_reset <= m_cg_i) {
36 m_cg_i = 0;
37 m_forceOld.setZero();
38 m_directionOld.setZero();
39 }
40 m_cg_i += 1;
41
42 return m_direction;
43}
Eigen::VectorXd m_force
Current force vector.
Eigen::VectorXd m_directionNorm
Normalised version of the current direction vector.
Eigen::VectorXd m_direction
Current step direction of the conjugate gradient.
const OptimizerConfig m_optConfig
Definition Optimizer.h:67
constexpr double safe_div(double num, double denom, double fallback=0.0)
Definition SafeMath.h:21

◆ line_search()

int ConjugateGradients::line_search ( double a_maxMove)
private

Steps the conjugate gradient.

Checks for convergence based on the ratio of the projected force along a line to the norm of the total force

Definition at line 57 of file ConjugateGradients.cpp.

57 {
58 Eigen::VectorXd pos;
59 Eigen::VectorXd posStep;
60 Eigen::VectorXd forceBeforeStep;
61 double stepSize;
62 double projectedForce;
63 double projectedForceBeforeStep;
64 double curvature;
65
66 forceBeforeStep = -m_objf->getGradient();
67 m_force = forceBeforeStep;
68 getStep();
69 pos = m_objf->getPositions();
70 projectedForceBeforeStep = m_force.dot(m_directionNorm);
71
72 // move system an infinitesimal step to determine the optimal step size along
73 // the search line
74 posStep = pos + m_directionNorm * m_optConfig.finiteDifference;
75 m_objf->setPositions(posStep);
76 m_force = -m_objf->getGradient(true);
77 projectedForce = m_force.dot(m_directionNorm);
78 stepSize = m_optConfig.finiteDifference;
79
80 int line_i = 0;
81 do {
82 // Determine curvature from last step (Secant method)
83 curvature = std::fabs(eonc::safemath::safe_div(
84 projectedForceBeforeStep - projectedForce, stepSize, 0.0));
85 stepSize = eonc::safemath::safe_div(projectedForce, curvature, a_maxMove);
86
87 if (a_maxMove < std::fabs(stepSize)) {
88 // first part get the sign of stepSize
89 stepSize = ((stepSize > 0) - (stepSize < 0)) * a_maxMove;
90 }
91
92 forceBeforeStep = m_force;
93 projectedForceBeforeStep = projectedForce;
94
95 pos += stepSize * m_directionNorm;
96 m_objf->setPositions(pos);
97 m_force = -m_objf->getGradient();
98 projectedForce = m_force.dot(m_directionNorm);
99
100 line_i += 1;
101
102 // Line search considered converged based in the ratio between the projected
103 // force and the norm of the true force
104 } while (m_optConfig.opts.cg.line_converged <
105 std::fabs(projectedForce) /
106 (std::sqrt(m_force.dot(m_force) +
107 m_optConfig.opts.cg.line_converged)) &&
108 line_i < m_optConfig.opts.cg.line_search_max_iter);
109 if (m_objf->isConverged())
110 return 1;
111 return 0;
112}
Eigen::VectorXd getStep()
Gets the direction of the next step.
std::shared_ptr< ObjectiveFunction > m_objf
Definition Optimizer.h:68

◆ run()

int ConjugateGradients::run ( size_t a_maxIterations,
double a_maxMove )
overridevirtual

Runs the conjugate gradient.

Todo
method should also return an error code and message if the algorithm errors out
Returns
algorithm convergence

Implements eonc::Optimizer.

Definition at line 192 of file ConjugateGradients.cpp.

192 {
193 size_t iterations = 0;
194 while (!m_objf->isConverged() && iterations <= a_maxIterations) {
195 step(a_maxMove);
196 iterations++;
197 }
198 return m_objf->isConverged() ? 1 : 0;
199}
int step(double a_maxMove) override
Calls the next step in the algorithm.

◆ single_step()

int ConjugateGradients::single_step ( double a_maxMove)
private

Steps the conjugate gradient.

Checks for convergence based on the change in displacement or direction

Definition at line 114 of file ConjugateGradients.cpp.

114 {
115 Eigen::VectorXd pos;
116 Eigen::VectorXd posStep;
117 Eigen::VectorXd forceAfterStep;
118
119 m_force = -m_objf->getGradient();
120 pos = m_objf->getPositions();
121 getStep();
122
123 // move system an infinitesimal step to determine the optimal step size along
124 // the search line
125 posStep = pos + m_directionNorm * m_optConfig.finiteDifference;
126 m_objf->setPositions(posStep);
127 forceAfterStep = -m_objf->getGradient(true);
128
129 // Determine curvature
130 double projectedForce1 = m_force.dot(m_directionNorm);
131 double projectedForce2 = forceAfterStep.dot(m_directionNorm);
132 double curvature = eonc::safemath::safe_div(
133 projectedForce1 - projectedForce2, m_optConfig.finiteDifference, 0.0);
134
135 double stepSize = a_maxMove;
136
137 if (curvature > 0.0) {
138 stepSize = eonc::safemath::safe_div(projectedForce1, curvature, a_maxMove);
139 }
140
141 if (m_optConfig.bowlBreakout && a_maxMove < 0.0) {
142 stepSize = -a_maxMove;
143 a_maxMove = -a_maxMove;
144 }
145
146 if (!m_optConfig.opts.cg.no_overshooting) {
147 if (m_optConfig.bowlBreakout) {
148 // max displacement is based on system not single atom
150 a_maxMove);
151 } else {
153 a_maxMove);
154 }
155 m_objf->setPositions(pos);
156 } else {
157 // negative if product of the projected forces before and after the step are
158 // in opposite directions
159 double passedMinimum = -1.;
160 double forceChange = 0.;
161 while (passedMinimum < 0.0 &&
162 0.1 * std::fabs(projectedForce1) < std::fabs(projectedForce2)) {
164 stepSize * m_directionNorm, a_maxMove);
165 m_objf->setPositions(posStep);
166 forceAfterStep = -m_objf->getGradient(true);
167 projectedForce2 = forceAfterStep.dot(m_directionNorm);
168
169 passedMinimum = projectedForce1 * projectedForce2;
170 if (passedMinimum < 0.0 &&
171 0.1 * std::fabs(projectedForce1) < std::fabs(projectedForce2)) {
172 forceChange = (projectedForce1 - projectedForce2);
173 stepSize = eonc::safemath::safe_div(projectedForce1, forceChange, 0.0) *
174 stepSize;
175 QUILL_LOG_DEBUG(m_log, "Force changed {}, step size adjusted to {}",
176 forceChange, stepSize);
177 }
178 }
179 }
180 if (m_optConfig.opts.cg.knock_out_max_move) {
181 if (stepSize >= a_maxMove) {
182 // knockout old search direction
183 m_directionOld.setZero();
184 m_forceOld.setZero();
185 QUILL_LOG_DEBUG(m_log, "Resetting the old search direction");
186 }
187 }
188
189 return m_objf->isConverged() ? 1 : 0;
190}
eonc::log::FileScoped m_log
VectorXd maxMotionAppliedV(const VectorXd v1, double maxMotion)
VectorXd maxAtomMotionAppliedV(const VectorXd v1, double maxMotion)

◆ step()

int ConjugateGradients::step ( double a_maxMove)
overridevirtual

Calls the next step in the algorithm.

Either calls the single_step or line_search method depending on the parameters

Returns
whether or not the algorithm has converged

Implements eonc::Optimizer.

Definition at line 45 of file ConjugateGradients.cpp.

45 {
46 bool converged;
47 if (m_optConfig.opts.cg.line_search) {
48 converged = line_search(a_maxMove);
49 } else {
50 converged = single_step(a_maxMove);
51 }
52 if (converged)
53 return 1;
54 return 0;
55}
int line_search(double a_maxMove)
Steps the conjugate gradient.
int single_step(double a_maxMove)
Steps the conjugate gradient.

Member Data Documentation

◆ m_cg_i

size_t eonc::ConjugateGradients::m_cg_i
private

Counts the number of discrete steps until algorithm convergence.

Definition at line 84 of file ConjugateGradients.h.

◆ m_direction

Eigen::VectorXd eonc::ConjugateGradients::m_direction
private

Current step direction of the conjugate gradient.

Definition at line 72 of file ConjugateGradients.h.

◆ m_directionNorm

Eigen::VectorXd eonc::ConjugateGradients::m_directionNorm
private

Normalised version of the current direction vector.

Definition at line 76 of file ConjugateGradients.h.

◆ m_directionOld

Eigen::VectorXd eonc::ConjugateGradients::m_directionOld
private

Algorithms previous step direction.

Definition at line 74 of file ConjugateGradients.h.

◆ m_force

Eigen::VectorXd eonc::ConjugateGradients::m_force
private

Current force vector.

Definition at line 78 of file ConjugateGradients.h.

◆ m_forceOld

Eigen::VectorXd eonc::ConjugateGradients::m_forceOld
private

Previous force vector.

Definition at line 80 of file ConjugateGradients.h.

◆ m_log

eonc::log::FileScoped eonc::ConjugateGradients::m_log {"cg", "_cg.log"}
private

Definition at line 81 of file ConjugateGradients.h.

81{"cg", "_cg.log"};

The documentation for this class was generated from the following files: