Skip to content
Snippets Groups Projects
Commit e7f2790b authored by michel's avatar michel
Browse files

new minimizer interface

parent ade26459
No related branches found
No related tags found
No related merge requests found
Showing
with 1204 additions and 39 deletions
......@@ -2,6 +2,8 @@
#include <iostream>
#include <cmath>
#include <sstream>
#include <vector>
#include <string>
// Boost header files go here
#include <boost/lexical_cast.hpp>
......@@ -20,53 +22,49 @@
* The main function.
*/
int main(int argc, char **argv){
double p0, p1, p2, p3, sigma_smear;
p0 = -10. ; p1 = 10.; p2 = 1. ; p3 = -0.01; sigma_smear = 3;
std::string whichMinimizer("all");
double p0=-10., p1=10., p2=1., p3=-0.01, sigma_smear=3;
// Generate data distribution
boost::shared_ptr<PolyFit> myFit(new PolyFit(p0, p1, p2, p3, sigma_smear));
boost::shared_ptr<MIData> myFit(new PolyFit(p0, p1, p2, p3, sigma_smear));
//--------------------------Minimizer IF --------------------------------------------------------
std::vector<boost::shared_ptr<MIBase> > myMinimizerList;
// Add minimizers
if (whichMinimizer=="Geneva") myMinimizerList.push_back(boost::shared_ptr<MIBase> (new MIGeneva(myFit)));
else if (whichMinimizer=="Minuit") myMinimizerList.push_back(boost::shared_ptr<MIBase> (new MIMinuit(myFit)));
else if (whichMinimizer=="all") {
myMinimizerList.push_back(boost::shared_ptr<MIBase> (new MIGeneva(myFit)));
myMinimizerList.push_back(boost::shared_ptr<MIBase> (new MIMinuit(myFit)));
}else{
std::cout << "Minimizer/t" << whichMinimizer << "\tdoesn't exist" << std::endl;
return 0;
}
// Initiate parameters
double val[4], min[4], max[4], err[4];
val[0] = -11; max[0] = 0; min[0] = -20; err[0] = 3;
val[1] = 9.8; max[1] = 15; min[1] = 5; err[1] = 2;
val[2] = 1.1; max[2] = 1.5; min[2] = 0.5; err[2] = 0.3;
val[3] = -0.008; max[3] = 0.; min[3] = -0.02; err[3] = 0.005;
MIGeneva genevaFit(myFit);
double genResult = genevaFit.exec(4, val, min, max, err);
std::cout << "final Geneva par :\t" << genResult << std::endl;
std::cout << "final a:\t" << val[0] << " +- " << err[0] << std::endl;
std::cout << "final b:\t" << val[1] << " +- " << err[1] << std::endl;
std::cout << "final c:\t" << val[2] << " +- " << err[2] << std::endl;
std::cout << "final d:\t" << val[3] << " +- " << err[3] << std::endl;
//myFit->drawGraph(resultPar.at(0),resultPar.at(1),resultPar.at(2),resultPar.at(3));
std::cout << "Done Geneva ..." << std::endl << std::endl;
//val[0] = -11; max[0] = 0; min[0] = -20; err[0] = 3;
//val[1] = 9.8; max[1] = 15; min[1] = 5; err[1] = 2;
//val[2] = 1.1; max[2] = 1.5; min[2] = 0.5; err[2] = 0.3;
//val[3] = -0.008; max[3] = 0.; min[3] = -0.02; err[3] = 0.005;
MIMinuit minuitFit(myFit);
double minResult = minuitFit.exec(4, val, min, max, err);
std::cout << "final Minuit par :\t" << minResult << std::endl;
std::cout << "final a:\t" << val[0] << " +- " << err[0] << std::endl;
std::cout << "final b:\t" << val[1] << " +- " << err[1] << std::endl;
std::cout << "final c:\t" << val[2] << " +- " << err[2] << std::endl;
std::cout << "final d:\t" << val[3] << " +- " << err[3] << std::endl;
std::cout<< "Done Minuit ..." << std::endl << std::endl;
myFit->drawGraph(val[0],val[1],val[2],val[3]);
// Loop over minimizers (at the moment this means: Geneva, Minuit or Geneva then Minuit)
for(unsigned int Nmin=0; Nmin<myMinimizerList.size(); Nmin++){
// Pointer to one ot the used minimizers
boost::shared_ptr<MIBase> minimizer = myMinimizerList[Nmin];
// Do the actual minimization
double genResult = minimizer->exec(4, val, min, max, err);
std::cout << "Minimizer " << Nmin << "\t final par :\t" << genResult << std::endl;
std::cout << "final a:\t" << val[0] << " +- " << err[0] << std::endl;
std::cout << "final b:\t" << val[1] << " +- " << err[1] << std::endl;
std::cout << "final c:\t" << val[2] << " +- " << err[2] << std::endl;
std::cout << "final d:\t" << val[3] << " +- " << err[3] << std::endl;
std::cout << "Done ..." << std::endl << std::endl;
}
// Plot results
//myFit->drawGraph(val[0],val[1],val[2],val[3]);
return 0;
}
#################################################################################
# Configuration file for GStartProject #
# #
# Authors: Ruediger Berlich #
#################################################################################
# The amount of random number producer threads
nProducerThreads = 8
# The amount of threads processing individuals simultaneously
nEvaluationThreads = 16
# The size of the population
populationSize = 100
# The number of parents in the population
nParents = 2
# Maximum number of iterations in the population
maxIterations = 10000
# The maximum number of minutes the optimization of the population should run (0 means "no limit")
maxMinutes = 0
# The number of iterations after which information should be emitted in the super-population
reportIteration = 1000
# The recombination scheme for the super-population (DEFAULTRECOMBINE=0, RANDOMRECOMBINE=1, VALUERECOMBINE=2)
rScheme = 2
# Determines whether sorting happens in MUPLUSNU (0), MUCOMMANU (1) or MUNU1PRETAIN (2) mode
sortingScheme = 1
# The size of the buffer with random arrays in the random factory
arraySize = 1000
# Determines whether information about configuration options should be emitted
verbose = 1
# The maximum number of cycles a client should perform adaptions before it returns without success
processingCycles = 1
# Specifies whether results should be returned even if they are not better than before
returnRegardless = 1
# Influences the maximum waiting time of the GBrokerEA after the arrival of the first evaluated individual
waitFactor = 2
\ No newline at end of file
build-project MinimizerInterface ;
project :
;
lib FitIF :
[ glob *.cc : *App.cc ]
MinimizerInterface//MinimizerInterface
$(TOP)/qft++//qft++
$(TOP)/ErrLogger//ErrLogger
: <use>$(TOP)//GemCom <use>$(TOP)//GemCou <use>$(TOP)//GemInd <use>$(TOP)//GemGen <use>$(TOP)//GemHap <use>$(TOP)//Minuit2
:
: <library>$(TOP)//GemCom <library>$(TOP)//GemCou <library>$(TOP)//GemInd <library>$(TOP)//GemGen <library>$(TOP)//GemHap <library>$(TOP)//Minuit2 ;
exe ToyFitApp : ToyFitApp.cc FitIF : ;
/**
* @file GArgumentParser.cpp
*/
/*
* Copyright (C) Gemfony scientific UG (haftungsbeschraenkt) and Karlsruhe
* Institute of Technology (University of the State of Baden-Wuerttemberg
* and National Laboratory of the Helmholtz Association).
*
* See the AUTHORS file in the top-level directory for a list of authors.
*
* Contact: contact [at] gemfony (dot) com
*
* This file is part of the Geneva library collection
*
* Geneva is free software: you can redistribute it and/or modify
* it under the terms of version 3 of the GNU Affero General Public License
* as published by the Free Software Foundation.
*
* Geneva 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 Affero General Public License for more details.
*
* You should have received a copy of the GNU Affero General Public License
* along with the Geneva library. If not, see <http://www.gnu.org/licenses/>.
*
* For further information on Gemfony scientific and Geneva, visit
* http://www.gemfony.com .
*/
#include "Examples/Tutorial/FitIF/MinimizerInterface/GArgumentParser.hh"
namespace Gem
{
namespace Geneva
{
/************************************************************************************************/
/**
* A function that parses the command line for all required parameters
*/
bool parseCommandLine(
int argc, char **argv
, std::string& configFile
, boost::uint16_t& parallelizationMode
, bool& serverMode
, std::string& ip
, unsigned short& port
, Gem::Common::serializationMode& serMode
) {
try{
// Check the command line options. Uses the Boost program options library.
po::options_description desc("Usage: evaluator [options]");
desc.add_options()
("help,h", "emit help message")
("configFile,c", po::value<std::string>(&configFile)->default_value(DEFAULTCONFIGFILE),
"The name of the configuration file holding further configuration options")
("parallelizationMode,p", po::value<boost::uint16_t>(&parallelizationMode)->default_value(DEFAULTPARALLELIZATIONMODE),
"Whether or not to run this optimization in serial mode (0), multi-threaded (1) or networked (2) mode")
("serverMode,s","Whether to run networked execution in server or client mode. The option only gets evaluated if \"--parallelizationMode=2\"")
("ip",po::value<std::string>(&ip)->default_value(DEFAULTIP), "The ip of the server")
("port",po::value<unsigned short>(&port)->default_value(DEFAULTPORT), "The port of the server")
("serMode", po::value<Gem::Common::serializationMode>(&serMode)->default_value(DEFAULTSERMODE),
"Specifies whether serialization shall be done in TEXTMODE (0), XMLMODE (1) or BINARYMODE (2)")
;
po::variables_map vm;
po::store(po::parse_command_line(argc, argv, desc), vm);
po::notify(vm);
// Emit a help message, if necessary
if (vm.count("help")) {
std::cerr << desc << std::endl;
return false;
}
serverMode=false;
if (vm.count("parallelizationMode")) {
if(parallelizationMode > 2) {
std::cout << "Error: the \"-p\" or \"--parallelizationMode\" option may only assume the"<< std::endl
<< "values 0 (serial), 1 (multi-threaded) or 2 (networked). Leaving ..." << std::endl;
return false;
}
if(parallelizationMode == 2) if(vm.count("serverMode")) serverMode = true;
}
if(parallelizationMode != DEFAULTPARALLELIZATIONMODE || ip != DEFAULTIP || port != DEFAULTPORT){
std::string parModeString;
switch(parallelizationMode) {
case 0:
parModeString = "serial";
break;
case 1:
parModeString = "multi-threaded";
break;
case 2:
parModeString = "networked";
break;
};
std::cout << std::endl
<< "Running with the following command line options:" << std::endl
<< "configFile = " << configFile << std::endl
<< "parallelizationMode = " << parModeString << std::endl
<< "serverMode = " << (serverMode?"true":"false") << std::endl
<< "ip = " << ip << std::endl
<< "port = " << port << std::endl
<< "serMode = " << serMode << std::endl
<< std::endl;
}
}
catch(...){
std::cout << "Error parsing the command line" << std::endl;
return false;
}
return true;
}
/************************************************************************************************/
/**
* A function that parses a config file for further parameters
*/
bool parseConfigFile(
const std::string& configFile
, boost::uint16_t& nProducerThreads
, boost::uint16_t& nEvaluationThreads
, std::size_t& populationSize
, std::size_t& nParents
, boost::uint32_t& maxIterations
, long& maxMinutes
, boost::uint32_t& reportIteration
, recoScheme& rScheme
, sortingMode& smode
, std::size_t& arraySize
, boost::uint32_t& processingCycles
, bool& returnRegardless
, boost::uint32_t& waitFactor
) {
boost::uint16_t recombinationScheme=0;
bool verbose;
// Check the name of the configuation file
if(configFile.empty() || configFile == "empty" || configFile == "unknown") {
std::cerr << "Error: Invalid configuration file name given: \"" << configFile << "\"" << std::endl;
return false;
}
try{
// Check the configuration file line options. Uses the Boost program options library.
po::options_description config("Allowed options");
config.add_options()
("nProducerThreads",po::value<boost::uint16_t>(&nProducerThreads)->default_value(DEFAULTNPRODUCERTHREADS),
"The amount of random number producer threads")
("nEvaluationThreads",po::value<boost::uint16_t>(&nEvaluationThreads)->default_value(DEFAULTNEVALUATIONTHREADS),
"The amount of threads processing individuals simultaneously")
("populationSize",po::value<std::size_t>(&populationSize)->default_value(DEFAULTPOPULATIONSIZE),
"The size of the super-population")
("nParents",po::value<std::size_t>(&nParents)->default_value(DEFAULTNPARENTS),
"The number of parents in the population") // Needs to be treated separately
("maxIterations", po::value<boost::uint32_t>(&maxIterations)->default_value(DEFAULTMAXITERATIONS),
"Maximum number of iterations in the population")
("maxMinutes", po::value<long>(&maxMinutes)->default_value(DEFAULTMAXMINUTES),
"The maximum number of minutes the optimization of the population should run")
("reportIteration",po::value<boost::uint32_t>(&reportIteration)->default_value(DEFAULTREPORTITERATION),
"The number of iterations after which information should be emitted in the super-population")
("rScheme",po::value<boost::uint16_t>(&recombinationScheme)->default_value(DEFAULTRSCHEME),
"The recombination scheme for the super-population")
("sortingScheme,o", po::value<sortingMode>(&smode)->default_value(DEFAULTSORTINGSCHEME),
"Determines whether sorting is done in MUCOMMANU (0), MUPLUSNU (1) or MUNU1PRETAIN (2) mode")
("arraySize", po::value<std::size_t>(&arraySize)->default_value(DEFAULTARRAYSIZE),
"The size of the buffer with random arrays in the random factory")
("verbose",po::value<bool>(&verbose)->default_value(DEFAULTVERBOSE),
"Whether additional information should be emitted")
("processingCycles", po::value<boost::uint32_t>(&processingCycles)->default_value(DEFAULTPROCESSINGCYCLES),
"The maximum number of cycles a client should perform adaptions before it returns without success")
("returnRegardless", po::value<bool>(&returnRegardless)->default_value(DEFAULTRETURNREGARDLESS),
"Specifies whether results should be returned even if they are not better than before")
("waitFactor", po::value<boost::uint32_t>(&waitFactor)->default_value(DEFAULTGBTCWAITFACTOR),
"Influences the maximum waiting time of the GBrokerEA after the arrival of the first evaluated individuum")
;
po::variables_map vm;
std::ifstream ifs(configFile.c_str());
if(!ifs.good()) {
std::cerr << "Error accessing configuration file " << configFile;
return false;
}
store(po::parse_config_file(ifs, config), vm);
notify(vm);
// Emit a help message, if necessary
if (vm.count("help")) {
std::cout << config << std::endl;
return false;
}
// Check the number of parents in the super-population
if(2*nParents > populationSize){
std::cout << "Error: Invalid number of parents inpopulation" << std::endl
<< "nParents = " << nParents << std::endl
<< "populationSize = " << populationSize << std::endl;
return false;
}
// Workaround for assigment problem with rScheme
if(recombinationScheme==(boost::uint16_t)VALUERECOMBINE)
rScheme=VALUERECOMBINE;
else if(recombinationScheme==(boost::uint16_t)RANDOMRECOMBINE)
rScheme=RANDOMRECOMBINE;
else if(recombinationScheme==(boost::uint16_t)DEFAULTRECOMBINE)
rScheme=DEFAULTRECOMBINE;
else {
std::cout << "Error: Invalid recombination scheme in population: " << recombinationScheme << std::endl;
return false;
}
if(verbose){
std::cout << std::endl
<< "Running with the following options from " << configFile << ":" << std::endl
<< "nProducerThreads = " << (boost::uint16_t)nProducerThreads << std::endl // boost::uint8_t not printable on gcc ???
<< "populationSize = " << populationSize << std::endl
<< "nParents = " << nParents << std::endl
<< "maxIterations = " << maxIterations << std::endl
<< "maxMinutes = " << maxMinutes << std::endl
<< "reportIteration = " << reportIteration << std::endl
<< "rScheme = " << (boost::uint16_t)rScheme << std::endl
<< "sortingScheme = " << smode << std::endl
<< "arraySize = " << arraySize << std::endl
<< "processingCycles = " << processingCycles << std::endl
<< "returnRegardless = " << (returnRegardless?"true":"false") << std::endl
<< "waitFactor = " << waitFactor << std::endl
<< std::endl;
}
}
catch(...){
std::cout << "Error parsing the configuration file " << configFile << std::endl;
return false;
}
return true;
}
/************************************************************************************************/
} /* namespace Geneva */
} /* namespace Gem */
/**
* @file GArgumentParser.hpp
*/
/*
* Copyright (C) Gemfony scientific UG (haftungsbeschraenkt) and Karlsruhe
* Institute of Technology (University of the State of Baden-Wuerttemberg
* and National Laboratory of the Helmholtz Association).
*
* See the AUTHORS file in the top-level directory for a list of authors.
*
* Contact: contact [at] gemfony (dot) com
*
* This file is part of the Geneva library collection
*
* Geneva is free software: you can redistribute it and/or modify
* it under the terms of version 3 of the GNU Affero General Public License
* as published by the Free Software Foundation.
*
* Geneva 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 Affero General Public License for more details.
*
* You should have received a copy of the GNU Affero General Public License
* along with the Geneva library. If not, see <http://www.gnu.org/licenses/>.
*
* For further information on Gemfony scientific and Geneva, visit
* http://www.gemfony.com .
*/
// Standard headers go here
#include <iostream>
#include <fstream>
#include <sstream>
#include <string>
#include <vector>
// Boost headers go here
#include <boost/version.hpp>
#if BOOST_VERSION < 103600
#error "Error: Boost should at least have version 1.36 !"
#endif /* BOOST_VERSION */
#include <boost/program_options.hpp>
#include <boost/filesystem.hpp>
#ifndef GARGUMENTPARSER_HPP_
#define GARGUMENTPARSER_HPP_
// For Microsoft-compatible compilers
#if defined(_MSC_VER) && (_MSC_VER >= 1020)
#pragma once
#endif
// Geneva headers go here
#include <common/GCommonEnums.hpp>
#include <common/GSerializationHelperFunctionsT.hpp>
#include <geneva/GOptimizationEnums.hpp>
namespace Gem
{
namespace Geneva
{
// Default settings
const boost::uint16_t DEFAULTNPRODUCERTHREADS=10;
const boost::uint16_t DEFAULTNEVALUATIONTHREADS=4;
const std::size_t DEFAULTPOPULATIONSIZE=100;
const std::size_t DEFAULTNPARENTS=5; // Allow to explore the parameter space from many starting points
const boost::uint32_t DEFAULTMAXITERATIONS=2000;
const long DEFAULTMAXMINUTES=10;
const boost::uint32_t DEFAULTREPORTITERATION=1;
const recoScheme DEFAULTRSCHEME=VALUERECOMBINE;
const bool DEFAULTVERBOSE=true;
const bool DEFAULTPARALLELIZATIONMODE=1;
const std::size_t DEFAULTARRAYSIZE=1000;
const bool DEFAULTPRODUCTIONPLACE=true; // local production
const bool DEFAULTUSECOMMONADAPTOR=false; // whether to use a common adaptor for all GParameterT objects
const unsigned short DEFAULTPORT=10000;
const std::string DEFAULTIP="localhost";
const std::string DEFAULTCONFIGFILE="./GStartProject.cfg";
const sortingMode DEFAULTSORTINGSCHEME=MUPLUSNU;
const boost::uint32_t DEFAULTSTARTITERATION=0;
const boost::uint32_t DEFAULTPROCESSINGCYCLES=1;
const bool DEFAULTRETURNREGARDLESS=true;
const std::size_t DEFAULTNBTCONSUMERTHREADS=2;
const boost::uint32_t DEFAULTGBTCWAITFACTOR=5;
const Gem::Common::serializationMode DEFAULTSERMODE=Gem::Common::SERIALIZATIONMODE_TEXT;
namespace po = boost::program_options;
bool parseCommandLine(
int argc, char **argv
, std::string& configFile
, boost::uint16_t& parallelizationMode
, bool& serverMode
, std::string& ip
, unsigned short& port
, Gem::Common::serializationMode& serMode
);
bool parseConfigFile(
const std::string& configFile
, boost::uint16_t& nProducerThreads
, boost::uint16_t& nEvaluationThreads
, std::size_t& populationSize
, std::size_t& nParents
, boost::uint32_t& maxIterations
, long& maxMinutes
, boost::uint32_t& reportIteration
, recoScheme& rScheme
, sortingMode& smode
, std::size_t& arraySize
, boost::uint32_t& processingCycles
, bool& returnRegardless
, boost::uint32_t& waitFactor
);
} /* namespace Geneva */
} /* namespace Gem */
#endif /* GARGUMENTPARSER_HPP_ */
//********************************************************************
// Holds parameter data in the form deemed necessary by the PWA team,
// plus information on how to evaluate the parameter set. It would be
// useful if this interface could be agreed upon by both the IHEP and
// PAWIAN teams in order to reduce the multiplicity. I would suggest
// to specify the parameter definition and access in the base class and
// to specify evaluators in derived classes. These can then be specific
// to PAWIAN and the IHEP code (i.e. interface to OpenCL/GPGPU).
class GenericInterface
{
public:
// Evaluate a given parameter set
virtual double evaluator() = 0;
// Retrieve a clone of a GenericInterface object-derivative
// Needed so that we can clone objects through a base-pointer
// without information on whether we are dealing with PAWIAN or GPGPU.
virtual boost::shared_ptr<GenericInterface> clone() = 0;
// some access code for adding, modifying and retrieving parameters
// ...
private:
std::vector<double> parameters_;
};
//********************************************************************
// GPGPU-Interface to PWA
class GPGPUPWA
:public GenericInterface
{
public:
GPGPUPWA(const GPGPUPWA&);
virtual double evaluator(); // GPGPU code for PWA analysis
virtual boost::shared_ptr<GenericInterface> clone(){
return boost::shared_ptr<GenericInterface>(new GPGPUPWA(*this));
}
};
//********************************************************************
// Pawian interface
class PawianInterface
:public GenericInterface
{
public:
PawianInterface(const PawianInterface&);
virtual double evaluator(); // Standard PAWIAN code for the evaluation
virtual boost::shared_ptr<GenericInterface> clone(){
return boost::shared_ptr<GenericInterface>(new PawianInterface(*this));
}
};
//********************************************************************
// This class creates an evaluator with the API required by MINUIT (with
// which I am not too familiar, so I leave it open for now). The constructor
// could e.g. copy the parameters contained in the GenericInterface object
// into the "TheMinuitEvaluatorBaseClassWhichIamNotFamiliarWith" object in
// the form Minuit expects and transform the data for the evaluator member
// of the GenericInterface. This implies some overhead. However, as the
// execution time of the evaluation will be large in comparison, this is not
// too much of a problem. The code would be generic and not even bound to
// the PWA use case, which is encoded only in GenericInterface::evaluator().
class MinuitInterface
:public TheMinuitEvaluatorBaseClassWhichIamNotFamiliarWith
{
public:
// Constructs the object from a generic evaluator
MinuitInterface(boost::shared_ptr<GenericInterface>);
// Gives access to a GenericInterface object from the converged solution
const GenericInterface& getResult() const;
private:
boost::shared_ptr<GenericInterface> genericEvaluator_;
};
//********************************************************************
// This class does the same as MinuitInterface, albeit for Geneva.
class GenevaInterface
:public Gem::Geneva::GParameterSet
{
public:
// Constructs the object from a generic evaluator
GenevaInterface(boost::shared_ptr<GenericInterface>);
// Gives access to a GenericInterface object from the converged solution
const GenericInterface& getResult() const;
private:
boost::shared_ptr<GenericInterface> genericEvaluator_;
};
//********************************************************************
// Same procedure as for the other interfaces. Not specified here
class ThirdPartyMinimizedInterface
{
public:
};
//********************************************************************
int main(int argc, char **argv) {
// Create a GenericInterface-derivative object and fill it with data
// ...
// Create objets suitable for minimization with Geneva, MINUIT2 or any
// Thirdparty-optimizer
// ...
// If required, extract the result of a first optimization e.g. with Geneva
// and create a new interface object for another optimizer, such as MINUIT2.
// Note that it would of course be possible to create wrapper code arount
// MINUIT2 and Geneva to give the *optimizer* a common interface. However,
// I believe that just creating copy constructors for the GenericInterface
// object outlined above would be more than sufficient in the beginning.
// In particular, creating wrapper code for the optimizers does not scale
// well if you want to add additional optimizers.
}
project :
;
lib MinimizerInterface :
[ glob *.cc ]
$(TOP)/qft++//qft++
$(TOP)/ErrLogger//ErrLogger
: <use>$(TOP)//Minuit2 <use>$(TOP)//Geneva
:
: <library>$(TOP)//Minuit2 <library>$(TOP)//Geneva ;
#ifndef MIGENERIC_HH_
#define MIGENERIC_HH_
#include <vector>
class MIGeneric
{
public:
MIGeneric()
{
}
virtual ~MIGeneric()
{ /* nothing */ }
// Evaluate a given parameter set
virtual double evaluator() = 0;
// Retrieve a clone of a GenericInterface object-derivative
// Needed so that we can clone objects through a base-pointer
// without information on whether we are dealing with PAWIAN or GPGPU.
virtual boost::shared_ptr<MIGeneric> clone() = 0;
// some access code for adding, modifying and retrieving parameters
virtual void setStartPar(int num, double *par, double* min, double* max, double* err) = 0;
virtual const double getPar(const unsigned int num) const {
if(num>=parValue_.size()) return 0;
return parValue_[num];
};
virtual const double getMin(const unsigned int num) const {
if(num>=minValue_.size()) return 0;
return minValue_[num];
};
virtual const double getMax(const unsigned int num) const {
if(num>=maxValue_.size()) return 0;
return maxValue_[num];
};
virtual const double getErr(const unsigned int num) const {
if(num>=parError_.size()) return 0;
return parError_[num];
};
virtual const unsigned int getNumPars() const {
return parValue_.size();
};
protected:
std::vector<double> parValue_;
std::vector<double> minValue_;
std::vector<double> maxValue_;
std::vector<double> parError_;
};
#endif
#include <vector>
#include <string>
#include <sstream>
#include <iostream>
#include "Examples/Tutorial/FitIF2/MinimizerInterface/MIGeneva.hh"
#include "Examples/Tutorial/FitIF2/MinimizerInterface/MIGeneric.hh"
using namespace std;
using namespace Gem::Geneva;
using namespace Gem::Courtier;
using namespace Gem::Hap;
MIGeneva::MIGeneva(boost::shared_ptr<MIGeneric> theData): GParameterSet(){
genericEvaluator_=theData;
// Add bounded double objects
for(std::size_t i=0; i<genericEvaluator_->getNumPars(); i++) {
boost::shared_ptr<GConstrainedDoubleObject> gbd_ptr(new GConstrainedDoubleObject(genericEvaluator_->getPar(i), genericEvaluator_->getMin(i), genericEvaluator_->getMax(i)) );
// Create a suitable adaptor (sigma=0.1, sigma-adaption=0.5, min sigma=0, max sigma=0,5)
boost::shared_ptr<GDoubleGaussAdaptor> gdga_ptr(new GDoubleGaussAdaptor(genericEvaluator_->getErr(i), 0.5, 0., 3*genericEvaluator_->getErr(i)));
gdga_ptr->setAdaptionThreshold(1); // Adaption parameters are modified after each adaption
gdga_ptr->setAdaptionProbability(0.05); // The likelihood for a parameter to be adapted
// Register the adaptor with GConstrainedDoubleObject objects
gbd_ptr->addAdaptor(gdga_ptr);
// Add a GConstrainedDoubleObject object to the collection
// gbdc_ptr->push_back(gbd_ptr);
// gpoc_ptr->push_back(gbd_ptr);
this->push_back(gbd_ptr);
}
}
MIGeneva::~MIGeneva()
{
//delete _myFcn;
}
double MIGeneva::fitnessCalculation(){
//double result = 0.;
// Extract the GConstrainedDoubleObjectCollection object. In a realistic scenario, you might want
// to add error checks here upon first invocation.
return genericEvaluator_->evaluator();
}
const MIGeneric& MIGeneva::getResult() const {
return *genericEvaluator_;// bestIndividual_ptr->fitnessCalculation();
}
#ifndef _MIGENEVA_H
#define _MIGENEVA_H
#include <vector>
#include <boost/shared_ptr.hpp>
#include "Examples/Tutorial/FitIF2/MinimizerInterface/MIGeneric.hh"
#include <boost/program_options.hpp>
#include <boost/filesystem.hpp>
// Geneva header files go here
#include <hap/GRandomT.hpp>
#include <common/GCommonEnums.hpp>
#include <common/GExceptions.hpp>
#include <geneva/GConstrainedDoubleObject.hpp>
#include <geneva/GConstrainedDoubleObjectCollection.hpp>
#include <geneva/GDoubleGaussAdaptor.hpp>
#include <geneva/GObjectExpectationChecksT.hpp>
#include <geneva/GParameterObjectCollection.hpp>
#include <geneva/GParameterSet.hpp>
#include <courtier/GAsioHelperFunctions.hpp>
#include <courtier/GAsioTCPClientT.hpp>
#include <courtier/GAsioTCPConsumerT.hpp>
#include <geneva/GParameterSet.hpp>
#include <geneva/GBrokerEA.hpp>
#include <geneva/GEvolutionaryAlgorithm.hpp>
#include <geneva/GIndividual.hpp>
#include <geneva/GMultiThreadedEA.hpp>
#include <common/GSerializationHelperFunctionsT.hpp>
#include <geneva/GOptimizationEnums.hpp>
using namespace std;
class MIGeneva : public Gem::Geneva::GParameterSet {
friend class boost::serialization::access;
template<typename Archive>
void serialize(Archive & ar, const unsigned int) {
using boost::serialization::make_nvp;
ar & BOOST_SERIALIZATION_BASE_OBJECT_NVP(GParameterSet);
/* Add your own class-variables here in the following way:
ar & BOOST_SERIALIZATION_NVP(myVar);
or
ar & make_nvp("myVar", myVar); // The latter form can be necessary when dealing with templates
*/
}
public:
// Constructs the object from a generic evaluator
MIGeneva(boost::shared_ptr<MIGeneric> theData);
// Gives access to a GenericInterface object from the converged solution
const MIGeneric& getResult() const;
/** Destructor */
virtual ~MIGeneva();
protected:
virtual double fitnessCalculation();
private:
boost::shared_ptr<MIGeneric> genericEvaluator_;
// boost::shared_ptr<MIData> _myData;
// vector<string> paramNames;
};
#endif /* _MIGeneva_H */
#include <getopt.h>
#include <fstream>
#include <sstream>
#include <iostream>
#include <string>
#include "Examples/Tutorial/FitIF2/PolyFit.hh"
#include "TFile.h"
#include "TMath.h"
#include "TGraph.h"
#include "TCanvas.h"
#include "TRandom.h"
#include "TF1.h"
#include "TGraphErrors.h"
#include "ErrLogger/ErrLogger.hh"
PolyFit::PolyFit(double p0, double p1, double p2, double p3, double sigma) :
_theTFile(0),
_sigma(sigma)
{
// Display parameters for test distribution
cout << endl;
Info <<"Set p0 as "<< p0 << endmsg;
Info <<"Set p1 as "<< p1 << endmsg;
Info <<"Set p2 as "<< p2 << endmsg;
Info <<"Set p3 as "<< p3 << endmsg;
Info <<"Set sigma as "<< sigma << endmsg;
// Generate test distribution and smear them with a gaussian
TRandom randomNumber;
for(int i=0; i<1000; i++) {
double tmpXvalue=static_cast<double>((rand() % 10000 + 1)) / 100;
double tmpYvalue=randomNumber.Gaus(p0 + p1 * tmpXvalue + p2 * tmpXvalue * tmpXvalue + p3 * tmpXvalue * tmpXvalue * tmpXvalue, _sigma);
_xValue.push_back(tmpXvalue);
_yValue.push_back(tmpYvalue);
}
}
PolyFit::PolyFit(const PolyFit& inFit) :
_theTFile(0)
{
_myGraph = inFit.getGraph();
_xValue = inFit.getX();
_yValue = inFit.getY();
_sigma = inFit.getSig();
}
double PolyFit::evaluator(){
// Calculate chi^2 for current set of fit parameters
double result=0.;
for (unsigned int i=0; i<_xValue.size(); i++){
double yValFit=parValue_[0]+parValue_[1]*_xValue[i]+parValue_[2]*_xValue[i]*_xValue[i]+parValue_[3]*_xValue[i]*_xValue[i]*_xValue[i];
double yValExp=_yValue[i];
double tmpChi=((yValExp-yValFit)*(yValExp-yValFit))/(_sigma*_sigma);
result+=tmpChi;
}
return result;
}
void PolyFit::setStartPar(int num, double *par, double* min, double* max, double* err){
for(int i=0; i<num ; i++){
parValue_.push_back(par[i]);
minValue_.push_back(min[i]);
maxValue_.push_back(max[i]);
parError_.push_back(err[i]);
}
}
void PolyFit::drawGraph(double a, double b, double c, double d){
// Create arrays (which will be needed for feeding the TGraph)
const unsigned int dataEntries=_xValue.size();
double x[dataEntries];
double y[dataEntries];
double xerr[dataEntries];
double yerr[dataEntries];
// Convert vectors to arrays
for(unsigned int i=0; i<dataEntries; i++) {
x[i]=_xValue[i];
y[i]=_yValue[i];
xerr[i] = 0;
yerr[i] = _sigma;
}
// Create ROOT file
_theTFile = new TFile("myFit.root","RECREATE");
// Create TGraph and write to ROOT file
TCanvas* c1=new TCanvas("c1","c1",1200,800);
c1->cd();
TGraphErrors* linDist = new TGraphErrors(100, x, y, xerr, yerr);
linDist->SetMarkerStyle(6);
linDist->Draw("AP");
linDist->Write("myGraph");
// Draw fit with calculated parameters
TF1* fit1 = new TF1("fit1","pol3",0.,100.);
fit1->SetParameters(a, b, c, d);
fit1->SetLineColor(46);
fit1->SetLineWidth(3);
fit1->Draw();
fit1->Write();
}
PolyFit::~PolyFit()
{
_theTFile->Write();
_theTFile->Close();
}
#ifndef _PolyFit_H
#define _PolyFit_H
#include <iostream>
#include <fstream>
#include <string>
#include <vector>
#include <map>
#include <cassert>
#include <boost/shared_ptr.hpp>
#include "TROOT.h"
#include "qft++/topincludes/relativistic-quantum-mechanics.hh"
#include "Examples/Tutorial/FitIF2/MinimizerInterface/MIGeneric.hh"
class TFile;
class TGraph;
class TCanvas;
class TRandom;
class PolyFit : public MIGeneric {
public:
// create/copy/destroy:
///Constructor
PolyFit(double p0, double p1, double p2, double p3, double sigma);
PolyFit(const PolyFit&);
virtual void setStartPar(int num, double *par, double* min, double* max, double* err);
virtual double evaluator(); // Standard PAWIAN code for the evaluation
virtual boost::shared_ptr<MIGeneric> clone(){
return boost::shared_ptr<MIGeneric>(new PolyFit(*this));
};
/** Destructor */
virtual ~PolyFit();
double controlParameter(const std::vector<double>& minPar);
void drawGraph(double a, double b, double c, double d);
// Getters:
std::map <unsigned int, TGraph* > getGraph() const { return _myGraph;};
std::vector< double > getX() const { return _xValue;};
std::vector< double > getY() const { return _yValue;};
const double getSig() const { return _sigma;};
protected:
private:
TFile* _theTFile;
std::map <unsigned int, TGraph* > _myGraph;
std::vector< double > _xValue;
std::vector< double > _yValue;
double _sigma;
};
#endif
// Standard header files go here
#include <iostream>
#include <cmath>
#include <sstream>
#include <vector>
#include <string>
// Boost header files go here
#include <boost/lexical_cast.hpp>
//#include "ErrLogger/ErrLogger.hh"
// Minimizer Interface header files go here
#include "Examples/Tutorial/FitIF2/MinimizerInterface/MIGeneric.hh"
#include "Examples/Tutorial/FitIF2/MinimizerInterface/MIGeneva.hh"
#include "Examples/Tutorial/FitIF2/MinimizerInterface/GArgumentParser.hh"
// Geneva header files go here
#include <hap/GRandomT.hpp>
#include <common/GCommonEnums.hpp>
#include <common/GExceptions.hpp>
#include <geneva/GConstrainedDoubleObject.hpp>
#include <geneva/GConstrainedDoubleObjectCollection.hpp>
#include <geneva/GDoubleGaussAdaptor.hpp>
#include <geneva/GObjectExpectationChecksT.hpp>
#include <geneva/GParameterObjectCollection.hpp>
#include <geneva/GParameterSet.hpp>
#include <courtier/GAsioHelperFunctions.hpp>
#include <courtier/GAsioTCPClientT.hpp>
#include <courtier/GAsioTCPConsumerT.hpp>
#include <geneva/GParameterSet.hpp>
#include <geneva/GBrokerEA.hpp>
#include <geneva/GEvolutionaryAlgorithm.hpp>
#include <geneva/GIndividual.hpp>
#include <geneva/GMultiThreadedEA.hpp>
#include <common/GSerializationHelperFunctionsT.hpp>
#include <geneva/GOptimizationEnums.hpp>
// The toy-data to fit to
#include "Examples/Tutorial/FitIF2/PolyFit.hh"
using namespace Gem::Geneva;
using namespace Gem::Courtier;
using namespace Gem::Hap;
/************************************************************************************************/
/**
* The main function.
*/
int main(int argc, char **argv){
//std::string whichMinimizer("Geneva");
double p0=-10., p1=10., p2=1., p3=-0.01, sigma_smear=3;
// Generate data distribution
boost::shared_ptr<MIGeneric> myFit(new PolyFit(p0, p1, p2, p3, sigma_smear));
// Initiate parameters
double val[4], min[4], max[4], err[4];
val[0] = -11; max[0] = 0; min[0] = -20; err[0] = 3;
val[1] = 9.8; max[1] = 15; min[1] = 5; err[1] = 2;
val[2] = 1.1; max[2] = 1.5; min[2] = 0.5; err[2] = 0.3;
val[3] = -0.008; max[3] = 0.; min[3] = -0.02; err[3] = 0.005;
myFit->setStartPar(4, val, min, max, err);
//--------------------------Minimizer IF --------------------------------------------------------
//boost::shared_ptr<MIGeneva> minGen(new MIGeneva(myFit));
std::string configFile;
boost::uint16_t parallelizationMode;
bool serverMode;
std::string ip;
unsigned short port;
Gem::Common::serializationMode serMode;
boost::uint16_t nProducerThreads;
boost::uint16_t nEvaluationThreads;
std::size_t populationSize;
std::size_t nParents;
boost::uint32_t maxIterations;
long maxMinutes;
boost::uint32_t reportIteration;
Gem::Geneva::recoScheme rScheme;
std::size_t arraySize;
Gem::Geneva::sortingMode smode;
boost::uint32_t processingCycles;
bool returnRegardless;
boost::uint32_t waitFactor;
configFile="./GStartProject.cfg"; parallelizationMode=1; serverMode=false; ip="localhost"; port=10000; serMode=Gem::Common::SERIALIZATIONMODE_TEXT;
bool parsedConfig = parseConfigFile(configFile,
nProducerThreads,
nEvaluationThreads,
populationSize,
nParents,
maxIterations,
maxMinutes,
reportIteration,
rScheme,
smode,
arraySize,
processingCycles,
returnRegardless,
waitFactor);
if(!parsedConfig) cout << "TODO" << endl;
////--------------------------------------------------------
// Random numbers are our most valuable good. Set the number of threads
GRANDOMFACTORY->setNProducerThreads(nProducerThreads);
GRANDOMFACTORY->setArraySize(arraySize);
//***************************************************************************
// If this is a client in networked mode, we can just start the listener and
// return when it has finished
if(parallelizationMode==2 && !serverMode) {
boost::shared_ptr<GAsioTCPClientT<GIndividual> > p(new GAsioTCPClientT<GIndividual>(ip, boost::lexical_cast<std::string>(port)));
p->setMaxStalls(0); // An infinite number of stalled data retrievals
p->setMaxConnectionAttempts(100); // Up to 100 failed connection attempts
// Prevent return of unsuccessful adaption attempts to the server
p->returnResultIfUnsuccessful(returnRegardless);
// Start the actual processing loop
p->run();
return 0;
}
//***************************************************************************
std::cout << "Create Parents ..." << std::endl << std::endl;
// Create the first set of parent individuals. Initialization of parameters is (should be) done randomly.
std::vector<boost::shared_ptr<MIGeneva> > parentIndividuals;
for(std::size_t p = 0 ; p<nParents; p++) {
//TODO: vary start parameters
boost::shared_ptr<MIGeneva> gdii_ptr(new MIGeneva(myFit));
gdii_ptr->setProcessingCycles(processingCycles);
parentIndividuals.push_back(gdii_ptr);
}
std::cout << "Parent created ..." << std::endl << std::endl;
//***************************************************************************
// We can now start creating populations. We refer to them through the base class
// This smart pointer will hold the different population types
boost::shared_ptr<GEvolutionaryAlgorithm> pop_ptr;
// Create the actual populations
switch (parallelizationMode) {
//-----------------------------------------------------------------------------------------------------
case 0: // Serial execution
// Create an empty population
pop_ptr = boost::shared_ptr<GEvolutionaryAlgorithm>(new GEvolutionaryAlgorithm());
break;
//-----------------------------------------------------------------------------------------------------
case 1: // Multi-threaded execution
{
// Create the multi-threaded population
boost::shared_ptr<GMultiThreadedEA> popPar_ptr(new GMultiThreadedEA());
// Population-specific settings
popPar_ptr->setNThreads(nEvaluationThreads);
// Assignment to the base pointer
pop_ptr = popPar_ptr;
}
break;
//-----------------------------------------------------------------------------------------------------
case 2: // Networked execution (server-side)
{
// Create a network consumer and enrol it with the broker
boost::shared_ptr<GAsioTCPConsumerT<GIndividual> > gatc(new GAsioTCPConsumerT<GIndividual>(port));
gatc->setSerializationMode(serMode);
GINDIVIDUALBROKER->enrol(gatc);
// Create the actual broker population
boost::shared_ptr<GBrokerEA> popBroker_ptr(new GBrokerEA());
popBroker_ptr->setWaitFactor(waitFactor);
// Assignment to the base pointer
pop_ptr = popBroker_ptr;
}
break;
}
///////////////////////////////////////////////////////////////////////////////////////////////////////////////
// Now we have suitable populations and can fill them with data
// Add individuals to the population
for(std::size_t p = 0 ; p<nParents; p++) {
pop_ptr->push_back(parentIndividuals[p]);
}
// Specify some general population settings
pop_ptr->setDefaultPopulationSize(populationSize,nParents);
pop_ptr->setMaxIteration(maxIterations);
pop_ptr->setMaxTime(boost::posix_time::minutes(maxMinutes));
pop_ptr->setReportIteration(reportIteration);
pop_ptr->setRecombinationMethod(rScheme);
pop_ptr->setSortingScheme(smode);
// Do the actual optimization
pop_ptr->optimize();
//--------------------------------------------------------------------------------------------
//boost::shared_ptr<MIGeneva> bestIndividual_ptr=pop_ptr->getBestIndividual<MIGeneva>();
////--------------------------------------------------------
//plot result
//const MIGeneric& bestResult= bestIndividual_ptr->getResult();
//std::cout << "Minimizer Geneva \t final par :\t" << bestResult.evaluator() << std::endl;
//for(int i=0; i<bestResult.getNumPars(); i++){
// std::cout << "final Par" << i <<": \t" << bestResult.getPar(i) << " +- " << bestResult.getErr(i) << std::endl;
//}
/*std::cout << "Minimizer " << Nmin << "\t final par :\t" << genResult << std::endl;
std::cout << "final a:\t" << val[0] << " +- " << err[0] << std::endl;
std::cout << "final b:\t" << val[1] << " +- " << err[1] << std::endl;
std::cout << "final c:\t" << val[2] << " +- " << err[2] << std::endl;
std::cout << "final d:\t" << val[3] << " +- " << err[3] << std::endl;
std::cout << "Done ..." << std::endl << std::endl;*/
//}
// Plot results
//myFit->drawGraph(val[0],val[1],val[2],val[3]);
return 0;
}
File added
build-project DfuncClebschG ;
build-project MinuitFit ;
build-project FitIF ;
\ No newline at end of file
build-project FitIF ;
build-project FitIF2 ;
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment