1 // This is core/vnl/algo/vnl_conjugate_gradient.cxx
2 //:
3 // \file
4 // \author Geoffrey Cross, Oxford RRG
5 // \date 15 Feb 99
6 //
7 //-----------------------------------------------------------------------------
8 #include <iostream>
9 #include "vnl_conjugate_gradient.h"
10
11 #include "vnl/vnl_cost_function.h"
12 #include "vnl/vnl_vector_ref.h"
13 #include <vnl/algo/vnl_netlib.h>
14
15 /////////////////////////////////////
16
17 vnl_conjugate_gradient::~vnl_conjugate_gradient() = default;
18
19 void
init(vnl_cost_function & f)20 vnl_conjugate_gradient::init(vnl_cost_function & f)
21 {
22 f_ = &f;
23 num_iterations_ = 0;
24 num_evaluations_ = 0;
25 start_error_ = 0;
26 end_error_ = 0;
27 }
28
29 ///////////////////////////////////////
30
31 double
valuecomputer_(double * x,void * userdata)32 vnl_conjugate_gradient::valuecomputer_(double * x, void * userdata)
33 {
34 auto * self = static_cast<vnl_conjugate_gradient *>(userdata);
35 vnl_cost_function * f = self->f_;
36 vnl_vector_ref<double> ref_x(f->get_number_of_unknowns(), x);
37
38 self->num_evaluations_++;
39
40 return f->f(ref_x);
41 }
42
43 void
gradientcomputer_(double * g,double * x,void * userdata)44 vnl_conjugate_gradient::gradientcomputer_(double * g, double * x, void * userdata)
45 {
46 auto * self = static_cast<vnl_conjugate_gradient *>(userdata);
47 vnl_cost_function * f = self->f_;
48 vnl_vector_ref<double> ref_x(f->get_number_of_unknowns(), x);
49 vnl_vector_ref<double> ref_g(f->get_number_of_unknowns(), g);
50
51 f->gradf(ref_x, ref_g);
52 }
53
54 void
valueandgradientcomputer_(double * v,double * g,double * x,void * userdata)55 vnl_conjugate_gradient::valueandgradientcomputer_(double * v, double * g, double * x, void * userdata)
56 {
57 auto * self = static_cast<vnl_conjugate_gradient *>(userdata);
58 vnl_cost_function * f = self->f_;
59 vnl_vector_ref<double> ref_x(f->get_number_of_unknowns(), x);
60 vnl_vector_ref<double> ref_g(f->get_number_of_unknowns(), g);
61
62 f->compute(ref_x, v, &ref_g);
63 }
64
65 void
preconditioner_(double * out,double * in,void * userdata)66 vnl_conjugate_gradient::preconditioner_(double * out, double * in, void * userdata)
67 {
68 // FIXME - there should be some way to set a preconditioner if you have one
69 // e.g. P = inv(diag(A'A)) for linear least squares systems.
70
71 auto * self = static_cast<vnl_conjugate_gradient *>(userdata);
72 vnl_cost_function * f = self->f_;
73
74 int n = f->get_number_of_unknowns();
75 for (int i = 0; i < n; ++i)
76 out[i] = in[i];
77 }
78
79 ///////////////////////////////////////
80 bool
minimize(vnl_vector<double> & x)81 vnl_conjugate_gradient::minimize(vnl_vector<double> & x)
82 {
83 double * xp = x.data_block();
84 double max_norm_of_gradient;
85 long number_of_iterations;
86 final_step_size_ = 0;
87 double gradient_tolerance = gtol;
88 vnl_vector<double> workspace(f_->get_number_of_unknowns() * 3);
89 long number_of_unknowns = f_->get_number_of_unknowns();
90 long error_code;
91
92 // Compute the initial value.
93 start_error_ = valuecomputer_(xp, this);
94 num_evaluations_ = 0;
95
96 // Run the conjugate gradient algorithm.
97 v3p_netlib_cg_(xp,
98 &max_norm_of_gradient,
99 &number_of_iterations,
100 &final_step_size_,
101 &gradient_tolerance,
102 &maxfev,
103 &number_of_unknowns,
104 &number_of_unknowns,
105 valuecomputer_,
106 gradientcomputer_,
107 valueandgradientcomputer_,
108 preconditioner_,
109 workspace.data_block(),
110 this,
111 &error_code);
112
113 // Check for an error condition.
114 if (error_code > 0)
115 {
116 failure_code_ = ERROR_DODGY_INPUT;
117 if (verbose_)
118 {
119 switch (error_code)
120 {
121 case 1:
122 std::cout << "UNABLE TO OBTAIN DESCENT DIRECTION\n";
123 break;
124 case 2:
125 std::cout << "THE FUNCTION DECREASES WITH NO MINIMUM\n";
126 break;
127 case 3:
128 std::cout << "PRECONDITIONER NOT POSITIVE DEFINITE\n";
129 break;
130 case 4:
131 std::cout << "UNABLE TO SATISFY ARMIJO CONDITION\n";
132 break;
133 default:
134 std::cout << "UNKNOWN ERROR CODE\n";
135 break;
136 }
137 }
138 }
139
140 // Compute the final value.
141 end_error_ = valuecomputer_(xp, this);
142 num_iterations_ = number_of_iterations;
143
144 return error_code == 0;
145 }
146
147
148 void
diagnose_outcome(std::ostream & os) const149 vnl_conjugate_gradient::diagnose_outcome(std::ostream & os) const
150 {
151 os << "vnl_conjugate_gradient: " << num_iterations_ << " iterations, " << num_evaluations_
152 << " evaluations. Cost function reported error" << f_->reported_error(start_error_) << '/'
153 << f_->reported_error(end_error_) << " . Final step size = " << final_step_size_ << std::endl;
154 }
155
156 void
diagnose_outcome() const157 vnl_conjugate_gradient::diagnose_outcome() const
158 {
159 diagnose_outcome(std::cout);
160 }
161