Ignore:
Timestamp:
Feb 27, 2013, 12:39:33 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:
c16cce
Parents:
d76ce74
git-author:
Frederik Heber <heber@…> (11/30/12 17:24:28)
git-committer:
Frederik Heber <heber@…> (02/27/13 12:39:33)
Message:

LevMartester now fits potentials given on the command-line.

  • instead of three we only fit just one potential.
  • potential is instantiated via a PotentialFactory.
File:
1 edited

Legend:

Unmodified
Added
Removed
  • src/LevMartester.cpp

    rd76ce74 r3e236e  
    6969#include "FunctionApproximation/writeDistanceEnergyTable.hpp"
    7070#include "Helpers/defs.hpp"
     71#include "Potentials/PotentialFactory.hpp"
     72#include "Potentials/PotentialRegistry.hpp"
    7173#include "Potentials/Specifics/PairPotential_Morse.hpp"
    7274#include "Potentials/Specifics/PairPotential_Angle.hpp"
     
    172174      ("help", "produce help message")
    173175      ("homology-file", po::value< boost::filesystem::path >(), "homology file to parse")
     176      ("fit-potential", po::value< std::string >(), "potential type to fit")
     177      ("charges", po::value< FunctionModel::charges_t >()->multitoken(), "charges specifying the fragment")
    174178  ;
    175179
     
    183187  }
    184188
     189  // homology-file
    185190  boost::filesystem::path homology_file;
    186191  if (vm.count("homology-file")) {
     
    188193    LOG(1, "INFO: Parsing " << homology_file.string() << ".");
    189194  } else {
    190     LOG(0, "homology-file level was not set.");
    191   }
     195    ELOG(0, "homology file (homology-file) was not set.");
     196    return 1;
     197  }
     198
     199  // type of potential to fit
     200  std::string potentialtype;
     201  if (vm.count("fit-potential")) {
     202    potentialtype = vm["fit-potential"].as<std::string>();
     203  } else {
     204    ELOG(0, "potential type to fit (fit-potential) was not set.");
     205    return 1;
     206  }
     207
     208  // charges
     209  FunctionModel::charges_t charges;
     210  if (vm.count("charges")) {
     211    charges = vm["charges"].as< FunctionModel::charges_t >();
     212  } else {
     213    ELOG(0, "Vector of charges specifying the fragment (charges) was not set.");
     214    return 1;
     215  }
     216
     217  // parse homologies into container
    192218  HomologyContainer homologies;
    193219  if (boost::filesystem::exists(homology_file)) {
     
    197223      ia >> homologies;
    198224    } else {
    199       ELOG(2, "Failed to parse from " << homology_file.string() << ".");
     225      ELOG(0, "Failed to parse from " << homology_file.string() << ".");
     226      return 1;
    200227    }
    201228    returnstream.close();
    202229  } else {
    203230    ELOG(0, homology_file << " does not exist.");
     231    return 1;
    204232  }
    205233
     
    212240  }
    213241
    214   /******************** Angle TRAINING ********************/
    215   FunctionModel::parameters_t angleparams(PairPotential_Angle::MAXPARAMS, 0.);
     242  LOG(0, "STATUS: I'm training now a " << potentialtype << " potential on charges " << charges << ".");
     243
     244  /******************** TRAINING ********************/
     245  // fit potential
     246  FunctionModel *model =
     247      PotentialFactory::getInstance().createInstance(
     248          potentialtype,
     249          charges);
     250  ASSERT( model != NULL,
     251      "main() - model returned from PotentialFactory is NULL.");
     252  FunctionModel::parameters_t params(model->getParameterDimension(), 0.);
    216253  {
    217254    // then we ought to pick the right HomologyGraph ...
    218255    const HomologyGraph graph = getFirstGraphwithTimesSpecificElement(homologies,8,1);
    219256    if (graph != HomologyGraph()) {
    220       LOG(1, "First representative graph containing three saturated carbons is " << graph << ".");
     257      LOG(1, "First representative graph containing charges is " << graph << ".");
    221258
    222259      // Afterwards we go through all of this type and gather the distance and the energy value
    223       PairPotential_Angle::ParticleTypes_t types =
    224             boost::assign::list_of<PairPotential_Angle::ParticleType_t>
    225             (8)(1)(1)
    226           ;
    227       PairPotential_Angle angle(types);
    228       Fragment::charges_t h2o;
    229       h2o += 8,1,1;
    230       TrainingData AngleData(angle.getFragmentSpecificExtractor(h2o));
    231       AngleData(homologies.getHomologousGraphs(graph));
    232       LOG(1, "INFO: I gathered the following training data: " << AngleData);
     260      TrainingData data(model->getFragmentSpecificExtractor(charges));
     261      data(homologies.getHomologousGraphs(graph));
     262      LOG(1, "INFO: I gathered the following training data:\n" <<
     263          _detail::writeDistanceEnergyTable(data.getDistanceEnergyTable()));
    233264      // NOTICE that distance are in bohrradi as they come from MPQC!
    234265
    235266      // now perform the function approximation by optimizing the model function
    236       FunctionModel &model = angle;
    237       FunctionApproximation approximator(AngleData, model);
    238       if (model.isBoxConstraint() && approximator.checkParameterDerivatives()) {
     267      FunctionApproximation approximator(data, *model);
     268      if (model->isBoxConstraint() && approximator.checkParameterDerivatives()) {
    239269        // we set parameters here because we want to test with default ones
    240270        srand((unsigned)time(0)); // seed with current time
    241         model.setParametersToRandomInitialValues(AngleData);
    242         LOG(0, "INFO: Initial parameters are " << angle.getParameters() << ".");
     271        model->setParametersToRandomInitialValues(data);
     272        LOG(0, "INFO: Initial parameters are " << model->getParameters() << ".");
    243273        approximator(FunctionApproximation::ParameterDerivative);
    244274      } else {
     
    247277      }
    248278
    249       LOG(0, "RESULT: " << angle << ".");
    250 
    251       angleparams = model.getParameters();
     279//      LOG(0, "RESULT: " << *model << ".");
     280//
     281      params = model->getParameters();
     282
    252283    }
    253284  }
    254 
    255   /******************** MORSE TRAINING ********************/
    256   FunctionModel::parameters_t morseparams(PairPotential_Morse::MAXPARAMS, 0.);
    257   {
    258     // then we ought to pick the right HomologyGraph ...
    259     const HomologyGraph graph = getFirstGraphwithTimesSpecificElement(homologies,8,1);
    260     if (graph != HomologyGraph()) {
    261       LOG(1, "First representative graph containing two saturated carbons is " << graph << ".");
    262 
    263       // Afterwards we go through all of this type and gather the distance and the energy value
    264       Fragment::charges_t h2o;
    265       h2o += 8,1;
    266       PairPotential_Morse::ParticleTypes_t types =
    267             boost::assign::list_of<PairPotential_Morse::ParticleType_t>
    268             (8)(1)
    269           ;
    270       PairPotential_Morse morse(types);
    271       TrainingData MorseData(morse.getFragmentSpecificExtractor(h2o));
    272       MorseData(homologies.getHomologousGraphs(graph));
    273       LOG(1, "INFO: I gathered the following training data: " << MorseData);
    274       // NOTICE that distance are in bohrradi as they come from MPQC!
    275 
    276       // now perform the function approximation by optimizing the model function
    277       FunctionModel &model = morse;
    278       FunctionApproximation approximator(MorseData, model); // we only give CC distance, hence 1 input dim
    279       if (model.isBoxConstraint() && approximator.checkParameterDerivatives()) {
    280         // we set parameters here because we want to test with default ones
    281         srand((unsigned)time(0)); // seed with current time
    282         model.setParametersToRandomInitialValues(MorseData);
    283         LOG(0, "INFO: Initial parameters are " << morse.getParameters() << ".");
    284         approximator(FunctionApproximation::ParameterDerivative);
    285       } else {
    286         ELOG(0, "We require parameter derivatives for a box constraint minimization.");
    287         return 1;
    288       }
    289 
    290       LOG(0, "RESULT: " << morse << ".");
    291 
    292       morseparams = model.getParameters();
    293     }
    294   }
    295 
    296   /******************* SATURATION TRAINING *******************/
    297   FunctionModel::parameters_t params(SaturationPotential::MAXPARAMS, 0.);
    298   {
    299     // then we ought to pick the right HomologyGraph ...
    300     const HomologyGraph graph = getFirstGraphwithTimesSpecificElement(homologies,8,1);
    301     if (graph != HomologyGraph()) {
    302       LOG(1, "First representative graph containing one saturated carbon is " << graph << ".");
    303 
    304       // Afterwards we go through all of this type and gather the distance and the energy value
    305       Fragment::charges_t h2o;
    306       h2o += 8,1,1;
    307       SaturationPotential::ParticleTypes_t types =
    308             boost::assign::list_of<SaturationPotential::ParticleType_t>
    309             (8)(1)
    310           ;
    311       SaturationPotential saturation(types);
    312       TrainingData TersoffData(saturation.getFragmentSpecificExtractor(h2o));
    313       boost::function< std::vector<FunctionModel::arguments_t>(const argument_t &, const double)> triplefunction =
    314           boost::bind(&getTripleFromArgument, boost::cref(TersoffData.getTrainingInputs()), _1, _2);
    315       saturation.setTriplefunction(triplefunction);
    316 
    317       TersoffData( homologies.getHomologousGraphs(graph) );
    318       LOG(1, "INFO: I gathered the following training data: " << TersoffData);
    319       // NOTICE that distance are in bohrradi as they come from MPQC!
    320 
    321       // now perform the function approximation by optimizing the model function
    322       FunctionModel &model = saturation;
    323       FunctionApproximation approximator(TersoffData, model); // CH4 has 5 atoms, hence 5*4/2 distances
    324       if (model.isBoxConstraint() && approximator.checkParameterDerivatives()) {
    325         // we set parameters here because we want to test with default ones
    326         srand((unsigned)time(0)); // seed with current time
    327         model.setParametersToRandomInitialValues(TersoffData);
    328         approximator(FunctionApproximation::ParameterDerivative);
    329       } else {
    330         ELOG(0, "We require parameter derivatives for a box constraint minimization.");
    331         return 1;
    332       }
    333 
    334       LOG(0, "RESULT: " << saturation << ".");
    335 
    336       // check L2 and Lmax error against training set
    337       LOG(1, "INFO: L2sum = " << TersoffData.getL2Error(model)
    338           << ", LMax = " << TersoffData.getLMaxError(model) << ".");
    339     }
    340 
    341   }
     285  delete model;
     286  // remove static instances
     287  PotentialFactory::purgeInstance();
    342288
    343289  return 0;
    344290}
     291
Note: See TracChangeset for help on using the changeset viewer.