Index: src/Actions/FragmentationAction/FitPotentialAction.cpp
===================================================================
--- src/Actions/FragmentationAction/FitPotentialAction.cpp	(revision 48d20d12c40193473e911ec536adb1dfb3d14bff)
+++ src/Actions/FragmentationAction/FitPotentialAction.cpp	(revision 48d20d12c40193473e911ec536adb1dfb3d14bff)
@@ -0,0 +1,285 @@
+/*
+ * Project: MoleCuilder
+ * Description: creates and alters molecular systems
+ * Copyright (C)  2013 University of Bonn. All rights reserved.
+ * 
+ *
+ *   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/>.
+ */
+
+/*
+ * FitPotentialAction.cpp
+ *
+ *  Created on: Apr 09, 2013
+ *      Author: heber
+ */
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+// needs to come before MemDebug due to placement new
+#include <boost/archive/text_iarchive.hpp>
+
+#include "CodePatterns/MemDebug.hpp"
+
+#include <algorithm>
+#include <boost/bind.hpp>
+#include <boost/filesystem.hpp>
+#include <boost/foreach.hpp>
+#include <map>
+#include <string>
+#include <sstream>
+
+#include "Actions/FragmentationAction/FitPotentialAction.hpp"
+
+#include "CodePatterns/Log.hpp"
+
+#include "Element/element.hpp"
+#include "Fragmentation/Homology/HomologyContainer.hpp"
+#include "Fragmentation/Homology/HomologyGraph.hpp"
+#include "Fragmentation/Summation/SetValues/Fragment.hpp"
+#include "FunctionApproximation/Extractors.hpp"
+#include "FunctionApproximation/FunctionApproximation.hpp"
+#include "FunctionApproximation/FunctionModel.hpp"
+#include "FunctionApproximation/TrainingData.hpp"
+#include "FunctionApproximation/writeDistanceEnergyTable.hpp"
+#include "Potentials/PotentialFactory.hpp"
+#include "Potentials/SerializablePotential.hpp"
+
+using namespace MoleCuilder;
+
+// and construct the stuff
+#include "FitPotentialAction.def"
+#include "Action_impl_pre.hpp"
+/** =========== define the function ====================== */
+
+HomologyGraph getFirstGraphwithSpecifiedElements(
+    const HomologyContainer &homologies,
+    const SerializablePotential::ParticleTypes_t &types)
+{
+  ASSERT( !types.empty(),
+      "getFirstGraphwithSpecifiedElements() - charges is empty?");
+  // create charges
+  Fragment::charges_t charges;
+  charges.resize(types.size());
+  std::transform(types.begin(), types.end(),
+      charges.begin(), boost::lambda::_1);
+  // convert into count map
+  Extractors::elementcounts_t counts_per_charge =
+      Extractors::_detail::getElementCounts(charges);
+  ASSERT( !counts_per_charge.empty(),
+      "getFirstGraphwithSpecifiedElements() - charge counts are empty?");
+  LOG(2, "DEBUG: counts_per_charge is " << counts_per_charge << ".");
+  // we want to check each (unique) key only once
+  HomologyContainer::const_key_iterator olditer = homologies.key_end();
+  for (HomologyContainer::const_key_iterator iter =
+      homologies.key_begin(); iter != homologies.key_end(); olditer = iter++) {
+    // if it's the same as the old one, skip it
+    if (*olditer == *iter)
+      continue;
+    // if it's a new key, check if every element has the right number of counts
+    Extractors::elementcounts_t::const_iterator countiter = counts_per_charge.begin();
+    for (; countiter != counts_per_charge.end(); ++countiter)
+      if (!(*iter).hasTimesAtomicNumber(countiter->first,countiter->second))
+        break;
+    if( countiter == counts_per_charge.end())
+      return *iter;
+  }
+  return HomologyGraph();
+}
+
+Action::state_ptr FragmentationFitPotentialAction::performCall() {
+  // charges specify the potential type
+  SerializablePotential::ParticleTypes_t chargenumbers;
+  {
+    const std::vector<const element *> &charges = params.charges.get();
+    std::transform(charges.begin(), charges.end(), std::back_inserter(chargenumbers),
+        boost::bind(&element::getAtomicNumber, _1));
+  }
+  // fragment specifies the homology fragment to use
+  SerializablePotential::ParticleTypes_t fragmentnumbers;
+  {
+    const std::vector<const element *> &fragment = params.fragment.get();
+    std::transform(fragment.begin(), fragment.end(), std::back_inserter(fragmentnumbers),
+        boost::bind(&element::getAtomicNumber, _1));
+  }
+
+  // parse homologies into container
+  HomologyContainer homologies;
+  if (boost::filesystem::exists(params.homology_file.get())) {
+    std::ifstream returnstream(params.homology_file.get().string().c_str());
+    if (returnstream.good()) {
+      boost::archive::text_iarchive ia(returnstream);
+      ia >> homologies;
+    } else {
+      ELOG(0, "Failed to parse from " << params.homology_file.get().string() << ".");
+      return Action::failure;
+    }
+    returnstream.close();
+  } else {
+    ELOG(0, params.homology_file.get() << " does not exist.");
+    return Action::failure;
+  }
+
+  // first we try to look into the HomologyContainer
+  LOG(1, "INFO: Listing all present homologies ...");
+  for (HomologyContainer::container_t::const_iterator iter =
+      homologies.begin(); iter != homologies.end(); ++iter) {
+    LOG(1, "INFO: graph " << iter->first << " has Fragment " << iter->second.first
+        << " and associated energy " << iter->second.second << ".");
+  }
+
+  LOG(0, "STATUS: I'm training now a " << params.potentialtype.get() << " potential on charges "
+      << chargenumbers << " on data from " << params.homology_file.get() << ".");
+
+  /******************** TRAINING ********************/
+  // fit potential
+  FunctionModel *model =
+      PotentialFactory::getInstance().createInstance(
+          params.potentialtype.get(),
+          chargenumbers);
+  ASSERT( model != NULL,
+      "main() - model returned from PotentialFactory is NULL.");
+  FunctionModel::parameters_t bestparams(model->getParameterDimension(), 0.);
+  {
+    // then we ought to pick the right HomologyGraph ...
+    const HomologyGraph graph = getFirstGraphwithSpecifiedElements(homologies,fragmentnumbers);
+    if (graph != HomologyGraph()) {
+      LOG(1, "First representative graph containing fragment "
+          << fragmentnumbers << " is " << graph << ".");
+
+      // Afterwards we go through all of this type and gather the distance and the energy value
+      TrainingData data(model->getFragmentSpecificExtractor());
+      data(homologies.getHomologousGraphs(graph));
+
+      // print distances and energies if desired for debugging
+      if (!data.getTrainingInputs().empty()) {
+        // print which distance is which
+        size_t counter=1;
+        if (DoLog(3)) {
+          const FunctionModel::arguments_t &inputs = data.getTrainingInputs()[0];
+          for (FunctionModel::arguments_t::const_iterator iter = inputs.begin();
+              iter != inputs.end(); ++iter) {
+            const argument_t &arg = *iter;
+            LOG(3, "DEBUG: distance " << counter++ << " is between (#"
+                << arg.indices.first << "c" << arg.types.first << ","
+                << arg.indices.second << "c" << arg.types.second << ").");
+          }
+        }
+
+        // print table
+        LOG(3, "DEBUG: I gathered the following training data:\n" <<
+            _detail::writeDistanceEnergyTable(data.getDistanceEnergyTable()));
+      }
+
+      // now perform the function approximation by optimizing the model function
+      FunctionApproximation approximator(data, *model);
+      if (model->isBoxConstraint() && approximator.checkParameterDerivatives()) {
+        double l2error = std::numeric_limits<double>::infinity();
+        // seed with current time
+        srand((unsigned)time(0));
+        for (unsigned int runs=0; runs < params.best_of_howmany.get(); ++runs) {
+          // generate new random initial parameter values
+          model->setParametersToRandomInitialValues(data);
+          LOG(1, "INFO: Initial parameters of run " << runs << " are "
+              << model->getParameters() << ".");
+          approximator(FunctionApproximation::ParameterDerivative);
+          LOG(1, "INFO: Final parameters of run " << runs << " are "
+              << model->getParameters() << ".");
+          const double new_l2error = data.getL2Error(*model);
+          if (new_l2error < l2error) {
+            // store currently best parameters
+            l2error = new_l2error;
+            bestparams = model->getParameters();
+            LOG(1, "STATUS: New fit from run " << runs
+                << " has better error of " << l2error << ".");
+          }
+        }
+        // reset parameters from best fit
+        model->setParameters(bestparams);
+        LOG(1, "INFO: Best parameters with L2 error of "
+            << l2error << " are " << model->getParameters() << ".");
+      } else {
+        ELOG(0, "We require parameter derivatives for a box constraint minimization.");
+        return Action::failure;
+      }
+
+      // create a map of each fragment with error.
+      typedef std::multimap< double, size_t > WorseFragmentMap_t;
+      WorseFragmentMap_t WorseFragmentMap;
+      HomologyContainer::range_t fragmentrange = homologies.getHomologousGraphs(graph);
+      // fragments make it into the container in reversed order, hence count from top down
+      size_t index= std::distance(fragmentrange.first, fragmentrange.second)-1;
+      for (HomologyContainer::const_iterator iter = fragmentrange.first;
+          iter != fragmentrange.second;
+          ++iter) {
+        const Fragment& fragment = iter->second.first;
+        const double &energy = iter->second.second;
+
+        // create arguments from the fragment
+        FunctionModel::extractor_t extractor = model->getFragmentSpecificExtractor();
+        FunctionModel::arguments_t args = extractor(fragment, 1);
+
+        // calculate value from potential
+        const double fitvalue = (*model)(args)[0];
+
+        // insert difference into map
+        const double error = fabs(energy - fitvalue);
+        WorseFragmentMap.insert( std::make_pair( error, index-- ) );
+
+        {
+          // give only the distances in the debugging text
+          std::stringstream streamargs;
+          BOOST_FOREACH (argument_t arg, args) {
+            streamargs << " " << arg.distance*AtomicLengthToAngstroem;
+          }
+          LOG(2, "DEBUG: frag.#" << index+1 << "'s error is |" << energy << " - " << fitvalue
+              << "| = " << error << " for args " << streamargs.str() << ".");
+        }
+      }
+      LOG(0, "RESULT: WorstFragmentMap " << WorseFragmentMap << ".");
+
+      SerializablePotential *potential = dynamic_cast<SerializablePotential *>(model);
+      if (potential != NULL) {
+        LOG(1, "STATUS: Resulting parameters are " << std::endl << *potential);
+      } else {
+        LOG(1, "INFO: FunctionModel is no serializable potential.");
+      }
+    }
+  }
+  delete model;
+
+  return Action::success;
+}
+
+Action::state_ptr FragmentationFitPotentialAction::performUndo(Action::state_ptr _state) {
+  return Action::success;
+}
+
+Action::state_ptr FragmentationFitPotentialAction::performRedo(Action::state_ptr _state){
+  return Action::success;
+}
+
+bool FragmentationFitPotentialAction::canUndo() {
+  return false;
+}
+
+bool FragmentationFitPotentialAction::shouldUndo() {
+  return false;
+}
+/** =========== end of function ====================== */
Index: src/Actions/FragmentationAction/FitPotentialAction.def
===================================================================
--- src/Actions/FragmentationAction/FitPotentialAction.def	(revision 48d20d12c40193473e911ec536adb1dfb3d14bff)
+++ src/Actions/FragmentationAction/FitPotentialAction.def	(revision 48d20d12c40193473e911ec536adb1dfb3d14bff)
@@ -0,0 +1,47 @@
+/*
+ * FitPotentialAction.def
+ *
+ *  Created on: Apr 08, 2013
+ *      Author: heber
+ */
+
+// all includes and forward declarations necessary for non-integral types below
+#include <boost/filesystem/path.hpp>
+#include <string>
+#include <vector>
+
+#include "Parameters/Validators/RangeValidator.hpp"
+#include "Parameters/Validators/STLVectorValidator.hpp"
+#include "Parameters/Validators/Specific/ElementValidator.hpp"
+#include "Parameters/Validators/Specific/FilePresentValidator.hpp"
+#include "Parameters/Validators/Specific/PotentialTypeValidator.hpp"
+
+// i.e. there is an integer with variable name Z that can be found in
+// ValueStorage by the token "Z" -> first column: int, Z, "Z"
+// "undefine" if no parameters are required, use (NOPARAM_DEFAULT) for each (undefined) default value
+#define paramtypes (boost::filesystem::path)(std::string)(std::vector<const element *>)(std::vector<const element *>)(unsigned int)
+#define paramtokens ("homology-file")("fit-potential")("potential-charges")("fragment-charges")("take-best-of")
+#define paramdescriptions ("homology file to parse")("potential type to fit")("charges specifying the potential")("charges specifying the fragment")("take the best among this many approximations")
+#define paramdefaults (NOPARAM_DEFAULT)(NOPARAM_DEFAULT)(NOPARAM_DEFAULT)(NOPARAM_DEFAULT)(PARAM_DEFAULT(3))
+#define paramreferences (homology_file)(potentialtype)(charges)(fragment)(best_of_howmany)
+#define paramvalids \
+(FilePresentValidator()) \
+(PotentialTypeValidator()) \
+(STLVectorValidator< std::vector<const element *> >(1,99, ElementValidator())) \
+(STLVectorValidator< std::vector<const element *> >(1,99, ElementValidator())) \
+(RangeValidator< unsigned int>(1,99))
+
+#undef statetypes
+#undef statereferences
+
+// some defines for all the names, you may use ACTION, STATE and PARAMS
+#define CATEGORY Fragmentation
+#define MENUNAME "fragmentation"
+#define MENUPOSITION 6
+#define ACTIONNAME FitPotential
+#define TOKEN "fit-potential"
+
+
+// finally the information stored in the ActionTrait specialization
+#define DESCRIPTION "fit an empirical potential to the given part of a specified fragment"
+#undef SHORTFORM
Index: src/Actions/FragmentationAction/FitPotentialAction.hpp
===================================================================
--- src/Actions/FragmentationAction/FitPotentialAction.hpp	(revision 48d20d12c40193473e911ec536adb1dfb3d14bff)
+++ src/Actions/FragmentationAction/FitPotentialAction.hpp	(revision 48d20d12c40193473e911ec536adb1dfb3d14bff)
@@ -0,0 +1,22 @@
+/*
+ * FitPotentialAction.hpp
+ *
+ *  Created on: Apr 09, 2013
+ *      Author: heber
+ */
+
+#ifndef FITPOTENTIALACTION_HPP_
+#define FITPOTENTIALACTION_HPP_
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+
+#include "Actions/Action.hpp"
+
+#include "FitPotentialAction.def"
+#include "Action_impl_header.hpp"
+
+#endif /* FITPOTENTIALACTION_HPP_ */
Index: src/Actions/GlobalListOfActions.hpp
===================================================================
--- src/Actions/GlobalListOfActions.hpp	(revision e2925fd7f1c9a7ab18fe71b5a40598b65ddb05dc)
+++ src/Actions/GlobalListOfActions.hpp	(revision 48d20d12c40193473e911ec536adb1dfb3d14bff)
@@ -122,6 +122,6 @@
   (SelectionNotAllShapes) \
   (SelectionNotShapeByName) \
+  (FragmentationAnalyseFragmentationResults) \
   (FragmentationFragmentation) \
-  (FragmentationAnalyseFragmentationResults) \
   (FragmentationStoreSaturatedFragment) \
   (FillRegularGrid) \
@@ -136,5 +136,5 @@
 // we need to append the automation action in case we have the JobMarket
 #ifdef HAVE_JOBMARKET
-#define GLOBALLISTOFACTIONS \
+#define GLOBALLISTOFACTIONS_JOBMARKET \
     BOOST_PP_SEQ_PUSH_BACK( \
         BOOST_PP_SEQ_PUSH_BACK( \
@@ -145,8 +145,22 @@
     )
 #else
-#define GLOBALLISTOFACTIONS \
+#define GLOBALLISTOFACTIONS_JOBMARKET \
     GLOBALLISTOFACTIONS_initial
 #endif /* HAVE_JOBMARKET */
 
+#ifdef HAVE_LEVMAR
+#define GLOBALLISTOFACTIONS_LEVMAR \
+    BOOST_PP_SEQ_PUSH_BACK( \
+        GLOBALLISTOFACTIONS_JOBMARKET, \
+        FragmentationFitPotential \
+    )
+#else
+#define GLOBALLISTOFACTIONS_LEVMAR \
+    GLOBALLISTOFACTIONS_JOBMARKET
+#endif /* HAVE_LEVMAR */
+
+// define final list
+#define GLOBALLISTOFACTIONS GLOBALLISTOFACTIONS_LEVMAR
+
 #endif /* GLOBALLISTOFACTIONS_HPP_ */
 
Index: src/Actions/Makefile.am
===================================================================
--- src/Actions/Makefile.am	(revision e2925fd7f1c9a7ab18fe71b5a40598b65ddb05dc)
+++ src/Actions/Makefile.am	(revision 48d20d12c40193473e911ec536adb1dfb3d14bff)
@@ -234,4 +234,14 @@
   Actions/FragmentationAction/ParseFragmentJobsAction.def
 endif
+
+if CONDLEVMAR
+FRAGMENTATIONACTIONSOURCE += \
+  Actions/FragmentationAction/FitPotentialAction.cpp
+FRAGMENTATIONACTIONHEADER += \
+  Actions/FragmentationAction/FitPotentialAction.hpp
+FRAGMENTATIONACTIONDEFS += \
+  Actions/FragmentationAction/FitPotentialAction.def
+endif
+
 
 GRAPHACTIONSOURCE = \
@@ -543,4 +553,7 @@
 libMolecuilderActionPrototypes_la_CPPFLAGS += $(JobMarket_CFLAGS)
 endif
+if CONDLEVMAR
+libMolecuilderActionPrototypes_la_CPPFLAGS += $(LEVMAR_CPPFLAGS)
+endif
 libMolecuilderActions_la_includedir = $(includedir)/MoleCuilder/
 libMolecuilderActionPrototypes_la_includedir = $(includedir)/MoleCuilder/
