[66cfc7] | 1 | /*
|
---|
| 2 | * Project: MoleCuilder
|
---|
| 3 | * Description: creates and alters molecular systems
|
---|
| 4 | * Copyright (C) 2012 University of Bonn. All rights reserved.
|
---|
| 5 | * Please see the COPYING file or "Copyright notice" in builder.cpp for details.
|
---|
| 6 | *
|
---|
| 7 | *
|
---|
| 8 | * This file is part of MoleCuilder.
|
---|
| 9 | *
|
---|
| 10 | * MoleCuilder is free software: you can redistribute it and/or modify
|
---|
| 11 | * it under the terms of the GNU General Public License as published by
|
---|
| 12 | * the Free Software Foundation, either version 2 of the License, or
|
---|
| 13 | * (at your option) any later version.
|
---|
| 14 | *
|
---|
| 15 | * MoleCuilder is distributed in the hope that it will be useful,
|
---|
| 16 | * but WITHOUT ANY WARRANTY; without even the implied warranty of
|
---|
| 17 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
---|
| 18 | * GNU General Public License for more details.
|
---|
| 19 | *
|
---|
| 20 | * You should have received a copy of the GNU General Public License
|
---|
| 21 | * along with MoleCuilder. If not, see <http://www.gnu.org/licenses/>.
|
---|
| 22 | */
|
---|
| 23 |
|
---|
| 24 | /*
|
---|
| 25 | * FunctionApproximation.cpp
|
---|
| 26 | *
|
---|
| 27 | * Created on: 02.10.2012
|
---|
| 28 | * Author: heber
|
---|
| 29 | */
|
---|
| 30 |
|
---|
| 31 | // include config.h
|
---|
| 32 | #ifdef HAVE_CONFIG_H
|
---|
| 33 | #include <config.h>
|
---|
| 34 | #endif
|
---|
| 35 |
|
---|
| 36 | #include "CodePatterns/MemDebug.hpp"
|
---|
| 37 |
|
---|
| 38 | #include "FunctionApproximation.hpp"
|
---|
| 39 |
|
---|
| 40 | #include <algorithm>
|
---|
| 41 | #include <boost/bind.hpp>
|
---|
| 42 | #include <boost/function.hpp>
|
---|
| 43 | #include <iostream>
|
---|
| 44 | #include <iterator>
|
---|
| 45 | #include <numeric>
|
---|
| 46 | #include <sstream>
|
---|
| 47 |
|
---|
| 48 | #include <levmar.h>
|
---|
| 49 |
|
---|
| 50 | #include "CodePatterns/Assert.hpp"
|
---|
| 51 | #include "CodePatterns/Log.hpp"
|
---|
| 52 |
|
---|
| 53 | #include "FunctionApproximation/FunctionModel.hpp"
|
---|
[69ab84] | 54 | #include "FunctionApproximation/TrainingData.hpp"
|
---|
| 55 |
|
---|
| 56 | FunctionApproximation::FunctionApproximation(
|
---|
| 57 | const TrainingData &_data,
|
---|
| 58 | FunctionModel &_model) :
|
---|
| 59 | input_dimension(_data.getTrainingInputs().size()),
|
---|
| 60 | output_dimension(_data.getTrainingOutputs().size()),
|
---|
| 61 | input_data(_data.getTrainingInputs()),
|
---|
| 62 | output_data(_data.getTrainingOutputs()),
|
---|
| 63 | model(_model)
|
---|
| 64 | {}
|
---|
[66cfc7] | 65 |
|
---|
[e1fe7e] | 66 | void FunctionApproximation::setTrainingData(const filtered_inputs_t &input, const outputs_t &output)
|
---|
[66cfc7] | 67 | {
|
---|
| 68 | ASSERT( input.size() == output.size(),
|
---|
| 69 | "FunctionApproximation::setTrainingData() - the number of input and output tuples differ: "+toString(input.size())+"!="
|
---|
| 70 | +toString(output.size())+".");
|
---|
| 71 | if (input.size() != 0) {
|
---|
| 72 | ASSERT( input[0].size() == input_dimension,
|
---|
| 73 | "FunctionApproximation::setTrainingData() - the dimension of the input tuples and input dimension differ: "+toString(input[0].size())+"!="
|
---|
| 74 | +toString(input_dimension)+".");
|
---|
| 75 | input_data = input;
|
---|
| 76 | ASSERT( output[0].size() == output_dimension,
|
---|
| 77 | "FunctionApproximation::setTrainingData() - the dimension of the output tuples and output dimension differ: "+toString(output[0].size())+"!="
|
---|
| 78 | +toString(output_dimension)+".");
|
---|
| 79 | output_data = output;
|
---|
| 80 | } else {
|
---|
| 81 | ELOG(2, "Given vectors of training data are empty, clearing internal vectors accordingly.");
|
---|
| 82 | input_data.clear();
|
---|
| 83 | output_data.clear();
|
---|
| 84 | }
|
---|
| 85 | }
|
---|
| 86 |
|
---|
| 87 | void FunctionApproximation::setModelFunction(FunctionModel &_model)
|
---|
| 88 | {
|
---|
| 89 | model= _model;
|
---|
| 90 | }
|
---|
| 91 |
|
---|
| 92 | /** Callback to circumvent boost::bind, boost::function and function pointer problem.
|
---|
| 93 | *
|
---|
| 94 | * See here (second answer!) to the nature of the problem:
|
---|
| 95 | * http://stackoverflow.com/questions/282372/demote-boostfunction-to-a-plain-function-pointer
|
---|
| 96 | *
|
---|
| 97 | * We cannot use a boost::bind bounded boost::function as a function pointer.
|
---|
| 98 | * boost::function::target() will just return NULL because the signature does not
|
---|
| 99 | * match. We have to use a C-style callback function and our luck is that
|
---|
| 100 | * the levmar signature provides for a void* additional data pointer which we
|
---|
| 101 | * can cast back to our FunctionApproximation class, as we need access to the
|
---|
| 102 | * data contained, e.g. the FunctionModel reference FunctionApproximation::model.
|
---|
| 103 | *
|
---|
| 104 | */
|
---|
| 105 | void FunctionApproximation::LevMarCallback(double *p, double *x, int m, int n, void *data)
|
---|
| 106 | {
|
---|
| 107 | FunctionApproximation *approximator = static_cast<FunctionApproximation *>(data);
|
---|
| 108 | ASSERT( approximator != NULL,
|
---|
| 109 | "LevMarCallback() - received data does not represent a FunctionApproximation object.");
|
---|
| 110 | boost::function<void(double*,double*,int,int,void*)> function =
|
---|
| 111 | boost::bind(&FunctionApproximation::evaluate, approximator, _1, _2, _3, _4, _5);
|
---|
| 112 | function(p,x,m,n,data);
|
---|
| 113 | }
|
---|
| 114 |
|
---|
[5b5724] | 115 | void FunctionApproximation::LevMarDerivativeCallback(double *p, double *x, int m, int n, void *data)
|
---|
| 116 | {
|
---|
| 117 | FunctionApproximation *approximator = static_cast<FunctionApproximation *>(data);
|
---|
| 118 | ASSERT( approximator != NULL,
|
---|
| 119 | "LevMarDerivativeCallback() - received data does not represent a FunctionApproximation object.");
|
---|
| 120 | boost::function<void(double*,double*,int,int,void*)> function =
|
---|
| 121 | boost::bind(&FunctionApproximation::evaluateDerivative, approximator, _1, _2, _3, _4, _5);
|
---|
| 122 | function(p,x,m,n,data);
|
---|
| 123 | }
|
---|
| 124 |
|
---|
[63b9f7] | 125 | void FunctionApproximation::prepareParameters(double *&p, int &m) const
|
---|
[66cfc7] | 126 | {
|
---|
| 127 | m = model.getParameterDimension();
|
---|
| 128 | const FunctionModel::parameters_t params = model.getParameters();
|
---|
[c62f96] | 129 | {
|
---|
| 130 | p = new double[m];
|
---|
| 131 | size_t index = 0;
|
---|
| 132 | for(FunctionModel::parameters_t::const_iterator paramiter = params.begin();
|
---|
| 133 | paramiter != params.end();
|
---|
| 134 | ++paramiter, ++index) {
|
---|
| 135 | p[index] = *paramiter;
|
---|
| 136 | }
|
---|
| 137 | }
|
---|
[63b9f7] | 138 | }
|
---|
| 139 |
|
---|
| 140 | void FunctionApproximation::prepareOutput(double *&x, int &n) const
|
---|
| 141 | {
|
---|
| 142 | n = output_data.size();
|
---|
[c62f96] | 143 | {
|
---|
| 144 | x = new double[n];
|
---|
| 145 | size_t index = 0;
|
---|
| 146 | for(outputs_t::const_iterator outiter = output_data.begin();
|
---|
| 147 | outiter != output_data.end();
|
---|
| 148 | ++outiter, ++index) {
|
---|
| 149 | x[index] = (*outiter)[0];
|
---|
| 150 | }
|
---|
[66cfc7] | 151 | }
|
---|
[63b9f7] | 152 | }
|
---|
| 153 |
|
---|
| 154 | void FunctionApproximation::operator()(const enum JacobianMode mode)
|
---|
| 155 | {
|
---|
| 156 | // let levmar optimize
|
---|
| 157 | register int i, j;
|
---|
| 158 | int ret;
|
---|
| 159 | double *p;
|
---|
| 160 | double *x;
|
---|
| 161 | int m, n;
|
---|
| 162 | double opts[LM_OPTS_SZ], info[LM_INFO_SZ];
|
---|
| 163 |
|
---|
| 164 | opts[0]=LM_INIT_MU; opts[1]=1E-15; opts[2]=1E-15; opts[3]=1E-20;
|
---|
| 165 | opts[4]= LM_DIFF_DELTA; // relevant only if the Jacobian is approximated using finite differences; specifies forward differencing
|
---|
| 166 | //opts[4]=-LM_DIFF_DELTA; // specifies central differencing to approximate Jacobian; more accurate but more expensive to compute!
|
---|
| 167 |
|
---|
| 168 | prepareParameters(p,m);
|
---|
| 169 | prepareOutput(x,n);
|
---|
[66cfc7] | 170 |
|
---|
| 171 | {
|
---|
| 172 | double *work, *covar;
|
---|
| 173 | work=(double *)malloc((LM_DIF_WORKSZ(m, n)+m*m)*sizeof(double));
|
---|
| 174 | if(!work){
|
---|
| 175 | ELOG(0, "FunctionApproximation::operator() - memory allocation request failed.");
|
---|
| 176 | return;
|
---|
| 177 | }
|
---|
| 178 | covar=work+LM_DIF_WORKSZ(m, n);
|
---|
| 179 |
|
---|
| 180 | // give this pointer as additional data to construct function pointer in
|
---|
| 181 | // LevMarCallback and call
|
---|
[d03292] | 182 | if (model.isBoxConstraint()) {
|
---|
| 183 | FunctionModel::parameters_t lowerbound = model.getLowerBoxConstraints();
|
---|
| 184 | FunctionModel::parameters_t upperbound = model.getUpperBoxConstraints();
|
---|
| 185 | double *lb = new double[m];
|
---|
| 186 | double *ub = new double[m];
|
---|
[a2a2f7] | 187 | for (size_t i=0;i<(size_t)m;++i) {
|
---|
[d03292] | 188 | lb[i] = lowerbound[i];
|
---|
| 189 | ub[i] = upperbound[i];
|
---|
| 190 | }
|
---|
| 191 | if (mode == FiniteDifferences) {
|
---|
| 192 | ret=dlevmar_bc_dif(
|
---|
| 193 | &FunctionApproximation::LevMarCallback,
|
---|
| 194 | p, x, m, n, lb, ub, NULL, 1000, opts, info, work, covar, this); // no Jacobian, caller allocates work memory, covariance estimated
|
---|
| 195 | } else if (mode == ParameterDerivative) {
|
---|
| 196 | ret=dlevmar_bc_der(
|
---|
| 197 | &FunctionApproximation::LevMarCallback,
|
---|
| 198 | &FunctionApproximation::LevMarDerivativeCallback,
|
---|
| 199 | p, x, m, n, lb, ub, NULL, 1000, opts, info, work, covar, this); // no Jacobian, caller allocates work memory, covariance estimated
|
---|
| 200 | } else {
|
---|
| 201 | ASSERT(0, "FunctionApproximation::operator() - Unknown jacobian method chosen.");
|
---|
| 202 | }
|
---|
| 203 | delete[] lb;
|
---|
| 204 | delete[] ub;
|
---|
[76e63d] | 205 | } else {
|
---|
| 206 | ASSERT(0, "FunctionApproximation::operator() - Unknown jacobian method chosen.");
|
---|
[d03292] | 207 | if (mode == FiniteDifferences) {
|
---|
| 208 | ret=dlevmar_dif(
|
---|
| 209 | &FunctionApproximation::LevMarCallback,
|
---|
| 210 | p, x, m, n, 1000, opts, info, work, covar, this); // no Jacobian, caller allocates work memory, covariance estimated
|
---|
| 211 | } else if (mode == ParameterDerivative) {
|
---|
| 212 | ret=dlevmar_der(
|
---|
| 213 | &FunctionApproximation::LevMarCallback,
|
---|
| 214 | &FunctionApproximation::LevMarDerivativeCallback,
|
---|
| 215 | p, x, m, n, 1000, opts, info, work, covar, this); // no Jacobian, caller allocates work memory, covariance estimated
|
---|
| 216 | } else {
|
---|
| 217 | ASSERT(0, "FunctionApproximation::operator() - Unknown jacobian method chosen.");
|
---|
| 218 | }
|
---|
[76e63d] | 219 | }
|
---|
[66cfc7] | 220 |
|
---|
| 221 | {
|
---|
| 222 | std::stringstream covar_msg;
|
---|
| 223 | covar_msg << "Covariance of the fit:\n";
|
---|
| 224 | for(i=0; i<m; ++i){
|
---|
| 225 | for(j=0; j<m; ++j)
|
---|
| 226 | covar_msg << covar[i*m+j] << " ";
|
---|
| 227 | covar_msg << std::endl;
|
---|
| 228 | }
|
---|
| 229 | covar_msg << std::endl;
|
---|
| 230 | LOG(1, "INFO: " << covar_msg.str());
|
---|
| 231 | }
|
---|
| 232 |
|
---|
| 233 | free(work);
|
---|
| 234 | }
|
---|
| 235 |
|
---|
| 236 | {
|
---|
| 237 | std::stringstream result_msg;
|
---|
| 238 | result_msg << "Levenberg-Marquardt returned " << ret << " in " << info[5] << " iter, reason " << info[6] << "\nSolution: ";
|
---|
| 239 | for(i=0; i<m; ++i)
|
---|
| 240 | result_msg << p[i] << " ";
|
---|
| 241 | result_msg << "\n\nMinimization info:\n";
|
---|
| 242 | std::vector<std::string> infonames(LM_INFO_SZ);
|
---|
| 243 | infonames[0] = std::string("||e||_2 at initial p");
|
---|
| 244 | infonames[1] = std::string("||e||_2");
|
---|
| 245 | infonames[2] = std::string("||J^T e||_inf");
|
---|
| 246 | infonames[3] = std::string("||Dp||_2");
|
---|
| 247 | infonames[4] = std::string("mu/max[J^T J]_ii");
|
---|
| 248 | infonames[5] = std::string("# iterations");
|
---|
| 249 | infonames[6] = std::string("reason for termination");
|
---|
| 250 | infonames[7] = std::string(" # function evaluations");
|
---|
| 251 | infonames[8] = std::string(" # Jacobian evaluations");
|
---|
| 252 | infonames[9] = std::string(" # linear systems solved");
|
---|
| 253 | for(i=0; i<LM_INFO_SZ; ++i)
|
---|
| 254 | result_msg << infonames[i] << ": " << info[i] << " ";
|
---|
| 255 | result_msg << std::endl;
|
---|
| 256 | LOG(1, "INFO: " << result_msg.str());
|
---|
| 257 | }
|
---|
| 258 |
|
---|
[c62f96] | 259 | delete[] p;
|
---|
| 260 | delete[] x;
|
---|
[66cfc7] | 261 | }
|
---|
| 262 |
|
---|
[371c8b] | 263 | bool FunctionApproximation::checkParameterDerivatives()
|
---|
| 264 | {
|
---|
| 265 | double *p;
|
---|
| 266 | int m;
|
---|
| 267 | const FunctionModel::parameters_t backupparams = model.getParameters();
|
---|
| 268 | prepareParameters(p,m);
|
---|
| 269 | int n = output_data.size();
|
---|
| 270 | double *err = new double[n];
|
---|
| 271 | dlevmar_chkjac(
|
---|
| 272 | &FunctionApproximation::LevMarCallback,
|
---|
| 273 | &FunctionApproximation::LevMarDerivativeCallback,
|
---|
| 274 | p, m, n, this, err);
|
---|
| 275 | int i;
|
---|
| 276 | for(i=0; i<n; ++i)
|
---|
| 277 | LOG(1, "INFO: gradient " << i << ", err " << err[i] << ".");
|
---|
| 278 | bool status = true;
|
---|
| 279 | for(i=0; i<n; ++i)
|
---|
| 280 | status &= err[i] > 0.5;
|
---|
| 281 |
|
---|
| 282 | if (!status) {
|
---|
| 283 | int faulty;
|
---|
| 284 | ELOG(0, "At least one of the parameter derivatives are incorrect.");
|
---|
| 285 | for (faulty=1; faulty<=m; ++faulty) {
|
---|
| 286 | LOG(1, "INFO: Trying with only the first " << faulty << " parameters...");
|
---|
| 287 | model.setParameters(backupparams);
|
---|
| 288 | dlevmar_chkjac(
|
---|
| 289 | &FunctionApproximation::LevMarCallback,
|
---|
| 290 | &FunctionApproximation::LevMarDerivativeCallback,
|
---|
| 291 | p, faulty, n, this, err);
|
---|
| 292 | bool status = true;
|
---|
| 293 | for(i=0; i<n; ++i)
|
---|
| 294 | status &= err[i] > 0.5;
|
---|
| 295 | for(i=0; i<n; ++i)
|
---|
| 296 | LOG(1, "INFO: gradient(" << faulty << ") " << i << ", err " << err[i] << ".");
|
---|
| 297 | if (!status)
|
---|
| 298 | break;
|
---|
| 299 | }
|
---|
| 300 | ELOG(0, "The faulty parameter derivative is with respect to the " << faulty << " parameter.");
|
---|
| 301 | } else
|
---|
| 302 | LOG(1, "INFO: parameter derivatives are ok.");
|
---|
| 303 |
|
---|
| 304 | delete[] err;
|
---|
| 305 | delete[] p;
|
---|
| 306 | model.setParameters(backupparams);
|
---|
| 307 |
|
---|
| 308 | return status;
|
---|
| 309 | }
|
---|
| 310 |
|
---|
[66cfc7] | 311 | double SquaredDifference(const double res1, const double res2)
|
---|
| 312 | {
|
---|
| 313 | return (res1-res2)*(res1-res2);
|
---|
| 314 | }
|
---|
| 315 |
|
---|
[5b5724] | 316 | void FunctionApproximation::prepareModel(double *p, int m)
|
---|
[66cfc7] | 317 | {
|
---|
[371c8b] | 318 | // ASSERT( (size_t)m == model.getParameterDimension(),
|
---|
| 319 | // "FunctionApproximation::prepareModel() - LevMar expects "+toString(m)
|
---|
| 320 | // +" parameters but the model function expects "+toString(model.getParameterDimension())+".");
|
---|
[66cfc7] | 321 | FunctionModel::parameters_t params(m, 0.);
|
---|
| 322 | std::copy(p, p+m, params.begin());
|
---|
| 323 | model.setParameters(params);
|
---|
[5b5724] | 324 | }
|
---|
| 325 |
|
---|
| 326 | void FunctionApproximation::evaluate(double *p, double *x, int m, int n, void *data)
|
---|
| 327 | {
|
---|
| 328 | // first set parameters
|
---|
| 329 | prepareModel(p,m);
|
---|
[66cfc7] | 330 |
|
---|
| 331 | // then evaluate
|
---|
[5b5724] | 332 | ASSERT( (size_t)n == output_data.size(),
|
---|
| 333 | "FunctionApproximation::evaluate() - LevMar expects "+toString(n)
|
---|
| 334 | +" outputs but we provide "+toString(output_data.size())+".");
|
---|
[c62f96] | 335 | if (!output_data.empty()) {
|
---|
[e1fe7e] | 336 | filtered_inputs_t::const_iterator initer = input_data.begin();
|
---|
[66cfc7] | 337 | outputs_t::const_iterator outiter = output_data.begin();
|
---|
| 338 | size_t index = 0;
|
---|
[5b5724] | 339 | for (; initer != input_data.end(); ++initer, ++outiter) {
|
---|
[66cfc7] | 340 | // result may be a vector, calculate L2 norm
|
---|
| 341 | const FunctionModel::results_t functionvalue =
|
---|
| 342 | model(*initer);
|
---|
[5b5724] | 343 | x[index++] = functionvalue[0];
|
---|
| 344 | // std::vector<double> differences(functionvalue.size(), 0.);
|
---|
| 345 | // std::transform(
|
---|
| 346 | // functionvalue.begin(), functionvalue.end(), outiter->begin(),
|
---|
| 347 | // differences.begin(),
|
---|
| 348 | // &SquaredDifference);
|
---|
| 349 | // x[index] = std::accumulate(differences.begin(), differences.end(), 0.);
|
---|
| 350 | }
|
---|
| 351 | }
|
---|
| 352 | }
|
---|
| 353 |
|
---|
| 354 | void FunctionApproximation::evaluateDerivative(double *p, double *jac, int m, int n, void *data)
|
---|
| 355 | {
|
---|
| 356 | // first set parameters
|
---|
| 357 | prepareModel(p,m);
|
---|
| 358 |
|
---|
| 359 | // then evaluate
|
---|
| 360 | ASSERT( (size_t)n == output_data.size(),
|
---|
| 361 | "FunctionApproximation::evaluateDerivative() - LevMar expects "+toString(n)
|
---|
| 362 | +" outputs but we provide "+toString(output_data.size())+".");
|
---|
| 363 | if (!output_data.empty()) {
|
---|
[e1fe7e] | 364 | filtered_inputs_t::const_iterator initer = input_data.begin();
|
---|
[5b5724] | 365 | outputs_t::const_iterator outiter = output_data.begin();
|
---|
| 366 | size_t index = 0;
|
---|
| 367 | for (; initer != input_data.end(); ++initer, ++outiter) {
|
---|
| 368 | // result may be a vector, calculate L2 norm
|
---|
| 369 | for (int paramindex = 0; paramindex < m; ++paramindex) {
|
---|
| 370 | const FunctionModel::results_t functionvalue =
|
---|
| 371 | model.parameter_derivative(*initer, paramindex);
|
---|
| 372 | jac[index++] = functionvalue[0];
|
---|
| 373 | }
|
---|
| 374 | // std::vector<double> differences(functionvalue.size(), 0.);
|
---|
| 375 | // std::transform(
|
---|
| 376 | // functionvalue.begin(), functionvalue.end(), outiter->begin(),
|
---|
| 377 | // differences.begin(),
|
---|
| 378 | // &SquaredDifference);
|
---|
| 379 | // x[index] = std::accumulate(differences.begin(), differences.end(), 0.);
|
---|
[66cfc7] | 380 | }
|
---|
| 381 | }
|
---|
| 382 | }
|
---|