Changes in / [c73e35:d8821e]


Ignore:
Files:
44 added
306 deleted
81 edited

Legend:

Unmodified
Added
Removed

  • Makefile.am

    rc73e35 rd8821e  
    1818        m4/cppunit.m4
    1919
    20 .PHONY: doc
     20.PHONY: doc extracheck installextracheck
    2121doc:
    2222        mkdir -p ${DX_DOCDIR}
    2323        cd doc && make doxygen-doc
    2424
     25extracheck:
     26        cd tests && make extracheck
     27installextracheck:
     28        cd tests && make installextracheck
     29
    2530unity:
    2631        cd src && make unity

  • configure.ac

    rc73e35 rd8821e  
    235235AC_MSG_RESULT($enable_jobmarket)
    236236AS_IF([test x"$enable_jobmarket" != x"no"],[
    237   # JobMarket library (needs priorizing of jobs)
    238   AM_PATH_JOBMARKET([1.1.4], $have_debug,[
     237  # JobMarket library (needs SystemCommandJob with suffix)
     238  AM_PATH_JOBMARKET([1.1.5], $have_debug,[
    239239    # the following is only required if we have JobMarket
    240240    BOOST_ASIO
    … …  
    409409        tests/Makefile])
    410410
     411AC_CONFIG_TESTDIR(tests/Calculations)
     412AC_CONFIG_FILES([
     413        tests/Calculations/atlocal
     414        tests/Calculations/Makefile])
     415AC_CONFIG_FILES([tests/Calculations/molecuilder], [chmod +x tests/Calculations/molecuilder])
     416
    411417AC_CONFIG_TESTDIR(tests/CodeChecks)
    412418AC_CONFIG_FILES([
    … …  
    418424        tests/Fragmentations/atlocal
    419425        tests/Fragmentations/Makefile])
    420 AC_CONFIG_FILES([tests/Fragmentations/analyzer], [chmod +x tests/Fragmentations/analyzer])
    421 AC_CONFIG_FILES([tests/Fragmentations/joiner], [chmod +x tests/Fragmentations/joiner])
    422426AC_CONFIG_FILES([tests/Fragmentations/molecuilder], [chmod +x tests/Fragmentations/molecuilder])
    423427

  • src/Actions/ActionQueue.cpp

    rc73e35 rd8821e  
    5858    AR(new ActionRegistry()),
    5959    history(new ActionHistory),
     60    CurrentAction(0),
    6061#ifndef HAVE_ACTION_THREAD
    61     CurrentAction(0)
     62    lastActionOk(true)
    6263#else
    63     CurrentAction(0),
     64    lastActionOk(true),
    6465    run_thread(boost::bind(&ActionQueue::run, this)),
    6566    run_thread_isIdle(true)
    … …  
    99100  try {
    100101    newaction->call();
     102    lastActionOk = true;
    101103  } catch(ActionFailureException &e) {
    102104    std::cerr << "Action " << *boost::get_error_info<ActionNameString>(e) << " has failed." << std::endl;
    103105    World::getInstance().setExitFlag(5);
     106    actionqueue.clear();
     107    tempqueue.clear();
     108    lastActionOk = false;
     109    std::cerr << "ActionQueue cleared." << std::endl;
    104110  }
    105111#else
    … …  
    157163        actionqueue[CurrentAction]->call();
    158164        pushStatus("SUCCESS: Action "+actionqueue[CurrentAction]->getName()+" successful.");
     165        lastActionOk = true;
    159166      } catch(ActionFailureException &e) {
    160167        pushStatus("FAIL: Action "+*boost::get_error_info<ActionNameString>(e)+" has failed.");
    161168        World::getInstance().setExitFlag(5);
     169        actionqueue.clear();
     170        tempqueue.clear();
     171        lastActionOk = false;
     172        std::cerr << "ActionQueue cleared." << std::endl;
     173        CurrentAction = (size_t)-1;
    162174      }
    163175      // access actionqueue, hence using mutex

  • src/Actions/ActionQueue.hpp

    rc73e35 rd8821e  
    123123  void redoLast();
    124124
     125  /** Return status of last executed action.
     126   *
     127   * \return true - action executed correctly, false - else
     128   */
     129  bool getLastActionOk() const
     130  {  return lastActionOk; }
     131
    125132#ifdef HAVE_ACTION_THREAD
    126133  boost::thread &getRunThread()
    … …  
    221228  ActionQueue_t tempqueue;
    222229
     230  //!> indicates that the last action has failed
     231  bool lastActionOk;
     232
    223233#ifdef HAVE_ACTION_THREAD
    224234  //!> internal thread to call Actions

  • src/Actions/CommandAction/ElementDbAction.cpp

    rc73e35 rd8821e  
    4141#include "config.hpp"
    4242#include "Element/element.hpp"    // we need element because of serialization!
     43#include "Element/ion.hpp"    // we need element because of serialization!
    4344#include "Element/periodentafel.hpp"
    4445#include "CodePatterns/Log.hpp"

  • src/Actions/FragmentationAction/AnalyseFragmentationResultsAction.cpp

    rc73e35 rd8821e  
    6363#include "Descriptors/AtomIdDescriptor.hpp"
    6464#include "Fragmentation/Summation/Containers/FragmentationChargeDensity.hpp"
    65 #include "Fragmentation/Summation/Containers/FragmentationLongRangeResults.hpp"
    6665#include "Fragmentation/Summation/Containers/FragmentationResultContainer.hpp"
    6766#include "Fragmentation/Summation/Containers/FragmentationShortRangeResults.hpp"
    … …  
    7877#include "Fragmentation/Summation/writeIndexedTable.hpp"
    7978#include "Fragmentation/Summation/writeTable.hpp"
    80 #ifdef HAVE_VMG
     79#if defined(HAVE_JOBMARKET) && defined(HAVE_VMG)
     80#include "Fragmentation/Summation/Containers/FragmentationLongRangeResults.hpp"
    8181#include "Fragmentation/Summation/Containers/VMGData.hpp"
    8282#include "Fragmentation/Summation/Containers/VMGDataFused.hpp"
    … …  
    232232
    233233
     234#if defined(HAVE_JOBMARKET) && defined(HAVE_VMG)
    234235/** Print long range energy from received results.
    235236 *
    … …  
    261262  }
    262263}
     264#endif
    263265
    264266void appendToHomologies(
    265267    const FragmentationShortRangeResults &shortrangeresults,
     268#if defined(HAVE_JOBMARKET) && defined(HAVE_VMG)
    266269    const FragmentationLongRangeResults &longrangeresults,
     270#endif
    267271    const bool storeGrids
    268272    )
    … …  
    311315    // only store sampled grids if desired
    312316    if (storeGrids) {
     317#if defined(HAVE_JOBMARKET) && defined(HAVE_VMG)
    313318      // obtain charge distribution
    314319      std::map<IndexSet::ptr, std::pair< MPQCDataGridMap_t, MPQCDataGridMap_t> >::const_iterator chargeiter
    … …  
    338343//          ++iter)
    339344//        *iter -= offset;
     345#else
     346      ELOG(1, "Long-range information in homology desired but long-range analysis capability not compiled in.");
     347#endif
    340348    }
    341349    values.insert( std::make_pair( graph, value) );
    … …  
    362370          << ", associated energy " << iter->second.energy;
    363371      if (iter->second.containsGrids)
     372#if defined(HAVE_JOBMARKET) && defined(HAVE_VMG)
    364373        output  << ", and sampled grid integral " << iter->second.charge_distribution.integral();
     374#else
     375        output << ", and there are sampled grids but capability not compiled in";
     376#endif
    365377      output << ".";
    366378      LOG(2, output.str());
    … …  
    610622  }
    611623
    612 #ifdef HAVE_VMG
     624#if defined(HAVE_JOBMARKET) && defined(HAVE_VMG)
    613625  if (DoLongrange) {
    614626    if ( World::getInstance().getAllAtoms().size() == 0) {
    … …  
    651663  }
    652664#else
    653   if (DoLongrange)
     665  if (DoLongrange) {
    654666    ELOG(2, "File contains long-range information but long-range analysis capability not compiled in.");
     667  }
    655668
    656669  // append all keysets to homology file with short-range info only (without grids)
    657   {
    658     std::map<JobId_t, VMGData> longrangeData;
    659     FragmentationLongRangeResults longrangeresults(
    660         shortrangedata,longrangeData,keysets, forcekeysets);
    661     appendToHomologies(shortrangeresults, longrangeresults, false);
    662   }
     670  appendToHomologies(shortrangeresults, false);
    663671#endif
    664672

  • src/Actions/FragmentationAction/FragmentationAction.cpp

    rc73e35 rd8821e  
    4141#include "Descriptors/AtomSelectionDescriptor.hpp"
    4242#include "Fragmentation/Exporters/ExportGraph_ToFiles.hpp"
    43 #ifdef HAVE_JOBMARKET
    4443#include "Fragmentation/Exporters/ExportGraph_ToJobs.hpp"
    45 #endif
    4644#include "Fragmentation/Fragmentation.hpp"
    4745#include "Fragmentation/Graph.hpp"
    … …  
    273271      exporter();
    274272    } else {
    275 #ifdef HAVE_JOBMARKET
    276273      // store molecule's fragment in FragmentJobQueue
    277274      ExportGraph_ToJobs exporter(TotalGraph, treatment, saturation);
    278275      exporter.setLevel(params.level.get());
    279276      exporter();
    280 #else
    281       STATUS("No output file types specified and JobMarket support is not compiled in.");
    282       return Action::failure;
    283 #endif
    284277    }
    285278  }

  • src/Actions/FragmentationAction/FragmentationAutomationAction.cpp

    rc73e35 rd8821e  
    5353#include "CodePatterns/Info.hpp"
    5454#include "CodePatterns/Log.hpp"
     55
     56#ifdef HAVE_JOBMARKET
    5557#include "JobMarket/Jobs/FragmentJob.hpp"
     58#else
     59#include "Jobs/JobMarket/FragmentJob.hpp"
     60#endif
    5661
    5762#include "Fragmentation/Automation/FragmentJobQueue.hpp"
     63#ifdef HAVE_JOBMARKET
    5864#include "Fragmentation/Automation/MPQCFragmentController.hpp"
     65#else
     66#include "Fragmentation/Automation/MPQCCommandFragmentController.hpp"
     67#endif
    5968#include "Fragmentation/Summation/Containers/FragmentationChargeDensity.hpp"
    6069#include "Fragmentation/Summation/Containers/FragmentationLongRangeResults.hpp"
    … …  
    6372#include "Fragmentation/Summation/Containers/MPQCData.hpp"
    6473#include "Fragmentation/KeySetsContainer.hpp"
    65 #ifdef HAVE_VMG
     74#if defined(HAVE_JOBMARKET) && defined(HAVE_VMG)
    6675#include "Fragmentation/Automation/VMGDebugGridFragmentController.hpp"
    6776#include "Fragmentation/Automation/VMGFragmentController.hpp"
    … …  
    99108}
    100109
     110#ifdef HAVE_JOBMARKET
    101111static void updateSteps(Process &p, const size_t step, const size_t total)
    102112{
    … …  
    104114  p.setCurrStep(step);
    105115}
     116#endif
    106117
    107118ActionState::ptr FragmentationFragmentationAutomationAction::performCall() {
    … …  
    117128  size_t Exitflag = 0;
    118129  std::map<JobId_t, MPQCData> shortrangedata;
     130#ifdef HAVE_JOBMARKET
    119131  {
    120132    const size_t NumberJobs = FragmentJobQueue::getInstance().size();
    … …  
    126138
    127139    // Phase Two: add MPQCJobs and send
    128     const size_t NoJobs = mpqccontroller.addJobsFromQueue(
     140    const bool AddJobsStatus = mpqccontroller.addJobsFromQueue(
    129141        params.DoLongrange.get() ? MPQCData::DoSampleDensity : MPQCData::DontSampleDensity,
    130142        params.DoValenceOnly.get() ? MPQCData::DoSampleValenceOnly : MPQCData::DontSampleValenceOnly
    131143        );
    132     LOG(1, "INFO: Added " << NoJobs << " from FragmentJobsQueue.");
     144    if (AddJobsStatus)
     145      LOG(1, "INFO: Added jobs from FragmentJobsQueue.");
     146    else {
     147      ELOG(1, "Adding jobs failed.");
     148      return Action::failure;
     149    }
    133150    mpqccontroller.run();
    134151
    … …  
    148165    Exitflag += mpqccontroller.getExitflag();
    149166  }
    150 
    151 #ifdef HAVE_VMG
     167#else
     168  {
     169    const size_t NumberJobs = FragmentJobQueue::getInstance().size();
     170    MPQCCommandFragmentController mpqccontroller;
     171    // Phase One: obtain ids: not needed, we have infinite pool
     172
     173    // Phase Two: add MPQCJobs and send
     174    const size_t NoJobs = mpqccontroller.addJobsFromQueue(
     175        params.DoLongrange.get() ? MPQCData::DoSampleDensity : MPQCData::DontSampleDensity,
     176        params.DoValenceOnly.get() ? MPQCData::DoSampleValenceOnly : MPQCData::DontSampleValenceOnly,
     177        params.executable.get().string()
     178        );
     179    LOG(1, "INFO: Added " << NoJobs << " from FragmentJobsQueue.");
     180    mpqccontroller.run();
     181
     182    // get back the results and place them in shortrangedata
     183    mpqccontroller.getResults(shortrangedata);
     184    ASSERT( shortrangedata.size() == NumberJobs,
     185        "FragmentationFragmentationAutomationAction::performCall() - number of converted results "
     186        +toString(shortrangedata.size())+" and number of jobs "+toString(NumberJobs)+ " differ.");
     187
     188    Exitflag += mpqccontroller.getExitflag();
     189  }
     190#endif
     191
     192#if defined(HAVE_JOBMARKET) && defined(HAVE_VMG)
    152193  if (params.DoLongrange.get()) {
    153194  if ( World::getInstance().getAllAtoms().size() == 0) {

  • src/Actions/GlobalListOfActions.hpp

    rc73e35 rd8821e  
    122122  (SelectionNotShapeByName) \
    123123  (FragmentationAnalyseFragmentationResults) \
     124  (FragmentationFragmentationAutomation) \
    124125  (FragmentationClearFragmentationResults) \
    125126  (FragmentationFragmentation) \
     127  (FragmentationMolecularDynamics) \
     128  (FragmentationParseFragmentJobs) \
    126129  (FragmentationStoreSaturatedFragment) \
    127130  (PotentialFitParticleCharges) \
    … …  
    137140  (ShapeTranslateShape)
    138141
    139 // we need to append the automation action in case we have the JobMarket
    140 #ifdef HAVE_JOBMARKET
    141 #define GLOBALLISTOFACTIONS_JOBMARKET \
    142     BOOST_PP_SEQ_PUSH_BACK( \
    143         BOOST_PP_SEQ_PUSH_BACK( \
    144             BOOST_PP_SEQ_PUSH_BACK( \
    145                 GLOBALLISTOFACTIONS_initial, \
    146                 FragmentationFragmentationAutomation \
    147             ), \
    148             FragmentationMolecularDynamics \
    149         ), \
    150         FragmentationParseFragmentJobs \
    151     )
    152 #else
    153 #define GLOBALLISTOFACTIONS_JOBMARKET \
    154     GLOBALLISTOFACTIONS_initial
    155 #endif /* HAVE_JOBMARKET */
    156 
     142// extend list of actions in case levmar is available
    157143#ifdef HAVE_LEVMAR
    158144#define GLOBALLISTOFACTIONS_LEVMAR \
    159145    BOOST_PP_SEQ_PUSH_BACK( \
    160         GLOBALLISTOFACTIONS_JOBMARKET, \
     146        GLOBALLISTOFACTIONS_initial, \
    161147        PotentialFitPotential \
    162148    )
    163149#else
    164150#define GLOBALLISTOFACTIONS_LEVMAR \
    165     GLOBALLISTOFACTIONS_JOBMARKET
     151    GLOBALLISTOFACTIONS_initial
    166152#endif /* HAVE_LEVMAR */
    167153

  • src/Actions/Makefile.am

    rc73e35 rd8821e  
    228228  Actions/FragmentationAction/ClearFragmentationResultsAction.cpp \
    229229  Actions/FragmentationAction/FragmentationAction.cpp \
     230  Actions/FragmentationAction/FragmentationAutomationAction.cpp \
     231  Actions/FragmentationAction/MolecularDynamicsAction.cpp \
     232  Actions/FragmentationAction/ParseFragmentJobsAction.cpp \
    230233  Actions/FragmentationAction/StoreSaturatedFragmentAction.cpp
    231234FRAGMENTATIONACTIONHEADER = \
    … …  
    233236  Actions/FragmentationAction/ClearFragmentationResultsAction.hpp \
    234237  Actions/FragmentationAction/FragmentationAction.hpp \
     238  Actions/FragmentationAction/FragmentationAutomationAction.hpp \
     239  Actions/FragmentationAction/MolecularDynamicsAction.hpp \
     240  Actions/FragmentationAction/ParseFragmentJobsAction.hpp \
    235241  Actions/FragmentationAction/StoreSaturatedFragmentAction.hpp
    236242FRAGMENTATIONACTIONDEFS = \
    … …  
    238244  Actions/FragmentationAction/ClearFragmentationResultsAction.def \
    239245  Actions/FragmentationAction/FragmentationAction.def \
    240   Actions/FragmentationAction/StoreSaturatedFragmentAction.def
    241 
    242 if CONDJOBMARKET
    243 FRAGMENTATIONACTIONSOURCE += \
    244   Actions/FragmentationAction/FragmentationAutomationAction.cpp \
    245   Actions/FragmentationAction/MolecularDynamicsAction.cpp \
    246   Actions/FragmentationAction/ParseFragmentJobsAction.cpp
    247 FRAGMENTATIONACTIONHEADER += \
    248   Actions/FragmentationAction/FragmentationAutomationAction.hpp \
    249   Actions/FragmentationAction/MolecularDynamicsAction.hpp \
    250   Actions/FragmentationAction/ParseFragmentJobsAction.hpp
    251 FRAGMENTATIONACTIONDEFS += \
    252246  Actions/FragmentationAction/FragmentationAutomationAction.def \
    253247  Actions/FragmentationAction/MolecularDynamicsAction.def \
    254   Actions/FragmentationAction/ParseFragmentJobsAction.def
    255 endif
     248  Actions/FragmentationAction/ParseFragmentJobsAction.def \
     249  Actions/FragmentationAction/StoreSaturatedFragmentAction.def
    256250
    257251GRAPHACTIONSOURCE = \

  • src/Actions/PotentialAction/FitParticleChargesAction.cpp

    rc73e35 rd8821e  
    4747#include <boost/foreach.hpp>
    4848#include <algorithm>
     49#include <functional>
    4950#include <iostream>
    5051#include <string>
    … …  
    128129    return Action::failure;
    129130  }
     131
     132  // average partial charges over all fragments
    130133  HomologyContainer::const_iterator iter = range.first;
    131134  if (!iter->second.containsGrids) {
    … …  
    133136    return Action::failure;
    134137  }
    135   const Fragment &fragment = iter->second.fragment;
    136 //  const double &energy = iter->second.energy;
    137 //  const SamplingGrid &charge = iter->second.charge_distribution;
    138   const SamplingGrid &potential = iter->second.potential_distribution;
    139 
    140   // then we extract positions from fragment
    141   PartialNucleiChargeFitter::positions_t positions;
    142   Fragment::positions_t fragmentpositions = fragment.getPositions();
    143   positions.reserve(fragmentpositions.size());
    144   BOOST_FOREACH( Fragment::position_t pos, fragmentpositions) {
    145     positions.push_back( Vector(pos[0], pos[1], pos[2]) );
    146   }
    147   PartialNucleiChargeFitter fitter(potential, positions, params.radius.get());
    148   fitter();
    149   PartialNucleiChargeFitter::charges_t return_charges = fitter.getSolutionAsCharges_t();
     138  PartialNucleiChargeFitter::charges_t averaged_charges;
     139  averaged_charges.resize(iter->second.fragment.getCharges().size(), 0.);
     140  size_t NoFragments = 0;
     141  for (;
     142      iter != range.second; ++iter, ++NoFragments) {
     143    if (!iter->second.containsGrids) {
     144      ELOG(2, "This HomologyGraph does not contain sampled grids,\ndid you forget to add '--store-grids 1' to AnalyseFragmentResults.");
     145      return Action::failure;
     146    }
     147    const Fragment &fragment = iter->second.fragment;
     148  //  const double &energy = iter->second.energy;
     149  //  const SamplingGrid &charge = iter->second.charge_distribution;
     150    const SamplingGrid &potential = iter->second.potential_distribution;
     151    if ((potential.level == 0)
     152        || ((potential.begin[0] == potential.end[0])
     153            && (potential.begin[1] == potential.end[1])
     154            && (potential.begin[2] == potential.end[2]))) {
     155      ELOG(1, "Sampled grid contains grid made of zero points.");
     156      return Action::failure;
     157    }
     158
     159    // then we extract positions from fragment
     160    PartialNucleiChargeFitter::positions_t positions;
     161    Fragment::positions_t fragmentpositions = fragment.getPositions();
     162    positions.reserve(fragmentpositions.size());
     163    BOOST_FOREACH( Fragment::position_t pos, fragmentpositions) {
     164      positions.push_back( Vector(pos[0], pos[1], pos[2]) );
     165    }
     166    PartialNucleiChargeFitter fitter(potential, positions, params.radius.get());
     167    fitter();
     168    PartialNucleiChargeFitter::charges_t return_charges = fitter.getSolutionAsCharges_t();
     169    std::transform(
     170        return_charges.begin(), return_charges.end(),
     171        averaged_charges.begin(),
     172        averaged_charges.begin(),
     173        std::plus<PartialNucleiChargeFitter::charge_t>());
     174  }
     175  std::transform(averaged_charges.begin(),averaged_charges.end(),
     176      averaged_charges.begin(),
     177      std::bind1st(std::multiplies<PartialNucleiChargeFitter::charge_t>(),1./NoFragments)
     178  );
     179
    150180
    151181  // output fitted charges
    152   LOG(0, "STATUS: We have fitted the following charges " << return_charges << ".");
     182  LOG(0, "STATUS: We have fitted the following charges " << averaged_charges
     183      << ", averaged over " << NoFragments << " fragments.");
    153184
    154185  return Action::success;

  • src/Actions/PotentialAction/FitPotentialAction.cpp

    rc73e35 rd8821e  
    113113}
    114114
     115SerializablePotential::ParticleTypes_t getNumbersFromElements(
     116    const std::vector<const element *> &fragment)
     117{
     118  SerializablePotential::ParticleTypes_t fragmentnumbers;
     119  std::transform(fragment.begin(), fragment.end(), std::back_inserter(fragmentnumbers),
     120      boost::bind(&element::getAtomicNumber, _1));
     121  return fragmentnumbers;
     122}
     123
     124
    115125ActionState::ptr PotentialFitPotentialAction::performCall() {
    116126  // fragment specifies the homology fragment to use
    117   SerializablePotential::ParticleTypes_t fragmentnumbers;
    118   {
    119     const std::vector<const element *> &fragment = params.fragment.get();
    120     std::transform(fragment.begin(), fragment.end(), std::back_inserter(fragmentnumbers),
    121         boost::bind(&element::getAtomicNumber, _1));
    122   }
     127  SerializablePotential::ParticleTypes_t fragmentnumbers =
     128      getNumbersFromElements(params.fragment.get());
    123129
    124130  // either charges and a potential is specified or a file
    … …  
    158164    } else {
    159165      // charges specify the potential type
    160       SerializablePotential::ParticleTypes_t chargenumbers;
    161       {
    162         const std::vector<const element *> &charges = params.charges.get();
    163         std::transform(charges.begin(), charges.end(), std::back_inserter(chargenumbers),
    164             boost::bind(&element::getAtomicNumber, _1));
    165       }
     166      SerializablePotential::ParticleTypes_t chargenumbers =
     167          getNumbersFromElements(params.charges.get());
    166168
    167169      LOG(0, "STATUS: I'm training now a " << params.potentialtype.get()
    … …  
    240242      size_t counter=1;
    241243      if (DoLog(3)) {
    242         const FunctionModel::arguments_t &inputs = data.getTrainingInputs()[0];
     244        const FunctionModel::arguments_t &inputs = data.getAllArguments()[0];
    243245        for (FunctionModel::arguments_t::const_iterator iter = inputs.begin();
    244246            iter != inputs.end(); ++iter) {
    … …  
    265267    }
    266268
     269    if ((params.threshold.get() < 1) && (params.best_of_howmany.isSet()))
     270      ELOG(2, "threshold parameter always overrules max_runs, both are specified.");
    267271    // now perform the function approximation by optimizing the model function
    268272    FunctionApproximation approximator(data, *model);
    … …  
    304308
    305309    // create a map of each fragment with error.
    306     typedef std::multimap< double, size_t > WorseFragmentMap_t;
    307     WorseFragmentMap_t WorseFragmentMap;
    308310    HomologyContainer::range_t fragmentrange = homologies.getHomologousGraphs(graph);
    309     // fragments make it into the container in reversed order, hence count from top down
    310     size_t index= std::distance(fragmentrange.first, fragmentrange.second)-1;
    311     for (HomologyContainer::const_iterator iter = fragmentrange.first;
    312         iter != fragmentrange.second;
    313         ++iter) {
    314       const Fragment& fragment = iter->second.fragment;
    315       const double &energy = iter->second.energy;
    316 
    317       // create arguments from the fragment
    318       FunctionModel::extractor_t extractor = model->getSpecificExtractor();
    319       FunctionModel::arguments_t args = extractor(fragment, 1);
    320 
    321       // calculate value from potential
    322       const double fitvalue = (*model)(args)[0];
    323 
    324       // insert difference into map
    325       const double error = fabs(energy - fitvalue);
    326       WorseFragmentMap.insert( std::make_pair( error, index-- ) );
    327 
    328       {
    329         // give only the distances in the debugging text
    330         std::stringstream streamargs;
    331         BOOST_FOREACH (argument_t arg, args) {
    332           streamargs << " " << arg.distance*AtomicLengthToAngstroem;
    333         }
    334         LOG(2, "DEBUG: frag.#" << index+1 << "'s error is |" << energy << " - " << fitvalue
    335             << "| = " << error << " for args " << streamargs.str() << ".");
    336       }
    337     }
     311    TrainingData::L2ErrorConfigurationIndexMap_t WorseFragmentMap =
     312        data.getWorstFragmentMap(*model, fragmentrange);
    338313    LOG(0, "RESULT: WorstFragmentMap " << WorseFragmentMap << ".");
    339314

  • src/Element/Makefile.am

    rc73e35 rd8821e  
    44ELEMENTSOURCE = \
    55        Element/element.cpp \
     6        Element/ion.cpp \
    67        Element/periodentafel.cpp \
    78        Element/elements_db.cpp
    … …  
    910ELEMENTHEADER = \
    1011        Element/element.hpp \
     12        Element/ion.hpp \
    1113        Element/periodentafel.hpp \
    1214        Element/elements_db.hpp

  • src/Element/element.hpp

    rc73e35 rd8821e  
    2525
    2626class periodentafel;
     27class IonTest;
    2728
    2829/********************************************** declarations *******************************/
    … …  
    3334class element {
    3435  friend class periodentafel;
     36  friend class IonTest;
    3537  public:
    3638    element();
    3739    element(const element&);
    38     ~element();
     40    virtual ~element();
    3941
    4042    element &operator=(const element&);
    … …  
    4749    double getVanDerWaalsRadius() const;
    4850    atomicNumber_t getAtomicNumber() const;
    49     double getValence() const;
    50     int getNoValenceOrbitals() const;
     51    virtual double getCharge() const
     52        { return 0.; }
     53    virtual double getValence() const;
     54    virtual int getNoValenceOrbitals() const;
    5155    double getHBondDistance(const size_t i) const;
    5256    double getHBondAngle(const size_t i) const;
    … …  
    8791    }
    8892
     93  protected:
     94
    8995    double mass;    //!< mass in g/mol
    9096    double CovalentRadius;  //!< covalent radius

  • src/Element/periodentafel.cpp

    rc73e35 rd8821e  
    4545#include "elements_db.hpp"
    4646#include "Helpers/defs.hpp"
     47#include "ion.hpp"
    4748#include "periodentafel.hpp"
    4849
    … …  
    120121  return res!=elements.end()?((*res).second):0;
    121122};
     123
     124/** Returns the desired ion to a specific element.
     125 *
     126 * If the respective element is not in the list, we return
     127 * NULL.
     128 * \return pointer to an element or NULL if not found
     129 */
     130const element * periodentafel::FindElement(atomicNumber_t Z, const int ionization)
     131{
     132  // if not ionization given, fall back to other function
     133  if (ionization == 0) {
     134    return FindElement(Z);
     135  }
     136  // element present?
     137  const_iterator elementiter = elements.find(Z);
     138  if (elementiter == elements.end())
     139    return NULL;
     140  const element & element_base = *(elementiter->second);
     141       
     142  // element has already got ions?
     143  IonsPerElement::iterator setiter =
     144      ions.find(Z);
     145  if (setiter != ions.end()) {
     146    // yes, found ion list
     147    ionSet::const_iterator res = setiter->second.find(ionization);
     148       
     149        if (res != setiter->second.end()) {
     150          // ion present already
     151          element * const _ion = res->second;
     152          return _ion;
     153        } else {
     154          // ion not present yet
     155          ion * const _ion = new ion(element_base, ionization);
     156          // insert ion
     157          setiter->second.insert( std::make_pair( ionization, _ion) );
     158          return _ion;
     159        }
     160  } else {
     161    // no ions yet, create map
     162    ion * const _ion = new ion(element_base, ionization);
     163    std::pair<IonsPerElement::iterator,bool> inserter =
     164        ions.insert( std::make_pair(Z, ionSet()) );
     165    // insert ion
     166    ASSERT( inserter.second,
     167        "periodentafel::FindElement() - could not insert new ionSet to element.");
     168    inserter.first->second.insert( std::make_pair( ionization, _ion) );
     169    return _ion;
     170  }
     171  return NULL;
     172}
    122173
    123174/** Finds an element by its atomic number.

  • src/Element/periodentafel.hpp

    rc73e35 rd8821e  
    2222
    2323class element;
     24class ion;
    2425
    2526/********************************************** declarations *******************************/
    … …  
    3334private:
    3435  typedef std::map<atomicNumber_t,element*> elementSet;
     36  typedef std::map<int,ion*> ionSet;
     37  typedef std::map<atomicNumber_t,ionSet> IonsPerElement;
    3538public:
    3639  typedef elementSet::iterator iterator;
    … …  
    4851  void CleanupPeriodtable();
    4952  const element * FindElement(atomicNumber_t) const;
     53  const element * FindElement(atomicNumber_t, const int);
    5054  const element * FindElement(const std::string &shorthand) const;
    5155  const element * AskElement() const;
    … …  
    7680    ar & boost::serialization::make_array<char>(header2, MAXSTRINGSIZE);
    7781    ar & elements;
     82    ar & ions;
    7883  }
    7984
    … …  
    9398
    9499  elementSet elements;
     100  IonsPerElement ions;
    95101};
    96102

  • src/Element/unittests/ElementUnitTest.cpp

    rc73e35 rd8821e  
    105105  CPPUNIT_ASSERT_EQUAL( 0., testelement->getElectronegativity() );
    106106  CPPUNIT_ASSERT_EQUAL( 0., testelement->getVanDerWaalsRadius() );
     107  CPPUNIT_ASSERT_EQUAL( 0., testelement->getCharge() );
    107108  CPPUNIT_ASSERT_EQUAL( 0., testelement->getValence() );
    108109  CPPUNIT_ASSERT_EQUAL( 0, testelement->getNoValenceOrbitals() );
    … …  
    172173  boost::archive::text_iarchive ia(stream);
    173174  // read class state from archive
    174   element * newelement;
     175  element * newelement = NULL;
    175176
    176177  ia >> newelement;

  • src/Element/unittests/Makefile.am

    rc73e35 rd8821e  
    44ELEMENTTESTSSOURCES = \
    55  ../Element/unittests/ElementUnitTest.cpp \
    6         ../Element/unittests/PeriodentafelUnitTest.cpp
     6  ../Element/unittests/IonUnitTest.cpp \
     7  ../Element/unittests/PeriodentafelUnitTest.cpp
    78
    89ELEMENTTESTSHEADERS = \
    9         ../Element/unittests/ElementUnitTest.hpp \
    10         ../Element/unittests/PeriodentafelUnitTest.hpp
     10  ../Element/unittests/ElementUnitTest.hpp \
     11  ../Element/unittests/IonUnitTest.hpp \
     12  ../Element/unittests/PeriodentafelUnitTest.hpp
    1113
    1214ELEMENTTESTS = \
    1315  ElementUnitTest \
     16  IonUnitTest \
    1417  PeriodentafelUnitTest
    1518
    … …  
    2427        ${CodePatterns_LIBS}
    2528
    26 ElementUnitTest_SOURCES = $(top_srcdir)/src/unittests/UnitTestMain.cpp \
    27         ../Element/unittests/ElementUnitTest.cpp \
    28         ../Element/unittests/ElementUnitTest.hpp
     29ElementUnitTest_SOURCES =   $(top_srcdir)/src/unittests/UnitTestMain.cpp \
     30  ../Element/unittests/ElementUnitTest.cpp \
     31  ../Element/unittests/ElementUnitTest.hpp
    2932ElementUnitTest_LDADD = $(ELEMENTLIBS)
     33
     34IonUnitTest_SOURCES =   $(top_srcdir)/src/unittests/UnitTestMain.cpp \
     35  ../Element/unittests/IonUnitTest.cpp \
     36  ../Element/unittests/IonUnitTest.hpp
     37IonUnitTest_LDADD = $(ELEMENTLIBS)
    3038
    3139PeriodentafelUnitTest_SOURCES = $(top_srcdir)/src/unittests/UnitTestMain.cpp \

  • src/Element/unittests/PeriodentafelUnitTest.cpp

    rc73e35 rd8821e  
    4545
    4646#include "Element/element.hpp"
     47#include "Element/ion.hpp"
    4748#include "Element/elements_db.hpp"
    4849#include "Element/periodentafel.hpp"

  • src/Fragmentation/Automation/MPQCFragmentController.cpp

    rc73e35 rd8821e  
    5757    )
    5858{
     59  bool status = true;
    5960  // give them all valid ids
    6061  std::vector<FragmentJob::ptr> newjobs = FragmentJobQueue::getInstance().getJobs();
    … …  
    6566    job->DoLongrange = _DoLongrange;
    6667    job->DoValenceOnly = _DoValenceOnly;
    67     changeJobId((*jobiter), getAvailableId());
     68    JobId_t id = getAvailableId();
     69    if (id == FragmentJob::IllegalJob) {
     70      status = false;
     71      break;
     72    }
     73    changeJobId((*jobiter), id);
    6874  }
    6975  // add the jobs
    7076  addJobs(newjobs);
     77  status &= (newjobs.size() != 0);
    7178
    72   return (newjobs.size() != 0);
     79  return status;
    7380}
    7481

  • src/Fragmentation/Automation/MPQCFragmentController.hpp

    rc73e35 rd8821e  
    7070private:
    7171  //!> type-specific result container
    72   SpecificFragmentController::ResultContainer<MPQCData> results;
     72  SpecificFragmentController::ReceiveResultContainer<MPQCData> results;
    7373};
    7474

  • src/Fragmentation/Automation/Makefile.am

    rc73e35 rd8821e  
    44FRAGMENTATIONAUTOMATIONSOURCE = \
    55        Fragmentation/Automation/FragmentJobQueue.cpp \
     6        Fragmentation/Automation/MPQCCommandFragmentController.cpp
     7
     8if CONDJOBMARKET
     9FRAGMENTATIONAUTOMATIONSOURCE += \
    610        Fragmentation/Automation/MPQCFragmentController.cpp \
    711        Fragmentation/Automation/SpecificFragmentController.cpp
    … …  
    1317endif
    1418
     19endif
     20
    1521FRAGMENTATIONAUTOMATIONHEADER = \
    1622        Fragmentation/Automation/FragmentJobQueue.hpp \
     23        Fragmentation/Automation/MPQCCommandFragmentController.hpp \
     24        Fragmentation/Automation/ResultContainer.hpp \
     25        Fragmentation/Automation/ResultContainer_impl.hpp
     26
     27if CONDJOBMARKET
     28FRAGMENTATIONAUTOMATIONHEADER += \
    1729        Fragmentation/Automation/MPQCFragmentController.hpp \
    1830        Fragmentation/Automation/SpecificFragmentController.hpp \
    19         Fragmentation/Automation/SpecificFragmentController_ResultContainer_impl.hpp
     31        Fragmentation/Automation/SpecificFragmentController_ReceiveResultContainer_impl.hpp
    2032       
    2133if CONDVMG
    … …  
    2335        Fragmentation/Automation/VMGDebugGridFragmentController.hpp \
    2436        Fragmentation/Automation/VMGFragmentController.hpp
     37endif
    2538endif
    2639

  • src/Fragmentation/Automation/SpecificFragmentController.hpp

    rc73e35 rd8821e  
    2323#include <string>
    2424#include <vector>
     25
     26#include "Fragmentation/Automation/ResultContainer.hpp"
    2527
    2628/** This class extends FragmentController by some commodity functions used
    … …  
    5860   * This struct takes of the waiting.
    5961   *
     62   * We extend ResultContainer by functions to receive results from a
     63   * JobMarket server.
     64   *
    6065   */
    6166  template <typename T>
    62   struct ResultContainer {
     67  struct ReceiveResultContainer : public ResultContainer<T> {
    6368    /** cycle to wait for results
    6469     *
    … …  
    7883     */
    7984    size_t receiveResults(SpecificFragmentController &callback);
    80 
    81     /** Extracts specific Data type from received vector of FragmentResults.
    82      *
    83      * @param results results to extract specific Data type from
    84      * @param fragmentData on return array filled with extracted specific Data type
    85      */
    86     void ConvertFragmentResultTo(
    87         const std::vector<FragmentResult::ptr> &results,
    88         std::vector<T> &fragmentData);
    89 
    90     /** Clears all internally stored results.
    91      *
    92      */
    93     void clear()
    94     { IdData.clear(); }
    95 
    96     //!> internal container for the received data
    97     std::map<JobId_t, T> IdData;
    9885  };
    9986
    … …  
    10996};
    11097
    111 #include "SpecificFragmentController_ResultContainer_impl.hpp"
     98#include "SpecificFragmentController_ReceiveResultContainer_impl.hpp"
    11299
    113100#endif /* SPECIFICFRAGMENTCONTROLLER_HPP_ */

  • src/Fragmentation/Automation/VMGDebugGridFragmentController.hpp

    rc73e35 rd8821e  
    5555private:
    5656  //!> type-specific result container
    57   SpecificFragmentController::ResultContainer<std::string> results;
     57  SpecificFragmentController::ReceiveResultContainer<std::string> results;
    5858};
    5959

  • src/Fragmentation/Automation/VMGFragmentController.hpp

    rc73e35 rd8821e  
    8989private:
    9090  //!> type-specific result container
    91   SpecificFragmentController::ResultContainer<VMGData> results;
     91  SpecificFragmentController::ReceiveResultContainer<VMGData> results;
    9292};
    9393

  • src/Fragmentation/Exporters/ExportGraph_ToJobs.cpp

    rc73e35 rd8821e  
    4646#include "Fragmentation/KeySet.hpp"
    4747#include "Fragmentation/Automation/FragmentJobQueue.hpp"
    48 #include "Fragmentation/Automation/MPQCFragmentController.hpp"
    4948#include "Helpers/defs.hpp"
     49#ifdef HAVE_JOBMARKET
    5050#include "Jobs/MPQCJob.hpp"
     51#else
     52#include "Jobs/MPQCCommandJob.hpp"
     53#endif
    5154#include "LinearAlgebra/RealSpaceMatrix.hpp"
    5255#include "Parser/FormatParserStorage.hpp"
    … …  
    9396      CurrentFragment->OutputConfig(output, jobtype);
    9497      // create job and insert
    95       FragmentJob::ptr testJob( new MPQCJob(JobId::IllegalJob, output.str(), begin, end, level) );
     98      FragmentJob::ptr testJob(
     99#ifdef HAVE_JOBMARKET
     100          new MPQCJob(JobId::IllegalJob, output.str(), begin, end, level)
     101#else
     102          new MPQCCommandJob(output.str(), JobId::IllegalJob)
     103#endif
     104      );
    96105      testJob->setPriority(CurrentFragment->getKeySet().size());
    97106      jobs.push_back(testJob);

  • src/Fragmentation/Makefile.am

    rc73e35 rd8821e  
    44FRAGMENTATIONSOURCE = \
    55        Fragmentation/Exporters/ExportGraph_ToFiles.cpp \
     6        Fragmentation/Exporters/ExportGraph_ToJobs.cpp \
    67        Fragmentation/Exporters/ExportGraph.cpp \
    78        Fragmentation/Exporters/HydrogenPool.cpp \
    … …  
    2930FRAGMENTATIONHEADER = \
    3031        Fragmentation/Exporters/ExportGraph_ToFiles.hpp \
     32        Fragmentation/Exporters/ExportGraph_ToJobs.hpp \
    3133        Fragmentation/Exporters/ExportGraph.hpp \
    3234        Fragmentation/Exporters/HydrogenPool.hpp \
    … …  
    5759        Fragmentation/SortIndex.hpp \
    5860        Fragmentation/UniqueFragments.hpp
    59 
    60 if CONDJOBMARKET
    61 FRAGMENTATIONSOURCE += \
    62         Fragmentation/Exporters/ExportGraph_ToJobs.cpp
    63 FRAGMENTATIONHEADER += \
    64         Fragmentation/Exporters/ExportGraph_ToJobs.hpp
    65 endif
    6661
    6762lib_LTLIBRARIES += \

  • src/Fragmentation/Summation/Containers/FragmentationResultContainer.cpp

    rc73e35 rd8821e  
    4141#include "CodePatterns/Singleton_impl.hpp"
    4242
    43 #ifdef HAVE_VMG
     43#if defined(HAVE_JOBMARKET) && defined(HAVE_VMG)
    4444void FragmentationResultContainer::addFullResults(
    4545    const KeySetsContainer &_keysets,
    … …  
    5151  OBSERVE;
    5252  ASSERT( ResultsType == BothRanges,
    53       "FragmentationResultContainer::getLongRangeResults() - mixing full and short range is not in the cards.");
     53      "FragmentationResultContainer::addFullResults() - mixing full and short range is not in the cards.");
    5454  ASSERT( _keysets.KeySets.size() == _forcekeysets.KeySets.size(),
    5555      "FragmentationResultContainer::addFullResults() - keysets ("
    … …  
    8080  OBSERVE;
    8181  ASSERT( (ResultsType == ShortRangeOnly) || (keysets.KeySets.empty()) ,
    82       "FragmentationResultContainer::getLongRangeResults() - mixing full and short range is not in the cards.");
     82      "FragmentationResultContainer::addFullResults() - mixing full and short range is not in the cards.");
    8383  ASSERT( _keysets.KeySets.size() == _forcekeysets.KeySets.size(),
    8484      "FragmentationResultContainer::addFullResults() - keysets ("
    … …  
    104104}
    105105
    106 #ifdef HAVE_VMG
     106#if defined(HAVE_JOBMARKET) && defined(HAVE_VMG)
    107107const FragmentationResultContainer::longrangedata_t&
    108108FragmentationResultContainer::getLongRangeResults() const

  • src/Fragmentation/Summation/Containers/FragmentationResultContainer.hpp

    rc73e35 rd8821e  
    3636#include "Fragmentation/Summation/Containers/MPQCData.hpp"
    3737#include "Fragmentation/Summation/Containers/MPQCDataMap.hpp"
    38 #ifdef HAVE_VMG
     38#if defined(HAVE_JOBMARKET) && defined(HAVE_VMG)
    3939#include "Fragmentation/Summation/Containers/VMGData.hpp"
    4040#endif
    … …  
    6363  //!> typedef for short range data container
    6464  typedef std::map<JobId_t, MPQCData> shortrangedata_t;
    65 #ifdef HAVE_VMG
     65#if defined(HAVE_JOBMARKET) && defined(HAVE_VMG)
    6666  //!> typedef for long range data container
    6767  typedef std::map<JobId_t, VMGData> longrangedata_t;
    6868#endif
    6969
    70 #ifdef HAVE_VMG
     70#if defined(HAVE_JOBMARKET) && defined(HAVE_VMG)
    7171  /** Adds a a set of both short and long range results.
    7272   *
    … …  
    128128    shortrangedata.clear();
    129129    summedshortrange.clear();
    130 #ifdef HAVE_VMG
     130#if defined(HAVE_JOBMARKET) && defined(HAVE_VMG)
    131131    longrangedata.clear();
    132132#endif
    … …  
    139139  const shortrangedata_t& getShortRangeResults() const  { return shortrangedata; }
    140140  const FragmentationShortRangeResults::summedshortrange_t& getShortRangeSummedResults() const  { return summedshortrange; }
    141 #ifdef HAVE_VMG
     141#if defined(HAVE_JOBMARKET) && defined(HAVE_VMG)
    142142  const longrangedata_t& getLongRangeResults() const;
    143143#endif
    … …  
    162162  FragmentationShortRangeResults::summedshortrange_t summedshortrange;
    163163
    164 #ifdef HAVE_VMG
     164#if defined(HAVE_JOBMARKET) && defined(HAVE_VMG)
    165165  //!> container for all long-range results
    166166  std::map<JobId_t, VMGData> longrangedata;
    … …  
    181181    if (version > 1)
    182182      ar & summedshortrange;
    183 #ifdef HAVE_VMG
     183#if defined(HAVE_JOBMARKET) && defined(HAVE_VMG)
    184184    ar & longrangedata;
    185185#endif
    … …  
    196196    if (version > 1)
    197197      ar & summedshortrange;
    198 #ifdef HAVE_VMG
     198#if defined(HAVE_JOBMARKET) && defined(HAVE_VMG)
    199199    ar & longrangedata;
    200200#endif

  • src/Fragmentation/Summation/Containers/Makefile.am

    rc73e35 rd8821e  
    55        Fragmentation/Summation/Containers/createMatrixNrLookup.cpp \
    66        Fragmentation/Summation/Containers/FragmentationChargeDensity.cpp \
    7         Fragmentation/Summation/Containers/FragmentationLongRangeResults.cpp \
    87        Fragmentation/Summation/Containers/FragmentationResultContainer.cpp \
    98        Fragmentation/Summation/Containers/FragmentationShortRangeResults.cpp \
    109        Fragmentation/Summation/Containers/MPQCData.cpp
    11        
     10
     11if CONDJOBMARKET       
    1212if CONDVMG
    1313FRAGMENTATIONCONTAINERSOURCE += \
     14        Fragmentation/Summation/Containers/FragmentationLongRangeResults.cpp \
    1415        Fragmentation/Summation/Containers/VMGData.cpp
     16endif
    1517endif
    1618
    … …  
    1921        Fragmentation/Summation/Containers/extractJobIds.hpp \
    2022        Fragmentation/Summation/Containers/FragmentationChargeDensity.hpp \
    21         Fragmentation/Summation/Containers/FragmentationLongRangeResults.hpp \
    2223        Fragmentation/Summation/Containers/FragmentationResultContainer.hpp \
    2324        Fragmentation/Summation/Containers/FragmentationShortRangeResults.hpp \
    … …  
    2627        Fragmentation/Summation/Containers/MPQCDataMap.hpp \
    2728        Fragmentation/Summation/Containers/MPQCData_printKeyNames.hpp
    28        
     29
     30if CONDJOBMARKET       
    2931if CONDVMG
    3032FRAGMENTATIONCONTAINERHEADER += \
     33        Fragmentation/Summation/Containers/FragmentationLongRangeResults.hpp \
    3134        Fragmentation/Summation/Containers/VMGData.hpp \
    3235        Fragmentation/Summation/Containers/VMGDataFused.hpp \
    3336        Fragmentation/Summation/Containers/VMGDataMap.hpp \
    3437        Fragmentation/Summation/Containers/VMGData_printKeyNames.hpp
     38endif
    3539endif
    3640

  • src/FunctionApproximation/Extractors.cpp

    rc73e35 rd8821e  
    532532}
    533533
    534 FunctionModel::arguments_t Extractors::reorderArgumentsByIncreasingDistance(
    535     const FunctionModel::arguments_t &args
    536     )
    537 {
    538   FunctionModel::arguments_t returnargs(args);
    539   std::sort(returnargs.begin(), returnargs.end(), argument_t::DistanceComparator);
     534FunctionModel::list_of_arguments_t Extractors::reorderArgumentsByIncreasingDistance(
     535    const FunctionModel::list_of_arguments_t &listargs
     536    )
     537{
     538  FunctionModel::list_of_arguments_t returnargs;
     539  for (FunctionModel::list_of_arguments_t::const_iterator iter = listargs.begin();
     540      iter != listargs.end(); ++iter) {
     541    const FunctionModel::arguments_t &args = *iter;
     542    FunctionModel::arguments_t sortedargs(args);
     543    std::sort(sortedargs.begin(), sortedargs.end(), argument_t::DistanceComparator);
     544    returnargs.push_back(sortedargs);
     545  }
    540546  return returnargs;
    541547}
    … …  
    587593}
    588594
    589 FunctionModel::arguments_t Extractors::reorderArgumentsByParticleTypes(
     595FunctionModel::list_of_arguments_t Extractors::reorderArgumentsByParticleTypes(
     596    const FunctionModel::list_of_arguments_t &listargs,
     597    const ParticleTypes_t &_types
     598    )
     599{
     600  FunctionModel::list_of_arguments_t returnargs;
     601  for (FunctionModel::list_of_arguments_t::const_iterator iter = listargs.begin();
     602      iter != listargs.end(); ++iter) {
     603    const FunctionModel::arguments_t &args = *iter;
     604    /// We  place all arguments into multimap according to particle type pair.
     605    // here, we need a special comparator such that types in key pair are always
     606    // properly ordered.
     607    typedef std::multimap<
     608        argument_t::types_t,
     609        argument_t,
     610        ParticleTypesComparator> TypePair_Argument_Map_t;
     611    TypePair_Argument_Map_t argument_map;
     612    for(FunctionModel::arguments_t::const_iterator iter = args.begin();
     613        iter != args.end(); ++iter) {
     614      argument_map.insert( std::make_pair(iter->types, *iter) );
     615    }
     616    LOG(4, "DEBUG: particle_type map is " << argument_map << ".");
     617
     618    /// Then, we create the desired unique keys
     619    typedef std::vector<argument_t::types_t> UniqueTypes_t;
     620    UniqueTypes_t UniqueTypes;
     621    for (ParticleTypes_t::const_iterator firstiter = _types.begin();
     622        firstiter != _types.end();
     623        ++firstiter) {
     624      for (ParticleTypes_t::const_iterator seconditer = firstiter;
     625          seconditer != _types.end();
     626          ++seconditer) {
     627        if (seconditer == firstiter)
     628          continue;
     629        UniqueTypes.push_back( std::make_pair(*firstiter, *seconditer) );
     630      }
     631    }
     632    LOG(4, "DEBUG: Created unique types as keys " << UniqueTypes << ".");
     633
     634    /// Finally, we use the unique key list to pick corresponding arguments from the map
     635    FunctionModel::arguments_t sortedargs;
     636    sortedargs.reserve(args.size());
     637    while (!argument_map.empty()) {
     638      // note that particle_types_t may be flipped, i.e. 1,8 is equal to 8,1, but we
     639      // must maintain the correct order in indices in accordance with the order
     640      // in _types, i.e. 1,8,1 must match with e.g. ids 1,0,2 where 1 has type 1,
     641      // 0 has type 8, and 2 has type 2.
     642      // In other words: We do not want to flip/modify arguments such that they match
     643      // with the specific type pair we seek but then this comes at the price that we
     644      // have flip indices when the types in a pair are flipped.
     645
     646      typedef std::vector<size_t> indices_t;
     647      //!> here, we gather the indices as we discover them
     648      indices_t indices;
     649      indices.resize(_types.size(), (size_t)-1);
     650
     651      // these are two iterators that create index pairs in the same way as we have
     652      // created type pairs. If a -1 is still present in indices, then the index is
     653      // still arbitrary but is then set by the next found index
     654      indices_t::iterator firstindex = indices.begin();
     655      indices_t::iterator secondindex = firstindex+1;
     656
     657      //!> here, we gather the current bunch of arguments as we find them
     658      FunctionModel::arguments_t argumentbunch;
     659      argumentbunch.reserve(UniqueTypes.size());
     660
     661      for (UniqueTypes_t::const_iterator typeiter = UniqueTypes.begin();
     662          typeiter != UniqueTypes.end(); ++typeiter) {
     663        // have all arguments to same type pair as list within the found range
     664        std::pair<
     665            TypePair_Argument_Map_t::iterator,
     666            TypePair_Argument_Map_t::iterator> range_t =
     667                argument_map.equal_range(*typeiter);
     668        LOG(4, "DEBUG: Set of arguments to current key [" << typeiter->first << ","
     669            << typeiter->second << "] is " << std::list<argument_t>(
     670                MapValueIterator<TypePair_Argument_Map_t::iterator>(range_t.first),
     671                MapValueIterator<TypePair_Argument_Map_t::iterator>(range_t.second)
     672                ) << ".");
     673        // the first key is always easy and is pivot which the rest has to be associated to
     674        if (typeiter == UniqueTypes.begin()) {
     675          const argument_t & arg = range_t.first->second;
     676          if ((typeiter->first == arg.types.first) && (typeiter->second == arg.types.second)) {
     677            // store in correct order
     678            *firstindex = arg.indices.first;
     679            *secondindex = arg.indices.second;
     680          } else {
     681            // store in flipped order
     682            *firstindex = arg.indices.second;
     683            *secondindex = arg.indices.first;
     684          }
     685          argumentbunch.push_back(arg);
     686          argument_map.erase(range_t.first);
     687          LOG(4, "DEBUG: Gathered first argument " << arg << ".");
     688        } else {
     689          // go through the range and pick the first argument matching the index constraints
     690          for (TypePair_Argument_Map_t::iterator argiter = range_t.first;
     691              argiter != range_t.second; ++argiter) {
     692            // seconditer may be -1 still
     693            const argument_t &arg = argiter->second;
     694            if (arg.indices.first == *firstindex) {
     695              if ((arg.indices.second == *secondindex) || (*secondindex == (size_t)-1)) {
     696                if (*secondindex == (size_t)-1)
     697                  *secondindex = arg.indices.second;
     698                argumentbunch.push_back(arg);
     699                argument_map.erase(argiter);
     700                LOG(4, "DEBUG: Gathered another argument " << arg << ".");
     701                break;
     702              }
     703            } else if ((arg.indices.first == *secondindex) || (*secondindex == (size_t)-1)) {
     704              if (arg.indices.second == *firstindex) {
     705                if (*secondindex == (size_t)-1)
     706                  *secondindex = arg.indices.first;
     707                argumentbunch.push_back(arg);
     708                argument_map.erase(argiter);
     709                LOG(4, "DEBUG: Gathered another (flipped) argument " << arg << ".");
     710                break;
     711              }
     712            }
     713          }
     714        }
     715        // move along in indices and check bounds
     716        ++secondindex;
     717        if (secondindex == indices.end()) {
     718          ++firstindex;
     719          if (firstindex != indices.end()-1)
     720            secondindex = firstindex+1;
     721        }
     722      }
     723      ASSERT( (firstindex == indices.end()-1) && (secondindex == indices.end()),
     724          "Extractors::reorderArgumentsByParticleTypes() - we have not gathered enough arguments.");
     725      ASSERT( argumentbunch.size() == UniqueTypes.size(),
     726          "Extractors::reorderArgumentsByParticleTypes() - we have not gathered enough arguments.");
     727      // place bunch of arguments in return args
     728      LOG(3, "DEBUG: Given types " << _types << " and found indices " << indices << ".");
     729      LOG(3, "DEBUG: Final bunch of arguments is " << argumentbunch << ".");
     730      sortedargs.insert(sortedargs.end(), argumentbunch.begin(), argumentbunch.end());
     731    }
     732    returnargs.push_back(sortedargs);
     733  }
     734
     735  return returnargs;
     736}
     737
     738FunctionModel::list_of_arguments_t Extractors::filterArgumentsByParticleTypes(
    590739    const FunctionModel::arguments_t &args,
    591740    const ParticleTypes_t &_types
    592741    )
    593742{
    594   /// We  place all arguments into multimap according to particle type pair.
    595   // here, we need a special comparator such that types in key pair are always
    596   // properly ordered.
    597   typedef std::multimap<
    598       argument_t::types_t,
    599       argument_t,
    600       ParticleTypesComparator> TypePair_Argument_Map_t;
    601   TypePair_Argument_Map_t argument_map;
    602   for(FunctionModel::arguments_t::const_iterator iter = args.begin();
    603       iter != args.end(); ++iter) {
    604     argument_map.insert( std::make_pair(iter->types, *iter) );
    605   }
    606   LOG(4, "DEBUG: particle_type map is " << argument_map << ".");
    607 
    608   /// Then, we create the desired unique keys
    609   typedef std::vector<argument_t::types_t> UniqueTypes_t;
    610   UniqueTypes_t UniqueTypes;
     743  typedef std::list< argument_t > ListArguments_t;
     744  ListArguments_t availableList(args.begin(), args.end());
     745  LOG(2, "DEBUG: Initial list of args is " << args << ".");
     746
     747
     748  // TODO: fill a lookup map such that we don't have O(M^3) scaling, if M is number
     749  // of types (and we always must have M(M-1)/2 args) but O(M^2 log(M)). However, as
     750  // M is very small (<=3), this is not necessary fruitful now.
     751//  typedef ParticleTypes_t firsttype;
     752//  typedef ParticleTypes_t secondtype;
     753//  typedef std::map< firsttype, std::map< secondtype, boost::ref(args) > > ArgsLookup_t;
     754//  ArgsLookup_t ArgsLookup;
     755
     756  // basically, we have two choose any two pairs out of types but only those
     757  // where the first is less than the latter. Hence, we start the second
     758  // iterator at the current position of the first one and skip the equal case.
     759  FunctionModel::arguments_t allargs;
     760  allargs.reserve(args.size());
    611761  for (ParticleTypes_t::const_iterator firstiter = _types.begin();
    612762      firstiter != _types.end();
    … …  
    617767      if (seconditer == firstiter)
    618768        continue;
    619       UniqueTypes.push_back( std::make_pair(*firstiter, *seconditer) );
    620     }
    621   }
    622   LOG(4, "DEBUG: Created unique types as keys " << UniqueTypes << ".");
    623 
    624   /// Finally, we use the unique key list to pick corresponding arguments from the map
    625   FunctionModel::arguments_t returnargs;
    626   returnargs.reserve(args.size());
    627   while (!argument_map.empty()) {
    628     // note that particle_types_t may be flipped, i.e. 1,8 is equal to 8,1, but we
    629     // must maintain the correct order in indices in accordance with the order
    630     // in _types, i.e. 1,8,1 must match with e.g. ids 1,0,2 where 1 has type 1,
    631     // 0 has type 8, and 2 has type 2.
    632     // In other words: We do not want to flip/modify arguments such that they match
    633     // with the specific type pair we seek but then this comes at the price that we
    634     // have flip indices when the types in a pair are flipped.
    635 
    636     typedef std::vector<size_t> indices_t;
    637     //!> here, we gather the indices as we discover them
    638     indices_t indices;
    639     indices.resize(_types.size(), (size_t)-1);
    640 
    641     // these are two iterators that create index pairs in the same way as we have
    642     // created type pairs. If a -1 is still present in indices, then the index is
    643     // still arbitrary but is then set by the next found index
    644     indices_t::iterator firstindex = indices.begin();
    645     indices_t::iterator secondindex = firstindex+1;
    646 
    647     //!> here, we gather the current bunch of arguments as we find them
    648     FunctionModel::arguments_t argumentbunch;
    649     argumentbunch.reserve(UniqueTypes.size());
    650 
    651     for (UniqueTypes_t::const_iterator typeiter = UniqueTypes.begin();
    652         typeiter != UniqueTypes.end(); ++typeiter) {
    653       // have all arguments to same type pair as list within the found range
    654       std::pair<
    655           TypePair_Argument_Map_t::iterator,
    656           TypePair_Argument_Map_t::iterator> range_t =
    657               argument_map.equal_range(*typeiter);
    658       LOG(4, "DEBUG: Set of arguments to current key [" << typeiter->first << ","
    659           << typeiter->second << "] is " << std::list<argument_t>(
    660               MapValueIterator<TypePair_Argument_Map_t::iterator>(range_t.first),
    661               MapValueIterator<TypePair_Argument_Map_t::iterator>(range_t.second)
    662               ) << ".");
    663       // the first key is always easy and is pivot which the rest has to be associated to
    664       if (typeiter == UniqueTypes.begin()) {
    665         const argument_t & arg = range_t.first->second;
    666         if ((typeiter->first == arg.types.first) && (typeiter->second == arg.types.second)) {
    667           // store in correct order
    668           *firstindex = arg.indices.first;
    669           *secondindex = arg.indices.second;
    670         } else {
    671           // store in flipped order
    672           *firstindex = arg.indices.second;
    673           *secondindex = arg.indices.first;
    674         }
    675         argumentbunch.push_back(arg);
    676         argument_map.erase(range_t.first);
    677         LOG(4, "DEBUG: Gathered first argument " << arg << ".");
    678       } else {
    679         // go through the range and pick the first argument matching the index constraints
    680         for (TypePair_Argument_Map_t::iterator argiter = range_t.first;
    681             argiter != range_t.second; ++argiter) {
    682           // seconditer may be -1 still
    683           const argument_t &arg = argiter->second;
    684           if (arg.indices.first == *firstindex) {
    685             if ((arg.indices.second == *secondindex) || (*secondindex == (size_t)-1)) {
    686               if (*secondindex == (size_t)-1)
    687                 *secondindex = arg.indices.second;
    688               argumentbunch.push_back(arg);
    689               argument_map.erase(argiter);
    690               LOG(4, "DEBUG: Gathered another argument " << arg << ".");
    691               break;
    692             }
    693           } else if ((arg.indices.first == *secondindex) || (*secondindex == (size_t)-1)) {
    694             if (arg.indices.second == *firstindex) {
    695               if (*secondindex == (size_t)-1)
    696                 *secondindex = arg.indices.first;
    697               argumentbunch.push_back(arg);
    698               argument_map.erase(argiter);
    699               LOG(4, "DEBUG: Gathered another (flipped) argument " << arg << ".");
    700               break;
    701             }
    702           }
    703         }
    704       }
    705       // move along in indices and check bounds
    706       ++secondindex;
    707       if (secondindex == indices.end()) {
    708         ++firstindex;
    709         if (firstindex != indices.end()-1)
    710           secondindex = firstindex+1;
    711       }
    712     }
    713     ASSERT( (firstindex == indices.end()-1) && (secondindex == indices.end()),
    714         "Extractors::reorderArgumentsByParticleTypes() - we have not gathered enough arguments.");
    715     ASSERT( argumentbunch.size() == UniqueTypes.size(),
    716         "Extractors::reorderArgumentsByParticleTypes() - we have not gathered enough arguments.");
    717     // place bunch of arguments in return args
    718     LOG(3, "DEBUG: Given types " << _types << " and found indices " << indices << ".");
    719     LOG(3, "DEBUG: Final bunch of arguments is " << argumentbunch << ".");
    720     returnargs.insert(returnargs.end(), argumentbunch.begin(), argumentbunch.end());
    721   }
    722 
    723   return returnargs;
    724 }
    725 
    726 FunctionModel::arguments_t Extractors::filterArgumentsByParticleTypes(
    727     const FunctionModel::arguments_t &args,
    728     const ParticleTypes_t &_types
    729     )
    730 {
    731   typedef std::list< argument_t > ListArguments_t;
    732   ListArguments_t availableList(args.begin(), args.end());
    733   LOG(2, "DEBUG: Initial list of args is " << args << ".");
    734 
    735 
    736   // TODO: fill a lookup map such that we don't have O(M^3) scaling, if M is number
    737   // of types (and we always must have M(M-1)/2 args) but O(M^2 log(M)). However, as
    738   // M is very small (<=3), this is not necessary fruitful now.
    739 //  typedef ParticleTypes_t firsttype;
    740 //  typedef ParticleTypes_t secondtype;
    741 //  typedef std::map< firsttype, std::map< secondtype, boost::ref(args) > > ArgsLookup_t;
    742 //  ArgsLookup_t ArgsLookup;
    743 
    744   // basically, we have two choose any two pairs out of types but only those
    745   // where the first is less than the latter. Hence, we start the second
    746   // iterator at the current position of the first one and skip the equal case.
    747   FunctionModel::arguments_t returnargs;
    748   returnargs.reserve(args.size());
    749   for (ParticleTypes_t::const_iterator firstiter = _types.begin();
    750       firstiter != _types.end();
    751       ++firstiter) {
    752     for (ParticleTypes_t::const_iterator seconditer = firstiter;
    753         seconditer != _types.end();
    754         ++seconditer) {
    755       if (seconditer == firstiter)
    756         continue;
     769      LOG(3, "DEBUG: Looking for (" << *firstiter << "," << *seconditer << ") in all args.");
    757770
    758771      // search the right one in _args (we might allow switching places of
    … …  
    760773      ListArguments_t::iterator iter = availableList.begin();
    761774      while (iter != availableList.end()) {
    762         LOG(3, "DEBUG: Current args is " << *iter << ".");
     775        LOG(4, "DEBUG: Current args is " << *iter << ".");
    763776        if ((iter->types.first == *firstiter)
    764777              && (iter->types.second == *seconditer)) {
    765           returnargs.push_back( *iter );
     778          allargs.push_back( *iter );
    766779          iter = availableList.erase(iter);
    767780          LOG(4, "DEBUG: Accepted argument.");
    768781        } else if ((iter->types.first == *seconditer)
    769782              && (iter->types.second == *firstiter)) {
    770           returnargs.push_back( *iter );
     783          allargs.push_back( *iter );
    771784          iter = availableList.erase(iter);
    772785          LOG(4, "DEBUG: Accepted (flipped) argument.");
    … …  
    778791    }
    779792  }
    780   LOG(2, "DEBUG: Final list of args is " << returnargs << ".");
     793  LOG(2, "DEBUG: Final list of args is " << allargs << ".");
     794
     795  // first, we bring together tuples of distances that belong together
     796  FunctionModel::list_of_arguments_t singlelist_allargs;
     797  singlelist_allargs.push_back(allargs);
     798  FunctionModel::list_of_arguments_t sortedargs =
     799      reorderArgumentsByParticleTypes(singlelist_allargs, _types);
     800  ASSERT( sortedargs.size() == (size_t)1,
     801      "Extractors::filterArgumentsByParticleTypes() - reordering did not generate a single list.");
     802  // then we split up the tuples of arguments and place each into single list
     803  FunctionModel::list_of_arguments_t returnargs;
     804  FunctionModel::arguments_t::const_iterator argiter = sortedargs.begin()->begin();
     805  const size_t num_types = _types.size();
     806  const size_t args_per_tuple = num_types * (num_types-1) / 2;
     807  while (argiter != sortedargs.begin()->end()) {
     808    FunctionModel::arguments_t currenttuple(args_per_tuple);
     809    const FunctionModel::arguments_t::const_iterator startiter = argiter;
     810    std::advance(argiter, args_per_tuple);
     811#ifndef NDEBUG
     812    FunctionModel::arguments_t::const_iterator endoutiter =
     813#endif
     814        std::copy(startiter, argiter, currenttuple.begin());
     815    ASSERT( endoutiter == currenttuple.end(),
     816        "Extractors::filterArgumentsByParticleTypes() - currenttuple not initialized to right size.");
     817    returnargs.push_back(currenttuple);
     818  }
     819
     820  LOG(2, "DEBUG: We have generated " << returnargs.size() << " tuples of distances.");
    781821
    782822  return returnargs;
    … …  
    807847}
    808848
     849FunctionModel::list_of_arguments_t Extractors::concatenateListOfArguments(
     850    const FunctionModel::list_of_arguments_t &firstlistargs,
     851    const FunctionModel::list_of_arguments_t &secondlistargs)
     852{
     853  FunctionModel::list_of_arguments_t listargs(firstlistargs);
     854  listargs.insert(listargs.end(), secondlistargs.begin(), secondlistargs.end());
     855  return listargs;
     856}

  • src/FunctionApproximation/Extractors.hpp

    rc73e35 rd8821e  
    2222
    2323/** Namespace containing all simple extractor functions.
     24 *
     25 * Extractor functions extract distances from a given fragment matching with
     26 * a given set of particle types (i.e. elements, e.h. H2O).
     27 * Filter functions extract a subset of distances from a given set of distances
     28 * to be used with a specific model.
     29 *
     30 * To this end, each FunctionModel has both a filter and an extractor function.
     31 *
     32 * The functions in this namespace act as helpers or basic building blocks in
     33 * constructing such filters and extractors.
    2434 *
    2535 */
    … …  
    294304  /** Reorder arguments by increasing distance.
    295305   *
    296    * \param args arguments to reorder
     306   * \param listargs list of arguments to reorder each
    297307   * \return reordered args
    298308   */
    299   FunctionModel::arguments_t reorderArgumentsByIncreasingDistance(
    300       const FunctionModel::arguments_t &args
     309  FunctionModel::list_of_arguments_t reorderArgumentsByIncreasingDistance(
     310      const FunctionModel::list_of_arguments_t &listargs
    301311      );
    302312
    … …  
    308318   * \sa filterArgumentsByParticleTypes()
    309319   *
     320   * \param listargs list of arguments to reorder each
     321   * \param _types particle type vector
     322   * \return reordered args
     323   */
     324  FunctionModel::list_of_arguments_t reorderArgumentsByParticleTypes(
     325      const FunctionModel::list_of_arguments_t &eachargs,
     326      const ParticleTypes_t &_types
     327      );
     328
     329  /** Filter arguments according to types, allowing multiples.
     330   *
     331   * If particle types is (0,1,2) and three arguments, each with a pair of types,
     332   * are given, then the alignment will be: (0,1), (0,2), and (1,2).
     333   *
    310334   * \param args arguments to reorder
    311335   * \param _types particle type vector
    312    * \return reordered args
    313    */
    314   FunctionModel::arguments_t reorderArgumentsByParticleTypes(
    315       const FunctionModel::arguments_t &args,
    316       const ParticleTypes_t &_types
    317       );
    318 
    319   /** Filter arguments according to types.
    320    *
    321    * If particle types is (0,1,2) and three arguments, each with a pair of types,
    322    * are given, then the alignment will be: (0,1), (0,2), and (1,2).
    323    *
    324    * \param args arguments to reorder
    325    * \param _types particle type vector
    326    * \return filtered args
    327    */
    328   FunctionModel::arguments_t filterArgumentsByParticleTypes(
     336   * \return filtered list of args
     337   */
     338  FunctionModel::list_of_arguments_t filterArgumentsByParticleTypes(
    329339      const FunctionModel::arguments_t &args,
    330340      const ParticleTypes_t &_types
    … …  
    351361      const FunctionModel::arguments_t &secondargs);
    352362
     363  /** Combines two argument lists by concatenation.
     364   *
     365   * @param firstlistargs first list of argument tuples
     366   * @param secondlistargs second list of argument tuples
     367   * @return concatenated lists
     368   */
     369  FunctionModel::list_of_arguments_t concatenateListOfArguments(
     370      const FunctionModel::list_of_arguments_t &firstlistargs,
     371      const FunctionModel::list_of_arguments_t &secondlistargs);
     372
    353373}; /* namespace Extractors */
    354374

  • src/FunctionApproximation/FunctionApproximation.cpp

    rc73e35 rd8821e  
    6464{}
    6565
    66 void FunctionApproximation::setTrainingData(const inputs_t &input, const outputs_t &output)
     66void FunctionApproximation::setTrainingData(const filtered_inputs_t &input, const outputs_t &output)
    6767{
    6868  ASSERT( input.size() == output.size(),
    … …  
    334334      +" outputs but we provide "+toString(output_data.size())+".");
    335335  if (!output_data.empty()) {
    336     inputs_t::const_iterator initer = input_data.begin();
     336    filtered_inputs_t::const_iterator initer = input_data.begin();
    337337    outputs_t::const_iterator outiter = output_data.begin();
    338338    size_t index = 0;
    … …  
    362362      +" outputs but we provide "+toString(output_data.size())+".");
    363363  if (!output_data.empty()) {
    364     inputs_t::const_iterator initer = input_data.begin();
     364    filtered_inputs_t::const_iterator initer = input_data.begin();
    365365    outputs_t::const_iterator outiter = output_data.begin();
    366366    size_t index = 0;

  • src/FunctionApproximation/FunctionApproximation.hpp

    rc73e35 rd8821e  
    4141 * package.
    4242 *
     43 * \section FunctionApproximation-details Details on the inner workings.
     44 *
     45 *  FunctionApproximation::operator() is the main function that performs the
     46 *  non-linear regression. It consists of the following steps:
     47 *  -# hand given (initial) parameters over to model.
     48 *  -# convert output vector to format suitable to levmar
     49 *  -# allocate memory for levmar to work in
     50 *  -# depending on whether the model is constrained or not and whether we
     51 *   have a derivative, we make use of various levmar functions with prepared
     52 *   parameters.
     53 *  -# memory is free'd and some final infos is given.
     54 *
     55 *  levmar needs to evaluate the model. To this end, FunctionApproximation has
     56 *  two functions whose signatures is such as to match with the one required
     57 *  by the levmar package. Hence,
     58 *  -# FunctionApproximation::LevMarCallback()
     59 *  -# FunctionApproximation::LevMarDerivativeCallback()
     60 *  are used as callbacks by levmar only.
     61 *  These hand over the current set of parameters to the model, then both bind
     62 *  FunctionApproximation::evaluate() and
     63 *  FunctionApproximation::evaluateDerivative(), respectively, and execute
     64 *  FunctionModel::operator() or FunctionModel::parameter_derivative(),
     65 *  respectively.
     66 *
    4367 */
    4468class FunctionApproximation
    … …  
    4771  //!> typedef for a vector of input arguments
    4872  typedef std::vector<FunctionModel::arguments_t> inputs_t;
     73  //!> typedef for a vector of input arguments
     74  typedef std::vector<FunctionModel::list_of_arguments_t> filtered_inputs_t;
    4975  //!> typedef for a vector of output values
    5076  typedef std::vector<FunctionModel::results_t> outputs_t;
    … …  
    86112   *        FunctionApproximation::output_dimension size
    87113   */
    88   void setTrainingData(const inputs_t &input, const outputs_t &output);
     114  void setTrainingData(const filtered_inputs_t &input, const outputs_t &output);
    89115
    90116  /** Setter for the model function to be used in the approximation.
    … …  
    167193
    168194  //!> current input set of training data
    169   inputs_t input_data;
     195  filtered_inputs_t input_data;
    170196  //!> current output set of training data
    171197  outputs_t output_data;

  • src/FunctionApproximation/FunctionModel.hpp

    rc73e35 rd8821e  
    1515
    1616#include <boost/function.hpp>
     17#include <list>
    1718#include <vector>
    1819
    … …  
    5859  //!> typedef for the whole set of parameters of the function
    5960  typedef std::vector<parameter_t> parameters_t;
    60   //!> typedef for the argument vector as input to the function
     61  //!> typedef for the argument vector as input to the function (subset of distances)
    6162  typedef std::vector<argument_t> arguments_t;
     63  //!> typedef for a list of argument vectors as input to the function (list of subsets)
     64  typedef std::list<arguments_t> list_of_arguments_t;
    6265  //!> typedef for a single result degree of freedom
    6366  typedef double result_t;
    … …  
    6568  typedef std::vector<result_t> results_t;
    6669  //!> typedef for a function containing how to extract required information from a Fragment.
    67   typedef boost::function< arguments_t (const Fragment &, const size_t)> extractor_t;
     70  typedef boost::function< list_of_arguments_t (const Fragment &, const size_t)> extractor_t;
    6871  //!> typedef for a function containing how to filter required distances from a full argument list.
    69   typedef boost::function< arguments_t (const arguments_t &)> filter_t;
     72  typedef boost::function< list_of_arguments_t (const arguments_t &)> filter_t;
    7073  //!> typedef for the magic triple function that gets the other two distances for a given argument
    7174  typedef boost::function< std::vector<arguments_t>(const argument_t &, const double)> triplefunction_t;
    … …  
    114117   * \return result of the function
    115118   */
    116   virtual results_t operator()(const arguments_t &arguments) const=0;
     119  virtual results_t operator()(const list_of_arguments_t &arguments) const=0;
    117120
    118121  /** Evaluates the derivative of the function with the given \a arguments
    … …  
    124127   *         input
    125128   */
    126   virtual results_t parameter_derivative(const arguments_t &arguments, const size_t index) const=0;
     129  virtual results_t parameter_derivative(const list_of_arguments_t &arguments, const size_t index) const=0;
    127130
    128131  /** States whether lower and upper boundaries should be used to constraint
    … …  
    152155   * \return bound function extracting distances from a fragment
    153156   */
    154   virtual extractor_t getSpecificExtractor() const=0;
    155 
    156   /** Returns a bound function to be used with TrainingData, extracting distances
    157    * from a Fragment.
    158    *
    159    * \return bound function extracting distances from a fragment
    160    */
    161157  virtual filter_t getSpecificFilter() const=0;
    162158

  • src/FunctionApproximation/TrainingData.cpp

    rc73e35 rd8821e  
    4141#include <algorithm>
    4242#include <boost/bind.hpp>
     43#include <boost/foreach.hpp>
    4344#include <boost/lambda/lambda.hpp>
    4445#include <iostream>
     46#include <sstream>
    4547
    4648#include "CodePatterns/Assert.hpp"
    … …  
    4951
    5052#include "Fragmentation/Summation/SetValues/Fragment.hpp"
     53#include "FunctionApproximation/FunctionArgument.hpp"
    5154#include "FunctionApproximation/FunctionModel.hpp"
    5255#include "FunctionApproximation/Extractors.hpp"
    … …  
    6568    EnergyVector.push_back( FunctionModel::results_t(1, energy) );
    6669    // filter distances out of list of all arguments
    67     FunctionModel::arguments_t args = filter(all_args);
     70    FunctionModel::list_of_arguments_t args = filter(all_args);
    6871    LOG(3, "DEBUG: Filtered arguments are " << args << ".");
    6972    ArgumentVector.push_back( args );
    … …  
    7578  double L2sum = 0.;
    7679
    77   FunctionApproximation::inputs_t::const_iterator initer = ArgumentVector.begin();
    78   FunctionApproximation::outputs_t::const_iterator outiter = EnergyVector.begin();
     80  FilteredInputVector_t::const_iterator initer = ArgumentVector.begin();
     81  OutputVector_t::const_iterator outiter = EnergyVector.begin();
    7982  for (; initer != ArgumentVector.end(); ++initer, ++outiter) {
    8083    const FunctionModel::results_t result = model((*initer));
    … …  
    8992  double Lmax = 0.;
    9093//  size_t maxindex = -1;
    91   FunctionApproximation::inputs_t::const_iterator initer = ArgumentVector.begin();
    92   FunctionApproximation::outputs_t::const_iterator outiter = EnergyVector.begin();
     94  FilteredInputVector_t::const_iterator initer = ArgumentVector.begin();
     95  OutputVector_t::const_iterator outiter = EnergyVector.begin();
    9396  for (; initer != ArgumentVector.end(); ++initer, ++outiter) {
    9497    const FunctionModel::results_t result = model((*initer));
    … …  
    105108}
    106109
     110const TrainingData::L2ErrorConfigurationIndexMap_t
     111TrainingData::getWorstFragmentMap(
     112      const FunctionModel &model,
     113      const range_t &range) const
     114{
     115  L2ErrorConfigurationIndexMap_t WorseFragmentMap;
     116  // fragments make it into the container in reversed order, hence count from top down
     117  size_t index= std::distance(range.first, range.second)-1;
     118  InputVector_t::const_iterator distanceiter = DistanceVector.begin();
     119  FilteredInputVector_t::const_iterator initer = ArgumentVector.begin();
     120  OutputVector_t::const_iterator outiter = EnergyVector.begin();
     121  for (; initer != ArgumentVector.end(); ++initer, ++outiter, ++distanceiter) {
     122    // calculate value from potential
     123    const FunctionModel::list_of_arguments_t &args = *initer;
     124    const FunctionModel::results_t result = model(args);
     125    const double energy = (*outiter)[0];
     126
     127    // insert difference into map
     128    const double error = fabs(energy - result[0]);
     129    WorseFragmentMap.insert( std::make_pair( error, index-- ) );
     130
     131    {
     132      // give only the distances in the debugging text
     133      std::stringstream streamargs;
     134      BOOST_FOREACH (argument_t arg, *distanceiter) {
     135        streamargs << " " << arg.distance;
     136      }
     137      LOG(2, "DEBUG: frag.#" << index+1 << "'s error is |" << energy << " - " << result[0]
     138          << "| = " << error << " for args " << streamargs.str() << ".");
     139    }
     140  }
     141
     142  return WorseFragmentMap;
     143}
     144
    107145const TrainingData::DistanceEnergyTable_t TrainingData::getDistanceEnergyTable() const
    108146{
    … …  
    111149  /// extract distance member variable from argument_t and first value from results_t
    112150  OutputVector_t::const_iterator ergiter = EnergyVector.begin();
    113   for (InputVector_t::const_iterator iter = ArgumentVector.begin();
    114       iter != ArgumentVector.end(); ++iter, ++ergiter) {
     151  for (InputVector_t::const_iterator iter = DistanceVector.begin();
     152      iter != DistanceVector.end(); ++iter, ++ergiter) {
    115153    ASSERT( ergiter != EnergyVector.end(),
    116154        "TrainingData::getDistanceEnergyTable() - less output than input values.");

  • src/FunctionApproximation/TrainingData.hpp

    rc73e35 rd8821e  
    2727 * Fragment.
    2828 *
    29  * In TrainingData::operator() we construct first all pair-wise distances as
    30  * list of all arguments. Then, these are filtered depending on the specific
    31  * FunctionModel's Filter and only these are handed to down to evaluate it.
     29 * TrainingData::operator() takes the set of all possible pair-wise  distances
     30 * (InputVector_t) and transforms it via the given filter into a list of subsets
     31 * of distances (FilteredInputVector_t) that is feedable to the model.
    3232 *
    3333 */
    … …  
    4141  //!> Training tuple input vector pair
    4242  typedef FunctionApproximation::inputs_t InputVector_t;
     43  //!> Training tuple modified input vector pair
     44  typedef FunctionApproximation::filtered_inputs_t FilteredInputVector_t;
    4345  //!> Training tuple output vector pair
    4446  typedef FunctionApproximation::outputs_t OutputVector_t;
    4547  //!> Typedef for a table with columns of all distances and the energy
    4648  typedef std::vector< std::vector<double> > DistanceEnergyTable_t;
     49  //!> Typedef for a map of each fragment with error.
     50  typedef std::multimap< double, size_t > L2ErrorConfigurationIndexMap_t;
    4751
    4852public:
    … …  
    7276   * \return const ref to training tuple of input vector
    7377   */
    74   const InputVector_t& getTrainingInputs() const {
     78  const FilteredInputVector_t& getTrainingInputs() const {
    7579    return ArgumentVector;
    7680  }
    … …  
    114118  const double getLMaxError(const FunctionModel &model) const;
    115119
     120  /** Calculate the Lmax error of a given \a model against the stored training data.
     121   *
     122   * \param model model whose Lmax error to calculate
     123   * \param range given range within a HomologyContainer of homologous fragments
     124   * \return map with L2 error per configuration
     125   */
     126  const L2ErrorConfigurationIndexMap_t getWorstFragmentMap(
     127      const FunctionModel &model,
     128      const range_t &range) const;
     129
    116130  /** Creates a table of columns with all distances and the energy.
    117131   *
    … …  
    129143  OutputVector_t EnergyVector;
    130144  //!> list of all filtered arguments over all tuples
    131   InputVector_t ArgumentVector;
     145  FilteredInputVector_t ArgumentVector;
    132146  //!> function to be used for training input data extraction from a fragment
    133147  const FunctionModel::filter_t filter;

  • src/FunctionApproximation/unittests/ExtractorsUnitTest.cpp

    rc73e35 rd8821e  
    182182  args.push_back(arg);
    183183  CPPUNIT_ASSERT_EQUAL( (size_t)3, args.size() );
     184  FunctionModel::list_of_arguments_t listargs(1, args);
    184185
    185186  // reorder
    186   FunctionModel::arguments_t args_sorted =
    187       Extractors::reorderArgumentsByIncreasingDistance(args);
     187  FunctionModel::list_of_arguments_t listargs_sorted =
     188      Extractors::reorderArgumentsByIncreasingDistance(listargs);
     189  CPPUNIT_ASSERT_EQUAL( (size_t)1, listargs_sorted.size() );
     190  FunctionModel::arguments_t args_sorted = *(listargs_sorted.begin());
    188191  CPPUNIT_ASSERT_EQUAL( (size_t)3, args_sorted.size() );
    189192  CPPUNIT_ASSERT_EQUAL( args[0].distance, args_sorted[0].distance );

  • src/Jobs/MPQCCommandJob.hpp

    rc73e35 rd8821e  
    1818#include <string>
    1919
     20#ifdef HAVE_JOBMARKET
    2021#include "JobMarket/Results/FragmentResult.hpp"
    2122#include "JobMarket/Jobs/SystemCommandJob.hpp"
     23#else
     24#include "Jobs/JobMarket/FragmentResult.hpp"
     25#include "Jobs/JobMarket/SystemCommandJob.hpp"
     26#endif
    2227
    2328#include "Fragmentation/Summation/Containers/MPQCData.hpp"
    2429
    2530class MPQCCommandJobTest;
     31class MPQCCommandFragmentController;
    2632
    2733/** This class calls mpqc for solving a specific Hartree Fock problem.
    … …  
    3238  //!> grant unit test access
    3339  friend class MPQCCommandJobTest;
     40  //!> grant access to controller
     41  friend class MPQCCommandFragmentController;
    3442public:
    3543  MPQCCommandJob(const std::string &_inputfile, const JobId_t _JobId, const std::string &_command = std::string("mpqc"));
    … …  
    4351
    4452private:
     53  //!> Allow controller access to changing the commands
     54  void setCommand(const std::string &_command)
     55  { command = _command; }
     56
     57  //!> Allow controller access to changing the suffix
     58  void setSuffix(const std::string &_suffix)
     59  { suffix = _suffix; }
     60
    4561  //!> private default cstor only for serializatio
    4662  MPQCCommandJob();

  • src/Jobs/MPQCJob.hpp

    rc73e35 rd8821e  
    1717#include "boost/serialization/export.hpp"
    1818
     19#ifdef HAVE_JOBMARKET
    1920#include "JobMarket/Jobs/FragmentJob.hpp"
     21#else
     22#include "Jobs/JobMarket/FragmentJob.hpp"
     23#endif
    2024
    2125#include "Fragmentation/Summation/SetValues/SamplingGridProperties.hpp"

  • src/Jobs/Makefile.am

    rc73e35 rd8821e  
    22# Also indentation by a single tab
    33
    4 JOBSSOURCE = \
     4JOBSSOURCE =
     5if CONDJOBMARKET
     6else
     7JOBSSOURCE += \
     8        Jobs/JobMarket/FragmentJob.cpp \
     9        Jobs/JobMarket/FragmentResult.cpp \
     10        Jobs/JobMarket/JobId.cpp \
     11        Jobs/JobMarket/SystemCommandJob.cpp
     12endif
     13JOBSSOURCE += \
    514        Jobs/MPQCCommandJob.cpp \
    615        Jobs/MPQCJob.cpp
     16if CONDJOBMARKET
    717if CONDVMG
    818JOBSSOURCE += \
    … …  
    1323        Jobs/WindowGrid_converter.cpp
    1424endif
     25endif
    1526
    16 JOBSHEADER = \
     27JOBSHEADER =
     28if CONDJOBMARKET
     29else
     30JOBSHEADER += \
     31        Jobs/JobMarket/FragmentJob.hpp \
     32        Jobs/JobMarket/FragmentResult.hpp \
     33        Jobs/JobMarket/JobId.hpp \
     34        Jobs/JobMarket/SystemCommandJob.hpp \
     35        Jobs/JobMarket/types.hpp
     36endif
     37JOBSHEADER += \
    1738        Jobs/MPQCCommandJob.hpp \
    1839        Jobs/MPQCCommandJob_binding.hpp \
    1940        Jobs/MPQCJob.hpp \
    2041        Jobs/MPQCJob_binding.hpp
     42if CONDJOBMARKET
    2143if CONDVMG
    2244JOBSHEADER += \
    … …  
    2951        Jobs/WindowGrid_converter.hpp
    3052endif
     53endif
    3154
    3255lib_LTLIBRARIES += libMolecuilderJobs.la
    … …  
    3558libMolecuilderJobs_la_CPPFLAGS = ${BOOST_CPPFLAGS} ${CodePatterns_CFLAGS} ${JobMarket_CFLAGS} -Dvmg_float=double -Dvmg_int=int $(VMG_CFLAGS)
    3659libMolecuilderJobs_la_LDFLAGS = $(AM_LDFLAGS) \
     60  $(BOOST_IOSTREAMS_LDFLAGS) \
    3761  $(BOOST_SERIALIZATION_LDFLAGS) \
    38   $(BOOST_SYSTEM_LDFLAGS) \
    39   $(JobMarket_LDFLAGS) \
     62  $(BOOST_SYSTEM_LDFLAGS)
     63if CONDJOBMARKET
     64libMolecuilderJobs_la_LDFLAGS += \
     65  $(JobMarket_LDFLAGS)
     66endif
     67
     68libMolecuilderJobs_la_LDFLAGS += \
    4069  $(CodePatterns_LDFLAGS)
    4170libMolecuilderJobs_la_LIBADD = \
    4271  libMolecuilderFragmentationSummation.la
     72if CONDJOBMARKET
    4373if CONDVMG
    4474libMolecuilderJobs_la_LIBADD += \
    4575  $(VMG_LIBS)
    4676endif
     77endif
     78
     79if CONDJOBMARKET
    4780libMolecuilderJobs_la_LIBADD += \
    48   $(JobMarket_LIBS) \
     81  $(JobMarket_LIBS)
     82endif
     83libMolecuilderJobs_la_LIBADD += \
     84  $(BOOST_IOSTREAMS_LIBS) \
    4985  $(BOOST_SERIALIZATION_LIBS) \
    5086  $(BOOST_SYSTEM_LIBS) \
    … …  
    91127#pkgconfigdir = $(libdir)/pkgconfig
    92128#pkgconfig_DATA = $(top_builddir)/MoleCuilder.pc
    93  
     129

  • src/Makefile.am

    rc73e35 rd8821e  
    1616include Filling/Makefile.am
    1717include Fragmentation/Makefile.am
     18include Fragmentation/Automation/Makefile.am
    1819include Fragmentation/Summation/Containers/Makefile.am
    1920include Fragmentation/Summation/Converter/Makefile.am
    … …  
    2324include Graph/Makefile.am
    2425include Helpers/Makefile.am
    25 
    26 if CONDJOBMARKET
    27 include Fragmentation/Automation/Makefile.am
    2826include Jobs/Makefile.am
    29 endif
    3027
    3128if CONDPYTHON
    … …  
    5754
    5855DESCRIPTORSOURCE = \
    59         Descriptors/AtomDescriptor.cpp \
     56  Descriptors/AtomDescriptor.cpp \
    6057  Descriptors/AtomIdDescriptor.cpp \
    6158  Descriptors/AtomOfMoleculeDescriptor.cpp \
    … …  
    288285        -I$(top_srcdir)/LinearAlgebra/src
    289286
    290 bin_PROGRAMS += molecuilder joiner analyzer
     287bin_PROGRAMS += molecuilder
    291288EXTRA_PROGRAMS = unity
    292289
    … …  
    313310        Actions/GlobalListOfActions.hpp \
    314311        Actions/ActionHistory.hpp
    315 pyMoleCuilder_la_CPPFLAGS = ${BOOST_CPPFLAGS} ${CodePatterns_CFLAGS} -I$(PYTHON_INCLUDE_DIR)
     312pyMoleCuilder_la_CPPFLAGS = ${BOOST_CPPFLAGS} ${CodePatterns_CFLAGS} $(JobMarket_CFLAGS) -I$(PYTHON_INCLUDE_DIR)
    316313pyMoleCuilder_la_LDFLAGS = -module -avoid-version -shared $(BOOST_PYTHON_LDFLAGS)
    317314pyMoleCuilder_la_LIBADD = \
    … …  
    403400endif
    404401
    405 joiner_SOURCES = Fragmentation/joiner.cpp Fragmentation/datacreator.cpp Fragmentation/datacreator.hpp
    406 joiner_CPPFLAGS = $(AM_CPPFLAGS)
    407 joiner_LDFLAGS = $(AM_LDFLAGS) $(BOOST_THREAD_LDFLAGS)
    408 joiner_LDADD = \
    409         libMolecuilderFragmentation.la \
    410         libMolecuilderFragmentation_KeysetsContainer.la \
    411         libMolecuilderHelpers.la \
    412         libMolecuilderElement.la \
    413         $(top_builddir)/LinearAlgebra/src/LinearAlgebra/libLinearAlgebra.la \
    414         ${CodePatterns_LIBS} \
    415         $(BOOST_THREAD_LIBS)
    416 
    417 analyzer_SOURCES = Fragmentation/analyzer.cpp Fragmentation/datacreator.cpp Fragmentation/datacreator.hpp
    418 analyzer_CPPFLAGS = $(AM_CPPFLAGS)
    419 analyzer_LDFLAGS = $(AM_LDFLAGS) $(BOOST_THREAD_LDFLAGS)
    420 analyzer_LDADD = \
    421         libMolecuilderFragmentation.la \
    422         libMolecuilderFragmentation_KeysetsContainer.la \
    423         libMolecuilderHelpers.la \
    424         libMolecuilderElement.la \
    425         $(top_builddir)/LinearAlgebra/src/LinearAlgebra/libLinearAlgebra.la \
    426         ${CodePatterns_LIBS} \
    427         $(BOOST_THREAD_LIBS)
    428 
    429402if CONDJOBMARKET
    430403CONTROLLERSOURCE = \

  • src/Potentials/CompoundPotential.cpp

    rc73e35 rd8821e  
    228228{
    229229  arguments_by_model_t partial_args;
    230   // go through each model and have it filter out its arguments
     230  // go through each model and have it filter out its arguments, this already
     231  // returns a list of tuples associated with the specific model
    231232  for(models_t::const_iterator modeliter = models.begin();
    232233      modeliter != models.end(); ++modeliter) {
    233234    FunctionModel::filter_t filterfunction = (*modeliter)->getSpecificFilter();
    234     arguments_t tempargs = filterfunction(arguments);
     235    list_of_arguments_t tempargs = filterfunction(arguments);
    235236    // then split up all the bunches, too.
    236     arguments_t::const_iterator advanceiter = tempargs.begin();
    237     for (arguments_t::const_iterator argiter = tempargs.begin();
    238         argiter != tempargs.end(); argiter = advanceiter) {
    239       advanceiter += (*modeliter)->getSpecificArgumentCount();
     237    for (list_of_arguments_t::const_iterator argiter = tempargs.begin();
     238        argiter != tempargs.end(); ++argiter) {
     239      const arguments_t &args = *argiter;
    240240      partial_args.push_back(
    241241          std::make_pair(
    242242              *modeliter,
    243               arguments_t(argiter, advanceiter)
     243              args
    244244            )
    245245          );
    … …  
    251251
    252252CompoundPotential::arguments_by_model_t CompoundPotential::splitUpArgumentsByModels(
    253     const arguments_t &arguments) const
     253    const list_of_arguments_t &listarguments) const
    254254{
    255255  arguments_by_model_t partial_args;
    256   arguments_t::const_iterator argiter = arguments.begin();
    257256  particletypes_per_model_t::const_iterator typesiter = particletypes_per_model.begin();
    258257  models_t::const_iterator modeliter = models.begin();
    259258
    260   // add constant model (which is always first model) with empty args if present
     259  /// add constant model (which is always first model) with empty args if present
    261260  if (typesiter->empty()) {
    262261    partial_args.push_back(
    … …  
    266265    ++typesiter;
    267266  }
     267
    268268  // then check other models
    269   while (argiter != arguments.end()) {
     269  /// we only have to check whether the current model still matches or whether
     270  /// have to use the next model.
     271  for (list_of_arguments_t::const_iterator argiter = listarguments.begin();
     272      argiter != listarguments.end(); ++argiter) {
     273    const arguments_t &arguments = *argiter;
    270274    if (typesiter+1 != particletypes_per_model.end()) {
    271275      // check whether next argument bunch is for same model or different one
    … …  
    275279
    276280      // we always expect N(N-1)/2 distances for N particle types
    277       arguments_t::const_iterator enditer = argiter+(types.size()*(types.size()-1)/2);
    278       arguments_t::const_iterator nextenditer = argiter+(nexttypes.size()*(nexttypes.size()-1)/2);
    279       arguments_t args(argiter, enditer);
    280       arguments_t nextargs(argiter, nextenditer);
    281       if (types.size() < nexttypes.size()) {
    282         if (areValidArguments(nexttypes, nextargs)) {
     281      // check first from sizes alone
     282      const size_t tuplesize = types.size()*(types.size()-1)/2;
     283      const size_t nexttuplesize = nexttypes.size()*(nexttypes.size()-1)/2;
     284      if ((tuplesize != nexttuplesize)) {
     285        if ((arguments.size() == tuplesize) &&  areValidArguments(types, arguments)) {
     286          // only former still matches, don't increment
     287          partial_args.push_back(
     288              std::make_pair(*modeliter, arguments)
     289              );
     290        } else if ((arguments.size() == nexttuplesize) &&  areValidArguments(nexttypes, arguments)) {
    283291          // latter matches, increment
    284292          ++typesiter;
    285293          partial_args.push_back(
    286               std::make_pair(*(++modeliter), arguments_t(argiter, nextenditer))
     294              std::make_pair(*(++modeliter), arguments)
    287295              );
    288           argiter = nextenditer;
    289         } else if (areValidArguments(types, args)) {
    290           // only former matches, don't increment
    291           partial_args.push_back(
    292               std::make_pair(*modeliter, arguments_t(argiter, enditer))
    293               );
    294           argiter = enditer;
    295         } else
     296        } else {
    296297          ASSERT(0,
    297               "CompoundPotential::splitUpArgumentsByModels() - neither type matches its argument bunch.");
    298       } else {
    299         if (areValidArguments(types, args)) {
    300           // only former matches, don't increment
    301           partial_args.push_back(
    302               std::make_pair(*modeliter, arguments_t(argiter, enditer))
    303               );
    304           argiter = enditer;
    305         } else if (areValidArguments(nexttypes, nextargs)) {
    306           // latter matches, increment
    307           ++typesiter;
    308           partial_args.push_back(
    309               std::make_pair(*(++modeliter), arguments_t(argiter, nextenditer))
    310               );
    311           argiter = nextenditer;
     298              "CompoundPotential::splitUpArgumentsByModels() - neither this model nor next model match (size) with current tuple.");
     299        }
     300      } else { // same size, now we have to check the types individually
     301        size_t encodeValidity = 0;
     302        encodeValidity += 1*areValidArguments(types, arguments);
     303        encodeValidity += 2*areValidArguments(nexttypes, arguments);
     304
     305        switch (encodeValidity) {
     306          case 1:
     307            // only former still matches, don't increment
     308            partial_args.push_back(
     309                std::make_pair(*modeliter, arguments)
     310                );
     311            break;
     312          case 2:
     313            ++typesiter;
     314            partial_args.push_back(
     315                std::make_pair(*(++modeliter), arguments)
     316                );
     317            break;
     318          case 0:
     319          case 3:
     320          default:
     321            ASSERT(0,
     322                "CompoundPotential::splitUpArgumentsByModels() - neither this model nor next model match (type) with current tuple.");
     323            break;
    312324        }
    313325      }
    314326    } else {
    315327      const SerializablePotential::ParticleTypes_t &types = *typesiter;
    316       // we always expect N(N-1)/2 distances for N particle types
    317       arguments_t::const_iterator enditer = argiter+(types.size()*(types.size()-1)/2);
    318       partial_args.push_back(
    319           std::make_pair(*modeliter, arguments_t(argiter, enditer))
    320           );
    321       argiter = enditer;
     328      if (areValidArguments(types, arguments)) {
     329        // only former matches, don't increment
     330        partial_args.push_back(
     331            std::make_pair(*modeliter, arguments)
     332            );
     333      } else {
     334        ASSERT(0,
     335            "CompoundPotential::splitUpArgumentsByModels() - last model does not match with current tuple.");
     336      }
    322337    }
    323338  }
    … …  
    326341}
    327342
    328 CompoundPotential::results_t CompoundPotential::operator()(const arguments_t &arguments) const
     343CompoundPotential::results_t CompoundPotential::operator()(
     344    const list_of_arguments_t &listarguments) const
    329345{
    330346  /// first, we have to split up the given arguments
    331347  arguments_by_model_t partial_args =
    332       splitUpArgumentsByModels(arguments);
     348      splitUpArgumentsByModels(listarguments);
    333349  // print split up argument list for debugging
    334350  if (DoLog(4)) {
    … …  
    348364      iter != partial_args.end(); ++iter) {
    349365    partial_results.push_back(
    350         (*iter->first)(iter->second)
     366        (*iter->first)(
     367            FunctionModel::list_of_arguments_t(1, iter->second))
    351368    );
    352369  }
    … …  
    366383}
    367384
    368 CompoundPotential::results_t CompoundPotential::parameter_derivative(const arguments_t &arguments, const size_t index) const
     385CompoundPotential::results_t CompoundPotential::parameter_derivative(
     386    const list_of_arguments_t &listarguments,
     387    const size_t index) const
    369388{
    370389  // first, we have to split up the given arguments
    371390  arguments_by_model_t partial_args =
    372       splitUpArgumentsByModels(arguments);
     391      splitUpArgumentsByModels(listarguments);
    373392  // then, with each bunch of arguments, we call the specific model
    374393  // get parameter dimensions per model
    … …  
    399418      const size_t indexbase = (iter == dimensions.begin()) ? 0 : *(iter-1);
    400419      CompoundPotential::results_t results =
    401           model->parameter_derivative(argiter->second, index-indexbase);
     420          model->parameter_derivative(
     421              FunctionModel::list_of_arguments_t(1, argiter->second), index-indexbase);
    402422
    403423      // either set results or add
    … …  
    457477}
    458478
    459 FunctionModel::extractor_t CompoundPotential::getSpecificExtractor() const
    460 {
    461   // create initial returnfunction
    462   FunctionModel::extractor_t returnfunction =
    463       boost::bind(&Helpers::returnEmptyArguments);
    464 
    465   // every following fragments combines its arguments with the initial function
    466   for (particletypes_per_model_t::const_iterator typesiter = particletypes_per_model.begin();
    467       typesiter != particletypes_per_model.end(); ++typesiter) {
    468     Fragment::charges_t charges;
    469     {
    470       charges.resize((*typesiter).size());
    471       std::transform((*typesiter).begin(), (*typesiter).end(),
    472           charges.begin(), boost::lambda::_1);
    473     }
    474 
    475     returnfunction =
    476           boost::bind(&Extractors::concatenateArguments,
    477               boost::bind(returnfunction, _1, _2),
    478               boost::bind(&Extractors::gatherAllDistancesFromFragment,
    479                   boost::bind(&Fragment::getPositions, _1),
    480                   boost::bind(&Fragment::getCharges, _1),
    481                   charges,  // no crefs here as are temporaries!
    482                   _2)
    483         );
    484   }
    485   return returnfunction;
    486 }
    487 
    488479FunctionModel::filter_t CompoundPotential::getSpecificFilter() const
    489480{
    … …  
    491482  // create initial returnfunction
    492483  FunctionModel::filter_t returnfunction =
    493       boost::bind(&Helpers::returnEmptyArguments);
     484      boost::bind(&Helpers::returnEmptyListArguments);
    494485
    495486  // every following fragments combines its arguments with the initial function
    … …  
    497488      modeliter != models.end(); ++modeliter) {
    498489    returnfunction =
    499           boost::bind(&Extractors::concatenateArguments,
     490          boost::bind(&Extractors::concatenateListOfArguments,
    500491              boost::bind(returnfunction, _1),
    501492              boost::bind((*modeliter)->getSpecificFilter(), _1)

  • src/Potentials/CompoundPotential.hpp

    rc73e35 rd8821e  
    2323class TrainingData;
    2424
    25 /** CompoundPotential combines several FunctionModel's as one to allow for
    26  * simultaneous FunctionApproximation to a single fragment/graph.
     25/** CompoundPotential combines several EmipiricalPotential's into one
     26 * FunctionModel to allow for simultaneous FunctionApproximation to a single
     27 * fragment/graph.
    2728 *
    2829 * The CompoundPotential obtains a Graph as parameter to cstor and looks through
    … …  
    3031 * matches. All of these are placed into an internal vector and used for
    3132 * fitting.
     33 *
     34 * The main work is then in operator() and parameter_derivative(), where the
     35 * contained set of models has to be evaluated one after the other:
     36 *
    3237 */
    3338class CompoundPotential :
    … …  
    8893  /** Evaluates the harmonic potential function for the given arguments.
    8994   *
    90    * @param arguments single distance
     95   * @param listarguments list of tuples of distances
    9196   * @return value of the potential function
    9297   */
    93   results_t operator()(const arguments_t &arguments) const;
     98  results_t operator()(const list_of_arguments_t &listarguments) const;
    9499
    95100  /** Evaluates the derivative of the function with the given \a arguments
    96101   * with respect to a specific parameter indicated by \a index.
    97102   *
    98    * \param arguments set of arguments as input variables to the function
     103   * \param arguments list of sets of arguments as input variables to the function
    99104   * \param index derivative of which parameter
    100105   * \return result vector containing the derivative with respect to the given
    101106   *         input
    102107   */
    103   results_t parameter_derivative(const arguments_t &arguments, const size_t index) const;
     108  results_t parameter_derivative(const list_of_arguments_t &listarguments, const size_t index) const;
    104109
    105110  /** States whether lower and upper boundaries should be used to constraint
    … …  
    123128   */
    124129  parameters_t getUpperBoxConstraints() const;
    125 
    126   /** Returns a bound function to be used with TrainingData, extracting distances
    127    * from a Fragment.
    128    *
    129    * Here, we simply concatenate the arguments delivered by the extractors of each model.
    130    *
    131    * \return bound function extracting distances from a fragment
    132    */
    133   FunctionModel::extractor_t getSpecificExtractor() const;
    134130
    135131  /** Returns a bound function to be used with TrainingData, extracting distances
    … …  
    157153  /** Helper function to split up concatenated arguments for operator() calls.
    158154   *
    159    * \param arguments arguments to split up
     155   * \param listarguments list of tuples of distances to assign to a model each
    160156   * \return vector of partial arguments with associated model
    161157   */
    162   arguments_by_model_t splitUpArgumentsByModels(const arguments_t &arguments) const;
     158  arguments_by_model_t splitUpArgumentsByModels(const list_of_arguments_t &listarguments) const;
    163159
    164160  /** Helper function to split up total list of arguments for operator() calls.

  • src/Potentials/EmpiricalPotential.hpp

    rc73e35 rd8821e  
    6060   * parameters.
    6161   *
    62    * \param arguments set of arguments as input variables to the function
     62   * \param listarguments list of sets of arguments as input variables to the function
    6363   * \return result of the function
    6464   */
    65   virtual results_t operator()(const arguments_t &arguments) const=0;
     65  virtual results_t operator()(const list_of_arguments_t &listarguments) const=0;
    6666
    6767  /** Evaluates the derivative of the function with the given \a arguments
    6868   * for each component.
    6969   *
    70    * \param arguments set of arguments as input variables to the function
     70   * \param listarguments list of sets of arguments as input variables to the function
    7171   * \return result vector containing the derivative with respect to each
    7272   *         input comonent of the function
    7373   */
    74   virtual derivative_components_t derivative(const arguments_t &arguments) const=0;
     74  virtual derivative_components_t derivative(const list_of_arguments_t &listarguments) const=0;
    7575
    7676  /** Returns the functor that converts argument_s into the

  • src/Potentials/Specifics/ConstantPotential.cpp

    rc73e35 rd8821e  
    112112ConstantPotential::results_t
    113113ConstantPotential::operator()(
    114     const arguments_t &arguments
     114    const list_of_arguments_t &listarguments
    115115    ) const
    116116{
    117   ASSERT( arguments.size() == 0,
    118       "ConstantPotential::operator() - requires no argument.");
    119   ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
    120       arguments, getParticleTypes()),
    121       "ConstantPotential::operator() - types don't match with ones in arguments.");
    122   const result_t result = params[energy_offset];
    123   return std::vector<result_t>(1, result);
     117  // directly set result as energy constant as independent of arg list
     118  result_t result = params[energy_offset];
     119  for(list_of_arguments_t::const_iterator iter = listarguments.begin();
     120      iter != listarguments.end(); ++iter) {
     121    const arguments_t &arguments = *iter;
     122    ASSERT( arguments.size() == 0,
     123        "ConstantPotential::operator() - requires no argument.");
     124    ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
     125        arguments, getParticleTypes()),
     126        "ConstantPotential::operator() - types don't match with ones in arguments.");
     127    result += 0.;
     128  }
     129  return results_t(1, result);
    124130}
    125131
    126132ConstantPotential::derivative_components_t
    127133ConstantPotential::derivative(
    128     const arguments_t &arguments
     134    const list_of_arguments_t &listarguments
    129135    ) const
    130136{
    131   ASSERT( arguments.size() == 0,
    132       "ConstantPotential::operator() - requires no argument.");
    133   ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
    134       arguments, getParticleTypes()),
    135       "ConstantPotential::operator() - types don't match with ones in arguments.");
    136   derivative_components_t result(1, 0.);
    137   return result;
     137  result_t result = 0.;
     138  for(list_of_arguments_t::const_iterator iter = listarguments.begin();
     139      iter != listarguments.end(); ++iter) {
     140    const arguments_t &arguments = *iter;
     141    ASSERT( arguments.size() == 0,
     142        "ConstantPotential::operator() - requires no argument.");
     143    ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
     144        arguments, getParticleTypes()),
     145        "ConstantPotential::operator() - types don't match with ones in arguments.");
     146    result += 0.;
     147  }
     148  return derivative_components_t(1, result);
    138149}
    139150
    140151ConstantPotential::results_t
    141152ConstantPotential::parameter_derivative(
    142     const arguments_t &arguments,
     153    const list_of_arguments_t &listarguments,
    143154    const size_t index
    144155    ) const
    145156{
    146   ASSERT( arguments.size() == 0,
    147       "ConstantPotential::parameter_derivative() - requires no argument.");
    148   ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
    149       arguments, getParticleTypes()),
    150       "ConstantPotential::operator() - types don't match with ones in arguments.");
    151   switch (index) {
    152     case energy_offset:
    153     {
    154       // Maple result: 1
    155       const result_t result = +1.;
    156       return std::vector<result_t>(1, result);
    157       break;
    158     }
    159     default:
    160       break;
     157  // Maple result: 1
     158  result_t result = 1.;   // energy_offset
     159  for(list_of_arguments_t::const_iterator iter = listarguments.begin();
     160      iter != listarguments.end(); ++iter) {
     161    const arguments_t &arguments = *iter;
     162    ASSERT( arguments.size() == 0,
     163        "ConstantPotential::parameter_derivative() - requires no argument.");
     164    ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
     165        arguments, getParticleTypes()),
     166        "ConstantPotential::operator() - types don't match with ones in arguments.");
     167//    switch (index) {
     168//      case energy_offset:
     169//      {
     170//        // Maple result: 1
     171//        result = +1.;
     172//        break;
     173//      }
     174//      default:
     175//        break;
     176//    }
    161177  }
    162   return std::vector<result_t>(1, 0.);
    163 }
    164 
    165 FunctionModel::extractor_t
    166 ConstantPotential::getSpecificExtractor() const
    167 {
    168   Fragment::charges_t charges;
    169   charges.resize(getParticleTypes().size());
    170   std::transform(getParticleTypes().begin(), getParticleTypes().end(),
    171       charges.begin(), boost::lambda::_1);
    172   FunctionModel::extractor_t returnfunction =
    173       boost::bind(&Extractors::gatherDistancesFromFragment,
    174           boost::bind(&Fragment::getPositions, _1),
    175           boost::bind(&Fragment::getCharges, _1),
    176           charges,
    177           _2);
    178   return returnfunction;
     178  return results_t(1, result);
    179179}
    180180

  • src/Potentials/Specifics/ConstantPotential.hpp

    rc73e35 rd8821e  
    3636  friend class PotentialFactory;
    3737  // some repeated typedefs to avoid ambiguities
     38  typedef FunctionModel::list_of_arguments_t list_of_arguments_t;
    3839  typedef FunctionModel::arguments_t arguments_t;
    3940  typedef FunctionModel::result_t result_t;
    … …  
    9091  /** Evaluates the harmonic potential function for the given arguments.
    9192   *
    92    * @param arguments single distance
     93   * @param listarguments empty list
    9394   * @return value of the potential function
    9495   */
    95   results_t operator()(const arguments_t &arguments) const;
     96  results_t operator()(const list_of_arguments_t &listarguments) const;
    9697
    9798  /** Evaluates the derivative of the potential function.
    9899   *
    99    * @param arguments single distance
     100   * @param listarguments empty list
    100101   * @return vector with derivative with respect to the input degrees of freedom
    101102   */
    102   derivative_components_t derivative(const arguments_t &arguments) const;
     103  derivative_components_t derivative(const list_of_arguments_t &listarguments) const;
    103104
    104105  /** Evaluates the derivative of the function with the given \a arguments
    105106   * with respect to a specific parameter indicated by \a index.
    106107   *
    107    * \param arguments set of arguments as input variables to the function
     108   * \param listarguments empty list
    108109   * \param index derivative of which parameter
    109110   * \return result vector containing the derivative with respect to the given
    110111   *         input
    111112   */
    112   results_t parameter_derivative(const arguments_t &arguments, const size_t index) const;
     113  results_t parameter_derivative(const list_of_arguments_t &listarguments, const size_t index) const;
    113114
    114115  /** Returns the functor that converts argument_s into the
    … …  
    163164    return parameters_t(getParameterDimension(), std::numeric_limits<double>::max());
    164165  }
    165 
    166   /** Returns a bound function to be used with TrainingData, extracting distances
    167    * from a Fragment.
    168    *
    169    * \return bound function extracting distances from a fragment
    170    */
    171   FunctionModel::extractor_t getSpecificExtractor() const;
    172166
    173167  /** Returns a bound function to be used with TrainingData, extracting distances

  • src/Potentials/Specifics/FourBodyPotential_Torsion.cpp

    rc73e35 rd8821e  
    140140FourBodyPotential_Torsion::results_t
    141141FourBodyPotential_Torsion::operator()(
    142     const arguments_t &arguments
     142    const list_of_arguments_t &listarguments
    143143    ) const
    144144{
    145   ASSERT( arguments.size() == getSpecificArgumentCount(),
    146       "FourBodyPotential_Torsion::operator() - requires exactly three arguments.");
    147   ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
    148       arguments, getParticleTypes()),
    149       "FourBodyPotential_Torsion::operator() - types don't match with ones in arguments.");
    150   const argument_t &r_ij = arguments[0]; // 01
    151   const argument_t &r_ik = arguments[1]; // 02
    152   const argument_t &r_il = arguments[2]; // 03
    153   const argument_t &r_jk = arguments[3]; // 12
    154   const argument_t &r_jl = arguments[4]; // 13
    155   const argument_t &r_kl = arguments[5]; // 23
    156   const result_t result =
    157       params[spring_constant]
    158              * Helpers::pow( function_theta(r_ij.distance, r_ik.distance, r_il.distance, r_jk.distance, r_jl.distance, r_kl.distance) - params[equilibrium_distance], 2 );
    159   return std::vector<result_t>(1, result);
     145  result_t result = 0.;
     146  for(list_of_arguments_t::const_iterator iter = listarguments.begin();
     147      iter != listarguments.end(); ++iter) {
     148    const arguments_t &arguments = *iter;
     149    ASSERT( arguments.size() == getSpecificArgumentCount(),
     150        "FourBodyPotential_Torsion::operator() - requires exactly three arguments.");
     151    ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
     152        arguments, getParticleTypes()),
     153        "FourBodyPotential_Torsion::operator() - types don't match with ones in arguments.");
     154    const argument_t &r_ij = arguments[0]; // 01
     155    const argument_t &r_ik = arguments[1]; // 02
     156    const argument_t &r_il = arguments[2]; // 03
     157    const argument_t &r_jk = arguments[3]; // 12
     158    const argument_t &r_jl = arguments[4]; // 13
     159    const argument_t &r_kl = arguments[5]; // 23
     160    result +=
     161        params[spring_constant]
     162               * Helpers::pow( function_theta(r_ij.distance, r_ik.distance, r_il.distance, r_jk.distance, r_jl.distance, r_kl.distance)
     163                   - params[equilibrium_distance], 2 );
     164  }
     165  return results_t(1, result);
    160166}
    161167
    162168FourBodyPotential_Torsion::derivative_components_t
    163169FourBodyPotential_Torsion::derivative(
    164     const arguments_t &arguments
     170    const list_of_arguments_t &listarguments
    165171    ) const
    166172{
    167   ASSERT( arguments.size() == getSpecificArgumentCount(),
    168       "FourBodyPotential_Torsion::operator() - requires exactly three arguments.");
    169   ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
    170       arguments, getParticleTypes()),
    171       "FourBodyPotential_Torsion::operator() - types don't match with ones in arguments.");
    172   derivative_components_t result;
    173   const argument_t &r_ij = arguments[0]; // 01
    174   const argument_t &r_ik = arguments[1]; // 02
    175   const argument_t &r_il = arguments[2]; // 03
    176   const argument_t &r_jk = arguments[3]; // 12
    177   const argument_t &r_jl = arguments[4]; // 13
    178   const argument_t &r_kl = arguments[5]; // 23
    179   result.push_back( 2. * params[spring_constant] * ( function_theta(r_ij.distance, r_ik.distance, r_il.distance, r_jk.distance, r_jl.distance, r_kl.distance) - params[equilibrium_distance]) );
    180   ASSERT( result.size() == 1,
    181       "FourBodyPotential_Torsion::operator() - we did not create exactly one component.");
    182   return result;
     173  result_t result = 0.;
     174  for(list_of_arguments_t::const_iterator iter = listarguments.begin();
     175      iter != listarguments.end(); ++iter) {
     176    const arguments_t &arguments = *iter;
     177    ASSERT( arguments.size() == getSpecificArgumentCount(),
     178        "FourBodyPotential_Torsion::operator() - requires exactly three arguments.");
     179    ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
     180        arguments, getParticleTypes()),
     181        "FourBodyPotential_Torsion::operator() - types don't match with ones in arguments.");
     182    const argument_t &r_ij = arguments[0]; // 01
     183    const argument_t &r_ik = arguments[1]; // 02
     184    const argument_t &r_il = arguments[2]; // 03
     185    const argument_t &r_jk = arguments[3]; // 12
     186    const argument_t &r_jl = arguments[4]; // 13
     187    const argument_t &r_kl = arguments[5]; // 23
     188    result +=
     189        2. * params[spring_constant] *
     190        ( function_theta(r_ij.distance, r_ik.distance, r_il.distance, r_jk.distance, r_jl.distance, r_kl.distance)
     191            - params[equilibrium_distance]);
     192  }
     193  return derivative_components_t(1, result);
    183194}
    184195
    185196FourBodyPotential_Torsion::results_t
    186197FourBodyPotential_Torsion::parameter_derivative(
    187     const arguments_t &arguments,
     198    const list_of_arguments_t &listarguments,
    188199    const size_t index
    189200    ) const
    190201{
    191   ASSERT( arguments.size() == getSpecificArgumentCount(),
    192       "FourBodyPotential_Torsion::parameter_derivative() - requires exactly three arguments.");
    193   ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
    194       arguments, getParticleTypes()),
    195       "FourBodyPotential_Torsion::operator() - types don't match with ones in arguments.");
    196   const argument_t &r_ij = arguments[0]; // 01
    197   const argument_t &r_ik = arguments[1]; // 02
    198   const argument_t &r_il = arguments[2]; // 03
    199   const argument_t &r_jk = arguments[3]; // 12
    200   const argument_t &r_jl = arguments[4]; // 13
    201   const argument_t &r_kl = arguments[5]; // 23
    202   switch (index) {
    203     case spring_constant:
    204     {
    205       const result_t result =
    206                  Helpers::pow( function_theta(r_ij.distance, r_ik.distance, r_il.distance, r_jk.distance, r_jl.distance, r_kl.distance) - params[equilibrium_distance], 2 );
    207       return std::vector<result_t>(1, result);
    208       break;
     202  result_t result = 0.;
     203  for(list_of_arguments_t::const_iterator iter = listarguments.begin();
     204      iter != listarguments.end(); ++iter) {
     205    const arguments_t &arguments = *iter;
     206    ASSERT( arguments.size() == getSpecificArgumentCount(),
     207        "FourBodyPotential_Torsion::parameter_derivative() - requires exactly three arguments.");
     208    ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
     209        arguments, getParticleTypes()),
     210        "FourBodyPotential_Torsion::operator() - types don't match with ones in arguments.");
     211    const argument_t &r_ij = arguments[0]; // 01
     212    const argument_t &r_ik = arguments[1]; // 02
     213    const argument_t &r_il = arguments[2]; // 03
     214    const argument_t &r_jk = arguments[3]; // 12
     215    const argument_t &r_jl = arguments[4]; // 13
     216    const argument_t &r_kl = arguments[5]; // 23
     217    switch (index) {
     218      case spring_constant:
     219      {
     220        result +=
     221            Helpers::pow( function_theta(r_ij.distance, r_ik.distance, r_il.distance, r_jk.distance, r_jl.distance, r_kl.distance) - params[equilibrium_distance], 2 );
     222        break;
     223      }
     224      case equilibrium_distance:
     225      {
     226        result +=
     227            -2. * params[spring_constant]
     228                   * ( function_theta(r_ij.distance, r_ik.distance, r_il.distance, r_jk.distance, r_jl.distance, r_kl.distance) - params[equilibrium_distance]);
     229        break;
     230      }
     231      default:
     232        ASSERT(0, "FourBodyPotential_Torsion::parameter_derivative() - derivative to unknown parameter desired.");
     233        break;
    209234    }
    210     case equilibrium_distance:
    211     {
    212       const result_t result =
    213           -2. * params[spring_constant]
    214                  * ( function_theta(r_ij.distance, r_ik.distance, r_il.distance, r_jk.distance, r_jl.distance, r_kl.distance) - params[equilibrium_distance]);
    215       return std::vector<result_t>(1, result);
    216       break;
    217     }
    218     default:
    219       ASSERT(0, "FourBodyPotential_Torsion::parameter_derivative() - derivative to unknown parameter desired.");
    220       break;
    221235  }
    222   return std::vector<result_t>(1);
    223 }
    224 
    225 FunctionModel::extractor_t
    226 FourBodyPotential_Torsion::getSpecificExtractor() const
    227 {
    228   Fragment::charges_t charges;
    229   charges.resize(getParticleTypes().size());
    230   std::transform(getParticleTypes().begin(), getParticleTypes().end(),
    231       charges.begin(), boost::lambda::_1);
    232   FunctionModel::extractor_t returnfunction =
    233       boost::bind(&Extractors::gatherDistancesFromFragment,
    234           boost::bind(&Fragment::getPositions, _1),
    235           boost::bind(&Fragment::getCharges, _1),
    236           charges,
    237           _2);
    238   return returnfunction;
     236  return results_t(1, result);
    239237}
    240238

  • src/Potentials/Specifics/FourBodyPotential_Torsion.hpp

    rc73e35 rd8821e  
    3737  friend class PotentialFactory;
    3838  // some repeated typedefs to avoid ambiguities
     39  typedef FunctionModel::list_of_arguments_t list_of_arguments_t;
    3940  typedef FunctionModel::arguments_t arguments_t;
    4041  typedef FunctionModel::result_t result_t;
    … …  
    9293  /** Evaluates the harmonic potential function for the given arguments.
    9394   *
    94    * @param arguments single distance
     95   * @param listarguments list of six distances
    9596   * @return value of the potential function
    9697   */
    97   results_t operator()(const arguments_t &arguments) const;
     98  results_t operator()(const list_of_arguments_t &listarguments) const;
    9899
    99100  /** Evaluates the derivative of the potential function.
    100101   *
    101    * @param arguments single distance
     102   * @param listarguments list of six distances
    102103   * @return vector with derivative with respect to the input degrees of freedom
    103104   */
    104   derivative_components_t derivative(const arguments_t &arguments) const;
     105  derivative_components_t derivative(const list_of_arguments_t &listarguments) const;
    105106
    106107  /** Evaluates the derivative of the function with the given \a arguments
    107108   * with respect to a specific parameter indicated by \a index.
    108109   *
    109    * \param arguments set of arguments as input variables to the function
     110   * \param listarguments list of six distances
    110111   * \param index derivative of which parameter
    111112   * \return result vector containing the derivative with respect to the given
    112113   *         input
    113114   */
    114   results_t parameter_derivative(const arguments_t &arguments, const size_t index) const;
     115  results_t parameter_derivative(const list_of_arguments_t &listarguments, const size_t index) const;
    115116
    116117  /** Returns the functor that converts argument_s into the
    … …  
    169170    return upperbounds;
    170171  }
    171 
    172   /** Returns a bound function to be used with TrainingData, extracting distances
    173    * from a Fragment.
    174    *
    175    * \return bound function extracting distances from a fragment
    176    */
    177   FunctionModel::extractor_t getSpecificExtractor() const;
    178172
    179173  /** Returns a bound function to be used with TrainingData, extracting distances

  • src/Potentials/Specifics/ManyBodyPotential_Tersoff.cpp

    rc73e35 rd8821e  
    184184ManyBodyPotential_Tersoff::results_t
    185185ManyBodyPotential_Tersoff::operator()(
    186     const arguments_t &arguments
     186    const list_of_arguments_t &listarguments
    187187    ) const
    188188{
    189189//  Info info(__func__);
    190   double result = 0.;
    191   for(arguments_t::const_iterator argiter = arguments.begin();
    192       argiter != arguments.end();
    193       ++argiter) {
    194     const argument_t &r_ij = *argiter;
    195     ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
    196         arguments_t(1, r_ij), getParticleTypes()),
    197         "ManyBodyPotential_Tersoff::operator() - types don't match with ones in arguments.");
    198  
    199     const double cutoff = function_cutoff(r_ij.distance);
    200     const double temp = (cutoff == 0.) ?
    201         0. :
    202         cutoff * (
    203             function_prefactor(
    204                 alpha,
    205                 function_eta(r_ij))
     190  result_t result = 0.;
     191  for(list_of_arguments_t::const_iterator iter = listarguments.begin();
     192      iter != listarguments.end(); ++iter) {
     193    const arguments_t &arguments = *iter;
     194    for(arguments_t::const_iterator argiter = arguments.begin();
     195        argiter != arguments.end();
     196        ++argiter) {
     197      const argument_t &r_ij = *argiter;
     198      ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
     199          arguments_t(1, r_ij), getParticleTypes()),
     200          "ManyBodyPotential_Tersoff::operator() - types don't match with ones in arguments.");
     201
     202      const double cutoff = function_cutoff(r_ij.distance);
     203      const double temp = (cutoff == 0.) ?
     204          0. :
     205          cutoff * (
     206              function_prefactor(
     207                  alpha,
     208                  function_eta(r_ij))
     209              * function_smoother(
     210                  params[A],
     211                  params[lambda],
     212                  r_ij.distance)
     213              +
     214              function_prefactor(
     215                  params[beta],
     216                  function_zeta(r_ij))
     217              * function_smoother(
     218                  -params[B],
     219                  params[mu],
     220                  r_ij.distance)
     221          );
     222      result += temp;
     223    }
     224  }
     225//  LOG(2, "DEBUG: operator()(" << r_ij.distance << ") = " << result);
     226  return results_t(1, result);
     227}
     228
     229ManyBodyPotential_Tersoff::derivative_components_t
     230ManyBodyPotential_Tersoff::derivative(
     231    const list_of_arguments_t &listarguments
     232    ) const
     233{
     234//  Info info(__func__);
     235  return derivative_components_t();
     236}
     237
     238ManyBodyPotential_Tersoff::results_t
     239ManyBodyPotential_Tersoff::parameter_derivative(
     240    const list_of_arguments_t &listarguments,
     241    const size_t index
     242    ) const
     243{
     244//  Info info(__func__);
     245//  ASSERT( arguments.size() == 1,
     246//      "ManyBodyPotential_Tersoff::parameter_derivative() - requires exactly one argument.");
     247
     248  result_t result = 0.;
     249  for(list_of_arguments_t::const_iterator iter = listarguments.begin();
     250      iter != listarguments.end(); ++iter) {
     251    const arguments_t &arguments = *iter;
     252    for(arguments_t::const_iterator argiter = arguments.begin();
     253        argiter != arguments.end();
     254        ++argiter) {
     255      const argument_t &r_ij = *argiter;
     256      ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
     257          arguments_t(1, r_ij), getParticleTypes()),
     258          "ManyBodyPotential_Tersoff::operator() - types don't match with ones in arguments.");
     259
     260    switch (index) {
     261//      case R:
     262//      {
     263//        result += 0.;
     264//        break;
     265//      }
     266//      case S:
     267//      {
     268//        result += 0.;
     269//        break;
     270//      }
     271      case A:
     272      {
     273        const double cutoff = function_cutoff(r_ij.distance);
     274        result += (cutoff == 0.) ?
     275            0. :
     276            cutoff *
     277                function_prefactor(
     278                    alpha,
     279                    function_eta(r_ij))
     280                * function_smoother(
     281                    1.,
     282                    params[lambda],
     283                    r_ij.distance);
     284//            cutoff * function_prefactor(
     285//                alpha,
     286//                function_eta(r_ij))
     287//            * function_smoother(
     288//                1.,
     289//                params[mu],
     290//                r_ij.distance);
     291        break;
     292      }
     293      case B:
     294      {
     295        const double cutoff = function_cutoff(r_ij.distance);
     296        result += (cutoff == 0.) ?
     297            0. :
     298            cutoff * function_prefactor(
     299                params[beta],
     300                function_zeta(r_ij))
    206301            * function_smoother(
    207                 params[A],
    208                 params[lambda],
    209                 r_ij.distance)
    210             +
     302                -1.,
     303                params[mu],
     304                r_ij.distance);
     305//            cutoff * function_prefactor(
     306//                beta,
     307//                function_zeta(r_ij))
     308//            * function_smoother(
     309//                -params[B],
     310//                params[mu],
     311//                r_ij.distance)/params[B];
     312        break;
     313      }
     314      case lambda:
     315      {
     316        const double cutoff = function_cutoff(r_ij.distance);
     317        result += (cutoff == 0.) ?
     318            0. :
     319            -r_ij.distance * cutoff *
     320                function_prefactor(
     321                    alpha,
     322                    function_eta(r_ij))
     323                * function_smoother(
     324                    params[A],
     325                    params[lambda],
     326                    r_ij.distance);
     327        break;
     328      }
     329      case mu:
     330      {
     331        const double cutoff = function_cutoff(r_ij.distance);
     332        result += (cutoff == 0.) ?
     333            0. :
     334            -r_ij.distance * cutoff *(
    211335            function_prefactor(
    212336                params[beta],
    … …  
    217341                r_ij.distance)
    218342        );
    219     result += temp;
    220   }
    221 //  LOG(2, "DEBUG: operator()(" << r_ij.distance << ") = " << result);
    222   return std::vector<result_t>(1, result);
    223 }
    224 
    225 ManyBodyPotential_Tersoff::derivative_components_t
    226 ManyBodyPotential_Tersoff::derivative(
    227     const arguments_t &arguments
    228     ) const
    229 {
    230 //  Info info(__func__);
    231   return ManyBodyPotential_Tersoff::derivative_components_t();
    232 }
    233 
    234 ManyBodyPotential_Tersoff::results_t
    235 ManyBodyPotential_Tersoff::parameter_derivative(
    236     const arguments_t &arguments,
    237     const size_t index
    238     ) const
    239 {
    240 //  Info info(__func__);
    241 //  ASSERT( arguments.size() == 1,
    242 //      "ManyBodyPotential_Tersoff::parameter_derivative() - requires exactly one argument.");
    243 
    244   double result = 0.;
    245   for(arguments_t::const_iterator argiter = arguments.begin();
    246       argiter != arguments.end();
    247       ++argiter) {
    248     const argument_t &r_ij = *argiter;
    249     ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
    250         arguments_t(1, r_ij), getParticleTypes()),
    251         "ManyBodyPotential_Tersoff::operator() - types don't match with ones in arguments.");
    252 
    253   switch (index) {
    254 //    case R:
    255 //    {
    256 //      result += 0.;
    257 //      break;
    258 //    }
    259 //    case S:
    260 //    {
    261 //      result += 0.;
    262 //      break;
    263 //    }
    264     case A:
    265     {
    266       const double cutoff = function_cutoff(r_ij.distance);
    267       result += (cutoff == 0.) ?
    268           0. :
    269           cutoff *
    270               function_prefactor(
    271                   alpha,
    272                   function_eta(r_ij))
    273               * function_smoother(
    274                   1.,
    275                   params[lambda],
    276                   r_ij.distance);
    277 //          cutoff * function_prefactor(
    278 //              alpha,
    279 //              function_eta(r_ij))
    280 //          * function_smoother(
    281 //              1.,
    282 //              params[mu],
    283 //              r_ij.distance);
    284       break;
     343        break;
     344      }
     345//      case lambda3:
     346//      {
     347//        result += 0.;
     348//        break;
     349//      }
     350//      case alpha:
     351//      {
     352//        const double temp =
     353//            pow(alpha*function_eta(r_ij), params[n]);
     354//        const double cutoff = function_cutoff(r_ij.distance);
     355//        result += (cutoff == 0.) || (alpha == 0. )?
     356//            0. :
     357//            function_smoother(
     358//                params[A],
     359//                params[lambda],
     360//                r_ij.distance)
     361//            * (-.5) * alpha * (temp/alpha)
     362//            / (1. + temp)
     363//            ;
     364//        break;
     365//      }
     366//      case chi:
     367//      {
     368//        result += 0.;
     369//        break;
     370//      }
     371//      case omega:
     372//      {
     373//        result += 0.;
     374//        break;
     375//      }
     376      case beta:
     377      {
     378        const double temp =
     379            pow(params[beta]*function_zeta(r_ij), params[n]);
     380        const double cutoff = function_cutoff(r_ij.distance);
     381        result += (cutoff == 0.) || (params[beta] == 0. )?
     382            0. : cutoff *
     383            function_smoother(
     384                -params[B],
     385                params[mu],
     386                r_ij.distance)
     387            * (-.5)
     388            * function_prefactor(
     389                params[beta],
     390                function_zeta(r_ij))
     391            * (temp/params[beta])
     392            / (1. + temp)
     393            ;
     394        break;
     395      }
     396      case n:
     397      {
     398        const double zeta = function_zeta(r_ij);
     399        const double temp = pow( params[beta]*zeta , params[n]);
     400        const double cutoff = function_cutoff(r_ij.distance);
     401        const double tempres = ((cutoff == 0.) || (zeta == 0.)) ? // zeta must be caught if zero due to log
     402            0. : .5 * cutoff *
     403            function_smoother(
     404                -params[B],
     405                params[mu],
     406                r_ij.distance)
     407            * function_prefactor(
     408                params[beta],
     409                function_zeta(r_ij))
     410            * ( log(1.+temp)/(params[n]*params[n]) - temp
     411                * (log(function_zeta(r_ij)) + log(params[beta]))
     412                /(params[n]*(1.+temp)))
     413            ;
     414//        if (tempres != tempres)
     415//    LOG(2, "DEBUG: tempres is NaN.");
     416//        LOG(2, "DEBUG: Adding " << tempres << " for p.d. w.r.t n, temp=" << temp << ", cutoff=" << cutoff);
     417        result += tempres;
     418        break;
     419      }
     420      case c:
     421      {
     422        const double zeta = function_zeta(r_ij);
     423        if (zeta == 0.)
     424          break;
     425        const double temp =
     426            pow(zeta, params[n]-1.) * pow(params[beta],params[n]);
     427        const double cutoff = function_cutoff(r_ij.distance);
     428        const double tempres = (cutoff == 0.) ?
     429            0. : cutoff *
     430            function_smoother(
     431                -params[B],
     432                params[mu],
     433                r_ij.distance)
     434            * function_prefactor(
     435                params[beta],
     436                zeta)
     437             * (-1.) * temp / (1.+temp*zeta);
     438        double factor = function_derivative_c(r_ij);
     439        result += tempres*factor;
     440        if (result != result)
     441          ELOG(1, "result is NaN, zeta=" << zeta << ", temp=" << temp << ", cutoff=" << cutoff << ", tempres=" << tempres << ", factor=" << factor);
     442        break;
     443      }
     444      case d:
     445      {
     446        const double zeta = function_zeta(r_ij);
     447        const double temp =
     448            pow(zeta, params[n]-1.) * pow(params[beta],params[n]);
     449        const double cutoff = function_cutoff(r_ij.distance);
     450        const double tempres = (cutoff == 0.) ?
     451            0. : cutoff *
     452            function_smoother(
     453                -params[B],
     454                params[mu],
     455                r_ij.distance)
     456            * function_prefactor(
     457                params[beta],
     458                zeta)
     459             * (-1.) * temp / (1.+temp*zeta);
     460        double factor = function_derivative_d(r_ij);
     461        result += tempres*factor;
     462        if (result != result)
     463          ELOG(1, "result is NaN, zeta=" << zeta << ", temp=" << temp << ", cutoff=" << cutoff << ", tempres=" << tempres << ", factor=" << factor);
     464        break;
     465      }
     466      case h:
     467      {
     468        const double zeta = function_zeta(r_ij);
     469        const double temp =
     470            pow(zeta, params[n]-1.) * pow(params[beta],params[n]);
     471        const double cutoff = function_cutoff(r_ij.distance);
     472        const double tempres = (cutoff == 0.) ?
     473            0. : cutoff *
     474            function_smoother(
     475                -params[B],
     476                params[mu],
     477                r_ij.distance)
     478            * function_prefactor(
     479                params[beta],
     480                zeta)
     481             * (-1.) * temp / (1.+temp*zeta);
     482        double factor = function_derivative_h(r_ij);
     483        result += tempres*factor;
     484        if (result != result)
     485          ELOG(1, "result is NaN, zeta=" << zeta << ", temp=" << temp << ", cutoff=" << cutoff << ", tempres=" << tempres << ", factor=" << factor);
     486        break;
     487      }
     488      default:
     489        ASSERT(0, "ManyBodyPotential_Tersoff::parameter_derivative() - derivative to unknown parameter desired.");
     490        break;
    285491    }
    286     case B:
    287     {
    288       const double cutoff = function_cutoff(r_ij.distance);
    289       result += (cutoff == 0.) ?
    290           0. :
    291           cutoff * function_prefactor(
    292               params[beta],
    293               function_zeta(r_ij))
    294           * function_smoother(
    295               -1.,
    296               params[mu],
    297               r_ij.distance);
    298 //          cutoff * function_prefactor(
    299 //              beta,
    300 //              function_zeta(r_ij))
    301 //          * function_smoother(
    302 //              -params[B],
    303 //              params[mu],
    304 //              r_ij.distance)/params[B];
    305       break;
     492    if (result != result)
     493      ELOG(1, "result is NaN.");
    306494    }
    307     case lambda:
    308     {
    309       const double cutoff = function_cutoff(r_ij.distance);
    310       result += (cutoff == 0.) ?
    311           0. :
    312           -r_ij.distance * cutoff *
    313               function_prefactor(
    314                   alpha,
    315                   function_eta(r_ij))
    316               * function_smoother(
    317                   params[A],
    318                   params[lambda],
    319                   r_ij.distance);
    320       break;
    321     }
    322     case mu:
    323     {
    324       const double cutoff = function_cutoff(r_ij.distance);
    325       result += (cutoff == 0.) ?
    326           0. :
    327           -r_ij.distance * cutoff *(
    328           function_prefactor(
    329               params[beta],
    330               function_zeta(r_ij))
    331           * function_smoother(
    332               -params[B],
    333               params[mu],
    334               r_ij.distance)
    335       );
    336       break;
    337     }
    338 //    case lambda3:
    339 //    {
    340 //      result += 0.;
    341 //      break;
    342 //    }
    343 //    case alpha:
    344 //    {
    345 //      const double temp =
    346 //          pow(alpha*function_eta(r_ij), params[n]);
    347 //      const double cutoff = function_cutoff(r_ij.distance);
    348 //      result += (cutoff == 0.) || (alpha == 0. )?
    349 //          0. :
    350 //          function_smoother(
    351 //              params[A],
    352 //              params[lambda],
    353 //              r_ij.distance)
    354 //          * (-.5) * alpha * (temp/alpha)
    355 //          / (1. + temp)
    356 //          ;
    357 //      break;
    358 //    }
    359 //    case chi:
    360 //    {
    361 //      result += 0.;
    362 //      break;
    363 //    }
    364 //    case omega:
    365 //    {
    366 //      result += 0.;
    367 //      break;
    368 //    }
    369     case beta:
    370     {
    371       const double temp =
    372           pow(params[beta]*function_zeta(r_ij), params[n]);
    373       const double cutoff = function_cutoff(r_ij.distance);
    374       result += (cutoff == 0.) || (params[beta] == 0. )?
    375           0. : cutoff *
    376           function_smoother(
    377               -params[B],
    378               params[mu],
    379               r_ij.distance)
    380           * (-.5)
    381           * function_prefactor(
    382               params[beta],
    383               function_zeta(r_ij))
    384           * (temp/params[beta])
    385           / (1. + temp)
    386           ;
    387       break;
    388     }
    389     case n:
    390     {
    391       const double zeta = function_zeta(r_ij);
    392       const double temp = pow( params[beta]*zeta , params[n]);
    393       const double cutoff = function_cutoff(r_ij.distance);
    394       const double tempres = ((cutoff == 0.) || (zeta == 0.)) ? // zeta must be caught if zero due to log
    395           0. : .5 * cutoff *
    396           function_smoother(
    397               -params[B],
    398               params[mu],
    399               r_ij.distance)
    400           * function_prefactor(
    401               params[beta],
    402               function_zeta(r_ij))
    403           * ( log(1.+temp)/(params[n]*params[n]) - temp
    404               * (log(function_zeta(r_ij)) + log(params[beta]))
    405               /(params[n]*(1.+temp)))
    406           ;
    407 //      if (tempres != tempres)
    408 //      LOG(2, "DEBUG: tempres is NaN.");
    409 //      LOG(2, "DEBUG: Adding " << tempres << " for p.d. w.r.t n, temp=" << temp << ", cutoff=" << cutoff);
    410       result += tempres;
    411       break;
    412     }
    413     case c:
    414     {
    415       const double zeta = function_zeta(r_ij);
    416       if (zeta == 0.)
    417         break;
    418       const double temp =
    419           pow(zeta, params[n]-1.) * pow(params[beta],params[n]);
    420       const double cutoff = function_cutoff(r_ij.distance);
    421       const double tempres = (cutoff == 0.) ?
    422           0. : cutoff *
    423           function_smoother(
    424               -params[B],
    425               params[mu],
    426               r_ij.distance)
    427           * function_prefactor(
    428               params[beta],
    429               zeta)
    430            * (-1.) * temp / (1.+temp*zeta);
    431       double factor = function_derivative_c(r_ij);
    432       result += tempres*factor;
    433       if (result != result)
    434         ELOG(1, "result is NaN, zeta=" << zeta << ", temp=" << temp << ", cutoff=" << cutoff << ", tempres=" << tempres << ", factor=" << factor);
    435       break;
    436     }
    437     case d:
    438     {
    439       const double zeta = function_zeta(r_ij);
    440       const double temp =
    441           pow(zeta, params[n]-1.) * pow(params[beta],params[n]);
    442       const double cutoff = function_cutoff(r_ij.distance);
    443       const double tempres = (cutoff == 0.) ?
    444           0. : cutoff *
    445           function_smoother(
    446               -params[B],
    447               params[mu],
    448               r_ij.distance)
    449           * function_prefactor(
    450               params[beta],
    451               zeta)
    452            * (-1.) * temp / (1.+temp*zeta);
    453       double factor = function_derivative_d(r_ij);
    454       result += tempres*factor;
    455       if (result != result)
    456         ELOG(1, "result is NaN, zeta=" << zeta << ", temp=" << temp << ", cutoff=" << cutoff << ", tempres=" << tempres << ", factor=" << factor);
    457       break;
    458     }
    459     case h:
    460     {
    461       const double zeta = function_zeta(r_ij);
    462       const double temp =
    463           pow(zeta, params[n]-1.) * pow(params[beta],params[n]);
    464       const double cutoff = function_cutoff(r_ij.distance);
    465       const double tempres = (cutoff == 0.) ?
    466           0. : cutoff *
    467           function_smoother(
    468               -params[B],
    469               params[mu],
    470               r_ij.distance)
    471           * function_prefactor(
    472               params[beta],
    473               zeta)
    474            * (-1.) * temp / (1.+temp*zeta);
    475       double factor = function_derivative_h(r_ij);
    476       result += tempres*factor;
    477       if (result != result)
    478         ELOG(1, "result is NaN, zeta=" << zeta << ", temp=" << temp << ", cutoff=" << cutoff << ", tempres=" << tempres << ", factor=" << factor);
    479       break;
    480     }
    481     default:
    482       ASSERT(0, "ManyBodyPotential_Tersoff::parameter_derivative() - derivative to unknown parameter desired.");
    483       break;
    484   }
    485   if (result != result)
    486     ELOG(1, "result is NaN.");
    487495  }
    488496  return results_t(1,-result);
    … …  
    689697//  LOG(2, "DEBUG: function_angle(" << r_ij << "," << r_ik << "," << r_jk << ") = " << result);
    690698  return result;
    691 }
    692 
    693 FunctionModel::extractor_t
    694 ManyBodyPotential_Tersoff::getSpecificExtractor() const
    695 {
    696   FunctionModel::extractor_t returnfunction =
    697       boost::bind(&Extractors::gatherAllDistances,
    698           boost::bind(&Fragment::getPositions, _1),
    699           boost::bind(&Fragment::getCharges, _1),
    700           _2);
    701   return returnfunction;
    702699}
    703700

  • src/Potentials/Specifics/ManyBodyPotential_Tersoff.hpp

    rc73e35 rd8821e  
    3838  friend class PotentialFactory;
    3939  // some repeated typedefs to avoid ambiguities
     40  typedef FunctionModel::list_of_arguments_t list_of_arguments_t;
    4041  typedef FunctionModel::arguments_t arguments_t;
    4142  typedef FunctionModel::result_t result_t;
    … …  
    106107  /** Evaluates the Tersoff potential for the given arguments.
    107108   *
    108    * @param arguments single distance
     109   * @param listarguments list of distances
    109110   * @return value of the potential function
    110111   */
    111   results_t operator()(const arguments_t &arguments) const;
     112  results_t operator()(const list_of_arguments_t &listarguments) const;
    112113
    113114  /** Evaluates the derivative of the Tersoff potential with respect to the
    114115   * input variables.
    115116   *
    116    * @param arguments single distance
     117   * @param listarguments list of distances
    117118   * @return vector with components of the derivative
    118119   */
    119   derivative_components_t derivative(const arguments_t &arguments) const;
     120  derivative_components_t derivative(const list_of_arguments_t &listarguments) const;
    120121
    121122  /** Evaluates the derivative of the function with the given \a arguments
    122123   * with respect to a specific parameter indicated by \a index.
    123124   *
    124    * \param arguments set of arguments as input variables to the function
     125   * \param listarguments list of distances
    125126   * \param index derivative of which parameter
    126127   * \return result vector containing the derivative with respect to the given
    127128   *         input
    128129   */
    129   results_t parameter_derivative(const arguments_t &arguments, const size_t index) const;
     130  results_t parameter_derivative(const list_of_arguments_t &listarguments, const size_t index) const;
    130131
    131132  /** Returns the functor that converts argument_s into the
    … …  
    188189    return parameters_t(getParameterDimension(), std::numeric_limits<double>::max());
    189190  }
    190 
    191   /** Returns a bound function to be used with TrainingData, extracting distances
    192    * from a Fragment.
    193    *
    194    * \return bound function extracting distances from a fragment
    195    */
    196   FunctionModel::extractor_t getSpecificExtractor() const;
    197191
    198192  /** Returns a bound function to be used with TrainingData, extracting distances

  • src/Potentials/Specifics/PairPotential_Harmonic.cpp

    rc73e35 rd8821e  
    117117PairPotential_Harmonic::results_t
    118118PairPotential_Harmonic::operator()(
    119     const arguments_t &arguments
     119    const list_of_arguments_t &listarguments
    120120    ) const
    121121{
    122   ASSERT( arguments.size() == 1,
    123       "PairPotential_Harmonic::operator() - requires exactly one argument.");
    124   ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
    125       arguments, getParticleTypes()),
    126       "PairPotential_Harmonic::operator() - types don't match with ones in arguments.");
    127   const argument_t &r_ij = arguments[0];
    128   const result_t result =
    129       params[spring_constant]
    130              * Helpers::pow( r_ij.distance - params[equilibrium_distance], 2 );
    131   return std::vector<result_t>(1, result);
     122  result_t result = 0.;
     123  for(list_of_arguments_t::const_iterator iter = listarguments.begin();
     124      iter != listarguments.end(); ++iter) {
     125    const arguments_t &arguments = *iter;
     126    ASSERT( arguments.size() == 1,
     127        "PairPotential_Harmonic::operator() - requires exactly one argument.");
     128    ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
     129        arguments, getParticleTypes()),
     130        "PairPotential_Harmonic::operator() - types don't match with ones in arguments.");
     131    const argument_t &r_ij = arguments[0];
     132    result +=
     133        params[spring_constant]
     134               * Helpers::pow( r_ij.distance - params[equilibrium_distance], 2 );
     135  }
     136  return results_t(1, result);
    132137}
    133138
    134139PairPotential_Harmonic::derivative_components_t
    135140PairPotential_Harmonic::derivative(
    136     const arguments_t &arguments
     141    const list_of_arguments_t &listarguments
    137142    ) const
    138143{
    139   ASSERT( arguments.size() == 1,
    140       "PairPotential_Harmonic::operator() - requires exactly one argument.");
    141   ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
    142       arguments, getParticleTypes()),
    143       "PairPotential_Harmonic::operator() - types don't match with ones in arguments.");
    144   derivative_components_t result;
    145   const argument_t &r_ij = arguments[0];
    146   result.push_back( 2. * params[spring_constant] * ( r_ij.distance - params[equilibrium_distance]) );
    147   ASSERT( result.size() == 1,
    148       "PairPotential_Harmonic::operator() - we did not create exactly one component.");
    149   return result;
     144  result_t result = 0.;
     145  for(list_of_arguments_t::const_iterator iter = listarguments.begin();
     146      iter != listarguments.end(); ++iter) {
     147    const arguments_t &arguments = *iter;
     148      ASSERT( arguments.size() == 1,
     149          "PairPotential_Harmonic::operator() - requires exactly one argument.");
     150      ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
     151          arguments, getParticleTypes()),
     152          "PairPotential_Harmonic::operator() - types don't match with ones in arguments.");
     153      const argument_t &r_ij = arguments[0];
     154      result +=
     155          2. * params[spring_constant] *
     156          ( r_ij.distance - params[equilibrium_distance]);
     157  }
     158  return derivative_components_t(1, result);
    150159}
    151160
    152161PairPotential_Harmonic::results_t
    153162PairPotential_Harmonic::parameter_derivative(
    154     const arguments_t &arguments,
     163    const list_of_arguments_t &listarguments,
    155164    const size_t index
    156165    ) const
    157166{
    158   ASSERT( arguments.size() == 1,
    159       "PairPotential_Harmonic::parameter_derivative() - requires exactly one argument.");
    160   ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
    161       arguments, getParticleTypes()),
    162       "PairPotential_Harmonic::operator() - types don't match with ones in arguments.");
    163   const argument_t &r_ij = arguments[0];
    164   switch (index) {
    165     case spring_constant:
    166     {
    167       const result_t result =
    168                  Helpers::pow( r_ij.distance - params[equilibrium_distance], 2 );
    169       return std::vector<result_t>(1, result);
    170       break;
     167  result_t result = 0.;
     168  for(list_of_arguments_t::const_iterator iter = listarguments.begin();
     169      iter != listarguments.end(); ++iter) {
     170    const arguments_t &arguments = *iter;
     171    ASSERT( arguments.size() == 1,
     172        "PairPotential_Harmonic::parameter_derivative() - requires exactly one argument.");
     173    ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
     174        arguments, getParticleTypes()),
     175        "PairPotential_Harmonic::operator() - types don't match with ones in arguments.");
     176    const argument_t &r_ij = arguments[0];
     177    switch (index) {
     178      case spring_constant:
     179      {
     180        result +=
     181            Helpers::pow( r_ij.distance - params[equilibrium_distance], 2 );
     182        break;
     183      }
     184      case equilibrium_distance:
     185      {
     186        result +=
     187            -2. * params[spring_constant]
     188                   * ( r_ij.distance - params[equilibrium_distance]);
     189        break;
     190      }
     191      default:
     192        ASSERT(0, "PairPotential_Harmonic::parameter_derivative() - derivative to unknown parameter desired.");
     193        break;
    171194    }
    172     case equilibrium_distance:
    173     {
    174       const result_t result =
    175           -2. * params[spring_constant]
    176                  * ( r_ij.distance - params[equilibrium_distance]);
    177       return std::vector<result_t>(1, result);
    178       break;
    179     }
    180     default:
    181       ASSERT(0, "PairPotential_Harmonic::parameter_derivative() - derivative to unknown parameter desired.");
    182       break;
    183195  }
    184 
    185   return PairPotential_Harmonic::results_t(1, 0.);
    186 }
    187 
    188 FunctionModel::extractor_t
    189 PairPotential_Harmonic::getSpecificExtractor() const
    190 {
    191   Fragment::charges_t charges;
    192   charges.resize(getParticleTypes().size());
    193   std::transform(getParticleTypes().begin(), getParticleTypes().end(),
    194       charges.begin(), boost::lambda::_1);
    195   FunctionModel::extractor_t returnfunction =
    196       boost::bind(&Extractors::gatherDistancesFromFragment,
    197           boost::bind(&Fragment::getPositions, _1),
    198           boost::bind(&Fragment::getCharges, _1),
    199           charges,
    200           _2);
    201   return returnfunction;
     196  return results_t(1, result);
    202197}
    203198

  • src/Potentials/Specifics/PairPotential_Harmonic.hpp

    rc73e35 rd8821e  
    3535  friend class PotentialFactory;
    3636  // some repeated typedefs to avoid ambiguities
     37  typedef FunctionModel::list_of_arguments_t list_of_arguments_t;
    3738  typedef FunctionModel::arguments_t arguments_t;
    3839  typedef FunctionModel::result_t result_t;
    … …  
    9091  /** Evaluates the harmonic potential function for the given arguments.
    9192   *
    92    * @param arguments single distance
     93   * @param listarguments list of single distances
    9394   * @return value of the potential function
    9495   */
    95   results_t operator()(const arguments_t &arguments) const;
     96  results_t operator()(const list_of_arguments_t &listarguments) const;
    9697
    9798  /** Evaluates the derivative of the potential function.
    9899   *
    99    * @param arguments single distance
     100   * @param listarguments list of single distances
    100101   * @return vector with derivative with respect to the input degrees of freedom
    101102   */
    102   derivative_components_t derivative(const arguments_t &arguments) const;
     103  derivative_components_t derivative(const list_of_arguments_t &listarguments) const;
    103104
    104105  /** Evaluates the derivative of the function with the given \a arguments
    105106   * with respect to a specific parameter indicated by \a index.
    106107   *
    107    * \param arguments set of arguments as input variables to the function
     108   * \param listarguments list of single distances
    108109   * \param index derivative of which parameter
    109110   * \return result vector containing the derivative with respect to the given
    110111   *         input
    111112   */
    112   results_t parameter_derivative(const arguments_t &arguments, const size_t index) const;
     113  results_t parameter_derivative(const list_of_arguments_t &listarguments, const size_t index) const;
    113114
    114115  /** Returns the functor that converts argument_s into the
    … …  
    148149    return parameters_t(getParameterDimension(), std::numeric_limits<double>::max());
    149150  }
    150 
    151   /** Returns a bound function to be used with TrainingData, extracting distances
    152    * from a Fragment.
    153    *
    154    * \return bound function extracting distances from a fragment
    155    */
    156   FunctionModel::extractor_t getSpecificExtractor() const;
    157151
    158152  /** Returns a bound function to be used with TrainingData, extracting distances

  • src/Potentials/Specifics/PairPotential_LennardJones.cpp

    rc73e35 rd8821e  
    123123PairPotential_LennardJones::results_t
    124124PairPotential_LennardJones::operator()(
    125     const arguments_t &arguments
     125    const list_of_arguments_t &listarguments
    126126    ) const
    127127{
    128   ASSERT( arguments.size() == 1,
    129       "PairPotential_LennardJones::operator() - requires one argument.");
    130   ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
    131       arguments, getParticleTypes()),
    132       "PairPotential_LennardJones::operator() - types don't match with ones in arguments.");
    133   const double &r = arguments[0].distance;
    134   const double temp = Helpers::pow(params[sigma]/r, 6);
    135   const result_t result = 4.*params[epsilon] * (temp*temp - temp);
    136   return std::vector<result_t>(1, result);
     128  result_t result = 0.;
     129  for(list_of_arguments_t::const_iterator iter = listarguments.begin();
     130      iter != listarguments.end(); ++iter) {
     131    const arguments_t &arguments = *iter;
     132    ASSERT( arguments.size() == 1,
     133        "PairPotential_LennardJones::operator() - requires one argument.");
     134    ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
     135        arguments, getParticleTypes()),
     136        "PairPotential_LennardJones::operator() - types don't match with ones in arguments.");
     137    const double &r = arguments[0].distance;
     138    const double temp = Helpers::pow(params[sigma]/r, 6);
     139    result += 4.*params[epsilon] * (temp*temp - temp);
     140  }
     141  return results_t(1, result);
    137142}
    138143
    139144PairPotential_LennardJones::derivative_components_t
    140145PairPotential_LennardJones::derivative(
    141     const arguments_t &arguments
     146    const list_of_arguments_t &listarguments
    142147    ) const
    143148{
    144   ASSERT( arguments.size() == 1,
    145       "PairPotential_LennardJones::operator() - requires no argument.");
    146   ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
    147       arguments, getParticleTypes()),
    148       "PairPotential_LennardJones::operator() - types don't match with ones in arguments.");
    149   const double &r = arguments[0].distance;
    150149  const double sigma6 = Helpers::pow(params[sigma], 6);
    151   const result_t result =
    152       4.*params[epsilon] * (
    153           sigma6*sigma6*(-12.) / Helpers::pow(r,13)
    154           - sigma6*(-6.) /Helpers::pow(r,7)
    155       );
    156   derivative_components_t results(1, result);
    157   return results;
     150  result_t result = 0.;
     151  for(list_of_arguments_t::const_iterator iter = listarguments.begin();
     152      iter != listarguments.end(); ++iter) {
     153    const arguments_t &arguments = *iter;
     154    ASSERT( arguments.size() == 1,
     155        "PairPotential_LennardJones::operator() - requires no argument.");
     156    ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
     157        arguments, getParticleTypes()),
     158        "PairPotential_LennardJones::operator() - types don't match with ones in arguments.");
     159    const double &r = arguments[0].distance;
     160    result +=
     161        4.*params[epsilon] * (
     162            sigma6*sigma6*(-12.) / Helpers::pow(r,13)
     163            - sigma6*(-6.) /Helpers::pow(r,7)
     164        );
     165  }
     166  return derivative_components_t(1, result);
    158167}
    159168
    160169PairPotential_LennardJones::results_t
    161170PairPotential_LennardJones::parameter_derivative(
    162     const arguments_t &arguments,
     171    const list_of_arguments_t &listarguments,
    163172    const size_t index
    164173    ) const
    165174{
    166   ASSERT( arguments.size() == 1,
    167       "PairPotential_LennardJones::parameter_derivative() - requires no argument.");
    168   ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
    169       arguments, getParticleTypes()),
    170       "PairPotential_LennardJones::operator() - types don't match with ones in arguments.");
    171   const double &r = arguments[0].distance;
    172   switch (index) {
    173     case epsilon:
    174     {
    175       const double temp = Helpers::pow(params[sigma]/r, 6);
    176       const result_t result = 4. * (temp*temp - temp);
    177       return std::vector<result_t>(1, result);
    178       break;
     175  result_t result = 0.;
     176  for(list_of_arguments_t::const_iterator iter = listarguments.begin();
     177      iter != listarguments.end(); ++iter) {
     178    const arguments_t &arguments = *iter;
     179    ASSERT( arguments.size() == 1,
     180        "PairPotential_LennardJones::parameter_derivative() - requires no argument.");
     181    ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
     182        arguments, getParticleTypes()),
     183        "PairPotential_LennardJones::operator() - types don't match with ones in arguments.");
     184    const double &r = arguments[0].distance;
     185    switch (index) {
     186      case epsilon:
     187      {
     188        const double temp = Helpers::pow(params[sigma]/r, 6);
     189        result += 4. * (temp*temp - temp);
     190        break;
     191      }
     192      case sigma:
     193      {
     194        const double r6 = Helpers::pow(r, 6);
     195        result +=
     196            4.*params[epsilon] * (
     197                12. * Helpers::pow(params[sigma],11)/(r6*r6)
     198                - 6. * Helpers::pow(params[sigma],5)/r6
     199            );
     200        break;
     201      }
     202      default:
     203        break;
    179204    }
    180     case sigma:
    181     {
    182       const double r6 = Helpers::pow(r, 6);
    183       const result_t result =
    184           4.*params[epsilon] * (
    185               12. * Helpers::pow(params[sigma],11)/(r6*r6)
    186               - 6. * Helpers::pow(params[sigma],5)/r6
    187           );
    188       return std::vector<result_t>(1, result);
    189       break;
    190     }
    191     default:
    192       break;
    193205  }
    194   return std::vector<result_t>(1, 0.);
    195 }
    196 
    197 FunctionModel::extractor_t
    198 PairPotential_LennardJones::getSpecificExtractor() const
    199 {
    200   Fragment::charges_t charges;
    201   charges.resize(getParticleTypes().size());
    202   std::transform(getParticleTypes().begin(), getParticleTypes().end(),
    203       charges.begin(), boost::lambda::_1);
    204   FunctionModel::extractor_t returnfunction =
    205       boost::bind(&Extractors::gatherDistancesFromFragment,
    206           boost::bind(&Fragment::getPositions, _1),
    207           boost::bind(&Fragment::getCharges, _1),
    208           charges,
    209           _2);
    210   return returnfunction;
     206  return results_t(1, result);
    211207}
    212208

  • src/Potentials/Specifics/PairPotential_LennardJones.hpp

    rc73e35 rd8821e  
    3636  friend class PotentialFactory;
    3737  // some repeated typedefs to avoid ambiguities
     38  typedef FunctionModel::list_of_arguments_t list_of_arguments_t;
    3839  typedef FunctionModel::arguments_t arguments_t;
    3940  typedef FunctionModel::result_t result_t;
    … …  
    9293  /** Evaluates the harmonic potential function for the given arguments.
    9394   *
    94    * @param arguments single distance
     95   * @param listarguments list of single distances
    9596   * @return value of the potential function
    9697   */
    97   results_t operator()(const arguments_t &arguments) const;
     98  results_t operator()(const list_of_arguments_t &listarguments) const;
    9899
    99100  /** Evaluates the derivative of the potential function.
    100101   *
    101    * @param arguments single distance
     102   * @param listarguments list of single distances
    102103   * @return vector with derivative with respect to the input degrees of freedom
    103104   */
    104   derivative_components_t derivative(const arguments_t &arguments) const;
     105  derivative_components_t derivative(const list_of_arguments_t &listarguments) const;
    105106
    106107  /** Evaluates the derivative of the function with the given \a arguments
    107108   * with respect to a specific parameter indicated by \a index.
    108109   *
    109    * \param arguments set of arguments as input variables to the function
     110   * \param listarguments list of single distances
    110111   * \param index derivative of which parameter
    111112   * \return result vector containing the derivative with respect to the given
    112113   *         input
    113114   */
    114   results_t parameter_derivative(const arguments_t &arguments, const size_t index) const;
     115  results_t parameter_derivative(const list_of_arguments_t &listarguments, const size_t index) const;
    115116
    116117  /** Returns the functor that converts argument_s into the
    … …  
    165166    return parameters_t(getParameterDimension(), std::numeric_limits<double>::max());
    166167  }
    167 
    168   /** Returns a bound function to be used with TrainingData, extracting distances
    169    * from a Fragment.
    170    *
    171    * \return bound function extracting distances from a fragment
    172    */
    173   FunctionModel::extractor_t getSpecificExtractor() const;
    174168
    175169  /** Returns a bound function to be used with TrainingData, extracting distances

  • src/Potentials/Specifics/PairPotential_Morse.cpp

    rc73e35 rd8821e  
    124124PairPotential_Morse::results_t
    125125PairPotential_Morse::operator()(
    126     const arguments_t &arguments
     126    const list_of_arguments_t &listarguments
    127127    ) const
    128128{
    129   ASSERT( arguments.size() == 1,
    130       "PairPotential_Morse::operator() - requires exactly one argument.");
    131   ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
    132       arguments, getParticleTypes()),
    133       "PairPotential_Morse::operator() - types don't match with ones in arguments.");
    134   const argument_t &r_ij = arguments[0];
    135   // Maple: f(r,D,k,R,c) := D * (1 - exp(-k*(r-R)))^(2)+c
    136   const result_t result =
    137       params[dissociation_energy] * Helpers::pow( 1.
    138           - exp( - params[spring_constant] * ( r_ij.distance - params[equilibrium_distance])),2);
     129  result_t result = 0.;
     130  for(list_of_arguments_t::const_iterator iter = listarguments.begin();
     131      iter != listarguments.end(); ++iter) {
     132    const arguments_t &arguments = *iter;
     133    ASSERT( arguments.size() == 1,
     134        "PairPotential_Morse::operator() - requires exactly one argument.");
     135    ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
     136        arguments, getParticleTypes()),
     137        "PairPotential_Morse::operator() - types don't match with ones in arguments.");
     138    const argument_t &r_ij = arguments[0];
     139    // Maple: f(r,D,k,R,c) := D * (1 - exp(-k*(r-R)))^(2)+c
     140    result +=
     141        params[dissociation_energy] * Helpers::pow( 1.
     142            - exp( - params[spring_constant] * ( r_ij.distance - params[equilibrium_distance])),2);
     143  }
    139144  return std::vector<result_t>(1, result);
    140145}
    … …  
    142147PairPotential_Morse::derivative_components_t
    143148PairPotential_Morse::derivative(
    144     const arguments_t &arguments
     149    const list_of_arguments_t &listarguments
    145150    ) const
    146151{
    147   ASSERT( arguments.size() == 1,
    148       "PairPotential_Morse::operator() - requires exactly one argument.");
    149   ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
    150       arguments, getParticleTypes()),
    151       "PairPotential_Morse::operator() - types don't match with ones in arguments.");
    152   derivative_components_t result;
    153   const argument_t &r_ij = arguments[0];
    154   // Maple result: 2*D*(1-exp(-k*(r-R)))*k*exp(-k*(r-R))
    155   result.push_back(
    156       2. * params[dissociation_energy]
    157       * ( 1. - exp( - params[spring_constant] * ( r_ij.distance - params[equilibrium_distance])))
    158       * (- params[spring_constant]) * exp( - params[spring_constant] * ( r_ij.distance - params[equilibrium_distance]))
    159   );
    160   ASSERT( result.size() == 1,
    161       "PairPotential_Morse::operator() - we did not create exactly one component.");
    162   return result;
     152  result_t result = 0.;
     153  for(list_of_arguments_t::const_iterator iter = listarguments.begin();
     154      iter != listarguments.end(); ++iter) {
     155    const arguments_t &arguments = *iter;
     156    ASSERT( arguments.size() == 1,
     157        "PairPotential_Morse::operator() - requires exactly one argument.");
     158    ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
     159        arguments, getParticleTypes()),
     160        "PairPotential_Morse::operator() - types don't match with ones in arguments.");
     161    const argument_t &r_ij = arguments[0];
     162    // Maple result: 2*D*(1-exp(-k*(r-R)))*k*exp(-k*(r-R))
     163    result +=
     164        2. * params[dissociation_energy]
     165        * ( 1. - exp( - params[spring_constant] * ( r_ij.distance - params[equilibrium_distance])))
     166        * (- params[spring_constant]) * exp( - params[spring_constant] * ( r_ij.distance - params[equilibrium_distance]));
     167  }
     168  return derivative_components_t(1, result);
    163169}
    164170
    165171PairPotential_Morse::results_t
    166172PairPotential_Morse::parameter_derivative(
    167     const arguments_t &arguments,
     173    const list_of_arguments_t &listarguments,
    168174    const size_t index
    169175    ) const
    170176{
    171   ASSERT( arguments.size() == 1,
    172       "PairPotential_Morse::parameter_derivative() - requires exactly one argument.");
    173   ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
    174       arguments, getParticleTypes()),
    175       "PairPotential_Morse::operator() - types don't match with ones in arguments.");
    176   const argument_t &r_ij = arguments[0];
    177   switch (index) {
    178     case spring_constant:
    179     {
    180       // Maple result: -2*D*(1-exp(-k*(r-R)))*(-r+R)*exp(-k*(r-R))
    181       const result_t result =
    182           - 2. * params[dissociation_energy]
    183           * ( 1. - exp( - params[spring_constant] * ( r_ij.distance - params[equilibrium_distance])))
    184           * (- r_ij.distance + params[equilibrium_distance])
    185           * exp( - params[spring_constant] * ( r_ij.distance - params[equilibrium_distance]))
    186           ;
    187       return std::vector<result_t>(1, result);
    188       break;
     177  result_t result = 0.;
     178  for(list_of_arguments_t::const_iterator iter = listarguments.begin();
     179      iter != listarguments.end(); ++iter) {
     180    const arguments_t &arguments = *iter;
     181    ASSERT( arguments.size() == 1,
     182        "PairPotential_Morse::parameter_derivative() - requires exactly one argument.");
     183    ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
     184        arguments, getParticleTypes()),
     185        "PairPotential_Morse::operator() - types don't match with ones in arguments.");
     186    const argument_t &r_ij = arguments[0];
     187    switch (index) {
     188      case spring_constant:
     189      {
     190        // Maple result: -2*D*(1-exp(-k*(r-R)))*(-r+R)*exp(-k*(r-R))
     191        result +=
     192            - 2. * params[dissociation_energy]
     193            * ( 1. - exp( - params[spring_constant] * ( r_ij.distance - params[equilibrium_distance])))
     194            * (- r_ij.distance + params[equilibrium_distance])
     195            * exp( - params[spring_constant] * ( r_ij.distance - params[equilibrium_distance]))
     196            ;
     197        break;
     198      }
     199      case equilibrium_distance:
     200      {
     201        // Maple result: -2*D*(1-exp(-k*(r-R)))*k*exp(-k*(r-R))
     202        result +=
     203            - 2. * params[dissociation_energy]
     204            * ( 1. - exp( - params[spring_constant] * ( r_ij.distance - params[equilibrium_distance])))
     205            * params[spring_constant] * exp( - params[spring_constant] * ( r_ij.distance - params[equilibrium_distance]))
     206            ;
     207        break;
     208      }
     209      case dissociation_energy:
     210      {
     211        // Maple result: (1-exp(-k*(r-R)))^2
     212        result +=
     213            Helpers::pow(1.
     214                - exp( - params[spring_constant] * ( r_ij.distance - params[equilibrium_distance])),2);
     215        break;
     216      }
     217      default:
     218        ASSERT(0, "PairPotential_Morse::parameter_derivative() - derivative to unknown parameter desired.");
     219        break;
    189220    }
    190     case equilibrium_distance:
    191     {
    192       // Maple result: -2*D*(1-exp(-k*(r-R)))*k*exp(-k*(r-R))
    193       const result_t result =
    194           - 2. * params[dissociation_energy]
    195           * ( 1. - exp( - params[spring_constant] * ( r_ij.distance - params[equilibrium_distance])))
    196           * params[spring_constant] * exp( - params[spring_constant] * ( r_ij.distance - params[equilibrium_distance]))
    197           ;
    198       return std::vector<result_t>(1, result);
    199       break;
    200     }
    201     case dissociation_energy:
    202     {
    203       // Maple result: (1-exp(-k*(r-R)))^2
    204       const result_t result =
    205           Helpers::pow(1.
    206               - exp( - params[spring_constant] * ( r_ij.distance - params[equilibrium_distance])),2);
    207       return std::vector<result_t>(1, result);
    208       break;
    209     }
    210     default:
    211       ASSERT(0, "PairPotential_Morse::parameter_derivative() - derivative to unknown parameter desired.");
    212       break;
    213221  }
    214   return std::vector<result_t>(1, 0.);
    215 }
    216 
    217 FunctionModel::extractor_t
    218 PairPotential_Morse::getSpecificExtractor() const
    219 {
    220   Fragment::charges_t charges;
    221   charges.resize(getParticleTypes().size());
    222   std::transform(getParticleTypes().begin(), getParticleTypes().end(),
    223       charges.begin(), boost::lambda::_1);
    224   FunctionModel::extractor_t returnfunction =
    225       boost::bind(&Extractors::gatherDistancesFromFragment,
    226           boost::bind(&Fragment::getPositions, _1),
    227           boost::bind(&Fragment::getCharges, _1),
    228           charges,
    229           _2);
    230   return returnfunction;
     222  return results_t(1, result);
    231223}
    232224

  • src/Potentials/Specifics/PairPotential_Morse.hpp

    rc73e35 rd8821e  
    3535  friend class PotentialFactory;
    3636  // some repeated typedefs to avoid ambiguities
     37  typedef FunctionModel::list_of_arguments_t list_of_arguments_t;
    3738  typedef FunctionModel::arguments_t arguments_t;
    3839  typedef FunctionModel::result_t result_t;
    … …  
    9192  /** Evaluates the harmonic potential function for the given arguments.
    9293   *
    93    * @param arguments single distance
     94   * @param listarguments list of single distances
    9495   * @return value of the potential function
    9596   */
    96   results_t operator()(const arguments_t &arguments) const;
     97  results_t operator()(const list_of_arguments_t &listarguments) const;
    9798
    9899  /** Evaluates the derivative of the potential function.
    99100   *
    100    * @param arguments single distance
     101   * @param listarguments list of single distances
    101102   * @return vector with derivative with respect to the input degrees of freedom
    102103   */
    103   derivative_components_t derivative(const arguments_t &arguments) const;
     104  derivative_components_t derivative(const list_of_arguments_t &listarguments) const;
    104105
    105106  /** Evaluates the derivative of the function with the given \a arguments
    106107   * with respect to a specific parameter indicated by \a index.
    107108   *
    108    * \param arguments set of arguments as input variables to the function
     109   * \param listarguments list of single distances
    109110   * \param index derivative of which parameter
    110111   * \return result vector containing the derivative with respect to the given
    111112   *         input
    112113   */
    113   results_t parameter_derivative(const arguments_t &arguments, const size_t index) const;
     114  results_t parameter_derivative(const list_of_arguments_t &listarguments, const size_t index) const;
    114115
    115116  /** Returns the functor that converts argument_s into the
    … …  
    165166    return parameters_t(getParameterDimension(), std::numeric_limits<double>::max());
    166167  }
    167 
    168   /** Returns a bound function to be used with TrainingData, extracting distances
    169    * from a Fragment.
    170    *
    171    * \return bound function extracting distances from a fragment
    172    */
    173   FunctionModel::extractor_t getSpecificExtractor() const;
    174168
    175169  /** Returns a bound function to be used with TrainingData, extracting distances

  • src/Potentials/Specifics/ThreeBodyPotential_Angle.cpp

    rc73e35 rd8821e  
    136136ThreeBodyPotential_Angle::results_t
    137137ThreeBodyPotential_Angle::operator()(
    138     const arguments_t &arguments
     138    const list_of_arguments_t &listarguments
    139139    ) const
    140140{
    141   ASSERT( arguments.size() == 3,
    142       "ThreeBodyPotential_Angle::operator() - requires exactly three arguments.");
    143   ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
    144       arguments, getParticleTypes()),
    145       "ThreeBodyPotential_Angle::operator() - types don't match with ones in arguments.");
    146   const argument_t &r_ij = arguments[0]; // 01
    147   const argument_t &r_jk = arguments[2]; // 12
    148   const argument_t &r_ik = arguments[1]; // 02
    149   const result_t result =
    150       params[spring_constant]
    151              * Helpers::pow( function_theta(r_ij.distance, r_jk.distance, r_ik.distance) - params[equilibrium_distance], 2 );
    152   return std::vector<result_t>(1, result);
     141  result_t result = 0.;
     142  for(list_of_arguments_t::const_iterator iter = listarguments.begin();
     143      iter != listarguments.end(); ++iter) {
     144    const arguments_t &arguments = *iter;
     145    ASSERT( arguments.size() == 3,
     146        "ThreeBodyPotential_Angle::operator() - requires exactly three arguments.");
     147    ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
     148        arguments, getParticleTypes()),
     149        "ThreeBodyPotential_Angle::operator() - types don't match with ones in arguments.");
     150    const argument_t &r_ij = arguments[0]; // 01
     151    const argument_t &r_jk = arguments[2]; // 12
     152    const argument_t &r_ik = arguments[1]; // 02
     153    result +=
     154        params[spring_constant]
     155               * Helpers::pow( function_theta(r_ij.distance, r_jk.distance, r_ik.distance)
     156                   - params[equilibrium_distance], 2 );
     157  }
     158  return results_t(1, result);
    153159}
    154160
    155161ThreeBodyPotential_Angle::derivative_components_t
    156162ThreeBodyPotential_Angle::derivative(
    157     const arguments_t &arguments
     163    const list_of_arguments_t &listarguments
    158164    ) const
    159165{
    160   ASSERT( arguments.size() == 3,
    161       "ThreeBodyPotential_Angle::operator() - requires exactly three arguments.");
    162   ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
    163       arguments, getParticleTypes()),
    164       "ThreeBodyPotential_Angle::operator() - types don't match with ones in arguments.");
    165   derivative_components_t result;
    166   const argument_t &r_ij = arguments[0]; //01
    167   const argument_t &r_jk = arguments[2]; //12
    168   const argument_t &r_ik = arguments[1]; //02
    169   result.push_back( 2. * params[spring_constant] * ( function_theta(r_ij.distance, r_jk.distance, r_ik.distance) - params[equilibrium_distance]) );
    170   ASSERT( result.size() == 1,
    171       "ThreeBodyPotential_Angle::operator() - we did not create exactly one component.");
    172   return result;
     166  result_t result = 0.;
     167  for(list_of_arguments_t::const_iterator iter = listarguments.begin();
     168      iter != listarguments.end(); ++iter) {
     169    const arguments_t &arguments = *iter;
     170    ASSERT( arguments.size() == 3,
     171        "ThreeBodyPotential_Angle::operator() - requires exactly three arguments.");
     172    ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
     173        arguments, getParticleTypes()),
     174        "ThreeBodyPotential_Angle::operator() - types don't match with ones in arguments.");
     175    const argument_t &r_ij = arguments[0]; //01
     176    const argument_t &r_jk = arguments[2]; //12
     177    const argument_t &r_ik = arguments[1]; //02
     178    result +=
     179        2. * params[spring_constant] *
     180        ( function_theta(r_ij.distance, r_jk.distance, r_ik.distance)
     181            - params[equilibrium_distance]);
     182  }
     183  return derivative_components_t(1, result);
    173184}
    174185
    175186ThreeBodyPotential_Angle::results_t
    176187ThreeBodyPotential_Angle::parameter_derivative(
    177     const arguments_t &arguments,
     188    const list_of_arguments_t &listarguments,
    178189    const size_t index
    179190    ) const
    180191{
    181   ASSERT( arguments.size() == 3,
    182       "ThreeBodyPotential_Angle::parameter_derivative() - requires exactly three arguments.");
    183   ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
    184       arguments, getParticleTypes()),
    185       "ThreeBodyPotential_Angle::operator() - types don't match with ones in arguments.");
    186   const argument_t &r_ij = arguments[0]; //01
    187   const argument_t &r_jk = arguments[2]; //12
    188   const argument_t &r_ik = arguments[1]; //02
    189   switch (index) {
    190     case spring_constant:
    191     {
    192       const result_t result =
    193                  Helpers::pow( function_theta(r_ij.distance, r_jk.distance, r_ik.distance) - params[equilibrium_distance], 2 );
    194       return std::vector<result_t>(1, result);
    195       break;
     192  result_t result = 0.;
     193  for(list_of_arguments_t::const_iterator iter = listarguments.begin();
     194      iter != listarguments.end(); ++iter) {
     195    const arguments_t &arguments = *iter;
     196    ASSERT( arguments.size() == 3,
     197        "ThreeBodyPotential_Angle::parameter_derivative() - requires exactly three arguments.");
     198    ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
     199        arguments, getParticleTypes()),
     200        "ThreeBodyPotential_Angle::operator() - types don't match with ones in arguments.");
     201    const argument_t &r_ij = arguments[0]; //01
     202    const argument_t &r_jk = arguments[2]; //12
     203    const argument_t &r_ik = arguments[1]; //02
     204    switch (index) {
     205      case spring_constant:
     206      {
     207        result +=
     208            Helpers::pow( function_theta(r_ij.distance, r_jk.distance, r_ik.distance) - params[equilibrium_distance], 2 );
     209        break;
     210      }
     211      case equilibrium_distance:
     212      {
     213        result +=
     214            -2. * params[spring_constant]
     215                   * ( function_theta(r_ij.distance, r_jk.distance, r_ik.distance) - params[equilibrium_distance]);
     216        break;
     217      }
     218      default:
     219        ASSERT(0, "ThreeBodyPotential_Angle::parameter_derivative() - derivative to unknown parameter desired.");
     220        break;
    196221    }
    197     case equilibrium_distance:
    198     {
    199       const result_t result =
    200           -2. * params[spring_constant]
    201                  * ( function_theta(r_ij.distance, r_jk.distance, r_ik.distance) - params[equilibrium_distance]);
    202       return std::vector<result_t>(1, result);
    203       break;
    204     }
    205     default:
    206       ASSERT(0, "ThreeBodyPotential_Angle::parameter_derivative() - derivative to unknown parameter desired.");
    207       break;
    208222  }
    209   return std::vector<result_t>(1);
    210 }
    211 
    212 FunctionModel::extractor_t
    213 ThreeBodyPotential_Angle::getSpecificExtractor() const
    214 {
    215   Fragment::charges_t charges;
    216   charges.resize(getParticleTypes().size());
    217   std::transform(getParticleTypes().begin(), getParticleTypes().end(),
    218       charges.begin(), boost::lambda::_1);
    219   FunctionModel::extractor_t returnfunction =
    220       boost::bind(&Extractors::gatherDistancesFromFragment,
    221           boost::bind(&Fragment::getPositions, _1),
    222           boost::bind(&Fragment::getCharges, _1),
    223           charges,
    224           _2);
    225   return returnfunction;
     223  return results_t(1, result);
    226224}
    227225

  • src/Potentials/Specifics/ThreeBodyPotential_Angle.hpp

    rc73e35 rd8821e  
    3535  friend class PotentialFactory;
    3636  // some repeated typedefs to avoid ambiguities
     37  typedef FunctionModel::list_of_arguments_t list_of_arguments_t;
    3738  typedef FunctionModel::arguments_t arguments_t;
    3839  typedef FunctionModel::result_t result_t;
    … …  
    9091  /** Evaluates the harmonic potential function for the given arguments.
    9192   *
    92    * @param arguments single distance
     93   * @param listarguments list of three distances
    9394   * @return value of the potential function
    9495   */
    95   results_t operator()(const arguments_t &arguments) const;
     96  results_t operator()(const list_of_arguments_t &listarguments) const;
    9697
    9798  /** Evaluates the derivative of the potential function.
    9899   *
    99    * @param arguments single distance
     100   * @param listarguments list of three distances
    100101   * @return vector with derivative with respect to the input degrees of freedom
    101102   */
    102   derivative_components_t derivative(const arguments_t &arguments) const;
     103  derivative_components_t derivative(const list_of_arguments_t &listarguments) const;
    103104
    104105  /** Evaluates the derivative of the function with the given \a arguments
    105106   * with respect to a specific parameter indicated by \a index.
    106107   *
    107    * \param arguments set of arguments as input variables to the function
     108   * \param listarguments list of three distances
    108109   * \param index derivative of which parameter
    109110   * \return result vector containing the derivative with respect to the given
    110111   *         input
    111112   */
    112   results_t parameter_derivative(const arguments_t &arguments, const size_t index) const;
     113  results_t parameter_derivative(const list_of_arguments_t &listarguments, const size_t index) const;
    113114
    114115  /** Returns the functor that converts argument_s into the
    … …  
    166167    return upperbounds;
    167168  }
    168 
    169   /** Returns a bound function to be used with TrainingData, extracting distances
    170    * from a Fragment.
    171    *
    172    * \return bound function extracting distances from a fragment
    173    */
    174   FunctionModel::extractor_t getSpecificExtractor() const;
    175169
    176170  /** Returns a bound function to be used with TrainingData, extracting distances

  • src/Potentials/Specifics/unittests/ConstantPotentialUnitTest.cpp

    rc73e35 rd8821e  
    8282      Helpers::isEqual(
    8383          offset,
    84           constant( FunctionModel::arguments_t() )[0]
     84          constant( FunctionModel::list_of_arguments_t() )[0]
    8585      )
    8686  );
    … …  
    9797          0.,
    9898          constant.derivative(
    99               FunctionModel::arguments_t()
     99              FunctionModel::list_of_arguments_t()
    100100          )[0]
    101101      )
    … …  
    114114          1.,
    115115          constant.parameter_derivative(
    116               FunctionModel::arguments_t(),
     116              FunctionModel::list_of_arguments_t(),
    117117              0
    118118          )[0]

  • src/Potentials/Specifics/unittests/FourBodyPotential_ImproperUnitTest.cpp

    rc73e35 rd8821e  
    121121  for (size_t index = 0; index < input.size(); ++index) {
    122122    const FourBodyPotential_Improper::results_t result =
    123         angle( input[index] );
     123        angle( FunctionModel::list_of_arguments_t(1, input[index]) );
    124124    CPPUNIT_ASSERT(
    125125        Helpers::isEqual(
    … …  
    144144          0.,
    145145          angle.derivative(
    146               input[5]
     146              FunctionModel::list_of_arguments_t(1, input[5])
    147147          )[0],
    148148          10.
    … …  
    165165          0.,
    166166          angle.parameter_derivative(
    167               input[5],
     167              FunctionModel::list_of_arguments_t(1, input[5]),
    168168              0
    169169          )[0],
    … …  
    175175          0.,
    176176          angle.parameter_derivative(
    177               input[5],
     177              FunctionModel::list_of_arguments_t(1, input[5]),
    178178              1
    179179          )[0],

  • src/Potentials/Specifics/unittests/FourBodyPotential_TorsionUnitTest.cpp

    rc73e35 rd8821e  
    121121  for (size_t index = 0; index < input.size(); ++index) {
    122122    const FourBodyPotential_Torsion::results_t result =
    123         angle( input[index] );
     123        angle( FunctionModel::list_of_arguments_t(1, input[index]) );
    124124    CPPUNIT_ASSERT(
    125125        Helpers::isEqual(
    … …  
    144144          0.,
    145145          angle.derivative(
    146               input[5]
     146              FunctionModel::list_of_arguments_t(1, input[5])
    147147          )[0],
    148148          10.
    … …  
    165165          0.,
    166166          angle.parameter_derivative(
    167               input[5],
     167              FunctionModel::list_of_arguments_t(1, input[5]),
    168168              0
    169169          )[0],
    … …  
    175175          0.,
    176176          angle.parameter_derivative(
    177               input[5],
     177              FunctionModel::list_of_arguments_t(1, input[5]),
    178178              1
    179179          )[0],

  • src/Potentials/Specifics/unittests/Makefile.am

    rc73e35 rd8821e  
    3838POTENTIALSSPECIFICSLIBS = \
    3939        ../libMolecuilderPotentials.la \
    40         ../libMolecuilderFragmentationSetValues.la \
    4140        ../libMolecuilderFragmentation.la \
    4241        ../libMolecuilderFunctionApproximation.la \
    4342        ../libMolecuilderRandomNumbers.la \
     43        ../libMolecuilderFragmentationSummation.la \
     44        ../libMolecuilderFragmentation_getFromKeysetStub.la \
    4445        $(top_builddir)/LinearAlgebra/src/LinearAlgebra/libLinearAlgebra.la \
    4546        ${CodePatterns_LIBS} \

  • src/Potentials/Specifics/unittests/ManyBodyPotential_TersoffUnitTest.cpp

    rc73e35 rd8821e  
    299299        arg.globalid = index; // this is needed for the triplefunction to the configuration
    300300        FunctionModel::arguments_t args(1,arg);
    301         const ManyBodyPotential_Tersoff::results_t res = tersoff(args);
     301        const ManyBodyPotential_Tersoff::results_t res =
     302            tersoff( FunctionModel::list_of_arguments_t(1, args) );
    302303        temp += res[0];
    303304      }
    … …  
    336337//          0.,
    337338//          tersoff.derivative(
    338 //              FunctionModel::arguments_t(1,argument_t(1.))
     339//              FunctionModel::list_of_arguments_t(1, FunctionModel::arguments_t(1,argument_t(1.)))
    339340//          )[0]
    340341//      )
    … …  
    361362//          0.,
    362363//          tersoff.parameter_derivative(
    363 //              FunctionModel::arguments_t(1,argument_t(1.)),
     364//              FunctionModel::list_of_arguments_t(1, FunctionModel::arguments_t(1,argument_t(1.))),
    364365//              0
    365366//          )[0]
    … …  
    370371//          0.,
    371372//          tersoff.parameter_derivative(
    372 //              FunctionModel::arguments_t(1,argument_t(1.)),
     373//              FunctionModel::list_of_arguments_t(1, FunctionModel::arguments_t(1,argument_t(1.))),
    373374//              1
    374375//          )[0]
    … …  
    379380//          1.,
    380381//          tersoff.parameter_derivative(
    381 //              FunctionModel::arguments_t(1,argument_t(1.)),
     382//              FunctionModel::list_of_arguments_t(1, FunctionModel::arguments_t(1,argument_t(1.))),
    382383//              2
    383384//          )[0]

  • src/Potentials/Specifics/unittests/PairPotential_HarmonicUnitTest.cpp

    rc73e35 rd8821e  
    102102        Helpers::isEqual(
    103103            output[index],
    104             harmonic( FunctionModel::arguments_t(1,arg) )[0]
     104            harmonic( FunctionModel::list_of_arguments_t(1, FunctionModel::arguments_t(1,arg)) )[0]
    105105        )
    106106    );
    … …  
    122122          0.,
    123123          harmonic.derivative(
    124               FunctionModel::arguments_t(1,arg)
     124              FunctionModel::list_of_arguments_t(1, FunctionModel::arguments_t(1,arg))
    125125          )[0]
    126126      )
    … …  
    143143          0.,
    144144          harmonic.parameter_derivative(
    145               FunctionModel::arguments_t(1,arg),
     145              FunctionModel::list_of_arguments_t(1, FunctionModel::arguments_t(1,arg)),
    146146              0
    147147          )[0]
    … …  
    152152          0.,
    153153          harmonic.parameter_derivative(
    154               FunctionModel::arguments_t(1,arg),
     154              FunctionModel::list_of_arguments_t(1, FunctionModel::arguments_t(1,arg)),
    155155              1
    156156          )[0]

  • src/Potentials/Specifics/unittests/PairPotential_LennardJonesUnitTest.cpp

    rc73e35 rd8821e  
    118118        Helpers::isEqual(
    119119            output[index],
    120             lj( FunctionModel::arguments_t(1,arg) )[0],
     120            lj( FunctionModel::list_of_arguments_t(1, FunctionModel::arguments_t(1,arg)) )[0],
    121121            1.e-4/std::numeric_limits<double>::epsilon() // only compare four digits
    122122        )
    … …  
    140140          0.,
    141141          lj.derivative(
    142               FunctionModel::arguments_t(1,arg)
     142              FunctionModel::list_of_arguments_t(1, FunctionModel::arguments_t(1,arg))
    143143          )[0],
    144144          1.e+6
    … …  
    162162          -1.,
    163163          lj.parameter_derivative(
    164               FunctionModel::arguments_t(1,arg),
     164              FunctionModel::list_of_arguments_t(1, FunctionModel::arguments_t(1,arg)),
    165165              0
    166166          )[0],
    … …  
    172172          0.,
    173173          lj.parameter_derivative(
    174               FunctionModel::arguments_t(1,arg),
     174              FunctionModel::list_of_arguments_t(1, FunctionModel::arguments_t(1,arg)),
    175175              1
    176176          )[0],

  • src/Potentials/Specifics/unittests/PairPotential_MorseUnitTest.cpp

    rc73e35 rd8821e  
    120120        Helpers::isEqual(
    121121            output[index],
    122             Morse( FunctionModel::arguments_t(1,arg) )[0],
     122            Morse( FunctionModel::list_of_arguments_t(1, FunctionModel::arguments_t(1,arg)) )[0],
    123123            1.e-4/std::numeric_limits<double>::epsilon() // only compare four digits
    124124        )
    … …  
    141141          0.,
    142142          Morse.derivative(
    143               FunctionModel::arguments_t(1,arg)
     143              FunctionModel::list_of_arguments_t(1, FunctionModel::arguments_t(1,arg))
    144144          )[0],
    145145          1.e+6
    … …  
    163163          0.,
    164164          Morse.parameter_derivative(
    165               FunctionModel::arguments_t(1,arg),
     165              FunctionModel::list_of_arguments_t(1, FunctionModel::arguments_t(1,arg)),
    166166              0
    167167          )[0],
    … …  
    173173          0.,
    174174          Morse.parameter_derivative(
    175               FunctionModel::arguments_t(1,arg),
     175              FunctionModel::list_of_arguments_t(1, FunctionModel::arguments_t(1,arg)),
    176176              1
    177177          )[0],
    … …  
    183183          0.,
    184184          Morse.parameter_derivative(
    185               FunctionModel::arguments_t(1,arg),
     185              FunctionModel::list_of_arguments_t(1, FunctionModel::arguments_t(1,arg)),
    186186              2
    187187          )[0],

  • src/Potentials/Specifics/unittests/ThreeBodyPotential_AngleUnitTest.cpp

    rc73e35 rd8821e  
    112112  for (size_t index = 0; index < input.size(); ++index) {
    113113    const ThreeBodyPotential_Angle::results_t result =
    114         angle( input[index] );
     114        angle( FunctionModel::list_of_arguments_t(1, input[index]) );
    115115    CPPUNIT_ASSERT(
    116116        Helpers::isEqual(
    … …  
    135135          0.,
    136136          angle.derivative(
    137               input[5]
     137              FunctionModel::list_of_arguments_t(1, input[5])
    138138          )[0],
    139139          10.
    … …  
    156156          0.,
    157157          angle.parameter_derivative(
    158               input[5],
     158              FunctionModel::list_of_arguments_t(1, input[5]),
    159159              0
    160160          )[0],
    … …  
    166166          0.,
    167167          angle.parameter_derivative(
    168               input[5],
     168              FunctionModel::list_of_arguments_t(1, input[5]),
    169169              1
    170170          )[0],

  • src/Potentials/helpers.hpp

    rc73e35 rd8821e  
    108108  }
    109109
     110  inline FunctionModel::list_of_arguments_t
     111  returnEmptyListArguments()
     112  {
     113    return FunctionModel::list_of_arguments_t();
     114  }
     115
    110116}; /* namespace Helpers */
    111117

  • src/Potentials/unittests/CompoundPotentialUnitTest.cpp

    rc73e35 rd8821e  
    159159    argument_t firstarg(argument_t::indices_t(0,1), argument_t::types_t(0,1), input[index]);
    160160    argument_t secondarg(argument_t::indices_t(0,2), argument_t::types_t(0,1), input[index]);
    161     FunctionModel::arguments_t args;
    162     args += firstarg,secondarg;
    163     const double result = compound( args )[0];
     161    FunctionModel::list_of_arguments_t listargs;
     162    listargs += FunctionModel::arguments_t(1,firstarg),FunctionModel::arguments_t(1,secondarg);
     163    const double result = compound( listargs )[0];
    164164    CPPUNIT_ASSERT(
    165165        Helpers::isEqual(
    … …  
    205205  argument_t firstarg(argument_t::indices_t(0,1), argument_t::types_t(0,1), c);
    206206  argument_t secondarg(argument_t::indices_t(0,2), argument_t::types_t(0,1), c);
    207   FunctionModel::arguments_t args;
    208   args += firstarg,secondarg;
    209   {
    210     const double result =
    211         compound.parameter_derivative(
    212             args,
     207  FunctionModel::list_of_arguments_t listargs;
     208  listargs += FunctionModel::arguments_t(1,firstarg),FunctionModel::arguments_t(1,secondarg);
     209  {
     210    const double result =
     211        compound.parameter_derivative(
     212            listargs,
    213213            0)[0];
    214214    CPPUNIT_ASSERT(
    … …  
    223223    const double result =
    224224        compound.parameter_derivative(
    225             args,
     225            listargs,
    226226            1)[0];
    227227    CPPUNIT_ASSERT(
    … …  
    236236    const double result =
    237237        compound.parameter_derivative(
    238             args,
     238            listargs,
    239239            2)[0];
    240240    CPPUNIT_ASSERT(
    … …  
    249249    const double result =
    250250        compound.parameter_derivative(
    251             args,
     251            listargs,
    252252            3)[0];
    253253    CPPUNIT_ASSERT(

  • src/Potentials/unittests/Makefile.am

    rc73e35 rd8821e  
    3434        ../libMolecuilderPotentials.la \
    3535        ../libMolecuilderFragmentation.la \
     36        ../libMolecuilderFunctionApproximation.la \
    3637        ../libMolecuilderFragmentationSummation.la \
    37         ../libMolecuilderFunctionApproximation.la \
    3838        ../libMolecuilderFragmentation_getFromKeysetStub.la \
    3939        ../libMolecuilderRandomNumbers.la \

  • src/UIElements/CommandLineUI/CommandLineWindow.cpp

    rc73e35 rd8821e  
    7373    if (AQ.isActionKnownByName(*CommandRunner)) {
    7474//      LOG1, "INFO: Checking presence of " << *CommandRunner << ": " << "queuing.");
    75       AQ.queueAction(*CommandRunner);
     75      if (AQ.getLastActionOk())
     76        AQ.queueAction(*CommandRunner);
     77      else
     78        LOG(1, "INFO: Skipping " << *CommandRunner << " as last action failed.");
    7679    } else {
    7780//      LOG(1, "INFO: Checking presence of " << *CommandRunner << ": " << "absent.");

  • src/UIElements/Makefile.am

    rc73e35 rd8821e  
    254254        libMolecuilderGraph.la \
    255255        libMolecuilderFilling.la \
    256         libMolecuilder.la
     256        libMolecuilder.la \
     257        libMolecuilderFragmentationAutomation.la \
     258        libMolecuilderFragmentation_getFromKeyset.la \
     259        libMolecuilderFragmentation.la \
     260        libMolecuilderJobs.la
    257261if CONDJOBMARKET
    258262libMolecuilderUI_la_LIBADD += \
    259         libMolecuilderFragmentationAutomation.la
    260 endif
    261 libMolecuilderUI_la_LIBADD += \
    262         libMolecuilderFragmentation_getFromKeyset.la \
    263         libMolecuilderFragmentation.la
    264 if CONDJOBMARKET
    265 libMolecuilderUI_la_LIBADD += \
    266         libMolecuilderJobs.la \
    267263        ${JobMarket_Controller_LIBS}
    268264endif

  • src/UIElements/Views/Qt4/QtHomologyList.cpp

    rc73e35 rd8821e  
    184184    if (!data.getTrainingInputs().empty()) {
    185185      // generate QSeisData
    186       const TrainingData::InputVector_t &inputs = data.getTrainingInputs();
     186      const TrainingData::InputVector_t &inputs = data.getAllArguments();
    187187      const TrainingData::OutputVector_t &outputs = data.getTrainingOutputs();
    188188      std::vector<double> xvalues;

  • src/cleanUp.cpp

    rc73e35 rd8821e  
    5151#include "Actions/OptionRegistry.hpp"
    5252
     53#include "Fragmentation/Summation/Containers/FragmentationResultContainer.hpp"
     54
    5355#include "RandomNumbers/RandomNumberDistributionFactory.hpp"
    5456#include "RandomNumbers/RandomNumberEngineFactory.hpp"
    … …  
    9395{
    9496  // make sure that ActionQueue is already purged!
    95 
     97  FragmentationResultContainer::purgeInstance();
    9698  Chronos::purgeInstance();
    9799  RandomNumberDistributionFactory::purgeInstance();

  • src/documentation/constructs/potentials.dox

    rc73e35 rd8821e  
    7272 *    perform the fit.
    7373 *
     74 * \section potentials-fit-potential-action What happens in FitPotentialAction.
     75 *
     76 *  First, either a potential file is parsed via PotentialDeserializer or charges
     77 *  and a potential type from the given options. This is used to instantiate
     78 *  EmpiricalPotentials via the PotentialFactory, stored within the
     79 *  PotentialRegistry. This is the available set of potentials (without requiring
     80 *  any knowledge as to the nature of the fragment employed in fitting).
     81 *
     82 *  Second, the given fragment is used to get a suitable HomologyGraph from
     83 *  the World's HomologyContainer. This is given to a CompoundPotential, that in
     84 *  turn browses through the PotentialRegistry, picking out those
     85 *  EmpiricalPotential's that match with a subset of the FragmentNode's of the
     86 *  given graph. These are stored as a list of FunctionModel's within the
     87 *  CompoundPotential instance. Here comes the specific fragment into play,
     88 *  picking a subset from the available potentials.
     89 *
     90 *  Third, we need to setup the training data. For this we need vectors of input
     91 *  and output data that are obtained from the HomologyContainer with the
     92 *  HomologyGraph as key. The output vector in our case is simply a number
     93 *  (although the interface allows for more). The input vector is the set of
     94 *  distances. In order to pre-process the input data for the specific model
     95 *  a filter is required in the TrainingData's constructor. The purpose of the
     96 *  filter is to pick out the subset of distance arguments for each model one
     97 *  after the other and concatenate them to a list. On evaluation of the model
     98 *  this concatenated list of distances is given to the model and it may easily
     99 *  dissect the list and hand over each contained potential its subset of
     100 *  arguments. See Extractors for more information.
     101 *
     102 *  Afterwards, training may commence: The goal is to find a set of parameters
     103 *  for the model such that it as good as possible reproduces the output vector
     104 *  for a given input vector. This non-linear regression is contained in the
     105 *  levmar package and its functionality is wrapped in the FunctionApproximation
     106 *  class. An instance is initialized with both the gathered training data and
     107 *  the model containing a list of potentials. See
     108 *  [FunctionApproximation-details] for more details.
     109 *
     110 *  During the fitting procedure, first the derivatives of the model is checked
     111 *  for consistency, then the model is initialized with a sensible guess of
     112 *  starting parameters, and afterwards the Levenberg-Marquardt algorithm
     113 *  commences that makes numerous calls to evaluate the model and its derivative
     114 *  to find the minimum in the L2-norm.
     115 *
     116 *  This is done more than once as high-dimensional regression is sensititive the
     117 *  the starting values as there are possible numerous local minima. The lowest
     118 *  of the found minima is taken, either via a given threshold or the best of a
     119 *  given number of attempts.
     120 *
     121 *  Eventually, these parameters of the best model are streamed via
     122 *  PotentialSerializer back into a potential file. Each EmpiricalPotential in
     123 *  the CompoundPotential making up the whole model is also a
     124 *  SerializablePotential. Hence, each in turn writes a single line with its
     125 *  respective subset of parameters and particle types, describing this
     126 *  particular fit function.
     127 *
     128 * \section potentials-function-evaluation How does the model evaluation work
     129 *
     130 *  We now come to the question of how the model and its derivative are actually
     131 *  evaluated. We have an input vector from the training data and we have the
     132 *  model with a specific set of parameters.
     133 *
     134 *  FunctionModel is just an abstract interface that is implemented by the
     135 *  potential functions, such as CompoundPotential, that combines multiple
     136 *  potentials into a single function for fitting, or PairPotential_Harmonic,
     137 *  that is a specific fit function for harmonic bonds.
     138 *
     139 *  The main issue with the evaluation is picking the right set of distances from
     140 *  ones given in the input vector and feed it to each potential contained in
     141 *  CompoundPotential. Note that the distances have already been prepared by
     142 *  the TrainingData instantiation.
     143 *
     144 *  Initially, the HomologyGraph only contains a list of configurations of a
     145 *  specific fragments (i.e. the position of each atom in the fragment) and an
     146 *  energy value. These first have to be converted into distances.
     147 *
     148 *
    74149 * \section potentials-howto-use Howto use the potentials
    75150 *

  • tests/Fragmentations/Makefile.am

    rc73e35 rd8821e  
    44        testsuite.at \
    55        $(TESTSUITE) \
    6         analyzer.in \
    76        atlocal.in \
    8         joiner.in \
    97        molecuilder.in \
    108        package.m4 \
    11         $(srcdir)/Analyzing/heptan \
    129        $(srcdir)/Fragmenting/1_2-dimethoxyethane \
    1310        $(srcdir)/Fragmenting/1_2-dimethylbenzene \
    … …  
    2825        $(srcdir)/Fragmenting/proline \
    2926        $(srcdir)/Fragmenting/putrescine \
    30         $(srcdir)/Fragmenting/tartaric_acid \
    31         $(srcdir)/Joining/heptan
     27        $(srcdir)/Fragmenting/tartaric_acid
    3228
    3329TESTSUITE = $(srcdir)/testsuite
    … …  
    139135        Fragmenting/tartaric_acid/testsuite-fragmenting-tartaric_acid-order4.at \
    140136        Fragmenting/tartaric_acid/testsuite-fragmenting-tartaric_acid-order5.at \
    141         Fragmenting/tartaric_acid/testsuite-fragmenting-tartaric_acid-order6.at \
    142         Joining/heptan/testsuite-joining-heptan.at \
    143         Analyzing/heptan/testsuite-analyzing-heptan.at
     137        Fragmenting/tartaric_acid/testsuite-fragmenting-tartaric_acid-order6.at
    144138
    145139max_jobs = 4

  • tests/Fragmentations/testsuite.at

    rc73e35 rd8821e  
    185185m4_include(Fragmenting/tartaric_acid/testsuite-fragmenting-tartaric_acid-order5.at)
    186186m4_include(Fragmenting/tartaric_acid/testsuite-fragmenting-tartaric_acid-order6.at)
    187 
    188 # Joining of heptan
    189 m4_include(Joining/heptan/testsuite-joining-heptan.at)
    190 
    191 # Analyzing of heptan
    192 m4_include(Analyzing/heptan/testsuite-analyzing-heptan.at)

  • tests/Makefile.am

    rc73e35 rd8821e  
    22
    33SUBDIRS += \
     4        Calculations \
    45        CodeChecks \
    56        regression \
    … …  
    89        Python \
    910        JobMarket
     11
     12.PHONY: extracheck installextracheck
     13
     14extracheck:
     15        cd Calculations && make extracheck
     16installextracheck:
     17        cd Calculations && make installextracheck

  • tests/regression/Fragmentation/testsuite-fragmentation.at

    rc73e35 rd8821e  
    3030# check storing saturated fragment
    3131m4_include([Fragmentation/StoreSaturatedFragment/testsuite-fragmentation-store-saturated-fragment.at])
     32
     33# check automation
     34m4_include([Fragmentation/FragmentationAutomation/testsuite-fragmentation-fragmentation-automation.at])

  • tests/regression/Makefile.am

    rc73e35 rd8821e  
    7272        $(srcdir)/Fragmentation/testsuite-fragmentation.at \
    7373        $(srcdir)/Fragmentation/AnalyseFragmentationResults/testsuite-fragmentation-analyse-fragment-results.at \
     74        $(srcdir)/Fragmentation/FragmentationAutomation/testsuite-fragmentation-fragmentation-automation.at \
    7475        $(srcdir)/Fragmentation/FragmentMolecule/testsuite-fragmentation-fragment-molecule.at \
    7576        $(srcdir)/Fragmentation/FragmentMolecule-MaxOrder/testsuite-fragmentation-fragment-molecule-maxorder.at \
Note: See TracChangeset for help on using the changeset viewer.