[98d166] | 1 | /*
|
---|
| 2 | * Project: MoleCuilder
|
---|
| 3 | * Description: creates and alters molecular systems
|
---|
| 4 | * Copyright (C) 2014 Frederik Heber. All rights reserved.
|
---|
| 5 | *
|
---|
| 6 | *
|
---|
| 7 | * This file is part of MoleCuilder.
|
---|
| 8 | *
|
---|
| 9 | * MoleCuilder is free software: you can redistribute it and/or modify
|
---|
| 10 | * it under the terms of the GNU General Public License as published by
|
---|
| 11 | * the Free Software Foundation, either version 2 of the License, or
|
---|
| 12 | * (at your option) any later version.
|
---|
| 13 | *
|
---|
| 14 | * MoleCuilder is distributed in the hope that it will be useful,
|
---|
| 15 | * but WITHOUT ANY WARRANTY; without even the implied warranty of
|
---|
| 16 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
---|
| 17 | * GNU General Public License for more details.
|
---|
| 18 | *
|
---|
| 19 | * You should have received a copy of the GNU General Public License
|
---|
| 20 | * along with MoleCuilder. If not, see <http://www.gnu.org/licenses/>.
|
---|
| 21 | */
|
---|
| 22 |
|
---|
| 23 | /*
|
---|
| 24 | * PotentialTrainer.cpp
|
---|
| 25 | *
|
---|
| 26 | * Created on: Sep 11, 2014
|
---|
| 27 | * Author: heber
|
---|
| 28 | */
|
---|
| 29 |
|
---|
| 30 | // include config.h
|
---|
| 31 | #ifdef HAVE_CONFIG_H
|
---|
| 32 | #include <config.h>
|
---|
| 33 | #endif
|
---|
| 34 |
|
---|
| 35 | // needs to come before MemDebug due to placement new
|
---|
| 36 | #include <boost/archive/text_iarchive.hpp>
|
---|
| 37 |
|
---|
| 38 | #include "CodePatterns/MemDebug.hpp"
|
---|
| 39 |
|
---|
| 40 | #include "PotentialTrainer.hpp"
|
---|
| 41 |
|
---|
| 42 | #include <algorithm>
|
---|
| 43 | #include <boost/lambda/lambda.hpp>
|
---|
| 44 | #include <boost/filesystem.hpp>
|
---|
| 45 | #include <fstream>
|
---|
| 46 | #include <sstream>
|
---|
| 47 |
|
---|
| 48 | #include "CodePatterns/Assert.hpp"
|
---|
| 49 | #include "CodePatterns/Log.hpp"
|
---|
| 50 |
|
---|
| 51 | #include "Element/element.hpp"
|
---|
| 52 | #include "Fragmentation/Homology/HomologyContainer.hpp"
|
---|
| 53 | #include "Fragmentation/Homology/HomologyGraph.hpp"
|
---|
| 54 | #include "FunctionApproximation/Extractors.hpp"
|
---|
| 55 | #include "FunctionApproximation/FunctionApproximation.hpp"
|
---|
| 56 | #include "FunctionApproximation/FunctionModel.hpp"
|
---|
| 57 | #include "FunctionApproximation/TrainingData.hpp"
|
---|
| 58 | #include "FunctionApproximation/writeDistanceEnergyTable.hpp"
|
---|
| 59 | #include "Potentials/CompoundPotential.hpp"
|
---|
[fde8e7] | 60 | #include "Potentials/RegistrySerializer.hpp"
|
---|
[98d166] | 61 | #include "Potentials/SerializablePotential.hpp"
|
---|
| 62 |
|
---|
| 63 | PotentialTrainer::PotentialTrainer()
|
---|
| 64 | {}
|
---|
| 65 |
|
---|
| 66 | PotentialTrainer::~PotentialTrainer()
|
---|
| 67 | {}
|
---|
| 68 |
|
---|
| 69 | bool PotentialTrainer::operator()(
|
---|
| 70 | const HomologyContainer &_homologies,
|
---|
| 71 | const HomologyGraph &_graph,
|
---|
| 72 | const boost::filesystem::path &_trainingfile,
|
---|
[b40690] | 73 | const unsigned int _maxiterations,
|
---|
[98d166] | 74 | const double _threshold,
|
---|
| 75 | const unsigned int _best_of_howmany) const
|
---|
| 76 | {
|
---|
| 77 | // fit potential
|
---|
| 78 | FunctionModel *model = new CompoundPotential(_graph);
|
---|
| 79 | ASSERT( model != NULL,
|
---|
| 80 | "PotentialTrainer::operator() - model is NULL.");
|
---|
| 81 |
|
---|
| 82 | /******************** TRAINING ********************/
|
---|
| 83 | // fit potential
|
---|
| 84 | FunctionModel::parameters_t bestparams(model->getParameterDimension(), 0.);
|
---|
| 85 | {
|
---|
| 86 | // Afterwards we go through all of this type and gather the distance and the energy value
|
---|
| 87 | TrainingData data(model->getSpecificFilter());
|
---|
| 88 | data(_homologies.getHomologousGraphs(_graph));
|
---|
| 89 |
|
---|
| 90 | // print distances and energies if desired for debugging
|
---|
| 91 | if (!data.getTrainingInputs().empty()) {
|
---|
| 92 | // print which distance is which
|
---|
| 93 | size_t counter=1;
|
---|
| 94 | if (DoLog(3)) {
|
---|
| 95 | const FunctionModel::arguments_t &inputs = data.getAllArguments()[0];
|
---|
| 96 | for (FunctionModel::arguments_t::const_iterator iter = inputs.begin();
|
---|
| 97 | iter != inputs.end(); ++iter) {
|
---|
| 98 | const argument_t &arg = *iter;
|
---|
| 99 | LOG(3, "DEBUG: distance " << counter++ << " is between (#"
|
---|
| 100 | << arg.indices.first << "c" << arg.types.first << ","
|
---|
| 101 | << arg.indices.second << "c" << arg.types.second << ").");
|
---|
| 102 | }
|
---|
| 103 | }
|
---|
| 104 |
|
---|
| 105 | // print table
|
---|
| 106 | if (_trainingfile.string().empty()) {
|
---|
| 107 | LOG(3, "DEBUG: I gathered the following training data:\n" <<
|
---|
| 108 | _detail::writeDistanceEnergyTable(data.getDistanceEnergyTable()));
|
---|
| 109 | } else {
|
---|
| 110 | std::ofstream trainingstream(_trainingfile.string().c_str());
|
---|
| 111 | if (trainingstream.good()) {
|
---|
| 112 | LOG(3, "DEBUG: Writing training data to file " <<
|
---|
| 113 | _trainingfile.string() << ".");
|
---|
| 114 | trainingstream << _detail::writeDistanceEnergyTable(data.getDistanceEnergyTable());
|
---|
| 115 | }
|
---|
| 116 | trainingstream.close();
|
---|
| 117 | }
|
---|
| 118 | }
|
---|
| 119 |
|
---|
| 120 | if ((_threshold < 1.) && (_best_of_howmany))
|
---|
| 121 | ELOG(2, "threshold parameter always overrules max_runs, both are specified.");
|
---|
| 122 | // now perform the function approximation by optimizing the model function
|
---|
[b40690] | 123 | FunctionApproximation approximator(data, *model, _threshold, _maxiterations);
|
---|
[98d166] | 124 | if (model->isBoxConstraint() && approximator.checkParameterDerivatives()) {
|
---|
| 125 | double l2error = std::numeric_limits<double>::max();
|
---|
| 126 | // seed with current time
|
---|
| 127 | srand((unsigned)time(0));
|
---|
| 128 | unsigned int runs=0;
|
---|
| 129 | // threshold overrules max_runs
|
---|
| 130 | const double threshold = _threshold;
|
---|
| 131 | const unsigned int max_runs = (threshold >= 1.) ? _best_of_howmany : 1;
|
---|
| 132 | LOG(1, "INFO: Maximum runs is " << max_runs << " and threshold set to " << threshold << ".");
|
---|
| 133 | do {
|
---|
| 134 | // generate new random initial parameter values
|
---|
| 135 | model->setParametersToRandomInitialValues(data);
|
---|
| 136 | LOG(1, "INFO: Initial parameters of run " << runs << " are "
|
---|
| 137 | << model->getParameters() << ".");
|
---|
| 138 | approximator(FunctionApproximation::ParameterDerivative);
|
---|
| 139 | LOG(1, "INFO: Final parameters of run " << runs << " are "
|
---|
| 140 | << model->getParameters() << ".");
|
---|
| 141 | const double new_l2error = data.getL2Error(*model);
|
---|
| 142 | if (new_l2error < l2error) {
|
---|
| 143 | // store currently best parameters
|
---|
| 144 | l2error = new_l2error;
|
---|
| 145 | bestparams = model->getParameters();
|
---|
| 146 | LOG(1, "STATUS: New fit from run " << runs
|
---|
| 147 | << " has better error of " << l2error << ".");
|
---|
| 148 | }
|
---|
| 149 | } while (( ++runs < max_runs) || (l2error > threshold));
|
---|
| 150 | // reset parameters from best fit
|
---|
| 151 | model->setParameters(bestparams);
|
---|
| 152 | LOG(1, "INFO: Best parameters with L2 error of "
|
---|
| 153 | << l2error << " are " << model->getParameters() << ".");
|
---|
| 154 | } else {
|
---|
| 155 | return false;
|
---|
| 156 | }
|
---|
| 157 |
|
---|
| 158 | // create a map of each fragment with error.
|
---|
| 159 | HomologyContainer::range_t fragmentrange = _homologies.getHomologousGraphs(_graph);
|
---|
| 160 | TrainingData::L2ErrorConfigurationIndexMap_t WorseFragmentMap =
|
---|
| 161 | data.getWorstFragmentMap(*model, fragmentrange);
|
---|
| 162 | LOG(0, "RESULT: WorstFragmentMap " << WorseFragmentMap << ".");
|
---|
| 163 |
|
---|
| 164 | }
|
---|
| 165 | delete model;
|
---|
| 166 |
|
---|
| 167 | return true;
|
---|
| 168 | }
|
---|
| 169 |
|
---|
| 170 | HomologyGraph PotentialTrainer::getFirstGraphwithSpecifiedElements(
|
---|
| 171 | const HomologyContainer &homologies,
|
---|
| 172 | const SerializablePotential::ParticleTypes_t &types)
|
---|
| 173 | {
|
---|
| 174 | ASSERT( !types.empty(),
|
---|
| 175 | "getFirstGraphwithSpecifiedElements() - charges is empty?");
|
---|
| 176 | // create charges
|
---|
| 177 | Fragment::charges_t charges;
|
---|
| 178 | charges.resize(types.size());
|
---|
| 179 | std::transform(types.begin(), types.end(),
|
---|
| 180 | charges.begin(), boost::lambda::_1);
|
---|
| 181 | // convert into count map
|
---|
| 182 | Extractors::elementcounts_t counts_per_charge =
|
---|
| 183 | Extractors::_detail::getElementCounts(charges);
|
---|
| 184 | ASSERT( !counts_per_charge.empty(),
|
---|
| 185 | "getFirstGraphwithSpecifiedElements() - charge counts are empty?");
|
---|
| 186 | LOG(2, "DEBUG: counts_per_charge is " << counts_per_charge << ".");
|
---|
| 187 | // we want to check each (unique) key only once
|
---|
| 188 | HomologyContainer::const_key_iterator olditer = homologies.key_end();
|
---|
| 189 | for (HomologyContainer::const_key_iterator iter =
|
---|
[e63edb] | 190 | homologies.key_begin(); iter != homologies.key_end();
|
---|
| 191 | iter = homologies.getNextKey(iter)) {
|
---|
[98d166] | 192 | // if it's the same as the old one, skip it
|
---|
[e63edb] | 193 | if (olditer == iter)
|
---|
[98d166] | 194 | continue;
|
---|
[e63edb] | 195 | else
|
---|
| 196 | olditer = iter;
|
---|
[98d166] | 197 | // if it's a new key, check if every element has the right number of counts
|
---|
| 198 | Extractors::elementcounts_t::const_iterator countiter = counts_per_charge.begin();
|
---|
| 199 | for (; countiter != counts_per_charge.end(); ++countiter)
|
---|
| 200 | if (!(*iter).hasTimesAtomicNumber(
|
---|
| 201 | static_cast<size_t>(countiter->first),
|
---|
| 202 | static_cast<size_t>(countiter->second))
|
---|
| 203 | )
|
---|
| 204 | break;
|
---|
| 205 | if( countiter == counts_per_charge.end())
|
---|
| 206 | return *iter;
|
---|
| 207 | }
|
---|
| 208 | return HomologyGraph();
|
---|
| 209 | }
|
---|
| 210 |
|
---|
| 211 | SerializablePotential::ParticleTypes_t PotentialTrainer::getNumbersFromElements(
|
---|
| 212 | const std::vector<const element *> &fragment)
|
---|
| 213 | {
|
---|
| 214 | SerializablePotential::ParticleTypes_t fragmentnumbers;
|
---|
| 215 | std::transform(fragment.begin(), fragment.end(), std::back_inserter(fragmentnumbers),
|
---|
| 216 | boost::bind(&element::getAtomicNumber, _1));
|
---|
| 217 | return fragmentnumbers;
|
---|
| 218 | }
|
---|