| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | |||
| 2 | #ifdef GEOGRAM_WITH_HLBFGS | ||
| 3 | |||
| 4 | #include <geogram/basic/command_line.h> | ||
| 5 | #include <geogram/basic/command_line_args.h> | ||
| 6 | #include <geogram/numerics/optimizer.h> | ||
| 7 | #include <geogram/numerics/lbfgs_optimizers.h> | ||
| 8 | #include <geogram/third_party/HLBFGS/Lite_Sparse_Matrix.h> | ||
| 9 | #include <geogram/third_party/HLBFGS/HLBFGS.h> | ||
| 10 | #include <iostream> | ||
| 11 | |||
| 12 | namespace { | ||
| 13 | Lite_Sparse_Matrix* m_sparse_matrix = nullptr; | ||
| 14 | |||
| 15 | /********************** Initial HLBFGS API ***********************************/ | ||
| 16 | |||
| 17 | void evalfunc_C(int N, double* x, double *prev_x, double* f, double* g) { | ||
| 18 | GEO::geo_argused(prev_x); | ||
| 19 | *f = 0; | ||
| 20 | for (int i = 0; i < N; i+=2) { | ||
| 21 | double T1 = 1 - x[i]; | ||
| 22 | double T2 = 10*(x[i+1]-x[i]*x[i]); | ||
| 23 | *f += T1*T1+T2*T2; | ||
| 24 | g[i+1] = 20*T2; | ||
| 25 | g[i] = -2*(x[i]*g[i+1]+T1); | ||
| 26 | } | ||
| 27 | } | ||
| 28 | |||
| 29 | void newiteration_C( | ||
| 30 | int iter, int call_iter, double *x, | ||
| 31 | double* f, double *g, double* gnorm | ||
| 32 | ) { | ||
| 33 | GEO::geo_argused(x); | ||
| 34 | GEO::geo_argused(g); | ||
| 35 | std::cout << iter <<": " << call_iter <<" " | ||
| 36 | << *f <<" " << *gnorm << std::endl; | ||
| 37 | } | ||
| 38 | |||
| 39 | void evalfunc_h_C( | ||
| 40 | int N, double *x, double *prev_x, double *f, double *g, | ||
| 41 | HESSIAN_MATRIX& hessian | ||
| 42 | ) { | ||
| 43 | GEO::geo_argused(prev_x); | ||
| 44 | |||
| 45 | //the following code is not optimal if the pattern of | ||
| 46 | // hessian matrix is fixed. | ||
| 47 | if (m_sparse_matrix) { | ||
| 48 | delete m_sparse_matrix; | ||
| 49 | } | ||
| 50 | |||
| 51 | m_sparse_matrix = new Lite_Sparse_Matrix( | ||
| 52 | (unsigned int)N, (unsigned int)N, SYM_LOWER, CCS, FORTRAN_TYPE, true | ||
| 53 | ); | ||
| 54 | |||
| 55 | m_sparse_matrix->begin_fill_entry(); | ||
| 56 | |||
| 57 | static bool first = true; | ||
| 58 | double *diag = m_sparse_matrix->get_diag(); | ||
| 59 | |||
| 60 | if (first) { | ||
| 61 | // you need to update f and g | ||
| 62 | *f = 0; | ||
| 63 | double tmp; | ||
| 64 | for (unsigned int i = 0; i < (unsigned int)N; i+=2) { | ||
| 65 | tmp = x[i]*x[i]; | ||
| 66 | double T1 = 1 - x[i]; | ||
| 67 | double T2 = 10*(x[i+1]-tmp); | ||
| 68 | *f += T1*T1+T2*T2; | ||
| 69 | g[i+1] = 20*T2; | ||
| 70 | g[i] = -2*(x[i]*g[i+1]+T1); | ||
| 71 | diag[i] = 2+1200*tmp-400*x[i+1]; | ||
| 72 | diag[i+1] = 200; | ||
| 73 | m_sparse_matrix->fill_entry(i, i+1, -400*x[i]); | ||
| 74 | } | ||
| 75 | } else { | ||
| 76 | for (unsigned int i = 0; i < (unsigned int)N; i+=2) { | ||
| 77 | diag[i] = 2+1200*x[i]*x[i]-400*x[i+1]; | ||
| 78 | diag[i+1] = 200; | ||
| 79 | m_sparse_matrix->fill_entry(i, i+1, -400*x[i]); | ||
| 80 | } | ||
| 81 | } | ||
| 82 | m_sparse_matrix->end_fill_entry(); | ||
| 83 | hessian.set_diag(m_sparse_matrix->get_diag()); | ||
| 84 | hessian.set_values(m_sparse_matrix->get_values()); | ||
| 85 | hessian.set_rowind((int*)m_sparse_matrix->get_rowind()); | ||
| 86 | hessian.set_colptr((int*)m_sparse_matrix->get_colptr()); | ||
| 87 | hessian.set_nonzeros((int)m_sparse_matrix->get_nonzero()); | ||
| 88 | first = false; | ||
| 89 | } | ||
| 90 | |||
| 91 | void Optimize_by_HLBFGS_C( | ||
| 92 | int N, double *init_x, int num_iter, int M, int T, bool with_hessian | ||
| 93 | ) { | ||
| 94 | double parameter[20]; | ||
| 95 | int info[20]; | ||
| 96 | //initialize | ||
| 97 | INIT_HLBFGS(parameter, info); | ||
| 98 | info[4] = num_iter; | ||
| 99 | info[6] = T; | ||
| 100 | info[7] = with_hessian?1:0; | ||
| 101 | info[10] = 0; | ||
| 102 | info[11] = 1; | ||
| 103 | |||
| 104 | if (with_hessian) { | ||
| 105 | HLBFGS( | ||
| 106 | N, M, init_x, | ||
| 107 | evalfunc_C, evalfunc_h_C, | ||
| 108 | HLBFGS_UPDATE_Hessian, newiteration_C, parameter, info | ||
| 109 | ); | ||
| 110 | } else { | ||
| 111 | HLBFGS( | ||
| 112 | N, M, init_x, | ||
| 113 | evalfunc_C, nullptr, | ||
| 114 | HLBFGS_UPDATE_Hessian, newiteration_C, parameter, info | ||
| 115 | ); | ||
| 116 | } | ||
| 117 | } | ||
| 118 | |||
| 119 | /********************** Geogram API **************************************/ | ||
| 120 | |||
| 121 | void evalfunc(GEO::index_t N, double* x, double& f, double* g) { | ||
| 122 | std::cerr << "eval func" << std::endl; | ||
| 123 | f = 0.0; | ||
| 124 | for (GEO::index_t i = 0; i < N; i+=2) { | ||
| 125 | double T1 = 1 - x[i]; | ||
| 126 | double T2 = 10*(x[i+1]-x[i]*x[i]); | ||
| 127 | f += T1*T1+T2*T2; | ||
| 128 | g[i+1] = 20*T2; | ||
| 129 | g[i] = -2*(x[i]*g[i+1]+T1); | ||
| 130 | } | ||
| 131 | } | ||
| 132 | |||
| 133 | void newiteration( | ||
| 134 | GEO::index_t iter, | ||
| 135 | const double *x, double f, const double *g, double gnorm | ||
| 136 | ) { | ||
| 137 | GEO::geo_argused(iter); | ||
| 138 | GEO::geo_argused(x); | ||
| 139 | GEO::geo_argused(g); | ||
| 140 | std::cout << " " << f <<" " << gnorm << std::endl; | ||
| 141 | } | ||
| 142 | |||
| 143 | void evalfunc_h( | ||
| 144 | GEO::index_t N, double *x, double& f, double *g, | ||
| 145 | HESSIAN_MATRIX& hessian | ||
| 146 | ) { | ||
| 147 | std::cerr << "eval func with Hessian" << std::endl; | ||
| 148 | |||
| 149 | //the following code is not optimal if the pattern of | ||
| 150 | // hessian matrix is fixed. | ||
| 151 | if (m_sparse_matrix) { | ||
| 152 | delete m_sparse_matrix; | ||
| 153 | } | ||
| 154 | |||
| 155 | m_sparse_matrix = new Lite_Sparse_Matrix( | ||
| 156 | N, N, SYM_LOWER, CCS, FORTRAN_TYPE, true | ||
| 157 | ); | ||
| 158 | |||
| 159 | m_sparse_matrix->begin_fill_entry(); | ||
| 160 | |||
| 161 | static bool first = true; | ||
| 162 | double *diag = m_sparse_matrix->get_diag(); | ||
| 163 | |||
| 164 | if (first) { | ||
| 165 | // you need to update f and g | ||
| 166 | f = 0.0; | ||
| 167 | double tmp; | ||
| 168 | for (unsigned int i = 0; i < N; i+=2) { | ||
| 169 | tmp = x[i]*x[i]; | ||
| 170 | double T1 = 1 - x[i]; | ||
| 171 | double T2 = 10*(x[i+1]-tmp); | ||
| 172 | f += T1*T1+T2*T2; | ||
| 173 | g[i+1] = 20*T2; | ||
| 174 | g[i] = -2*(x[i]*g[i+1]+T1); | ||
| 175 | diag[i] = 2+1200*tmp-400*x[i+1]; | ||
| 176 | diag[i+1] = 200; | ||
| 177 | m_sparse_matrix->fill_entry(i, i+1, -400*x[i]); | ||
| 178 | } | ||
| 179 | } else { | ||
| 180 | for (unsigned int i = 0; i < N; i+=2) { | ||
| 181 | diag[i] = 2+1200*x[i]*x[i]-400*x[i+1]; | ||
| 182 | diag[i+1] = 200; | ||
| 183 | m_sparse_matrix->fill_entry(i, i+1, -400*x[i]); | ||
| 184 | } | ||
| 185 | } | ||
| 186 | m_sparse_matrix->end_fill_entry(); | ||
| 187 | hessian.set_diag(m_sparse_matrix->get_diag()); | ||
| 188 | hessian.set_values(m_sparse_matrix->get_values()); | ||
| 189 | hessian.set_rowind((int*)m_sparse_matrix->get_rowind()); | ||
| 190 | hessian.set_colptr((int*)m_sparse_matrix->get_colptr()); | ||
| 191 | hessian.set_nonzeros((int)m_sparse_matrix->get_nonzero()); | ||
| 192 | first = false; | ||
| 193 | } | ||
| 194 | |||
| 195 | void Optimize_by_HLBFGS( | ||
| 196 | int N, double *init_x, int num_iter, int M, int T, bool with_hessian | ||
| 197 | ) { | ||
| 198 | GEO::Optimizer_var optimizer = with_hessian ? | ||
| 199 | GEO::Optimizer::create("HLBFGS_HESS") : | ||
| 200 | GEO::Optimizer::create("HLBFGS"); | ||
| 201 | |||
| 202 | |||
| 203 | optimizer->set_newiteration_callback(newiteration); | ||
| 204 | if(with_hessian) { | ||
| 205 | optimizer->set_evalhessian_callback(evalfunc_h); | ||
| 206 | } | ||
| 207 | optimizer->set_funcgrad_callback(evalfunc); | ||
| 208 | |||
| 209 | optimizer->set_N((unsigned int)N); | ||
| 210 | optimizer->set_M((unsigned int)M); | ||
| 211 | optimizer->set_max_iter((unsigned int)num_iter); | ||
| 212 | |||
| 213 | GEO::HLBFGS_HessOptimizer* hess = | ||
| 214 | dynamic_cast<GEO::HLBFGS_HessOptimizer*>( | ||
| 215 | (GEO::Optimizer*)(optimizer) | ||
| 216 | ); | ||
| 217 | |||
| 218 | if(hess != nullptr) { | ||
| 219 | hess->set_T((unsigned int)T); | ||
| 220 | } | ||
| 221 | optimizer->optimize(init_x); | ||
| 222 | } | ||
| 223 | } | ||
| 224 | |||
| 225 | /****************************************************************************/ | ||
| 226 | |||
| 227 | int main(int argc, char** argv) { | ||
| 228 | |||
| 229 | GEO::initialize(GEO::GEOGRAM_INSTALL_ALL); | ||
| 230 | GEO::CmdLine::import_arg_group("standard"); | ||
| 231 | GEO::CmdLine::declare_arg("Newton",false,"Use Newton solver"); | ||
| 232 | GEO::CmdLine::declare_arg("C_api",false,"Use HLBFGS C api"); | ||
| 233 | GEO::CmdLine::declare_arg("N", 1000, "Nb variables"); | ||
| 234 | |||
| 235 | if(!GEO::CmdLine::parse(argc, argv)) { | ||
| 236 | return 1; | ||
| 237 | } | ||
| 238 | |||
| 239 | std::cout.precision(16); | ||
| 240 | std::cout << std::scientific; | ||
| 241 | |||
| 242 | int N = GEO::CmdLine::get_arg_int("N"); | ||
| 243 | std::vector<double> x((unsigned int)N); | ||
| 244 | |||
| 245 | for (unsigned int i = 0; i < (unsigned int)(N/2); i++) { | ||
| 246 | x[2*i] = -1.2; | ||
| 247 | x[2*i+1] = 1.0; | ||
| 248 | } | ||
| 249 | |||
| 250 | int M = 7; | ||
| 251 | int T = 0; | ||
| 252 | |||
| 253 | if(GEO::CmdLine::get_arg_bool("C_api")) { | ||
| 254 | if(GEO::CmdLine::get_arg_bool("Newton")) { | ||
| 255 | //use Hessian | ||
| 256 | // if M = 0, T = 0, it is Newton | ||
| 257 | Optimize_by_HLBFGS_C(N, &x[0], 1000, M, T, true); | ||
| 258 | } else { | ||
| 259 | //without Hessian | ||
| 260 | // it is LBFGS(M) actually, T is not used | ||
| 261 | Optimize_by_HLBFGS_C(N, &x[0], 1000, M, T, false); | ||
| 262 | } | ||
| 263 | } else { | ||
| 264 | if(GEO::CmdLine::get_arg_bool("Newton")) { | ||
| 265 | //use Hessian | ||
| 266 | // if M = 0, T = 0, it is Newton | ||
| 267 | Optimize_by_HLBFGS(N, &x[0], 1000, M, T, true); | ||
| 268 | } else { | ||
| 269 | //without Hessian | ||
| 270 | // it is LBFGS(M) actually, T is not used | ||
| 271 | Optimize_by_HLBFGS(N, &x[0], 1000, M, T, false); | ||
| 272 | } | ||
| 273 | } | ||
| 274 | |||
| 275 | if (m_sparse_matrix) { | ||
| 276 | delete m_sparse_matrix; | ||
| 277 | } | ||
| 278 | |||
| 279 | return 0; | ||
| 280 | } | ||
| 281 | |||
| 282 | #else | ||
| 283 | |||
| 284 | #include <iostream> | ||
| 285 | |||
| 286 | ✗ | int main() { | |
| 287 | std::cout << "This geogram was not compiled with HLBFGS support" | ||
| 288 | << std::endl; | ||
| 289 | return 0; | ||
| 290 | } | ||
| 291 | |||
| 292 | #endif | ||
| 293 |