Changeset c62f96

Timestamp:

Dec 19, 2012, 3:25:30 PM (12 years ago)

Author:

Frederik Heber <heber@…>

Branches:

Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Added_ParseSaveFragmentResults, AddingActions_SaveParseParticleParameters, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_ParticleName_to_Atom, Adding_StructOpt_integration_tests, AtomFragments, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, CombiningParticlePotentialParsing, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_BoundInBox_CenterInBox_MoleculeActions, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_PopActions, Fix_QtFragmentList_sorted_selection, Fix_Restrictedkeyset_FragmentMolecule, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, Fix_fitting_potentials, Fixes, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, FragmentAction_writes_AtomFragments, FragmentMolecule_checks_bonddegrees, GeometryObjects, Gui_Fixes, Gui_displays_atomic_force_velocity, ImplicitCharges, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, MoreRobust_FragmentAutomation, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PdbParser_setsAtomName, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, Rewrite_FitPartialCharges, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, SaturateAtoms_singleDegree, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, Switchable_LogView, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, TremoloParser_setsAtomName, Ubuntu_1604_changes, stable

Children:

65c42c

Parents:

4f82f8

git-author:

Frederik Heber <heber@…> (10/03/12 10:47:16)

git-committer:

Frederik Heber <heber@…> (12/19/12 15:25:30)

Message:

LevMartester uses now FunctionApproximation.

Location:

src

Files:

• FunctionApproximation/FunctionApproximation.cpp (modified) (5 diffs)
• LevMartester.cpp (modified) (5 diffs)
• Makefile.am (modified) (1 diff)

src/FunctionApproximation/FunctionApproximation.cpp

@@ -145 +145 @@
 
   m = model.getParameterDimension();
-  n = 1;
+  n = output_data.size();
   const FunctionModel::parameters_t params = model.getParameters();
-  p = new double(m);
-  x = new double(n);
-  size_t index = 0;
-  for(FunctionModel::parameters_t::const_iterator paramiter = params.begin();
-      paramiter != params.end();
-      ++paramiter, ++index) {
-    p[index] = params[index];
-  }
-  x[0] = 0.;
+  {
+    p = new double[m];
+    size_t index = 0;
+    for(FunctionModel::parameters_t::const_iterator paramiter = params.begin();
+        paramiter != params.end();
+        ++paramiter, ++index) {
+      p[index] = *paramiter;
+    }
+  }
+  {
+    x = new double[n];
+    size_t index = 0;
+    for(outputs_t::const_iterator outiter = output_data.begin();
+        outiter != output_data.end();
+        ++outiter, ++index) {
+      x[index] = (*outiter)[0];
+    }
+  }
 
   {
@@ -214 +223 @@
   }
 
-  delete p;
-  delete x;
+  delete[] p;
+  delete[] x;
 }
 
@@ -229 +238 @@
       "FunctionApproximation::evaluate() - LevMar expects "+toString(m)
       +" parameters but the model function expects "+toString(model.getParameterDimension())+".");
-  ASSERT( n == 1,
-      "FunctionApproximation::evaluate() - we return precisely a single value, LevMar expects something else.");
+  ASSERT( (size_t)n == output_data.size(),
+      "FunctionApproximation::evaluate() - LevMar expects "+toString(n)
+      +" outputs but we provide "+toString(output_data.size())+".");
   FunctionModel::parameters_t params(m, 0.);
   std::copy(p, p+m, params.begin());
@@ -236 +246 @@
 
   // then evaluate
-  if (!input_data.empty()) {
+  if (!output_data.empty()) {
     inputs_t::const_iterator initer = input_data.begin();
     outputs_t::const_iterator outiter = output_data.begin();
@@ -251 +261 @@
       x[index] = std::accumulate(differences.begin(), differences.end(), 0.);
     }
-  } else {
-    for (register int i=0;i<n;++i)
-      x[i] = 0.;
-  }
-}
+  }
+}

src/LevMartester.cpp

@@ -55 +55 @@
 #include "Fragmentation/Homology/HomologyContainer.hpp"
 #include "Fragmentation/SetValues/Fragment.hpp"
+#include "FunctionApproximation/FunctionApproximation.hpp"
+#include "FunctionApproximation/FunctionModel.hpp"
+#include "Potentials/Specifics/PairPotential_Harmonic.hpp"
 
 namespace po = boost::program_options;
@@ -67 +70 @@
   }
   return HomologyGraph();
-}
-
-/* Meyer's (reformulated) problem, minimum at (2.48, 6.18, 3.45) */
-void meyer(double *p, double *x, int m, int n, void *data)
-{
-register int i;
-double ui;
-
-        for(i=0; i<n; ++i){
-                ui=0.45+0.05*i;
-                x[i]=p[0]*exp(10.0*p[1]/(ui+p[2]) - 13.0);
-        }
-}
-
-void jacmeyer(double *p, double *jac, int m, int n, void *data)
-{
-register int i, j;
-double ui, tmp;
-
-  for(i=j=0; i<n; ++i){
-          ui=0.45+0.05*i;
-          tmp=exp(10.0*p[1]/(ui+p[2]) - 13.0);
-
-          jac[j++]=tmp;
-          jac[j++]=10.0*p[0]*tmp/(ui+p[2]);
-          jac[j++]=-10.0*p[0]*p[1]*tmp/((ui+p[2])*(ui+p[2]));
-  }
 }
 
@@ -151 +127 @@
 
   // Afterwards we go through all of this type and gather the distance and the energy value
-  typedef std::vector< std::pair<double, double> > DistanceEnergyVector_t;
-  DistanceEnergyVector_t DistanceEnergyVector;
+  typedef std::pair<
+      FunctionApproximation::inputs_t,
+      FunctionApproximation::outputs_t> InputOutputVector_t;
+  InputOutputVector_t DistanceEnergyVector;
   std::pair<HomologyContainer::const_iterator, HomologyContainer::const_iterator> range =
       homologies.getHomologousGraphs(graph);
@@ -160 +138 @@
     const Fragment::charges_t charges = fragment.getCharges();
     const Fragment::positions_t positions = fragment.getPositions();
-    std::vector< Vector > DistanceVectors;
+    std::vector< std::pair<Vector, size_t> > DistanceVectors;
     for (Fragment::charges_t::const_iterator chargeiter = charges.begin();
         chargeiter != charges.end(); ++chargeiter) {
@@ -167 +145 @@
         const size_t steps = std::distance(charges.begin(), chargeiter);
         std::advance(positer, steps);
-        DistanceVectors.push_back(Vector((*positer)[0], (*positer)[1], (*positer)[2]));
+        DistanceVectors.push_back(
+            std::make_pair(Vector((*positer)[0], (*positer)[1], (*positer)[2]),
+                steps));
       }
     }
-    ASSERT( DistanceVectors.size() == (size_t)2,
-        "main() - found not exactly two carbon atoms in fragment.");
-    const double distance = DistanceVectors[0].distance(DistanceVectors[1]);
-    const double energy = iter->second.second;
-    DistanceEnergyVector.push_back( std::make_pair(distance, energy) );
+    if (DistanceVectors.size() == (size_t)2) {
+      argument_t arg;
+      arg.indices.first = DistanceVectors[0].second;
+      arg.indices.second = DistanceVectors[1].second;
+      arg.distance = DistanceVectors[0].first.distance(DistanceVectors[1].first);
+      const double energy = iter->second.second;
+      DistanceEnergyVector.first.push_back( FunctionModel::arguments_t(1,arg) );
+      DistanceEnergyVector.second.push_back( FunctionModel::results_t(1,energy) );
+    } else {
+      ELOG(2, "main() - found not exactly two carbon atoms in fragment "
+          << fragment << ".");
+    }
   }
-  LOG(1, "INFO: I gathered the following " << DistanceEnergyVector.size() << " data pairs: ");
-  for (DistanceEnergyVector_t::const_iterator dataiter = DistanceEnergyVector.begin();
-      dataiter != DistanceEnergyVector.end(); ++dataiter) {
-    LOG(1, "INFO: " << dataiter->first << " with energy " << dataiter->second);
+  // print training data for debugging
+  {
+    LOG(1, "INFO: I gathered the following (" << DistanceEnergyVector.first.size()
+        << "," << DistanceEnergyVector.second.size() << ") data pairs: ");
+    FunctionApproximation::inputs_t::const_iterator initer = DistanceEnergyVector.first.begin();
+    FunctionApproximation::outputs_t::const_iterator outiter = DistanceEnergyVector.second.begin();
+    for (; initer != DistanceEnergyVector.first.end(); ++initer, ++outiter) {
+      LOG(1, "INFO: (" << (*initer)[0].indices.first << "," << (*initer)[0].indices.second
+          << ") " << (*initer)[0].distance << " with energy " << *outiter);
+    }
   }
   // NOTICE that distance are in bohrradi as they come from MPQC!
 
-  // let levmar optimize
-  register int i, j;
-  int ret;
-  double p[5], // 5 is max(2, 3, 5)
-             x[16]; // 16 is max(2, 3, 5, 6, 16)
-  int m, n;
-  double opts[LM_OPTS_SZ], info[LM_INFO_SZ];
+  PairPotential_Harmonic harmonic(1., 1.);
+  FunctionApproximation approximator(1, 1, harmonic);
+  approximator();
+  const FunctionModel::parameters_t params = harmonic.getParameters();
 
-  opts[0]=LM_INIT_MU; opts[1]=1E-15; opts[2]=1E-15; opts[3]=1E-20;
-  opts[4]= LM_DIFF_DELTA; // relevant only if the Jacobian is approximated using finite differences; specifies forward differencing
-  //opts[4]=-LM_DIFF_DELTA; // specifies central differencing to approximate Jacobian; more accurate but more expensive to compute!
-
-  m=3; n=16;
-  p[0]=8.85; p[1]=4.0; p[2]=2.5;
-  x[0]=34.780;        x[1]=28.610; x[2]=23.650; x[3]=19.630;
-  x[4]=16.370;        x[5]=13.720; x[6]=11.540; x[7]=9.744;
-  x[8]=8.261; x[9]=7.030; x[10]=6.005; x[11]=5.147;
-  x[12]=4.427;        x[13]=3.820; x[14]=3.307; x[15]=2.872;
-  ret=dlevmar_der(meyer, jacmeyer, p, x, m, n, 1000, opts, info, NULL, NULL, NULL); // with analytic Jacobian
-
- { double *work, *covar;
-  work=(double *)malloc((LM_DIF_WORKSZ(m, n)+m*m)*sizeof(double));
-  if(!work){
-      fprintf(stderr, "memory allocation request failed in main()\n");
-    exit(1);
-  }
-  covar=work+LM_DIF_WORKSZ(m, n);
-
-  ret=dlevmar_dif(meyer, p, x, m, n, 1000, opts, info, work, covar, NULL); // no Jacobian, caller allocates work memory, covariance estimated
-
-  printf("Covariance of the fit:\n");
-  for(i=0; i<m; ++i){
-    for(j=0; j<m; ++j)
-      printf("%g ", covar[i*m+j]);
-    printf("\n");
-  }
-  printf("\n");
-
-  free(work);
- }
-
-//  {
-//   double err[16];
-//   dlevmar_chkjac(meyer, jacmeyer, p, m, n, NULL, err);
-//   for(i=0; i<n; ++i) printf("gradient %d, err %g\n", i, err[i]);
-//  }
-
-  printf("Levenberg-Marquardt returned %d in %g iter, reason %g\nSolution: ", ret, info[5], info[6]);
-  for(i=0; i<m; ++i)
-    printf("%.7g ", p[i]);
-  printf("\n\nMinimization info:\n");
-  for(i=0; i<LM_INFO_SZ; ++i)
-    printf("%g ", info[i]);
-  printf("\n");
+  LOG(0, "RESULT: Best parameters are " << params[0] << " and " << params[1] << ".");
 
   return 0;

src/Makefile.am

@@ -482 +482 @@
 LevMartester_LDFLAGS = $(AM_LDFLAGS) $(BOOST_SERIALIZATION_LDFLAGS) $(BOOST_PROGRAM_OPTIONS_LDFLAGS) $(BOOST_FILESYSTEM_LDFLAGS) $(BOOST_SYSTEM_LDFLAGS) ${CodePatterns_LDFLAGS} $(LEVMAR_LDFLAGS)
 LevMartester_LDADD = \
+        libMolecuilderPotentials.la \
+        libMolecuilderFunctionApproximation.la \
         libMolecuilderFragmentation.la \
         libMolecuilderFragmentation_getFromKeysetStub.la \