From e255b9e2c97827ad1e731d973f37f04983b6dfdb Mon Sep 17 00:00:00 2001 From: Anne de Jong Date: Thu, 5 Jan 2017 16:25:16 +0100 Subject: [PATCH] Updated all GUI stuff to new protobuf config of model containing a system, solver and solution --- src/CMakeLists.txt | 150 +++++++------- src/ductbc/pressurebc.h | 5 +- src/gui/mainwindow.cpp | 18 +- src/gui/solver_dialog.cpp | 20 +- src/gui/solver_dialog.h | 8 +- src/gui/solver_worker.cpp | 12 +- src/gui/solver_worker.h | 4 +- src/main.cpp | 170 ++++++++-------- src/protobuf/message_tools.cpp | 10 +- src/protobuf/message_tools.h | 2 +- src/solver/bracket_root.h | 194 +++++++++--------- src/solver/brent.cpp | 360 ++++++++++++++++----------------- src/sys/tasystem.cpp | 18 +- src/tasmet_constants.h | 2 +- src/tasmet_solvemodel.cpp | 140 ++++++------- src/tasmet_types.h | 172 ++++++++-------- 16 files changed, 650 insertions(+), 635 deletions(-) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index aa82546..e5f709f 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,75 +1,75 @@ -#================ -# The TaSMET code -#================ - -include_directories( - ${PYTHON_INCLUDE_PATH} - ${ARMADILLO_INCLUDE_DIRS} - . - protobuf - duct - ductbc - # duct/cell - # duct/connectors - # duct/eq - # duct/geom - # duct/geom/grid - # duct/drag - # duct/heat - # mech - # python - # seg - # sol - sys - solver - # var - # volume - ) - -add_library(tasmet_src - tasmet_tracer.cpp - tasmet_exception.cpp - tasmet_assert.cpp - tasmet_pyeval.cpp - - protobuf/message_tools.cpp - - duct/grid.cpp - duct/geom.cpp - duct/duct.cpp - - ductbc/ductbc.cpp - ductbc/pressurebc.cpp - ductbc/adiabaticwall.cpp - - funcs/bessel.cpp - funcs/cbessj.cpp - funcs/rottfuncs.cpp - funcs/skewsine.cpp - funcs/sph_j.cpp - - solver/solver.cpp - solver/system.cpp - solver/newton_raphson.cpp - solver/brent.cpp - - material/gas.cpp - material/air.cpp - material/helium.cpp - material/nitrogen.cpp - material/solid.cpp - material/adiabatictemp.cpp - - # sys/enginesystem.cpp - sys/globalconf.cpp - sys/tasmet_variable.cpp - sys/jaccol.cpp - sys/jacobian.cpp - sys/jacrow.cpp - sys/tasystem.cpp - sys/triplets.cpp -) - -add_subdirectory(gui) -add_subdirectory(protobuf) -target_link_libraries(tasmet_src Qt5::Widgets) +#================ +# The TaSMET code +#================ + +include_directories( + ${PYTHON_INCLUDE_PATH} + ${ARMADILLO_INCLUDE_DIRS} + . + protobuf + duct + ductbc + # duct/cell + # duct/connectors + # duct/eq + # duct/geom + # duct/geom/grid + # duct/drag + # duct/heat + # mech + # python + # seg + # sol + sys + solver + # var + # volume + ) + +add_library(tasmet_src + tasmet_tracer.cpp + tasmet_exception.cpp + tasmet_assert.cpp + tasmet_pyeval.cpp + + protobuf/message_tools.cpp + + duct/grid.cpp + duct/geom.cpp + duct/duct.cpp + + ductbc/ductbc.cpp + ductbc/pressurebc.cpp + ductbc/adiabaticwall.cpp + + funcs/bessel.cpp + funcs/cbessj.cpp + funcs/rottfuncs.cpp + funcs/skewsine.cpp + funcs/sph_j.cpp + + solver/solver.cpp + solver/system.cpp + solver/newton_raphson.cpp + solver/brent.cpp + + material/gas.cpp + material/air.cpp + material/helium.cpp + material/nitrogen.cpp + material/solid.cpp + material/adiabatictemp.cpp + + # sys/enginesystem.cpp + sys/globalconf.cpp + sys/tasmet_variable.cpp + sys/jaccol.cpp + sys/jacobian.cpp + sys/jacrow.cpp + sys/tasystem.cpp + sys/triplets.cpp +) + +add_subdirectory(gui) +add_subdirectory(protobuf) +target_link_libraries(tasmet_src Qt5::Widgets) diff --git a/src/ductbc/pressurebc.h b/src/ductbc/pressurebc.h index a7c37fc..4dccb2f 100644 --- a/src/ductbc/pressurebc.h +++ b/src/ductbc/pressurebc.h @@ -9,11 +9,8 @@ #ifndef PRESSUREBC_H #define PRESSUREBC_H #include "segment.h" - +#include "ductbc.pb.h" class TaSystem; -namespace pb{ - class DuctBc; -} class Variable; class PressureBc: public Segment { diff --git a/src/gui/mainwindow.cpp b/src/gui/mainwindow.cpp index 11e6879..ed532fb 100644 --- a/src/gui/mainwindow.cpp +++ b/src/gui/mainwindow.cpp @@ -54,8 +54,8 @@ namespace { TaSMETMainWindow::TaSMETMainWindow(): _window(new Ui::MainWindow()), - _model(pb::Model::default_instance()) - _system(*_model.mutable_sytem()) + _model(pb::Model::default_instance()), + _system(*_model.mutable_system()) { if(!_window) throw TaSMETBadAlloc(); @@ -224,7 +224,7 @@ void TaSMETMainWindow::changed() { windowtitle += QString::fromStdString(_filepath); } else { - windowtitle += default_system_name; + windowtitle += default_model_name; } setWindowTitle(windowtitle); @@ -254,10 +254,11 @@ void TaSMETMainWindow::changed() { } -void TaSMETMainWindow::set(const pb::System& sys) { +void TaSMETMainWindow::set(const pb::Model& model) { TRACE(15,"set()"); _init = true; + const pb::System& sys = model.system(); _window->nf->setValue(sys.nf()); _window->freq->setText(QString::number(sys.freq())); _window->p0->setText(QString::number(sys.p0())); @@ -265,7 +266,8 @@ void TaSMETMainWindow::set(const pb::System& sys) { _window->systemtype->setCurrentIndex((int) sys.systemtype()); _window->gastype->setCurrentIndex((int) sys.gastype()); - _system = sys; + _model = model; + _system = *_model.mutable_system(); _init = false; changed(); @@ -406,7 +408,7 @@ void TaSMETMainWindow::on_actionSolve_triggered() { SolverDialog *d; try { - d = new SolverDialog(this,_system); + d = new SolverDialog(this,_system,*_model.mutable_sparams()); } catch(TaSMETError &e) { e.show_user("Solver failed to initialize"); @@ -425,8 +427,8 @@ bool TaSMETMainWindow::isDirty() const { if(_filepath.size()==0) return true; try { - pb::System filesys = loadSystem(_filepath); - bool dirty = !compareSys(filesys,_system); + pb::Model filemodel = loadMessage(_filepath); + bool dirty = !compareMessage(filemodel,_model); VARTRACE(15,dirty); return dirty; } diff --git a/src/gui/solver_dialog.cpp b/src/gui/solver_dialog.cpp index 68dc031..1e15de6 100644 --- a/src/gui/solver_dialog.cpp +++ b/src/gui/solver_dialog.cpp @@ -9,6 +9,8 @@ #include "solver_dialog.h" #include "ui_solver_dialog.h" #include "system.pb.h" +#include "solver.pb.h" + #include "tasmet_exception.h" #include "tasmet_constants.h" #include "tasmet_tracer.h" @@ -20,9 +22,11 @@ #include SolverDialog::SolverDialog(QWidget* parent, - pb::System& sys): + pb::System& sys, + pb::SolverParams& sparams): QDialog(parent), _sys(sys), + _sparams(sparams), _dialog(new Ui::solver_dialog()) { @@ -85,7 +89,7 @@ SolverDialog::SolverDialog(QWidget* parent, _plot->legend->setVisible(true); _plot->replot(); - set(_sys.solverparams()); + set(_sparams); setEnabled(true); @@ -113,9 +117,9 @@ void SolverDialog::set(const pb::SolverParams& sol) { void SolverDialog::changed() { TRACE(15,"changed"); if(_init) return; - _sys.mutable_solverparams()->set_funtol(_dialog->funtol->text().toDouble()); - _sys.mutable_solverparams()->set_reltol(_dialog->reltol->text().toDouble()); - _sys.mutable_solverparams()->set_solvertype((pb::SolverType) _dialog->solvertype->currentIndex()); + _sparams.set_funtol(_dialog->funtol->text().toDouble()); + _sparams.set_reltol(_dialog->reltol->text().toDouble()); + _sparams.set_solvertype((pb::SolverType) _dialog->solvertype->currentIndex()); } void SolverDialog::solver_progress(const SolverProgress& progress){ @@ -123,8 +127,8 @@ void SolverDialog::solver_progress(const SolverProgress& progress){ // VARTRACE(15,progress.fun_err); // VARTRACE(15,progress.iteration); - d funtol = _sys.solverparams().funtol(); - d reltol = _sys.solverparams().reltol(); + d funtol = _sparams.funtol(); + d reltol = _sparams.reltol(); _funer->addData(progress.iteration,progress.fun_err); _reler->addData(progress.iteration,progress.rel_err); @@ -160,7 +164,7 @@ void SolverDialog::on_solve_clicked() { qRegisterMetaType(); - _solver_worker = new SolverWorker(_sys); + _solver_worker = new SolverWorker(_sys,_sparams); QThread* thread = new QThread; _solver_worker->moveToThread(thread); diff --git a/src/gui/solver_dialog.h b/src/gui/solver_dialog.h index fe1b4c1..87b57f1 100644 --- a/src/gui/solver_dialog.h +++ b/src/gui/solver_dialog.h @@ -15,10 +15,12 @@ namespace pb { class System; + class SolverParams; } namespace Ui { class solver_dialog; } + class QCustomPlot; class QCPGraph; class SolverWorker; @@ -28,6 +30,7 @@ class SolverDialog: public QDialog { Q_OBJECT pb::System& _sys; // Reference to system + pb::SolverParams& _sparams; Ui::solver_dialog* _dialog; @@ -41,11 +44,11 @@ class SolverDialog: public QDialog { public: SolverDialog(QWidget* parent, - pb::System& sys); + pb::System& sys, + pb::SolverParams& sparams); ~SolverDialog(); - void set(const pb::SolverParams&); public slots: void solver_progress(const SolverProgress&); private slots: @@ -64,6 +67,7 @@ private: // Called whenever the user changes input values void changed(); void setEnabled(bool); + void set(const pb::SolverParams&); }; diff --git a/src/gui/solver_worker.cpp b/src/gui/solver_worker.cpp index 575c623..ba5f3f7 100644 --- a/src/gui/solver_worker.cpp +++ b/src/gui/solver_worker.cpp @@ -9,13 +9,17 @@ #include "solver_worker.h" #include #include "tasmet_tracer.h" + #include "system.pb.h" +#include "solver.pb.h" + #include -#include -SolverWorker::SolverWorker(pb::System& sys): + + +SolverWorker::SolverWorker(const pb::System& sys,const pb::SolverParams& sparams): _run(false), - _reltol(sys.solverparams().reltol()), - _funtol(sys.solverparams().funtol()) + _reltol(sparams.reltol()), + _funtol(sparams.funtol()) { } diff --git a/src/gui/solver_worker.h b/src/gui/solver_worker.h index e907081..0a9a5fb 100644 --- a/src/gui/solver_worker.h +++ b/src/gui/solver_worker.h @@ -16,6 +16,7 @@ #include namespace pb{ class System; + class SolverParams; } Q_DECLARE_METATYPE(SolverProgress); @@ -28,8 +29,7 @@ class SolverWorker: public QObject { bool _converged = false; d _funtol,_reltol; public: - - SolverWorker(pb::System& sys); + SolverWorker(const pb::System& sys,const pb::SolverParams& sparams); ~SolverWorker(); void solver_stop(); diff --git a/src/main.cpp b/src/main.cpp index 6e842c5..9e84c06 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -1,85 +1,85 @@ -// main.cpp -// -// last-edit-by: J.A. de Jong -// -// Description: -// Main program to run -////////////////////////////////////////////////////////////////////// -#include -#include -#include - -// For using python from within Qt -#include - - -#include "tasmet_config.h" -#include "tasmet_tracer.h" -#include "gui/mainwindow.h" - -#include -#include - -void catchUnixSignals(const std::vector& quitSignals, - const std::vector& ignoreSignals = std::vector()) { - - auto handler = [](int sig) ->void { - printf("\nquit the application (user request signal = %d).\n", sig); - QCoreApplication::quit(); - }; - - // all these signals will be ignored. - for ( int sig : ignoreSignals ) - signal(sig, SIG_IGN); - - // each of these signals calls the handler (quits the QCoreApplication). - for ( int sig : quitSignals ) - signal(sig, handler); -} - - -int main(int argc, char *argv[]) { - - // This can be used to store files in the application. Very useful - // when our application grows bigger - // Q_INIT_RESOURCE(application); - INITTRACE(15); - - // Initialize PythonQt - // PythonQt::init(PythonQt::IgnoreSiteModule | PythonQt::RedirectStdOut); - PythonQt::init(); - - PythonQt* pyqt = PythonQt::self(); - PythonQtObjectPtr context = pyqt->getMainModule(); - - QVariant rv = context.evalScript("from math import *\n"); - if(pyqt->hadError()) { - return -1; - } - - QApplication app(argc, argv); - - catchUnixSignals({SIGQUIT, SIGINT, SIGTERM, SIGHUP}); - - // QCoreApplication::setOrganizationName("QtProject"); - QCoreApplication::setApplicationName(appname); - QCoreApplication::setApplicationVersion(appversion); - // QCommandLineParser parser; - // parser.setApplicationDescription(QCoreApplication::applicationName()); - // parser.addHelpOption(); - // parser.addVersionOption(); - // parser.addPositionalArgument("file", "The file to open."); - // parser.process(app); - - // if (!parser.positionalArguments().isEmpty()) - // mainWin.loadFile(parser.positionalArguments().first()); - - - - TaSMETMainWindow win; - - win.show(); - - return app.exec(); -} -////////////////////////////////////////////////////////////////////// +// main.cpp +// +// last-edit-by: J.A. de Jong +// +// Description: +// Main program to run +////////////////////////////////////////////////////////////////////// +#include +#include +#include + +// For using python from within Qt +#include + + +#include "tasmet_config.h" +#include "tasmet_tracer.h" +#include "gui/mainwindow.h" + +#include +#include + +void catchUnixSignals(const std::vector& quitSignals, + const std::vector& ignoreSignals = std::vector()) { + + auto handler = [](int sig) ->void { + printf("\nquit the application (user request signal = %d).\n", sig); + QCoreApplication::quit(); + }; + + // all these signals will be ignored. + for ( int sig : ignoreSignals ) + signal(sig, SIG_IGN); + + // each of these signals calls the handler (quits the QCoreApplication). + for ( int sig : quitSignals ) + signal(sig, handler); +} + + +int main(int argc, char *argv[]) { + + // This can be used to store files in the application. Very useful + // when our application grows bigger + // Q_INIT_RESOURCE(application); + INITTRACE(15); + + // Initialize PythonQt + // PythonQt::init(PythonQt::IgnoreSiteModule | PythonQt::RedirectStdOut); + PythonQt::init(); + + PythonQt* pyqt = PythonQt::self(); + PythonQtObjectPtr context = pyqt->getMainModule(); + + QVariant rv = context.evalScript("from math import *\n"); + if(pyqt->hadError()) { + return -1; + } + + QApplication app(argc, argv); + + catchUnixSignals({SIGQUIT, SIGINT, SIGTERM, SIGHUP}); + + // QCoreApplication::setOrganizationName("QtProject"); + QCoreApplication::setApplicationName(appname); + QCoreApplication::setApplicationVersion(appversion); + // QCommandLineParser parser; + // parser.setApplicationDescription(QCoreApplication::applicationName()); + // parser.addHelpOption(); + // parser.addVersionOption(); + // parser.addPositionalArgument("file", "The file to open."); + // parser.process(app); + + // if (!parser.positionalArguments().isEmpty()) + // mainWin.loadFile(parser.positionalArguments().first()); + + + + TaSMETMainWindow win; + + win.show(); + + return app.exec(); +} +////////////////////////////////////////////////////////////////////// diff --git a/src/protobuf/message_tools.cpp b/src/protobuf/message_tools.cpp index 772415b..2265ca9 100644 --- a/src/protobuf/message_tools.cpp +++ b/src/protobuf/message_tools.cpp @@ -7,7 +7,6 @@ ////////////////////////////////////////////////////////////////////// #include "message_tools.h" #include -#include #include "tasmet_types.h" #include "tasmet_io.h" #include "tasmet_exception.h" @@ -77,9 +76,14 @@ bool compareMessage(const T& s1,const T& s2) { return (s1.SerializeAsString()==s2.SerializeAsString()); } -template <> +// Explicit instantiation for pb::Model +template bool compareMessage(const pb::Model& s1,const pb::Model& s2); -template <> + +template pb::Model loadMessage(const string& filepath); + +template +void saveMessage(const string& filepath,const pb::Model& model); ////////////////////////////////////////////////////////////////////// diff --git a/src/protobuf/message_tools.h b/src/protobuf/message_tools.h index 083f7f8..9aeeeea 100644 --- a/src/protobuf/message_tools.h +++ b/src/protobuf/message_tools.h @@ -21,7 +21,7 @@ void saveMessage(const string& filepath,const T& sys); // Returns true when the two systems are equal template -bool compareSys(const T& s1,const T& s2); +bool compareMessage(const T& s1,const T& s2); diff --git a/src/solver/bracket_root.h b/src/solver/bracket_root.h index 4ff86a6..72273c1 100644 --- a/src/solver/bracket_root.h +++ b/src/solver/bracket_root.h @@ -1,97 +1,97 @@ -// bracket_lin.h -// -// Author: J.A. de Jong -// -// Description: Bracket a function's root given a guess value, returns -// another guess where the solution has the opposite sign -////////////////////////////////////////////////////////////////////// -#pragma once -#ifndef BRACKET_LIN_H -#define BRACKET_LIN_H -#include "system.h" -#include "tasmet_tracer.h" -#include - -inline int sign(d& x){ - return x < 0? -1: 1; -} -const us maxiter = 100; - -inline std::pair bracket_root(NoGradientNonlinearSystem& sys,const d guess) { - - TRACE(10,"bracket_root()"); - - d xa = guess; - d fa = sys.residual(xa); - - if(fa == 0) { - // Guess was exactly right - return std::make_pair(xa,xa); - } - - d fb = fa; - - // We add here a small amount, such that, in case the first guess - // is 0, we - d xb = xa+std::numeric_limits::epsilon(); - - d factor = 1.001; - d delta = xa/10 + std::numeric_limits::epsilon(); - - us iter = 0; - - int step = 64; - - bool switched = false; - - while (iter < maxiter) { - - VARTRACE(5,xb); - - xb += delta*factor; - - fb = sys.residual(xb); - - if(sign(fa) != sign(fb)) { - // Found it! - VARTRACE(5,xa); - VARTRACE(5,xb); - VARTRACE(5,fa); - VARTRACE(5,fb); - return std::make_pair(xa,xb); - } - - if (abs(fb) < abs(fa)){ - // We found a point which is closer to the real root - xa = xb; - fa = fb; - } - else if(!switched) { - TRACE(5,"Switching search direction"); - switched = true; - // We searched in the wrong direction - factor*=-1; - } - - // Heuristic: normally it's best not to increase the step sizes as we'll just end up - // with a really wide range to search for the root. However, if the initial guess was *really* - // bad then we need to speed up the search otherwise we'll take forever if we're orders of - // magnitude out. This happens most often if the guess is a small value (say 1) and the result - // we're looking for is close to std::numeric_limits::min(). - // - - if((10000 - iter) % step == 0) { - factor *=2; - if(step > 1) step/=2; - } - - iter++; - } - - WARN("Failed to bracket root"); - - return std::make_pair(xa,xb); - -} -#endif // BRACKET_LIN_H -////////////////////////////////////////////////////////////////////// +// bracket_lin.h +// +// Author: J.A. de Jong +// +// Description: Bracket a function's root given a guess value, returns +// another guess where the solution has the opposite sign +////////////////////////////////////////////////////////////////////// +#pragma once +#ifndef BRACKET_LIN_H +#define BRACKET_LIN_H +#include "system.h" +#include "tasmet_tracer.h" +#include + +inline int sign(d& x){ + return x < 0? -1: 1; +} +const us maxiter = 100; + +inline std::pair bracket_root(NoGradientNonlinearSystem& sys,const d guess) { + + TRACE(10,"bracket_root()"); + + d xa = guess; + d fa = sys.residual(xa); + + if(fa == 0) { + // Guess was exactly right + return std::make_pair(xa,xa); + } + + d fb = fa; + + // We add here a small amount, such that, in case the first guess + // is 0, we + d xb = xa+std::numeric_limits::epsilon(); + + d factor = 1.001; + d delta = xa/10 + std::numeric_limits::epsilon(); + + us iter = 0; + + int step = 64; + + bool switched = false; + + while (iter < maxiter) { + + VARTRACE(5,xb); + + xb += delta*factor; + + fb = sys.residual(xb); + + if(sign(fa) != sign(fb)) { + // Found it! + VARTRACE(5,xa); + VARTRACE(5,xb); + VARTRACE(5,fa); + VARTRACE(5,fb); + return std::make_pair(xa,xb); + } + + if (abs(fb) < abs(fa)){ + // We found a point which is closer to the real root + xa = xb; + fa = fb; + } + else if(!switched) { + TRACE(5,"Switching search direction"); + switched = true; + // We searched in the wrong direction + factor*=-1; + } + + // Heuristic: normally it's best not to increase the step sizes as we'll just end up + // with a really wide range to search for the root. However, if the initial guess was *really* + // bad then we need to speed up the search otherwise we'll take forever if we're orders of + // magnitude out. This happens most often if the guess is a small value (say 1) and the result + // we're looking for is close to std::numeric_limits::min(). + // + + if((10000 - iter) % step == 0) { + factor *=2; + if(step > 1) step/=2; + } + + iter++; + } + + WARN("Failed to bracket root"); + + return std::make_pair(xa,xb); + +} +#endif // BRACKET_LIN_H +////////////////////////////////////////////////////////////////////// diff --git a/src/solver/brent.cpp b/src/solver/brent.cpp index 8c1a5eb..d17fc44 100644 --- a/src/solver/brent.cpp +++ b/src/solver/brent.cpp @@ -1,180 +1,180 @@ -// brent.cpp -// -// last-edit-by: J.A. de Jong -// -// Description: -// Brent's root finding algorithm as implemented from Wikipedia: -// https://en.wikipedia.org/wiki/Brent's_method -////////////////////////////////////////////////////////////////////// -#define TRACERPLUS (10) -#include "brent.h" -#include "tasmet_tracer.h" -#include "tasmet_exception.h" -#include -#include "tasmet_constants.h" -#include // std::swap -#include "bracket_root.h" - -const d eps = std::numeric_limits::epsilon(); - -Brent::Brent(const NoGradientNonlinearSystem& sys,d reltol): - Solver(sys), - _reltol(reltol) -{ - - TRACE(15,"Brent::Brent"); - #ifdef TASMET_DEBUG - bool ok=true; - - ok&=(reltol >= eps); - - if(!ok){ - throw TaSMETError("Brent: invalid arguments"); - } - #endif -} -namespace { - inline bool is_between(d x,d a,d b) { - const d& mi = std::min(a,b); - const d& ma = std::max(a,b); - return ((mi <= x) && (x <=ma)); - } -} - -void Brent::start_implementation(NoGradientNonlinearSystem& system, - progress_callback* callback) { - - TRACE(15,"Brent::start_implementation"); - - us niter = 0; - - d a = system.getSolution(); - - std::pair brackets = bracket_root(system,a); - a = brackets.first; - d b = brackets.second; - - d fa = system.residual(a); - if((fa) == 0) { - TRACE(15,"Found root during bracketing"); - return; - } - - d fb = system.residual(b); - if((fb) == 0) { - TRACE(15,"Found root during bracketing"); - return; - } - d c = a; - d fc = fa; - - d s; // Test point - d fs; // Function at test point - - d t,u,v; - d d_; - - if(fa*fb>=0){ - WARN("Brent solver failed: root is not bracketed"); - return; - } - - if(abs(fa) < abs(fb)) { - std::swap(a,b); - std::swap(fa,fb); - } - - bool mflag = true; - bool bisec_flag; - SolverProgress progress; - - while(true) { - - if((fa!=fc) && (fb!=fc)){ - // Inverse quadratic interpolation - TRACE(5,"IQI"); - t = a*fb*fc/((fa-fb)*(fa-fc)); - u = b*fa*fc/((fb-fa)*(fb-fc)); - v = c*fa*fb/((fc-fa)*(fc-fb)); - - s = t+u+v; - } - else { - // Secant method - TRACE(5,"Secant"); - s = b-fb*(b-a)/(fb-fa); - } - - bisec_flag = false; - - d abssmb = abs(s-b); - d absbmc = abs(b-c); - d abscmd = abs(c-d_); - - VARTRACE(5,s); - - if((bisec_flag |= (!is_between(s,(3*a+b)/4,b)))) goto bflag; - TRACE(5,"Survived 1"); - if(bisec_flag |= (mflag && (abssmb >= absbmc/2))) goto bflag; - TRACE(5,"Survived 2"); - if(bisec_flag |= ((!mflag) && (abssmb >= abscmd/2))) goto bflag; - TRACE(5,"Survived 3"); - if(bisec_flag |= (mflag && (absbmc < abs(_reltol)))) goto bflag;; - TRACE(5,"Survived 4"); - bflag: - if(bisec_flag || ((!mflag) && (abscmd < abs(_reltol)))) { - TRACE(5,"Bisection"); - s = (a+b)/2; - mflag = true; - } - else { - mflag = false; - } - - fs = system.residual(s); - d_ = c; - - if(fa*fs < 0) { - c = b; - fc = fb; - b=s; - fb = fs; - } - else { - c = a; - fc = fa; - a = s; - fa = fs; - } - if(abs(fa) +#include "tasmet_constants.h" +#include // std::swap +#include "bracket_root.h" + +const d eps = std::numeric_limits::epsilon(); + +Brent::Brent(const NoGradientNonlinearSystem& sys,d reltol): + Solver(sys), + _reltol(reltol) +{ + + TRACE(15,"Brent::Brent"); + #ifdef TASMET_DEBUG + bool ok=true; + + ok&=(reltol >= eps); + + if(!ok){ + throw TaSMETError("Brent: invalid arguments"); + } + #endif +} +namespace { + inline bool is_between(d x,d a,d b) { + const d& mi = std::min(a,b); + const d& ma = std::max(a,b); + return ((mi <= x) && (x <=ma)); + } +} + +void Brent::start_implementation(NoGradientNonlinearSystem& system, + progress_callback* callback) { + + TRACE(15,"Brent::start_implementation"); + + us niter = 0; + + d a = system.getSolution(); + + std::pair brackets = bracket_root(system,a); + a = brackets.first; + d b = brackets.second; + + d fa = system.residual(a); + if((fa) == 0) { + TRACE(15,"Found root during bracketing"); + return; + } + + d fb = system.residual(b); + if((fb) == 0) { + TRACE(15,"Found root during bracketing"); + return; + } + d c = a; + d fc = fa; + + d s; // Test point + d fs; // Function at test point + + d t,u,v; + d d_; + + if(fa*fb>=0){ + WARN("Brent solver failed: root is not bracketed"); + return; + } + + if(abs(fa) < abs(fb)) { + std::swap(a,b); + std::swap(fa,fb); + } + + bool mflag = true; + bool bisec_flag; + SolverProgress progress; + + while(true) { + + if((fa!=fc) && (fb!=fc)){ + // Inverse quadratic interpolation + TRACE(5,"IQI"); + t = a*fb*fc/((fa-fb)*(fa-fc)); + u = b*fa*fc/((fb-fa)*(fb-fc)); + v = c*fa*fb/((fc-fa)*(fc-fb)); + + s = t+u+v; + } + else { + // Secant method + TRACE(5,"Secant"); + s = b-fb*(b-a)/(fb-fa); + } + + bisec_flag = false; + + d abssmb = abs(s-b); + d absbmc = abs(b-c); + d abscmd = abs(c-d_); + + VARTRACE(5,s); + + if((bisec_flag |= (!is_between(s,(3*a+b)/4,b)))) goto bflag; + TRACE(5,"Survived 1"); + if(bisec_flag |= (mflag && (abssmb >= absbmc/2))) goto bflag; + TRACE(5,"Survived 2"); + if(bisec_flag |= ((!mflag) && (abssmb >= abscmd/2))) goto bflag; + TRACE(5,"Survived 3"); + if(bisec_flag |= (mflag && (absbmc < abs(_reltol)))) goto bflag;; + TRACE(5,"Survived 4"); + bflag: + if(bisec_flag || ((!mflag) && (abscmd < abs(_reltol)))) { + TRACE(5,"Bisection"); + s = (a+b)/2; + mflag = true; + } + else { + mflag = false; + } + + fs = system.residual(s); + d_ = c; + + if(fa*fs < 0) { + c = b; + fc = fb; + b=s; + fb = fs; + } + else { + c = a; + fc = fa; + a = s; + fa = fs; + } + if(abs(fa)0) { - _solution = vd(size); - for(auto& val: sol) { - _solution(i) = val; - i++; - } - } + // const auto& sol = sys.solution(); + // us size = sol.size(), i=0; + // if(size>0) { + // _solution = vd(size); + // for(auto& val: sol) { + // _solution(i) = val; + // i++; + // } + // } } TaSystem::TaSystem(const TaSystem& o): diff --git a/src/tasmet_constants.h b/src/tasmet_constants.h index cce9f28..5cbe2b8 100644 --- a/src/tasmet_constants.h +++ b/src/tasmet_constants.h @@ -105,7 +105,7 @@ namespace constants { const int nvars_reserve=7; const int neqs_reserve=7; - const char* const system_fileext = ".tasmet"; + const char* const model_fileext = ".tasmet"; } // namespace constants diff --git a/src/tasmet_solvemodel.cpp b/src/tasmet_solvemodel.cpp index 723e0d6..c079be5 100644 --- a/src/tasmet_solvemodel.cpp +++ b/src/tasmet_solvemodel.cpp @@ -1,70 +1,70 @@ -// tasmet_solvemodel.cpp -// -// last-edit-by: J.A. de Jong -// -// Description: -// Program which can be used to solve a TaSMET Model from the CLI -////////////////////////////////////////////////////////////////////// -#include "protobuf/system_tools.h" // loadSystem, saveSystem -#include "tasmet_io.h" -#include "tasmet_exception.h" -#include "sys/tasystem.h" - -// For using python from within Qt -#include - -void usage(const char* fn) { - cout << "Usage: " << endl; - cout << fn << " [Input file] " << endl; -} - -int main(int argc, char *argv[]) { - - INITTRACE(10); - - // Initialize PythonQt - // PythonQt::init(PythonQt::IgnoreSiteModule | PythonQt::RedirectStdOut); - PythonQt::init(); - - PythonQt* pyqt = PythonQt::self(); - PythonQtObjectPtr context = pyqt->getMainModule(); - - QVariant rv = context.evalScript("from math import *\n"); - if(pyqt->hadError()) { - return -1; - } - - if(argc != 2) { - usage(argv[0]); - return -1; - } - - pb::System sys; - - string filename = argv[1]; - try { - sys = loadSystem(filename); - } - catch(TaSMETError &e) { - cerr << e.what() << endl; - return -1; - } - - cout << "Loaded model: " << endl; - cout << sys.DebugString(); - - TaSystem* system; - - try { - system = new TaSystem(sys); - } - catch(TaSMETError &e) { - cerr << "Model initialization error: " << endl; - cerr << e.what() << endl; - return -1; - } - -} - - -////////////////////////////////////////////////////////////////////// +// tasmet_solvemodel.cpp +// +// last-edit-by: J.A. de Jong +// +// Description: +// Program which can be used to solve a TaSMET Model from the CLI +////////////////////////////////////////////////////////////////////// +#include "protobuf/message_tools.h" // loadSystem, saveSystem +#include "tasmet_io.h" +#include "tasmet_exception.h" +#include "sys/tasystem.h" +#include "protobuf/model.pb.h" +// For using python from within Qt +#include + +void usage(const char* fn) { + cout << "Usage: " << endl; + cout << fn << " [Input file] " << endl; +} + +int main(int argc, char *argv[]) { + + INITTRACE(10); + + // Initialize PythonQt + // PythonQt::init(PythonQt::IgnoreSiteModule | PythonQt::RedirectStdOut); + PythonQt::init(); + + PythonQt* pyqt = PythonQt::self(); + PythonQtObjectPtr context = pyqt->getMainModule(); + + QVariant rv = context.evalScript("from math import *\n"); + if(pyqt->hadError()) { + return -1; + } + + if(argc != 2) { + usage(argv[0]); + return -1; + } + + pb::Model model; + + string filename = argv[1]; + try { + model = loadMessage(filename); + } + catch(TaSMETError &e) { + cerr << e.what() << endl; + return -1; + } + + cout << "Loaded model: " << endl; + cout << model.system().DebugString(); + + TaSystem* system; + + try { + system = new TaSystem(model.system()); + } + catch(TaSMETError &e) { + cerr << "Model initialization error: " << endl; + cerr << e.what() << endl; + return -1; + } + +} + + +////////////////////////////////////////////////////////////////////// diff --git a/src/tasmet_types.h b/src/tasmet_types.h index fdb21bf..3fe1277 100644 --- a/src/tasmet_types.h +++ b/src/tasmet_types.h @@ -1,86 +1,86 @@ -// tasmet_types.h -// -// Author: J.A. de Jong -// -// Description: -// Typedefs and namespace pollution for stuff that is almost always -// needed. -////////////////////////////////////////////////////////////////////// -#pragma once -#ifndef TASMET_TYPES_H -#define TASMET_TYPES_H -#include - -// We require C++ -#ifndef __cplusplus -#error The c++ compiler should be used. -#endif - -// Some header files I (almost) always need -#include // C++ std vector -#include // C++ string class -#include -#include // Our linear algebra package - -typedef size_t us; /* Size type I always use */ - -#define fillwith arma::fill // Fill vector or matrix with ... - -// To change the whole code to 32-bit floating points, change this to -// float. -#if TASMET_FLOAT == 32 -typedef float d; /* Shortcut for double */ -#elif TASMET_FLOAT == 64 -typedef double d; /* Shortcut for double */ -#else -static_assert(false,"TASMET_FLOAT should be either 32 or 64"); -#endif - -typedef std::complex c; /* Shortcut for c++ complex number */ - -using arma::abs; -using std::string; /* C++ Strings */ -using std::vector; // C++ Vectors - -/* Armadillo functions and objects we often need */ -using arma::sin; -using arma::cos; -using arma::exp; -using arma::ones; -using arma::zeros; -using arma::eye; -using arma::inv; -using arma::det; -using arma::log_det; -using arma::norm; -using arma::span; -using arma::linspace; -using arma::diagmat; - -// Constants -// I need these numbers so often, that they can be in the global namespace - -const c I(0,1); //Complex unity -const c sqI=sqrt(I); -const d sq2=sqrt(2.0); -const c sqmI=sqrt(-1.0*I); -const c minI(0,-1); //Minus complex unity -const d number_pi=arma::datum::pi; // The number 3.14159265359.. - -template -using vdi = arma::Col::fixed; -template -using vci = arma::Col::fixed; - -typedef arma::Col vd; /* Column vector of doubles */ -typedef arma::Col vc; /* Column vector of complex numbers */ -typedef arma::Col vus; /* Column vector of unsigned integers */ -typedef arma::Mat dmat; /* (Dense) Matrix of doubles */ -typedef arma::Mat cmat; /* (Dense) matrix of complex numbers */ - -typedef arma::SpMat sdmat; /* Sparse matrix of doubles */ - - -#endif // TASMET_TYPES_H -////////////////////////////////////////////////////////////////////// - +// tasmet_types.h +// +// Author: J.A. de Jong +// +// Description: +// Typedefs and namespace pollution for stuff that is almost always +// needed. +////////////////////////////////////////////////////////////////////// +#pragma once +#ifndef TASMET_TYPES_H +#define TASMET_TYPES_H +#include + +// We require C++ +#ifndef __cplusplus +#error The c++ compiler should be used. +#endif + +// Some header files I (almost) always need +#include // C++ std vector +#include // C++ string class +#include +#include // Our linear algebra package + +typedef size_t us; /* Size type I always use */ + +#define fillwith arma::fill // Fill vector or matrix with ... + +// To change the whole code to 32-bit floating points, change this to +// float. +#if TASMET_FLOAT == 32 +typedef float d; /* Shortcut for double */ +#elif TASMET_FLOAT == 64 +typedef double d; /* Shortcut for double */ +#else +static_assert(false,"TASMET_FLOAT should be either 32 or 64"); +#endif + +typedef std::complex c; /* Shortcut for c++ complex number */ + +using arma::abs; +using std::string; /* C++ Strings */ +using std::vector; // C++ Vectors + +/* Armadillo functions and objects we often need */ +using arma::sin; +using arma::cos; +using arma::exp; +using arma::ones; +using arma::zeros; +using arma::eye; +using arma::inv; +using arma::det; +using arma::log_det; +using arma::norm; +using arma::span; +using arma::linspace; +using arma::diagmat; + +// Constants +// I need these numbers so often, that they can be in the global namespace + +const c I(0,1); //Complex unity +const c sqI=sqrt(I); +const d sq2=sqrt(2.0); +const c sqmI=sqrt(-1.0*I); +const c minI(0,-1); //Minus complex unity +const d number_pi=arma::datum::pi; // The number 3.14159265359.. + +template +using vdi = arma::Col::fixed; +template +using vci = arma::Col::fixed; + +typedef arma::Col vd; /* Column vector of doubles */ +typedef arma::Col vc; /* Column vector of complex numbers */ +typedef arma::Col vus; /* Column vector of unsigned integers */ +typedef arma::Mat dmat; /* (Dense) Matrix of doubles */ +typedef arma::Mat cmat; /* (Dense) matrix of complex numbers */ + +typedef arma::SpMat sdmat; /* Sparse matrix of doubles */ + + +#endif // TASMET_TYPES_H +////////////////////////////////////////////////////////////////////// +