/* * Project: MoleCuilder * Description: creates and alters molecular systems * Copyright (C) 2013 University of Bonn. All rights reserved. * Copyright (C) 2013 Frederik Heber. All rights reserved. * * * This file is part of MoleCuilder. * * MoleCuilder is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 2 of the License, or * (at your option) any later version. * * MoleCuilder is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with MoleCuilder. If not, see . */ /* * FitPotentialAction.cpp * * Created on: Apr 09, 2013 * Author: heber */ // include config.h #ifdef HAVE_CONFIG_H #include #endif // needs to come before MemDebug due to placement new #include #include "CodePatterns/MemDebug.hpp" #include #include #include #include #include #include #include #include "Actions/FragmentationAction/FitPotentialAction.hpp" #include "CodePatterns/Log.hpp" #include "Element/element.hpp" #include "Fragmentation/Homology/HomologyContainer.hpp" #include "Fragmentation/Homology/HomologyGraph.hpp" #include "Fragmentation/Summation/SetValues/Fragment.hpp" #include "FunctionApproximation/Extractors.hpp" #include "FunctionApproximation/FunctionApproximation.hpp" #include "FunctionApproximation/FunctionModel.hpp" #include "FunctionApproximation/TrainingData.hpp" #include "FunctionApproximation/writeDistanceEnergyTable.hpp" #include "Potentials/CompoundPotential.hpp" #include "Potentials/PotentialDeserializer.hpp" #include "Potentials/PotentialFactory.hpp" #include "Potentials/PotentialRegistry.hpp" #include "Potentials/PotentialSerializer.hpp" #include "Potentials/SerializablePotential.hpp" using namespace MoleCuilder; // and construct the stuff #include "FitPotentialAction.def" #include "Action_impl_pre.hpp" /** =========== define the function ====================== */ HomologyGraph getFirstGraphwithSpecifiedElements( const HomologyContainer &homologies, const SerializablePotential::ParticleTypes_t &types) { ASSERT( !types.empty(), "getFirstGraphwithSpecifiedElements() - charges is empty?"); // create charges Fragment::charges_t charges; charges.resize(types.size()); std::transform(types.begin(), types.end(), charges.begin(), boost::lambda::_1); // convert into count map Extractors::elementcounts_t counts_per_charge = Extractors::_detail::getElementCounts(charges); ASSERT( !counts_per_charge.empty(), "getFirstGraphwithSpecifiedElements() - charge counts are empty?"); LOG(2, "DEBUG: counts_per_charge is " << counts_per_charge << "."); // we want to check each (unique) key only once HomologyContainer::const_key_iterator olditer = homologies.key_end(); for (HomologyContainer::const_key_iterator iter = homologies.key_begin(); iter != homologies.key_end(); olditer = iter++) { // if it's the same as the old one, skip it if (*olditer == *iter) continue; // if it's a new key, check if every element has the right number of counts Extractors::elementcounts_t::const_iterator countiter = counts_per_charge.begin(); for (; countiter != counts_per_charge.end(); ++countiter) if (!(*iter).hasTimesAtomicNumber( static_cast(countiter->first), static_cast(countiter->second)) ) break; if( countiter == counts_per_charge.end()) return *iter; } return HomologyGraph(); } Action::state_ptr FragmentationFitPotentialAction::performCall() { // fragment specifies the homology fragment to use SerializablePotential::ParticleTypes_t fragmentnumbers; { const std::vector &fragment = params.fragment.get(); std::transform(fragment.begin(), fragment.end(), std::back_inserter(fragmentnumbers), boost::bind(&element::getAtomicNumber, _1)); } // either charges and a potential is specified or a file if (boost::filesystem::exists(params.potential_file.get())) { std::ifstream returnstream(params.potential_file.get().string().c_str()); if (returnstream.good()) { PotentialDeserializer deserialize(returnstream); deserialize(); } else { ELOG(0, "Failed to parse from " << params.potential_file.get().string() << "."); return Action::failure; } returnstream.close(); LOG(0, "STATUS: I'm training now a set of potentials parsed from " << params.potential_file.get().string() << " on a fragment " << fragmentnumbers << " on data from " << params.homology_file.get() << "."); } else { if (params.charges.get().empty()) { ELOG(1, "Neither charges nor potential file given!"); return Action::failure; } else { // charges specify the potential type SerializablePotential::ParticleTypes_t chargenumbers; { const std::vector &charges = params.charges.get(); std::transform(charges.begin(), charges.end(), std::back_inserter(chargenumbers), boost::bind(&element::getAtomicNumber, _1)); } LOG(0, "STATUS: I'm training now a " << params.potentialtype.get() << " potential on charges " << chargenumbers << " on data from " << params.homology_file.get() << "."); // register desired potential and an additional constant one EmpiricalPotential *potential = PotentialFactory::getInstance().createInstance( params.potentialtype.get(), chargenumbers); PotentialRegistry::getInstance().registerInstance(potential); EmpiricalPotential *constant = PotentialFactory::getInstance().createInstance( std::string("constant"), SerializablePotential::ParticleTypes_t()); PotentialRegistry::getInstance().registerInstance(constant); } } // parse homologies into container HomologyContainer homologies; if (boost::filesystem::exists(params.homology_file.get())) { std::ifstream returnstream(params.homology_file.get().string().c_str()); if (returnstream.good()) { boost::archive::text_iarchive ia(returnstream); ia >> homologies; } else { ELOG(0, "Failed to parse from " << params.homology_file.get().string() << "."); return Action::failure; } returnstream.close(); } else { ELOG(0, params.homology_file.get() << " does not exist."); return Action::failure; } // first we try to look into the HomologyContainer LOG(1, "INFO: Listing all present homologies ..."); for (HomologyContainer::container_t::const_iterator iter = homologies.begin(); iter != homologies.end(); ++iter) { LOG(1, "INFO: graph " << iter->first << " has Fragment " << iter->second.first << " and associated energy " << iter->second.second << "."); } // then we ought to pick the right HomologyGraph ... const HomologyGraph graph = getFirstGraphwithSpecifiedElements(homologies,fragmentnumbers); if (graph != HomologyGraph()) { LOG(1, "First representative graph containing fragment " << fragmentnumbers << " is " << graph << "."); } else { ELOG(1, "Specific fragment " << fragmentnumbers << " not found in homologies!"); return Action::failure; } // fit potential FunctionModel *model = new CompoundPotential(graph); ASSERT( model != NULL, "FragmentationFitPotentialAction::performCall() - model is NULL."); /******************** TRAINING ********************/ // fit potential FunctionModel::parameters_t bestparams(model->getParameterDimension(), 0.); { // Afterwards we go through all of this type and gather the distance and the energy value TrainingData data(model->getFragmentSpecificExtractor()); data(homologies.getHomologousGraphs(graph)); // print distances and energies if desired for debugging if (!data.getTrainingInputs().empty()) { // print which distance is which size_t counter=1; if (DoLog(3)) { const FunctionModel::arguments_t &inputs = data.getTrainingInputs()[0]; for (FunctionModel::arguments_t::const_iterator iter = inputs.begin(); iter != inputs.end(); ++iter) { const argument_t &arg = *iter; LOG(3, "DEBUG: distance " << counter++ << " is between (#" << arg.indices.first << "c" << arg.types.first << "," << arg.indices.second << "c" << arg.types.second << ")."); } } // print table LOG(3, "DEBUG: I gathered the following training data:\n" << _detail::writeDistanceEnergyTable(data.getDistanceEnergyTable())); } // now perform the function approximation by optimizing the model function FunctionApproximation approximator(data, *model); if (model->isBoxConstraint() && approximator.checkParameterDerivatives()) { double l2error = std::numeric_limits::infinity(); // seed with current time srand((unsigned)time(0)); for (unsigned int runs=0; runs < params.best_of_howmany.get(); ++runs) { // generate new random initial parameter values model->setParametersToRandomInitialValues(data); LOG(1, "INFO: Initial parameters of run " << runs << " are " << model->getParameters() << "."); approximator(FunctionApproximation::ParameterDerivative); LOG(1, "INFO: Final parameters of run " << runs << " are " << model->getParameters() << "."); const double new_l2error = data.getL2Error(*model); if (new_l2error < l2error) { // store currently best parameters l2error = new_l2error; bestparams = model->getParameters(); LOG(1, "STATUS: New fit from run " << runs << " has better error of " << l2error << "."); } } // reset parameters from best fit model->setParameters(bestparams); LOG(1, "INFO: Best parameters with L2 error of " << l2error << " are " << model->getParameters() << "."); } else { ELOG(0, "We require parameter derivatives for a box constraint minimization."); return Action::failure; } // create a map of each fragment with error. typedef std::multimap< double, size_t > WorseFragmentMap_t; WorseFragmentMap_t WorseFragmentMap; HomologyContainer::range_t fragmentrange = homologies.getHomologousGraphs(graph); // fragments make it into the container in reversed order, hence count from top down size_t index= std::distance(fragmentrange.first, fragmentrange.second)-1; for (HomologyContainer::const_iterator iter = fragmentrange.first; iter != fragmentrange.second; ++iter) { const Fragment& fragment = iter->second.first; const double &energy = iter->second.second; // create arguments from the fragment FunctionModel::extractor_t extractor = model->getFragmentSpecificExtractor(); FunctionModel::arguments_t args = extractor(fragment, 1); // calculate value from potential const double fitvalue = (*model)(args)[0]; // insert difference into map const double error = fabs(energy - fitvalue); WorseFragmentMap.insert( std::make_pair( error, index-- ) ); { // give only the distances in the debugging text std::stringstream streamargs; BOOST_FOREACH (argument_t arg, args) { streamargs << " " << arg.distance*AtomicLengthToAngstroem; } LOG(2, "DEBUG: frag.#" << index+1 << "'s error is |" << energy << " - " << fitvalue << "| = " << error << " for args " << streamargs.str() << "."); } } LOG(0, "RESULT: WorstFragmentMap " << WorseFragmentMap << "."); // print fitted potentials std::stringstream potentials; PotentialSerializer serialize(potentials); serialize(); LOG(1, "STATUS: Resulting parameters are " << std::endl << potentials.str()); } delete model; return Action::success; } Action::state_ptr FragmentationFitPotentialAction::performUndo(Action::state_ptr _state) { return Action::success; } Action::state_ptr FragmentationFitPotentialAction::performRedo(Action::state_ptr _state){ return Action::success; } bool FragmentationFitPotentialAction::canUndo() { return false; } bool FragmentationFitPotentialAction::shouldUndo() { return false; } /** =========== end of function ====================== */