Index: src/FunctionApproximation/FunctionApproximation.cpp
===================================================================
--- src/FunctionApproximation/FunctionApproximation.cpp	(revision 66cfc791a31134a511bad433ec78b3e5532d5106)
+++ src/FunctionApproximation/FunctionApproximation.cpp	(revision 66cfc791a31134a511bad433ec78b3e5532d5106)
@@ -0,0 +1,257 @@
+/*
+ * Project: MoleCuilder
+ * Description: creates and alters molecular systems
+ * Copyright (C)  2012 University of Bonn. All rights reserved.
+ * Please see the COPYING file or "Copyright notice" in builder.cpp for details.
+ *
+ *
+ *   This file is part of MoleCuilder.
+ *
+ *    MoleCuilder is free software: you can redistribute it and/or modify
+ *    it under the terms of the GNU General Public License as published by
+ *    the Free Software Foundation, either version 2 of the License, or
+ *    (at your option) any later version.
+ *
+ *    MoleCuilder is distributed in the hope that it will be useful,
+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *    GNU General Public License for more details.
+ *
+ *    You should have received a copy of the GNU General Public License
+ *    along with MoleCuilder.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/*
+ * FunctionApproximation.cpp
+ *
+ *  Created on: 02.10.2012
+ *      Author: heber
+ */
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "CodePatterns/MemDebug.hpp"
+
+#include "FunctionApproximation.hpp"
+
+#include <algorithm>
+#include <boost/bind.hpp>
+#include <boost/function.hpp>
+#include <iostream>
+#include <iterator>
+#include <numeric>
+#include <sstream>
+
+#include <levmar.h>
+
+#include "CodePatterns/Assert.hpp"
+#include "CodePatterns/Log.hpp"
+
+#include "FunctionApproximation/FunctionModel.hpp"
+
+void FunctionApproximation::setTrainingData(const inputs_t &input, const outputs_t &output)
+{
+  ASSERT( input.size() == output.size(),
+      "FunctionApproximation::setTrainingData() - the number of input and output tuples differ: "+toString(input.size())+"!="
+      +toString(output.size())+".");
+  if (input.size() != 0) {
+    ASSERT( input[0].size() == input_dimension,
+        "FunctionApproximation::setTrainingData() - the dimension of the input tuples and input dimension differ: "+toString(input[0].size())+"!="
+        +toString(input_dimension)+".");
+    input_data = input;
+    ASSERT( output[0].size() == output_dimension,
+        "FunctionApproximation::setTrainingData() - the dimension of the output tuples and output dimension differ: "+toString(output[0].size())+"!="
+        +toString(output_dimension)+".");
+    output_data = output;
+  } else {
+    ELOG(2, "Given vectors of training data are empty, clearing internal vectors accordingly.");
+    input_data.clear();
+    output_data.clear();
+  }
+}
+
+void FunctionApproximation::setModelFunction(FunctionModel &_model)
+{
+  model= _model;
+}
+
+
+/* Meyer's (reformulated) problem, minimum at (2.48, 6.18, 3.45) */
+void meyer(double *p, double *x, int m, int n, void *data)
+{
+register int i;
+double ui;
+
+        for(i=0; i<n; ++i){
+                ui=0.45+0.05*i;
+                x[i]=p[0]*exp(10.0*p[1]/(ui+p[2]) - 13.0);
+        }
+}
+
+void jacmeyer(double *p, double *jac, int m, int n, void *data)
+{
+register int i, j;
+double ui, tmp;
+
+  for(i=j=0; i<n; ++i){
+          ui=0.45+0.05*i;
+          tmp=exp(10.0*p[1]/(ui+p[2]) - 13.0);
+
+          jac[j++]=tmp;
+          jac[j++]=10.0*p[0]*tmp/(ui+p[2]);
+          jac[j++]=-10.0*p[0]*p[1]*tmp/((ui+p[2])*(ui+p[2]));
+  }
+}
+
+/** Callback to circumvent boost::bind, boost::function and function pointer problem.
+ *
+ * See here (second answer!) to the nature of the problem:
+ * http://stackoverflow.com/questions/282372/demote-boostfunction-to-a-plain-function-pointer
+ *
+ * We cannot use a boost::bind bounded boost::function as a function pointer.
+ * boost::function::target() will just return NULL because the signature does not
+ * match. We have to use a C-style callback function and our luck is that
+ * the levmar signature provides for a void* additional data pointer which we
+ * can cast back to our FunctionApproximation class, as we need access to the
+ * data contained, e.g. the FunctionModel reference FunctionApproximation::model.
+ *
+ */
+void FunctionApproximation::LevMarCallback(double *p, double *x, int m, int n, void *data)
+{
+  FunctionApproximation *approximator = static_cast<FunctionApproximation *>(data);
+  ASSERT( approximator != NULL,
+      "LevMarCallback() - received data does not represent a FunctionApproximation object.");
+  boost::function<void(double*,double*,int,int,void*)> function =
+      boost::bind(&FunctionApproximation::evaluate, approximator, _1, _2, _3, _4, _5);
+  function(p,x,m,n,data);
+}
+
+void FunctionApproximation::operator()()
+{
+  // let levmar optimize
+  register int i, j;
+  int ret;
+  double *p;
+  double *x; // we set zero vector by giving NULL
+  int m, n;
+  double opts[LM_OPTS_SZ], info[LM_INFO_SZ];
+
+  opts[0]=LM_INIT_MU; opts[1]=1E-15; opts[2]=1E-15; opts[3]=1E-20;
+  opts[4]= LM_DIFF_DELTA; // relevant only if the Jacobian is approximated using finite differences; specifies forward differencing
+  //opts[4]=-LM_DIFF_DELTA; // specifies central differencing to approximate Jacobian; more accurate but more expensive to compute!
+
+  m = model.getParameterDimension();
+  n = 1;
+  const FunctionModel::parameters_t params = model.getParameters();
+  p = new double(m);
+  x = new double(n);
+  size_t index = 0;
+  for(FunctionModel::parameters_t::const_iterator paramiter = params.begin();
+      paramiter != params.end();
+      ++paramiter, ++index) {
+    p[index] = params[index];
+  }
+  x[0] = 0.;
+
+  {
+    double *work, *covar;
+    work=(double *)malloc((LM_DIF_WORKSZ(m, n)+m*m)*sizeof(double));
+    if(!work){
+      ELOG(0, "FunctionApproximation::operator() - memory allocation request failed.");
+      return;
+    }
+    covar=work+LM_DIF_WORKSZ(m, n);
+
+    // give this pointer as additional data to construct function pointer in
+    // LevMarCallback and call
+    ret=dlevmar_dif(&FunctionApproximation::LevMarCallback, p, x, m, n, 1000, opts, info, work, covar, this); // no Jacobian, caller allocates work memory, covariance estimated
+
+    {
+      std::stringstream covar_msg;
+      covar_msg << "Covariance of the fit:\n";
+      for(i=0; i<m; ++i){
+        for(j=0; j<m; ++j)
+          covar_msg << covar[i*m+j] << " ";
+        covar_msg << std::endl;
+      }
+      covar_msg << std::endl;
+      LOG(1, "INFO: " << covar_msg.str());
+    }
+
+    free(work);
+  }
+
+//  {
+//   double err[16];
+//   dlevmar_chkjac(meyer, jacmeyer, p, m, n, NULL, err);
+//   for(i=0; i<n; ++i) printf("gradient %d, err %g\n", i, err[i]);
+//  }
+
+  {
+   std::stringstream result_msg;
+   result_msg << "Levenberg-Marquardt returned " << ret << " in " << info[5] << " iter, reason " << info[6] << "\nSolution: ";
+   for(i=0; i<m; ++i)
+     result_msg << p[i] << " ";
+   result_msg << "\n\nMinimization info:\n";
+   std::vector<std::string> infonames(LM_INFO_SZ);
+   infonames[0] = std::string("||e||_2 at initial p");
+   infonames[1] = std::string("||e||_2");
+   infonames[2] = std::string("||J^T e||_inf");
+   infonames[3] = std::string("||Dp||_2");
+   infonames[4] = std::string("mu/max[J^T J]_ii");
+   infonames[5] = std::string("# iterations");
+   infonames[6] = std::string("reason for termination");
+   infonames[7] = std::string(" # function evaluations");
+   infonames[8] = std::string(" # Jacobian evaluations");
+   infonames[9] = std::string(" # linear systems solved");
+   for(i=0; i<LM_INFO_SZ; ++i)
+      result_msg << infonames[i] << ": " << info[i] << " ";
+    result_msg << std::endl;
+    LOG(1, "INFO: " << result_msg.str());
+  }
+
+  delete p;
+  delete x;
+}
+
+double SquaredDifference(const double res1, const double res2)
+{
+  return (res1-res2)*(res1-res2);
+}
+
+void FunctionApproximation::evaluate(double *p, double *x, int m, int n, void *data)
+{
+  // first set parameters
+  ASSERT( (size_t)m == model.getParameterDimension(),
+      "FunctionApproximation::evaluate() - LevMar expects "+toString(m)
+      +" parameters but the model function expects "+toString(model.getParameterDimension())+".");
+  ASSERT( n == 1,
+      "FunctionApproximation::evaluate() - we return precisely a single value, LevMar expects something else.");
+  FunctionModel::parameters_t params(m, 0.);
+  std::copy(p, p+m, params.begin());
+  model.setParameters(params);
+
+  // then evaluate
+  if (!input_data.empty()) {
+    inputs_t::const_iterator initer = input_data.begin();
+    outputs_t::const_iterator outiter = output_data.begin();
+    size_t index = 0;
+    for (; initer != input_data.end(); ++initer, ++outiter, ++index) {
+      // result may be a vector, calculate L2 norm
+      const FunctionModel::results_t functionvalue =
+          model(*initer);
+      std::vector<double> differences(functionvalue.size(), 0.);
+      std::transform(
+          functionvalue.begin(), functionvalue.end(), outiter->begin(),
+          differences.begin(),
+          &SquaredDifference);
+      x[index] = std::accumulate(differences.begin(), differences.end(), 0.);
+    }
+  } else {
+    for (register int i=0;i<n;++i)
+      x[i] = 0.;
+  }
+}
Index: src/FunctionApproximation/FunctionApproximation.hpp
===================================================================
--- src/FunctionApproximation/FunctionApproximation.hpp	(revision 66cfc791a31134a511bad433ec78b3e5532d5106)
+++ src/FunctionApproximation/FunctionApproximation.hpp	(revision 66cfc791a31134a511bad433ec78b3e5532d5106)
@@ -0,0 +1,122 @@
+/*
+ * FunctionApproximation.hpp
+ *
+ *  Created on: 02.10.2012
+ *      Author: heber
+ */
+
+#ifndef FUNCTIONAPPROXIMATION_HPP_
+#define FUNCTIONAPPROXIMATION_HPP_
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <vector>
+
+#include "FunctionApproximation/FunctionModel.hpp"
+
+/** This class encapsulates the solution to approximating a high-dimensional
+ * function represented by two vectors of tuples, being input variables and
+ * output of the function via a model function, manipulated by a set of
+ * parameters.
+ *
+ * \note For this reason the input and output dimension has to be given in
+ * the constructor since these are fixed parameters to the problem as a
+ * whole and usually: a different input dimension means we have a completely
+ * different problem (and hence we may as well construct and new instance of
+ * this class).
+ *
+ * The "training data", i.e. the two sets of input and output values, is
+ * given extra.
+ *
+ * The problem is then that a given high-dimensional function is supplied,
+ * the "model", and we have to fit this function via its set of variable
+ * parameters. This fitting procedure is executed via a Levenberg-Marquardt
+ * algorithm as implemented in the
+ * <a href="http://www.ics.forth.gr/~lourakis/levmar/index.html">LevMar</a>
+ * package.
+ *
+ */
+class FunctionApproximation
+{
+public:
+  //!> typedef for a vector of input arguments
+  typedef std::vector<FunctionModel::arguments_t> inputs_t;
+  //!> typedef for a vector of output values
+  typedef std::vector<FunctionModel::results_t> outputs_t;
+public:
+  /** Constructor of the class FunctionApproximation.
+   *
+   * \param _input_dimension input dimension for this function approximation
+   * \param _output_dimension output dimension for this function approximation
+   */
+  FunctionApproximation(
+      const size_t &_input_dimension,
+      const size_t &_output_dimension,
+      FunctionModel &_model) :
+    input_dimension(_input_dimension),
+    output_dimension(_output_dimension),
+    model(_model)
+  {}
+  /** Destructor for class FunctionApproximation.
+   *
+   */
+  ~FunctionApproximation()
+  {}
+
+  /** Setter for the training data to be used.
+   *
+   * \param input vector of input tuples, needs to be of
+   *        FunctionApproximation::input_dimension size
+   * \param output vector of output tuples, needs to be of
+   *        FunctionApproximation::output_dimension size
+   */
+  void setTrainingData(const inputs_t &input, const outputs_t &output);
+
+  /** Setter for the model function to be used in the approximation.
+   *
+   */
+  void setModelFunction(FunctionModel &_model);
+
+  /** This starts the fitting process, resulting in the parameters to
+   * the model function being optimized with respect to the given training
+   * data.
+   */
+  void operator()();
+
+  /** Evaluates the model function for each pair of training tuple and returns
+   * the error between the output of the function and the expected output as a
+   * vector.
+   *
+   * This function as a signature compatible to the one required by the
+   * LevMar package (with double precision).
+   *
+   * \param *p array of parameters for the model function of dimension \a m
+   * \param *x array of result values of dimension \a n
+   * \param m parameter dimension
+   * \param n output dimension
+   * \param *data additional data, unused here
+   */
+  void evaluate(double *p, double *x, int m, int n, void *data);
+
+private:
+  static void LevMarCallback(double *p, double *x, int m, int n, void *data);
+
+private:
+  //!> input dimension (is fixed from construction)
+  const size_t input_dimension;
+  //!> output dimension (is fixed from construction)
+  const size_t output_dimension;
+
+  //!> current input set of training data
+  inputs_t input_data;
+  //!> current output set of training data
+  outputs_t output_data;
+
+  //!> the model function to be used in the high-dimensional approximation
+  FunctionModel &model;
+};
+
+#endif /* FUNCTIONAPPROXIMATION_HPP_ */
Index: src/FunctionApproximation/FunctionArgument.hpp
===================================================================
--- src/FunctionApproximation/FunctionArgument.hpp	(revision 66cfc791a31134a511bad433ec78b3e5532d5106)
+++ src/FunctionApproximation/FunctionArgument.hpp	(revision 66cfc791a31134a511bad433ec78b3e5532d5106)
@@ -0,0 +1,34 @@
+/*
+ * FunctionArgument.hpp
+ *
+ *  Created on: 02.10.2012
+ *      Author: heber
+ */
+
+#ifndef FUNCTIONARGUMENT_HPP_
+#define FUNCTIONARGUMENT_HPP_
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <utility>
+
+/** This class encapsulates all information with respect to a single argument
+ *  for a high-dimensional model function.
+ *
+ *  We restrict ourselves here to a function that dependent on a set of
+ *  three-dimensional vectors, i.e. a set of positions in space. And for
+ *  the moment to distances in between these sets.
+ *
+ */
+struct argument_t
+{
+  //!> indices between which the distance is given
+  std::pair<size_t, size_t> indices;
+  //!> distance
+  double distance;
+};
+
+#endif /* FUNCTIONARGUMENT_HPP_ */
Index: src/FunctionApproximation/FunctionModel.hpp
===================================================================
--- src/FunctionApproximation/FunctionModel.hpp	(revision 66cfc791a31134a511bad433ec78b3e5532d5106)
+++ src/FunctionApproximation/FunctionModel.hpp	(revision 66cfc791a31134a511bad433ec78b3e5532d5106)
@@ -0,0 +1,83 @@
+/*
+ * FunctionModel.hpp
+ *
+ *  Created on: 02.10.2012
+ *      Author: heber
+ */
+
+#ifndef FUNCTIONMODEL_HPP_
+#define FUNCTIONMODEL_HPP_
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <vector>
+
+#include "FunctionApproximation/FunctionArgument.hpp"
+
+/** This class represents the interface for a given function to model a
+ * high-dimensional data set in FunctionApproximation.
+ *
+ * As the parameters may be stored differently, the interface functions for
+ * getting and setting them are as light-weight (and not speed-optimized)
+ * as possible.
+ *
+ */
+class FunctionModel
+{
+public:
+  //!> typedef for a single parameter degree of freedom of the function
+  typedef double parameter_t;
+  //!> typedef for the whole set of parameters of the function
+  typedef std::vector<parameter_t> parameters_t;
+  //!> typedef for the argument vector as input to the function
+  typedef std::vector<argument_t> arguments_t;
+  //!> typedef for a single result degree of freedom
+  typedef double result_t;
+  //!> typedef for the result vector as returned by the function
+  typedef std::vector<result_t> results_t;
+  //!> typedef for the components of the derivative
+  typedef std::vector<result_t> derivative_components_t;
+public:
+  FunctionModel() {}
+  virtual ~FunctionModel() {}
+
+  /** Setter for the parameters of the model function.
+   *
+   * \param params set of parameters to set
+   */
+  virtual void setParameters(const parameters_t &params)=0;
+
+  /** Getter for the parameters of this model function.
+   *
+   * \return current set of parameters of the model function
+   */
+  virtual parameters_t getParameters() const=0;
+
+  /** Getter for the number of parameters of this model function.
+   *
+   * \return number of parameters
+   */
+  virtual size_t getParameterDimension() const=0;
+
+  /** Evaluates the function with the given \a arguments and the current set of
+   * parameters.
+   *
+   * \param arguments set of arguments as input variables to the function
+   * \return result of the function
+   */
+  virtual results_t operator()(const arguments_t &arguments) const=0;
+
+  /** Evaluates the derivative of the function with the given \a arguments
+   * for each component.
+   *
+   * \param arguments set of arguments as input variables to the function
+   * \return result vector containing the derivative with respect to each
+   *         input comonent of the function
+   */
+  virtual derivative_components_t derivative(const arguments_t &arguments) const=0;
+};
+
+#endif /* FUNCTIONMODEL_HPP_ */
Index: src/FunctionApproximation/Makefile.am
===================================================================
--- src/FunctionApproximation/Makefile.am	(revision 66cfc791a31134a511bad433ec78b3e5532d5106)
+++ src/FunctionApproximation/Makefile.am	(revision 66cfc791a31134a511bad433ec78b3e5532d5106)
@@ -0,0 +1,65 @@
+# PLEASE adhere to the alphabetical ordering in this Makefile!
+# Also indentation by a single tab
+
+FUNCTIONAPPROXIMATIONSOURCE =
+	
+FUNCTIONAPPROXIMATIONHEADER = \
+	FunctionApproximation/FunctionModel.hpp \
+	FunctionApproximation/FunctionArgument.hpp
+
+if CONDLEVMAR
+FUNCTIONAPPROXIMATIONSOURCE += \
+	FunctionApproximation/FunctionApproximation.cpp
+FUNCTIONAPPROXIMATIONHEADER += \
+	FunctionApproximation/FunctionApproximation.hpp
+endif
+
+noinst_LTLIBRARIES += libMolecuilderFunctionApproximation.la
+libMolecuilderFunctionApproximation_la_includedir = $(includedir)/MoleCuilder/FunctionApproximation/
+libMolecuilderFunctionApproximation_la_CPPFLAGS = $(AM_CPPFLAGS)
+libMolecuilderFunctionApproximation_la_LDFLAGS = $(AM_LDFLAGS)
+libMolecuilderFunctionApproximation_la_LIBADD =
+if CONDLEVMAR
+libMolecuilderFunctionApproximation_la_CPPFLAGS += $(AM_CPPFLAGS) $(LEVMAR_CPPFLAGS)
+libMolecuilderFunctionApproximation_la_LDFLAGS += $(LEVMAR_LDFLAGS)
+libMolecuilderFunctionApproximation_la_LIBADD += \
+	$(LEVMAR_LIBS)
+endif
+
+nobase_libMolecuilderFunctionApproximation_la_include_HEADERS = ${FUNCTIONAPPROXIMATIONHEADER}
+
+## Define the source file list for the "libexample-@MOLECUILDER_API_VERSION@.la"
+## target.  Note that @MOLECUILDER_API_VERSION@ is not interpreted by Automake and
+## will therefore be treated as if it were literally part of the target name,
+## and the variable name derived from that.
+## The file extension .cc is recognized by Automake, and makes it produce
+## rules which invoke the C++ compiler to produce a libtool object file (.lo)
+## from each source file.  Note that it is not necessary to list header files
+## which are already listed elsewhere in a _HEADERS variable assignment.
+libMolecuilderFunctionApproximation_la_SOURCES = ${FUNCTIONAPPROXIMATIONSOURCE}
+
+## Instruct libtool to include ABI version information in the generated shared
+## library file (.so).  The library ABI version is defined in configure.ac, so
+## that all version information is kept in one place.
+#libMolecuilderFunctionApproximation_la_LDFLAGS = -version-info $(MOLECUILDER_SO_VERSION)
+
+## The generated configuration header is installed in its own subdirectory of
+## $(libdir).  The reason for this is that the configuration information put
+## into this header file describes the target platform the installed library
+## has been built for.  Thus the file must not be installed into a location
+## intended for architecture-independent files, as defined by the Filesystem
+## Hierarchy Standard (FHS).
+## The nodist_ prefix instructs Automake to not generate rules for including
+## the listed files in the distribution on 'make dist'.  Files that are listed
+## in _HEADERS variables are normally included in the distribution, but the
+## configuration header file is generated at configure time and should not be
+## shipped with the source tarball.
+#libMolecuilderFunctionApproximation_libincludedir = $(libdir)/MoleCuilder/include
+#nodist_libMolecuilderFunctionApproximation_libinclude_HEADERS = $(top_builddir)/libmolecuilder_config.h
+
+## Install the generated pkg-config file (.pc) into the expected location for
+## architecture-dependent package configuration information.  Occasionally,
+## pkg-config files are also used for architecture-independent data packages,
+## in which case the correct install location would be $(datadir)/pkgconfig.
+#pkgconfigdir = $(libdir)/pkgconfig
+#pkgconfig_DATA = $(top_builddir)/MoleCuilder.pc
Index: src/Makefile.am
===================================================================
--- src/Makefile.am	(revision dda3e8e0ed52a27bea7146869d59f76442212c40)
+++ src/Makefile.am	(revision 66cfc791a31134a511bad433ec78b3e5532d5106)
@@ -19,4 +19,5 @@
 include Fragmentation/Summation/Makefile.am
 include Fragmentation/SetValues/Makefile.am
+include FunctionApproximation/Makefile.am
 include Graph/Makefile.am
 include Helpers/Makefile.am
