Go to the documentation of this file.00001
00002 #ifdef VCL_NEEDS_PRAGMA_INTERFACE
00003 #pragma implementation
00004 #endif
00005
00006
00007
00008
00009
00010
00011
00012
00013 #include "vnl_lbfgsb.h"
00014 #include <vcl_cstring.h>
00015 #include <vcl_iostream.h>
00016
00017 #include <vnl/algo/vnl_netlib.h>
00018
00019
00020 vnl_lbfgsb::vnl_lbfgsb(): f_(0)
00021 {
00022 init_parameters();
00023 }
00024
00025
00026 vnl_lbfgsb::vnl_lbfgsb(vnl_cost_function& f): f_(&f)
00027 {
00028 init_parameters();
00029 }
00030
00031
00032 void vnl_lbfgsb::init_parameters()
00033 {
00034 long n = this->f_->get_number_of_unknowns();
00035 this->bound_selection_.set_size(n);
00036 this->bound_selection_.fill(0);
00037 this->max_corrections_ = 5;
00038 this->convergence_factor_ = 1e+7;
00039 this->projected_gradient_tolerance_ = 1e-5;
00040 }
00041
00042
00043 bool vnl_lbfgsb::minimize(vnl_vector<double>& x)
00044 {
00045
00046 long n = this->f_->get_number_of_unknowns();
00047 long m = this->max_corrections_;
00048
00049
00050 double f = 0;
00051 vnl_vector<double> gradient(n);
00052
00053
00054 vnl_vector<double> wa((2*m+4)*n + 12*m*m + 12*m);
00055 vnl_vector<long> iwa(3*n);
00056 char csave[60];
00057 long lsave[4];
00058 long isave[44];
00059 double dsave[29];
00060
00061
00062 char task[61]="START ";
00063
00064
00065
00066 long const iprint = trace ? 1 : -1;
00067
00068
00069 this->num_evaluations_ = 0;
00070 this->num_iterations_ = 0;
00071
00072
00073
00074
00075
00076 vnl_vector<double> x_best(x);
00077
00078 bool ok = true;
00079 for (;;)
00080 {
00081
00082 v3p_netlib_setulb_(
00083 &n,
00084 &m,
00085 x.data_block(),
00086 this->lower_bound_.data_block(),
00087 this->upper_bound_.data_block(),
00088 this->bound_selection_.data_block(),
00089 &f, gradient.data_block(),
00090 &this->convergence_factor_,
00091 &this->projected_gradient_tolerance_,
00092 wa.data_block(),
00093 iwa.data_block(),
00094 task,
00095 &iprint,
00096 csave, lsave, isave, dsave
00097 );
00098
00099
00100 if (vcl_strncmp("FG", task, 2) == 0)
00101 {
00102
00103 this->f_->compute(x, &f, &gradient);
00104
00105 if (this->num_evaluations_ == 0)
00106 {
00107 x_best = x;
00108 this->start_error_ = f;
00109 this->end_error_ = f;
00110 }
00111 else if (f < this->end_error_)
00112 {
00113 x_best = x;
00114 this->end_error_ = f;
00115 }
00116 this->report_eval(f);
00117 }
00118 else if (vcl_strncmp("NEW_X", task, 5) == 0)
00119 {
00120
00121 this->inf_norm_projected_gradient_ = dsave[12];
00122
00123
00124 if (this->report_iter())
00125 {
00126 this->failure_code_ = FAILED_USER_REQUEST;
00127 ok = false;
00128 break;
00129 }
00130 }
00131 else if (vcl_strncmp("ERROR", task, 5) == 0)
00132 {
00133
00134 this->failure_code_ = ERROR_FAILURE;
00135 ok = false;
00136 break;
00137 }
00138 else if (vcl_strncmp("CONVERGENCE", task, 11) == 0)
00139 {
00140
00141 if (f < this->end_error_)
00142 {
00143 x_best = x;
00144 this->end_error_ = f;
00145 }
00146
00147 if (vcl_strncmp("CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH",
00148 task, 47) == 0)
00149 {
00150
00151 this->failure_code_ = CONVERGED_FTOL;
00152 }
00153 else if (vcl_strncmp("CONVERGENCE: NORM OF PROJECTED GRADIENT <= PGTOL",
00154 task, 48) == 0)
00155 {
00156
00157 this->failure_code_ = CONVERGED_GTOL;
00158 }
00159 else
00160 {
00161 this->failure_code_ = ERROR_FAILURE;
00162 if (trace)
00163 {
00164 vcl_cerr << "Unknown convergence type: " << task << vcl_endl;
00165 }
00166 }
00167 break;
00168 }
00169 else
00170 {
00171
00172 this->failure_code_ = ERROR_FAILURE;
00173 if (trace)
00174 {
00175 vcl_cerr << "Unknown failure with task: " << task << vcl_endl;
00176 }
00177 ok = false;
00178 break;
00179 }
00180
00181 if (this->num_evaluations_ > this->get_max_function_evals())
00182 {
00183 this->failure_code_ = TOO_MANY_ITERATIONS;
00184 ok = false;
00185 break;
00186 }
00187 }
00188
00189
00190 x = x_best;
00191
00192 return ok;
00193 }