Commit 6b3cbdf8 authored by Dmitry Bagaev's avatar Dmitry Bagaev
Browse files

TTSP Optimizer alpha M.1

parent 49bc48d4
......@@ -28,15 +28,15 @@ int main(int argc, char **argv) {
{
std::string seriesFileName = "";
std::string seriesDirectory = "";
std::string seriesFileName = "";
std::string seriesDirectory = "";
std::string parametersFileName = "";
std::string solverName = "fcbiilu2";
std::string optimizerType = "bruteforce";
std::string solverName = "fcbiilu2";
std::string optimizerType = "bruteforce";
bool seriesFound = false;
bool seriesFound = false;
bool parametersFound = false;
bool waitNext = false;
bool waitNext = false;
//Parse argv parameters
if (argc == 1) goto helpMessage;
......@@ -60,12 +60,12 @@ int main(int argc, char **argv) {
std::cout << " Available solvers:" << std::endl;
Solver::Initialize(NULL, NULL, NULL);
std::vector<std::string> availableSolvers = Solver::getAvailableSolvers();
for (auto it = availableSolvers.begin(); it != availableSolvers.end(); ++it) {
for (auto it = availableSolvers.begin(); it != availableSolvers.end(); ++it) {
std::cout << " " << *it << std::endl;
}
std::cout << " Available optimizers:" << std::endl;
std::vector<std::string> availableOptimizers = TTSP::OptimizerInterface::getAvailableOptimizers();
for (auto it = availableOptimizers.begin(); it != availableOptimizers.end(); ++it) {
for (auto it = availableOptimizers.begin(); it != availableOptimizers.end(); ++it) {
std::cout << " " << *it << std::endl;
}
Solver::Finalize();
......@@ -74,8 +74,8 @@ int main(int argc, char **argv) {
}
//Series file name found with -s or --series options
if (strcmp(argv[i], "-s") == 0 || strcmp(argv[i], "--series") == 0) {
seriesFound = true;
seriesFileName = std::string(argv[i + 1]);
seriesFound = true;
seriesFileName = std::string(argv[i + 1]);
FILE *seriesFile = fopen(seriesFileName.c_str(), "r");
if (seriesFile == NULL) {
if (rank == 0) {
......@@ -105,7 +105,7 @@ int main(int argc, char **argv) {
if (rank == 0) {
std::cout << "Solver parameters file found: " << argv[i + 1] << std::endl;
}
parametersFound = true;
parametersFound = true;
parametersFileName = std::string(argv[i + 1]);
i++;
continue;
......@@ -165,26 +165,33 @@ int main(int argc, char **argv) {
Solver solver = Solver(solverName, "test");
solver.SetVerbosityLevel(SolverVerbosityLevel::Level0);
if (rank == 0) std::cout << "Solving with " << solverName << std::endl;
TTSP::OptimizationParameter tau("tau", {3e-2, 5e-2, 6e-2, 7e-2, 8e-2, 9e-2, 1e-1, 2e-1, 3e-1, 5e-1, 7e-1, 9e-1}, 1e-3);
TTSP::OptimizationParameter tau("tau", {3e-2, 5e-2, 6e-2, 7e-2, 8e-2, 9e-2, 1e-1, 2e-1, 3e-1, 5e-1, 7e-1, 9e-1}, 1e-3);
TTSP::OptimizationParameters parameters;
parameters.push_back(std::make_pair(tau, 1e-3));
TTSP::OptimizerProperties properties;
properties["tau:use_closest"] = "false";
TTSP::OptimizationParametersSpace space(solverName, "test", parameters);
TTSP::OptimizerInterface *optimizer = TTSP::OptimizerInterface::getOptimizer(optimizerType, space);
TTSP::OptimizerInterface *optimizer = TTSP::OptimizerInterface::getOptimizer(optimizerType, space, properties, 10);
if (optimizer == nullptr) {
if (rank == 0) {
std::cout << "Optimizer " << optimizerType << " not found" << std::endl;
std::cout << " Available optimizers:" << std::endl;
std::vector<std::string> availableOptimizers = TTSP::OptimizerInterface::getAvailableOptimizers();
for (auto it = availableOptimizers.begin(); it != availableOptimizers.end(); ++it) {
for (auto it = availableOptimizers.begin(); it != availableOptimizers.end(); ++it) {
std::cout << " " << *it << std::endl;
}
}
std::exit(0);
}
optimizer->SetVerbosityLevel(TTSP::OptimizerVerbosityLevel::Level1);
while (!series.end()) {
......@@ -211,43 +218,10 @@ int main(int argc, char **argv) {
matrix.Load(next.first);
if (rank == 0) std::cout << "Solving with A = " << next.first << " and b = " << next.second << std::endl;
// if (rank == 0) std::cout << "Solving with A = " << next.first << " and b = " << next.second << std::endl;
optimizer->Solve(solver, matrix, rhs, x);
if (rank == 0) {
std::cout << std::endl << "Next optimization parameters found for current iteration:" << std::endl;
const TTSP::OptimizationParameterPoints &best = optimizer->GetSpace().GetPoints();
std::for_each(best.begin(), best.end(), [](const TTSP::OptimizationParameterPoint &p) {
std::cout << "\t" << p.first << " = " << p.second << std::endl;
});
std::cout << std::endl << "Optimization results buffer output:" << std::endl;
const TTSP::OptimizationParameterResultsBuffer &results = optimizer->GetResults();
int index = 1;
std::for_each(results.begin(), results.end(), [&index](const TTSP::OptimizationParameterResult &result) {
std::cout << "\t" << index++ << "\t" << " [";
const TTSP::OptimizationParameterPoints &points = result.GetPoints();
std::for_each(points.begin(), points.end(), [](const TTSP::OptimizationParameterPoint &point) {
std::cout << " " << point.first << "=" << point.second << " ";
});
std::cout << "] " << result.GetPreconditionerTime() << "\t" << result.GetSolveTime() << "\t" << result.GetTime() << std::endl;
});
}
if (rank == 0) {
std::cout << std::endl
<< "Solved with " << solver.SolverName()
<< " on " << solver.Iterations()
<< " iterations and " << solver.Residual()
<< " residual. Reason: " << solver.ReturnReason()
<< std::endl;
}
if (rank == 0 && waitNext) {
std::cin.get();
}
......
......@@ -11,6 +11,12 @@ namespace INMOST
class SolverInterface;
class SolverParameters;
enum SolverVerbosityLevel {
Level0 = 0,
Level1 = 1,
Level2 = 2
};
/// Main class to set and solve linear system.
/// Solver class is used to set the coefficient Matrix, the right-hand side Vector
/// and the initial guess Vector, construct the preconditioner and Solve
......@@ -45,7 +51,7 @@ namespace INMOST
/// @param value Value for the parameter.
/// @see Solver::SetParameter
void SetParameterEnum(std::string name, INMOST_DATA_ENUM_TYPE value);
/// Backward-compatibility access to real-typed parameters.
/// Backward-compatibility access to real-typed parameters.
/// Check Solver::SetParameter for availible parameter names.
/// @param name Name of the parameter.
/// @param value Value for the parameter.
......@@ -62,6 +68,10 @@ namespace INMOST
static bool is_finalized; ///< Indicator that solvers were finalized for MPI operations.
SolverInterface *solver; ///< Pointer to abstract linear solver.
std::string prefix; ///< Prescribed solver name that is used to assign solver parameters.
SolverVerbosityLevel versbosity;
double preconditioner_time;
double iterations_time;
bool is_solved;
/// Reads the parameters from database file stored in xml format.
/// @param xml_database Path to xml file with solver parameters.
static void parseXMLDatabase(const char* xml_database);
......@@ -83,7 +93,7 @@ namespace INMOST
std::vector<INMOST_DATA_REAL_TYPE> recv_storage; ///< Storage used to receive data of the vector.
std::vector<INMOST_MPI_Request> send_requests; ///< Sturctures used to wait complition of send operations.
std::vector<INMOST_MPI_Request> recv_requests; ///< Sturctures used to detect complition of receive operations.
std::vector<INMOST_DATA_ENUM_TYPE> extended_indexes; ///< All the indices that appear outside of the local
std::vector<INMOST_DATA_ENUM_TYPE> extended_indexes; ///< All the indices that appear outside of the local
INMOST_DATA_ENUM_TYPE local_vector_begin; ///< Remember the first index in expanded vector.
INMOST_DATA_ENUM_TYPE local_vector_end; ///< Remember the last index in expanded vector.
INMOST_DATA_ENUM_TYPE initial_matrix_begin; ///< Initial the first index of the matrix before overlapping was computed.
......@@ -148,10 +158,10 @@ namespace INMOST
INMOST_DATA_ENUM_TYPE GetSize() const { return size; }
/// Update the shared data in parallel vector.
/// @param x Vector for which the shared data should be updated.
void Update(Sparse::Vector &x);
void Update(Sparse::Vector &x);
/// Sum shared values in parallel vector.
/// @param x Vector for which the shared data should be accumulated.
void Accumulate(Sparse::Vector &x);
void Accumulate(Sparse::Vector &x);
/// Get the sum of num elements of real array on all processes.
/// @param inout Data that should be integrated.
/// @param num Number of entries in inout array.
......@@ -198,12 +208,17 @@ namespace INMOST
/// Assign a solver.
/// \warning Not all solvers support assignment. This operation may be very expensive.
Solver &operator=(const Solver &other);
/// Return the solver name
/// Return the solver name
/// @see Sparse::Solve
std::string SolverName() const;
/// Return the solver user specified name of the current solver
/// @see Sparse::Solve
std::string SolverPrefix() const;
/// Return the solver verbosity specified level of the current solver
/// @see Sparse::Solve
SolverVerbosityLevel VerbosityLevel() const;
///
void SetVerbosityLevel(SolverVerbosityLevel level);
/// Initialize the stage of parallel solution.
/// If MPI is not initialized yet, then it will be initialized.
///
......@@ -238,11 +253,11 @@ namespace INMOST
static bool isFinalized();
/// Set the matrix and construct the preconditioner.
/// @param A Matrix A in linear problem Ax = b
/// @param ModifiedPattern Indicates whether the structure of the matrix have
/// @param ModifiedPattern Indicates whether the structure of the matrix have
/// changed since last call to Solver::SetMatrix.
/// @param OldPreconditioner If this parameter is set to true,
/// then the previous preconditioner will be used,
/// otherwise the new preconditioner will be constructed.
/// otherwise the new preconditioner will be constructed.
///
/// Preconditioner will be constructed on call to this function
/// - for INNER_*, PETSc and ANI packages
......@@ -348,6 +363,18 @@ namespace INMOST
/// Get the reason of convergence or divergence of the last solution.
/// @see Sparse::Solve
const std::string ReturnReason() const;
/// Return the time of preconditioner evaluation by the last solution
/// @see Sparse::Solve
double PreconditionerTime() const;
/// Return the time of iteration stage by the last solution
/// @see Sparse::Solve
double IterationsTime() const;
/// Return the last solution successfullness
/// @see Sparse::Solve
bool IsSolved() const;
/// Get the useful info of the last solution.
/// @see Sparse::Solve
const std::string SolutionMetadataLine(const std::string &delimiter) const;
/// Computes the smallest and the largest eigenvalue with the power method.
/// Requires SetMatrix to be called to compute the preconditioner.
/// Currently works for internal methods only, since it uses internal matrix-vector multiplication.
......
......@@ -4,27 +4,23 @@
namespace INMOST {
template<>
int from_string<int>(std::string str)
{
int from_string<int>(std::string str) {
return atoi(str.c_str());
}
template<>
double from_string<double>(std::string str)
{
double from_string<double>(std::string str) {
return atof(str.c_str());
}
template<>
std::string from_string<std::string>(std::string str)
{
std::string from_string<std::string>(std::string str) {
return str;
}
std::string string_to_lower(const std::string &str)
{
std::string string_to_lower(const std::string &str) {
std::string lower = std::string(str);
for (int i = 0; i < lower.length(); i++)
for (int i = 0; i < lower.length(); i++)
lower[i] = (char) tolower(lower[i]);
return lower;
}
......@@ -35,4 +31,12 @@ namespace INMOST {
#endif
}
int MPIGetRank() {
int rank = 0;
#if defined(USE_MPI)
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
#endif
return rank;
}
}
......@@ -28,6 +28,8 @@ namespace INMOST
std::string string_to_lower(const std::string &str);
void MPIBarrier();
int MPIGetRank();
class BinaryHeap
{
......
#include <inmost.h>
#include "SolverMaster.h"
#include "../Misc/utils.h"
#if defined(USE_SOLVER)
namespace INMOST {
int *Solver::argc = NULL;
char ***Solver::argv = NULL;
int *Solver::argc = NULL;
char ***Solver::argv = NULL;
bool Solver::is_initialized = false;
bool Solver::is_finalized = false;
bool Solver::is_finalized = false;
std::vector<SolverParameters> Solver::parameters = std::vector<SolverParameters>();
void Solver::SetParameterEnum(std::string name, INMOST_DATA_ENUM_TYPE value) {SetParameter(name,to_string(value));}
void Solver::SetParameterReal(std::string name, INMOST_DATA_REAL_TYPE value) {SetParameter(name,to_string(value));}
Solver::Solver(std::string solverName, std::string prefix, INMOST_MPI_Comm _comm) {
void Solver::SetParameterEnum(std::string name, INMOST_DATA_ENUM_TYPE value) { SetParameter(name, to_string(value)); }
void Solver::SetParameterReal(std::string name, INMOST_DATA_REAL_TYPE value) { SetParameter(name, to_string(value)); }
Solver::Solver(std::string solverName, std::string prefix, INMOST_MPI_Comm _comm) :
versbosity(SolverVerbosityLevel::Level0), preconditioner_time(0), iterations_time(0), is_solved(false) {
std::string lowerName = string_to_lower(solverName);
this->solver = SolverMaster::getSolver(solverName);
this->prefix = string_to_lower(prefix);
......@@ -37,7 +39,8 @@ namespace INMOST {
}
Solver::Solver(const Solver &other) {
Solver::Solver(const Solver &other) :
versbosity(other.versbosity), preconditioner_time(other.preconditioner_time), iterations_time(other.iterations_time), is_solved(other.is_solved) {
this->solver = SolverMaster::getSolver(other.solver->SolverName());
this->solver->Copy(other.solver);
this->prefix = other.prefix;
......@@ -62,8 +65,13 @@ namespace INMOST {
Solver &Solver::operator=(const Solver &other) {
if (this != &other) {
this->solver->SetCommunicator(other.solver->GetCommunicator());
this->prefix = other.prefix;
this->solver->Assign(other.solver);
this->prefix = other.prefix;
this->versbosity = other.versbosity;
this->preconditioner_time = other.preconditioner_time;
this->iterations_time = other.iterations_time;
this->is_solved = other.is_solved;
}
return *this;
}
......@@ -76,11 +84,19 @@ namespace INMOST {
return prefix;
}
SolverVerbosityLevel Solver::VerbosityLevel() const {
return versbosity;
}
void Solver::SetVerbosityLevel(INMOST::SolverVerbosityLevel level) {
versbosity = level;
}
void Solver::Initialize(int *argc, char ***argv, const char *database) {
Solver::argc = argc;
Solver::argv = argv;
Solver::argc = argc;
Solver::argv = argv;
Solver::is_initialized = true;
Solver::is_finalized = false;
Solver::is_finalized = false;
#if defined(USE_MPI)
{
int flag = 0;
......@@ -109,7 +125,7 @@ namespace INMOST {
}
}
#endif
Solver::is_finalized = true;
Solver::is_finalized = true;
Solver::is_initialized = false;
}
......@@ -122,7 +138,10 @@ namespace INMOST {
}
void Solver::SetMatrix(Sparse::Matrix &A, bool ModifiedPattern, bool OldPreconditioner) {
double preconditioner_timer_start = Timer();
solver->SetMatrix(A, ModifiedPattern, OldPreconditioner);
INMOST::MPIBarrier();
preconditioner_time = Timer() - preconditioner_timer_start;
}
bool Solver::Solve(INMOST::Sparse::Vector &RHS, INMOST::Sparse::Vector &SOL) {
......@@ -138,7 +157,24 @@ namespace INMOST {
for (Sparse::Vector::iterator ri = SOL.Begin(); ri != SOL.End(); ++ri) *ri = 0.0;
} else throw InconsistentSizesInSolver;
}
return solver->Solve(RHS, SOL);
double iterations_timer_start = Timer();
is_solved = solver->Solve(RHS, SOL);
INMOST::MPIBarrier();
iterations_time = Timer() - iterations_timer_start;
if (versbosity > SolverVerbosityLevel::Level1) {
#if defined(USE_MPI)
int rank;
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
if (rank == 0) {
std::cout << SolutionMetadataLine("\t") << std::endl;
}
#endif
}
return is_solved;
}
bool Solver::Clear() {
......@@ -165,6 +201,29 @@ namespace INMOST {
return solver->ReturnReason();
}
double Solver::PreconditionerTime() const {
return preconditioner_time;
}
double Solver::IterationsTime() const {
return iterations_time;
}
bool Solver::IsSolved() const {
return is_solved;
}
const std::string Solver::SolutionMetadataLine(const std::string &delimiter) const {
return this->SolverName() + delimiter +
this->SolverPrefix() + delimiter +
INMOST::to_string(this->IsSolved()) + delimiter +
INMOST::to_string(this->Iterations()) + delimiter +
INMOST::to_string(this->Residual()) + delimiter +
INMOST::to_string(this->PreconditionerTime()) + delimiter +
INMOST::to_string(this->IterationsTime()) + delimiter +
INMOST::to_string(this->PreconditionerTime() + this->IterationsTime());
}
INMOST_DATA_REAL_TYPE Solver::Condest(INMOST_DATA_REAL_TYPE tol, unsigned int maxits) {
if (!solver->isMatrixSet()) throw MatrixNotSetInSolver;
return solver->Condest(tol, maxits);
......@@ -200,12 +259,12 @@ namespace INMOST {
XMLReader::XMLTree root = reader.ReadXML();
if (root.tag.name != "SolverParameters") {
std::cout << __FILE__ << ":" << __LINE__ << ": Bad XML database file " << xml_database << "!" << std::endl;
std::cout << __FILE__ << ":" << __LINE__ << ": Bad XML database file " << xml_database << "!" << std::endl;
return;
}
for (xml_reader_tree_iterator_t solver = root.children.begin(); solver < root.children.end(); solver++) {
std::string solverName = string_to_lower((*solver).tag.name);
std::string solverName = string_to_lower((*solver).tag.name);
std::string internalFile = "";
if ((*solver).tag.attributes.size() != 0) {
......
......@@ -64,7 +64,8 @@ namespace TTSP {
current_index = index;
}
AlternatingOptimizer::AlternatingOptimizer(const OptimizationParametersSpace &space) : OptimizerInterface(space, 10), current_handler_index(0) {
AlternatingOptimizer::AlternatingOptimizer(const OptimizationParametersSpace &space, const OptimizerProperties &properties, std::size_t buffer_capacity) :
OptimizerInterface(space, properties, buffer_capacity), current_handler_index(0) {
const OptimizationParameters &parameters = space.GetParameters();
handlers.reserve(parameters.size());
std::for_each(parameters.begin(), parameters.end(), [this](const OptimizationParametersEntry &entry) {
......
......@@ -43,7 +43,7 @@ namespace TTSP {
std::size_t current_handler_index;
std::vector<AlternatingParameterHandler> handlers;
public:
AlternatingOptimizer(const OptimizationParametersSpace &space);
AlternatingOptimizer(const OptimizationParametersSpace &space, const OptimizerProperties &properties, std::size_t buffer_capacity);
OptimizationParameterPoints MakeOptimizationIteration(INMOST::Solver &solver, INMOST::Sparse::Matrix &matrix,
INMOST::Sparse::Vector &RHS) override;
......
......@@ -6,6 +6,12 @@
namespace TTSP {
template<typename T>
int sgn(T val) {
return (T(0) < val) - (val < T(0));
}
AnnealingUniformDistribution::AnnealingUniformDistribution() : distribution(0.0, 1.0) {
int rank, size;
......@@ -39,7 +45,7 @@ namespace TTSP {
bool AnnealingParameterHandler::DEFAULT_USE_CLOSEST = true;
AnnealingParameterHandler::AnnealingParameterHandler(const OptimizationParameter &parameter, const OptimizerInterface &optimizer) :
parameter(parameter),
count(0), parameter(parameter), current_value(parameter.GetDefaultValue()),
optimization_type(AnnealingParameterHandler::DEFAULT_OPTIMIZATION_TYPE),
temp0(AnnealingParameterHandler::DEFAULT_TEMP0),
decrement(AnnealingParameterHandler::DEFAULT_DECREMENT),
......@@ -105,6 +111,7 @@ namespace TTSP {
}
AnnealingParameterHandler::AnnealingParameterHandler(const AnnealingParameterHandler &other) :
count(other.count), current_value(other.current_value),
parameter(other.parameter), optimization_type(other.optimization_type), temp0(other.temp0), decrement(other.decrement),
allow_oscillation(other.allow_oscillation), oscillation_temp(other.oscillation_temp), strict_bound(other.strict_bound),
use_closest(other.use_closest) {}
......@@ -121,6 +128,10 @@ namespace TTSP {
return temp0;
}
double AnnealingParameterHandler::GetCurrentTemp() const {
return temp0 * std::exp(-decrement * std::pow(count, 1.0 / 2.0));
}
double AnnealingParameterHandler::GetDecrement() const {
return decrement;
}
......@@ -141,17 +152,113 @@ namespace TTSP {
return use_closest;
}
AnnealingOptimizer::AnnealingOptimizer(const OptimizationParametersSpace &space) : OptimizerInterface(space, 10), current_handler_index(0) {
double AnnealingParameterHandler::GetNextValue() const {
double current = GetCurrentValue();
double temp = GetCurrentTemp();
double alpha = random.next();
double z = sgn(alpha - 1.0 / 2.0) * temp * (std::pow(1.0 + 1.0 / temp, std::abs(2 * alpha - 1)) - 1);
double a = parameter.GetMinimalValue();
double b = parameter.GetMaximumValue();
double bound = (b - a);
if (strict_bound) {
double aa = current - a;
double bb = b - current;
if (bb < aa) {
bound = 2 * bb;
} else {
bound = 2 * aa;
}
}
double next = current + z * bound;
while (!((next >= a) && (next <= b))) {
alpha = random.next();
z = sgn(alpha - 1.0 / 2.0) * temp * (std::pow(1.0 + 1.0 / temp, std::abs(2 * alpha - 1)) - 1);
next = current + z * bound;
}
if (use_closest) {
next = parameter.GetClosestTo(next);
}
if (optimization_type == EXPONENT) {
return std::pow(10, -next);
} else {
return next;
}
}
double AnnealingParameterHandler::GetRandom() const {
return random.next();
}
void AnnealingParameterHandler::SetValue(double value) {
current_value = value;
count += 1;
}
double AnnealingParameterHandler::GetCurrentValue() const {
return current_value;
}
AnnealingOptimizer::AnnealingOptimizer(const OptimizationParametersSpace &space, const OptimizerProperties &properties, std::size_t buffer_capacity) :
OptimizerInterface(space, properties, buffer_capacity), current_handler_index(0), values(space.GetParameters().size()) {
const OptimizationParameters &parameters = space.GetParameters();
handlers.reserve(parameters.size());
std::for_each(parameters.begin(), parameters.end(), [this](const OptimizationParametersEntry &entry) {
handlers.emplace_back(AnnealingParameterHandler(entry.first, *this));
});
std::transform(parameters.begin(), parameters.end(), values.begin(), [](const OptimizationParametersEntry &entry) {
return entry.first.GetDefaultValue();
});
}
OptimizationParameterPoints AnnealingOptimizer::MakeOptimizationIteration(INMOST::Solver &solver, INMOST::Sparse::Matrix &matrix,
INMOST::Sparse::Vector &RHS) {
return TTSP::OptimizationParameterPoints();