diff --git a/include/lmscpp/nummatrix.hpp b/include/lmscpp/nummatrix.hpp index efa5f95..765302e 100644 --- a/include/lmscpp/nummatrix.hpp +++ b/include/lmscpp/nummatrix.hpp @@ -8,6 +8,8 @@ #include #include +#include + namespace NumMatrix { using namespace std; @@ -18,10 +20,13 @@ namespace NumMatrix vector> m_; public: NumMatrix(unsigned int rows, unsigned int cols): m_{rows,vector(cols, 0)} {} + NumMatrix(initializer_list> l): m_{l} {}; + + inline unsigned int nrows() const {return m_.size();} inline unsigned int ncols() const {return m_[0].size();} - vector& operator[](unsigned int idx) {return m_[idx];} - const vector& operator[](unsigned int idx) const {return m_[idx];} + vector& operator[](unsigned int idx) {return m_.at(idx);} + const vector& operator[](unsigned int idx) const {return m_.at(idx);} T get(unsigned int i, unsigned int j) const {return (*this)[i][j];} void set(unsigned int i, unsigned int j, T t) {(*this)[i][j] = t;} @@ -96,6 +101,10 @@ namespace NumMatrix return os; } + using nmatrix = NumMatrix; + + nmatrix operator*(double, const nmatrix &); + double max_eigen_value(const nmatrix &); }; diff --git a/include/lmscpp/stochastic.hpp b/include/lmscpp/stochastic.hpp index c31bc3b..07093ac 100644 --- a/include/lmscpp/stochastic.hpp +++ b/include/lmscpp/stochastic.hpp @@ -65,9 +65,9 @@ namespace stochastic{ takes only a function that check if two random variables are independent. */ class ExpectedOperator{ // The function responsible for checking if two random variables are independent. - function is_independent_; + function is_independent_; public: - ExpectedOperator(function is_ind_func); + ExpectedOperator(function is_ind_func); // Returns E(arg) with the following simplifications: // 1. If 'arg' is a number, then just return 'arg'. @@ -76,18 +76,18 @@ namespace stochastic{ // a FunctionSymbol whose name is 'E' and whose argument is 'arg'. RCP operator()(const RCP &arg) const; // Basic can carry a FunctionSymbol or a Pow. It returns whether the two random variables are independent or not. - bool is_independent(const Basic& b1, const Basic& b2) const; + bool is_independent(const Basic& b1, const Basic& b2, int high) const; // NOTE: This should be turned into a generator when C++20 becomes available. // Takes a expression and breaks it into a list of independent random variables. // The last element is the part of the expression that cannot be breaked into independent parts. - vec_basic split_independents(const RCP &expr) const ; + vec_basic split_independents(const RCP &expr, int high) const ; // This method takes a FunctionSymbol in the form E(...) and try to expand it accordingly with // the properties of the expected operator. Example: if it receives E(u(x)**2 + 2*v(x+1) + n(x-1)) and // the second moment of u(x) is γ2 , the first moment of v(x+1) is β and n(x-1) doesn't have a moment then the result of this // method will be γ2 + 2*β + E(n(x-1)). - RCP expand(const RCP &term) const; + RCP expand(const RCP &term, int high) const; }; /* This class is a set in which each element represents a state variable's update equation. @@ -159,6 +159,9 @@ namespace stochastic{ /* This is a help class that can compute the equations, mount both symbolic and numerical matrices, save/load these matrices into/from a cache, etc ...*/ using fallback_func_type = function (const Symbol &x)>; + + // Convert the symbolic DenseMatrix in the numerical nmatrix. + nmatrix sym2num(const DenseMatrix &, const map_basic_basic &inivalsmap, const fallback_func_type fallback); class Experiment { DenseMatrix A_,B_,Yk_; @@ -170,8 +173,6 @@ namespace stochastic{ const fallback_func_type fallback_; // Why I can't put a & here ? public: - // Convert the symbolic DenseMatrix in the numerical nmatrix. - nmatrix sym2num(const DenseMatrix&) const; /* 'inieqs' is the initial equations from the model, 'E' is the expected operator, 'inivalsmap' is a map associating each constant with its value (which can be symbolic or numeric). @@ -192,10 +193,10 @@ namespace stochastic{ inline const DenseMatrix & get_A() const {return A_;}; inline const DenseMatrix & get_Yk() const {return Yk_;}; inline const DenseMatrix & get_B() const {return B_;}; - nmatrix get_num_A() const {return sym2num(A_);}; + nmatrix get_num_A() const {return sym2num(A_, inivalsmap_, fallback_);}; DenseMatrix get_sym_Y0() const; - nmatrix get_num_Y0() const {return sym2num(get_sym_Y0());} - nmatrix get_num_B() const {return sym2num(B_);} + nmatrix get_num_Y0() const {return sym2num(get_sym_Y0(), inivalsmap_, fallback_);} + nmatrix get_num_B() const {return sym2num(B_, inivalsmap_, fallback_);} inline size_t get_number_of_eqs() const {return number_of_eqs_;}; void write_skewness(int niter,RCP randexpr, ofstream & os) const; diff --git a/lib/CMakeLists.txt b/lib/CMakeLists.txt index 41e8841..d1fcd46 100644 --- a/lib/CMakeLists.txt +++ b/lib/CMakeLists.txt @@ -22,7 +22,7 @@ FetchContent_MakeAvailable(bitsery) find_package(SymEngine) find_package(Threads) -add_library(stochastic stochastic.cpp utils.cpp) +add_library(stochastic stochastic.cpp utils.cpp nummatrix.cpp) target_compile_features(stochastic PUBLIC cxx_std_17) target_include_directories(stochastic PUBLIC ../include ${SYMENGINE_INCLUDE_DIRS}) target_link_libraries(stochastic PUBLIC ${SYMENGINE_LIBRARIES} Bitsery::bitsery ${CMAKE_THREAD_LIBS_INIT}) diff --git a/lib/nummatrix.cpp b/lib/nummatrix.cpp new file mode 100644 index 0000000..05b67b9 --- /dev/null +++ b/lib/nummatrix.cpp @@ -0,0 +1,59 @@ +#include + + +namespace NumMatrix { + +nmatrix operator*(double v, const nmatrix &a) { + nmatrix m{a.nrows(), a.ncols()}; + for (auto i = 0; i < a.nrows(); i++) + for (auto j = 0; j < a.ncols(); j++) + m.set(i, j, a.get(i, j) * v); + + return m; +} + +// TODO: Move it to a separate file! +double max_eigen_value(const nmatrix &A) { + if (A.ncols() != A.nrows()) + throw "This is not a square matrix"; + + auto mnorm = [](const auto &z) { + double res{-1.0}; + for (auto r = 0; r < z.nrows(); r++) + res = max(res, abs(z.get(r, 0))); + return res; + }; + + auto inner = [](const auto &a, const auto &b) { + double res{0.0}; + auto n{a.nrows()}; + for (auto i = 0; i < n; i++) + res += a.get(i,0) * b.get(i,0); + return res; + }; + + constexpr double precision{pow(10.0, -5.0)}; + const int n = A.ncols(); + nmatrix y(n, 1); + for (auto it = 0; it < n; it++) + y.set(it, 0, 1); + + double l{0.0}, ll{0.0}; + bool flag; + + do { + + auto z = A * y; + + ll = l; + l = inner(y, z)/inner(y,y); + + y = (1/mnorm(z))*z; + + } while (abs(l - ll) >= precision*abs(ll)); + + + return l; +} + +} // namespace NumMatrix diff --git a/lib/stochastic.cpp b/lib/stochastic.cpp index 010bd90..a5b094b 100644 --- a/lib/stochastic.cpp +++ b/lib/stochastic.cpp @@ -58,7 +58,7 @@ namespace stochastic{ return false; } - ExpectedOperator::ExpectedOperator(function is_ind_func): + ExpectedOperator::ExpectedOperator(function is_ind_func): is_independent_{is_ind_func} {} RCP ExpectedOperator::operator()(const RCP &arg) const{ @@ -74,7 +74,7 @@ namespace stochastic{ return make_rcp("E", arg); } - bool ExpectedOperator::is_independent(const Basic& b1, const Basic& b2) const { + bool ExpectedOperator::is_independent(const Basic& b1, const Basic& b2, int high) const { // NOTE: This will throw a exception if something diferent of a Pow and FunctionSymbol // is passed in. If the Pow's base is not a FunctionSymbol, this will raise a exception too. auto convert = [](const Basic &arg) -> const FunctionSymbol& { @@ -88,10 +88,10 @@ namespace stochastic{ const auto& f1{convert(b1)}; const auto& f2{convert(b2)}; - return is_independent_(f1, f2); + return is_independent_(f1, f2, high); }; - vec_basic ExpectedOperator::split_independents(const RCP &expr) const + vec_basic ExpectedOperator::split_independents(const RCP &expr, int high) const { if (not is_a(*expr)) return {expr}; @@ -110,7 +110,7 @@ namespace stochastic{ bool ind = true; for (auto &y : rvs) { - if (x != y and not is_independent(*x,*y)) + if (x != y and not is_independent(*x,*y, high)) { not_ind = mul(not_ind, x); ind = false; @@ -126,7 +126,7 @@ namespace stochastic{ return res; } - RCP ExpectedOperator::expand(const RCP &term) const + RCP ExpectedOperator::expand(const RCP &term, int high) const { vec_basic addtuple{get_addtuple(SymEngine::expand(term->get_args()[0]))}; auto s{addtuple[0]}; @@ -148,7 +148,7 @@ namespace stochastic{ } vec_basic e{}; - for (auto &t : split_independents(not_cnts)) + for (auto &t : split_independents(not_cnts, high)) e.emplace_back((*this)(t)); auto mul_e = reduce(begin(e), end(e), rcp_static_cast(one), [](auto &a1, auto &a2) @@ -247,13 +247,17 @@ namespace stochastic{ void Experiment::compute(const vec_func & seed) { + using rcp_function_symbol = RCP; EquationSet eqs{inieqs_.get_var()}; - list> s{seed.begin(), seed.end()}; - vector> Y; + list> s; + for (auto &sd:seed) + s.push_back(make_pair(sd, 0)); + + vector Y; while (not s.empty()) { - auto t = s.front(); + auto [t, h] = s.front(); s.pop_front(); if (not eqs.contains(t)) { @@ -262,10 +266,12 @@ namespace stochastic{ // This should not result in a number. Maybe I should do the expansion inside // the operator () of E instead of calling it separately. auto eee = E_(expand(t->get_args()[0], inieqs_)); - auto u = E_.expand(rcp_dynamic_cast(eee)); + auto u = E_.expand(rcp_dynamic_cast(eee), h); auto stvec = states_vars(u); - move(stvec.begin(), stvec.end(), back_inserter(s)); + for (auto &tt: stvec) + s.push_back(move(make_pair(tt, h + 1))); + eqs.setitem(t, u); } } @@ -367,13 +373,15 @@ namespace stochastic{ return s.apply(x); } - nmatrix Experiment::sym2num(const DenseMatrix & m) const + nmatrix sym2num(const DenseMatrix & m, + const map_basic_basic& inivalsmap, + const fallback_func_type fallback) { nmatrix tmp{m.nrows(), m.ncols()}; for (auto i = 0; i < m.nrows(); i++) for (auto j = 0; j < m.ncols(); j++) { - RCP term{ fallxreplace(m.get(i, j),inivalsmap_, fallback_) }; + RCP term{ fallxreplace(m.get(i, j),inivalsmap, fallback) }; if (auto dptr = dynamic_cast(term.get()); dptr != nullptr) tmp[i][j] = dptr->as_double(); diff --git a/models/CMakeLists.txt b/models/CMakeLists.txt index 1f4991e..d813618 100644 --- a/models/CMakeLists.txt +++ b/models/CMakeLists.txt @@ -13,3 +13,7 @@ install(TARGETS classical RUNTIME DESTINATION bin COMPONENT lmscpp) add_executable(skewness skewness.cpp common.cpp) target_link_libraries(skewness PRIVATE stochastic argparse::argparse) install(TARGETS skewness RUNTIME DESTINATION bin COMPONENT lmscpp) + +add_executable(intermediary intermediary.cpp common.cpp) +target_link_libraries(intermediary PRIVATE stochastic argparse::argparse) +install(TARGETS intermediary RUNTIME DESTINATION bin COMPONENT lmscpp) diff --git a/models/classical.cpp b/models/classical.cpp index a2c3b16..5720c8d 100644 --- a/models/classical.cpp +++ b/models/classical.cpp @@ -105,8 +105,8 @@ int main(int argc, char ** argv) return null; }); - ExpectedOperator E{[](const FunctionSymbol& x, const FunctionSymbol& y){ - auto xy = [](const FunctionSymbol &x, const FunctionSymbol &y) + ExpectedOperator E{[](const FunctionSymbol& x, const FunctionSymbol& y, int high){ + auto xy = [](const FunctionSymbol &x, const FunctionSymbol &y, int high) { auto xname{x.get_name()}; auto yname{y.get_name()}; @@ -125,7 +125,7 @@ int main(int argc, char ** argv) return false; }; - return xy(x,y) or xy(y,x); + return xy(x,y, high) or xy(y,x, high); }}; @@ -162,7 +162,8 @@ int main(int argc, char ** argv) for (auto i = 0; i < L; i++) eqs.setitem(V[i](k+one), V[i](k) - step_size*x(k - integer(i)) * inn + step_size*n(k)*x(k - integer(i))); - auto epislonk = E.expand(rcp_dynamic_cast(E(expand(pow(inn, 2_i))))); + // TODO: Default value of 0 should come here. + auto epislonk = E.expand(rcp_dynamic_cast(E(expand(pow(inn, 2_i)))), 0); MeasureDuration duration; Experiment todo{eqs, E, cache, [](const Symbol & x) -> RCP diff --git a/models/intermediary.cpp b/models/intermediary.cpp new file mode 100644 index 0000000..828e7f8 --- /dev/null +++ b/models/intermediary.cpp @@ -0,0 +1,457 @@ +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include + +#include "common.hpp" + +using namespace argparse; +using namespace SymEngine; +using namespace SymEngine::OverloadedOperators; +using namespace std; +using namespace stochastic; + +using namespace chrono; + + +// NOTE: This a redefininition. Maybe I should create a header file with definitions. +using vec_func = vector>; + +enum cachemode + { + IGNORE=0, + READ, + WRITE, + }; +// NOTE: This should be shared with classical.cpp +void read_write_compute(const string& cachename, cachemode cmode, Experiment& todo, const vec_func & seeds) +{ + MeasureDuration duration{}; + + if (cmode == cachemode::READ) + { + fstream fs{cachename, fstream::binary | fstream::in}; + if (fs.is_open()) + { + cout << "Using cached data..." << endl; + duration.reset(); + todo.load(fs); + duration.show("todo.load()"); + } + else + throw runtime_error{"It was not possible to open cache " + cachename + " for reading."}; + } + else + { + duration.reset(); + todo.compute(seeds); + duration.show("todo.compute()"); + + if (cmode == cachemode::WRITE) + { + cout << "Writing cache..." << endl; + fstream fs{cachename, fstream::binary | fstream::out}; + + duration.reset(); + todo.save(fs); + duration.show("todo.save()"); + } + } +} + +enum indmode + { + IA = 0, + EEA, + }; + +enum outmode + { + SK = 0, + MSE, + }; + +enum distmode + { + UNSPECIFIED = 0, + GAUSS, + LAP, + }; + + + + +string latex(const DenseMatrix& d) { + string out{}; + + auto n{d.nrows()}, m{d.ncols()}; + for (auto i = 0; i < n; i++) { + for (auto j = 0; j < m; j++){ + out += latex(*d.get(i, j)) + " &"; + } + out += " \\\\\n"; + } + + return out; +} + +const double binary_precision{pow(10,-3)}; +double binary(function f, double lo, double hi, double ml) { + double m{(lo + hi)/2.0}; + if (abs(m - ml) < binary_precision*abs(ml)) + return m; + + return f(m) ? binary(f, m, hi, m) : binary(f, lo, m, m); +} + +const auto sigmav{ symbol("σ_v") }, step_size{ symbol("β") }, k{ symbol("k") }; + +double get_max_beta(double lo, double hi, const DenseMatrix & m, + map_basic_basic ini, + const fallback_func_type fallback) { + + return binary([&](double beta) { + ini[step_size] = number(beta); + return NumMatrix::max_eigen_value(sym2num(m, ini, fallback)) < 1; + }, lo, hi, lo); +} + +int main(int argc, char ** argv) +{ + argparse::ArgumentParser program(argv[0]); + + program.add_argument("--readcache") + .help("The cache file to use.") + .default_value(""s) + .action([](const string &val){return val;}); + + program.add_argument("--writecache") + .help("The cache file to write.") + .default_value(""s) + .action([](const string &val){return val;}); + + program.add_argument("-N").help("filter length").required().action([](const string &val){return stoi(val);}); + + program.add_argument("-M").help("data length").required().action([](const string &val){return stoi(val);}); + + program.add_argument("--indmode").help("ia or eea").required() + .action([](const string &val) + { + static const map choices {{"ia", indmode::IA}, {"eea", indmode::EEA}}; + auto it = choices.find(val); + if (it == choices.end()) + throw runtime_error{"Invalid value for --iamode"}; + + return it->second; + }); + + program.add_argument("--outmode").help("sk or mse").required() + .action([](const string &val) + { + static const map choices {{"sk", outmode::SK}, {"mse", outmode::MSE}}; + auto it = choices.find(val); + if (it == choices.end()) + throw runtime_error{"Invalid value for --outmode"}; + + return it->second; + }); + + program.add_argument("-b","--beta").help("β").default_value(-1.0) + .action([](const string &val){return stod(val);}); + + program.add_argument("--sv2","--sigmav2").help("variance (σᵥ²)").default_value(-1.0) + .action([](const string &val){return stod(val);}); + + program.add_argument("-n","--niter").help("Number of iterations.").default_value(-1) + .action([](const string &val){return stoi(val);}); + + program.add_argument("-o","--output").help("The file where the output will be stored.").default_value(""s) + .action([](const string &val){return val;}); + + program.add_argument("-d", "--dist").help("gauss or lap") + .default_value(distmode::UNSPECIFIED) + .action([](const string &val) + { + static const map choices {{"gauss", distmode::GAUSS}, {"lap", distmode::LAP}}; + auto it = choices.find(val); + if (it == choices.end()) + throw runtime_error{"Invalid value for --dist"}; + + return it->second; + }); + + program.add_argument("--steady-state") + .help("The steady-state value. Currently only support the skewness.") + .implicit_value(true) + .default_value(false); + + program.add_argument("--latex") + .help("print the matrices to the standard output ") + .implicit_value(true) + .default_value(false); + + program.add_argument("--max-beta").help("get the maximum beta").default_value(0.0) + .action([](const string &val){return stod(val);}); + + // TODO: Maybe instead of an int we should return a uintmax_t ? + program.add_argument("--max-high").help("Max high/depth of the recursion!").default_value(numeric_limits::max()) + .action([](const string &val){return stoi(val);}); + + try { + program.parse_args(argc, argv); + } catch (const runtime_error &err) { + cout << err.what() << endl; + cout << program; + exit(1); + } + + const auto readcache{program.get("--readcache")}; + const auto writecache{program.get("--writecache")}; + const auto N{program.get("-N")}; + const auto M{program.get("-M")}; + const auto beta{program.get("--beta")}; + const auto sigmav2{program.get("--sigmav2")}; + const auto niter{program.get("--niter")}; + const auto ofilename {program.get("--output")}; + const auto ind_mode = program.get("--indmode"); + const auto out_mode = program.get("--outmode"); + const auto dist_mode = program.get("--dist"); + const auto steady_state = program.get("--steady-state"); + const auto print_latex = program.get("--latex"); + const auto max_beta = program.get("--max-beta"); + const auto max_high = program.get("--max-high"); + + + if (not ofilename.empty()) + if (beta < 0 or sigmav2 < 0 or niter < 0 or dist_mode == distmode::UNSPECIFIED) + throw runtime_error{"Missing some of the following: --beta, --sv2, --niter, --dist "}; + + + + const string cachename{ (not readcache.empty()) ? readcache : writecache }; + const auto cache_mode = [&]() + { + if (not readcache.empty()) + return cachemode::READ; + if (not writecache.empty()) + return cachemode::WRITE; + + return cachemode::IGNORE; + }(); + + map_basic_basic cache{ {sigmav, number(sqrt(sigmav2))}, {step_size, number(beta)} }; + + StochasticProcess::clear(); + + StochasticProcess u{"u", [](const vec_basic&, const RCP& n) -> RCP + { + if (n->as_int() % 2 == 0) + return symbol("γ_" + to_string(n->as_int())); + + return zero; + }}; + + StochasticProcess v{"v", [](const vec_basic&, const RCP& n) -> RCP + { + if (even(*n)) + return pow(sigmav, n); + + return zero; + }}; + + vector wtil; + for (auto i = 0; i < N; i++) + wtil.emplace_back("ẅ_" + to_string(i), [](const vec_basic& vec, const RCP& n) -> RCP + { + if (eq(*vec[0], *zero)) + return one; + + return null; + }); + + int recursion_depth = 0; + ExpectedOperator E{[ind_mode, max_high, &recursion_depth](const FunctionSymbol& x, const FunctionSymbol& y, int high){ + recursion_depth = max(recursion_depth, high); + auto local_ind_mode = high <= max_high ? ind_mode : indmode::IA; + auto xy = [local_ind_mode, max_high](const FunctionSymbol &x, const FunctionSymbol &y) + { + auto xname{x.get_name()}; + auto yname{y.get_name()}; + + if (xname == "u" and yname == "u") + return not eq(x,y); + + // NOTE: This is extremely ugly. Maybe there is some better way to do this + // comparasion. + if (xname.find("ẅ") != string::npos and yname == "u") + { + if (local_ind_mode == indmode::EEA) + return not rcp_dynamic_cast(expand(x.get_args()[0] - y.get_args()[0]))->is_positive(); + else if (local_ind_mode == indmode::IA) + return true; + } + + if (xname == "v") + return true; + + return false; + }; + return xy(x,y) or xy(y,x); + }}; + + array bvalues{1., 0.8, 0.8, -0.7, 0.6, -0.5, 0.4, -0.3, 0.2, -0.1}; + if (bvalues.size() < M) + { + cerr << "There size of 'bvalues' is less than " << M << endl; + return 1; + } + + vector> b; + for (auto i = 0; i < M; i++) + { + auto bi { function_symbol("b", integer(i)) }; + b.push_back( rcp_static_cast(bi) ); + cache[bi] = number(bvalues[i]); + } + + auto x = [&](const auto & k) + { + RCP sum{zero}; + // NOTE: This integer could be memorised. + for (auto i = 0; i < M; i++) + sum = sum + b[i] * u(k - integer(i)); + + return sum; + }; + + EquationSet eqs{k}; + + RCP inn{zero}; + for (auto j = 0; j < N; j++) + inn = inn + x(k - integer(j))*wtil[j](k); + + for (auto i = 0; i < N; i++) + eqs.setitem(wtil[i](k+one), wtil[i](k) - step_size*x(k - integer(i)) * inn - step_size*v(k)*x(k - integer(i))); + + const auto fallback = [dist_mode](const Symbol & x) -> RCP + { + auto name = x.__str__(); + if (name.size() >= 4 and name.substr(0, 2) == "γ") + { + auto p{ stoll(name.substr(3)) }; + if (dist_mode == distmode::LAP) + { + static const double scale{ 1.0/sqrt(2) }; + auto res = numfactorial(p)*pow(scale, p); + return number(res); + } + else if (dist_mode == distmode::GAUSS) + { + static constexpr double sd{1.0}; + auto res = pow(sd, p) * num_doublefactorial(p - 1); + return number(res); + } + } + + return null; + }; + + Experiment todo{eqs, E, cache, fallback}; + + RCP mse_expanded{null}; + vec_func seeds{}; + + if (out_mode == outmode::SK) + { + seeds.push_back(rcp_dynamic_cast(E(wtil[0](k)))); + seeds.push_back(rcp_dynamic_cast(E(pow(wtil[0](k), 2_i)))); + seeds.push_back(rcp_dynamic_cast(E(pow(wtil[0](k), 3_i)))); + } + else if (out_mode == outmode::MSE) + { + auto ee = rcp_dynamic_cast( E(expand(pow(inn + v(k), 2_i ))) ); + // TODO: The expand should have a default value of 0! + mse_expanded = E.expand(ee, 0); + seeds = states_vars(mse_expanded); + } + + read_write_compute(cachename, cache_mode, todo, seeds); + cout << "NUMBER_OF_EQUATIONS: " << todo.get_number_of_eqs() << endl; + cout << "RECURSION_DEPTH: " << recursion_depth << endl; + + MeasureDuration duration{}; + if (not ofilename.empty()) + { + ofstream os{ofilename, ofstream::out | ofstream::trunc}; + if (!os) + { + cerr << "It was not possible to open file " << ofilename << endl; + return 1; + } + + if (out_mode == outmode::SK) + { + duration.reset(); + todo.write_skewness(niter, wtil[0](k), os); + duration.show("todo.write_skewness()"); + } + else if (out_mode == outmode::MSE) + { + duration.reset(); + todo.write_expression(niter, mse_expanded, os); + duration.show("todo.write_expression()"); + } + } + + if (print_latex) + { + cout << "Printing latex...\n" << endl; + cout << "Matrix A: " << endl; + cout << latex(todo.get_A()) << endl; + + cout << "Matrix Yk: " << endl; + cout << latex(todo.get_Yk()) << endl; + + cout << "Matrix B: " << endl; + cout << latex(todo.get_B()) << endl; + + cout << "Printing numeric matrix...\n" << endl; + + cout << "Matrix A: " << endl; + cout << todo.get_num_A() << endl; + cout << "Matrix B: " << endl; + cout << todo.get_num_B() << endl; + } + + if (steady_state) + { + if (out_mode == outmode::SK) { + auto p = todo.skewness_steady_state(wtil[0](k)); + cout << "steady-state: " << p.second << endl; + cerr << "Steady-state found with " << p.first << " iterations." << endl; + } + else + throw runtime_error{"The steady-state value of MSE is currently not supported!"}; + } + + if (max_beta) { + cout << "beta max = " << get_max_beta(0.0, max_beta, todo.get_A(), cache, fallback) << endl; + } + + return 0; +} diff --git a/models/skewness.cpp b/models/skewness.cpp index 577229e..d73262a 100644 --- a/models/skewness.cpp +++ b/models/skewness.cpp @@ -254,8 +254,8 @@ int main(int argc, char ** argv) return null; }); - ExpectedOperator E{[ind_mode](const FunctionSymbol& x, const FunctionSymbol& y){ - auto xy = [ind_mode](const FunctionSymbol &x, const FunctionSymbol &y) + ExpectedOperator E{[ind_mode](const FunctionSymbol& x, const FunctionSymbol& y, int high){ + auto xy = [ind_mode](const FunctionSymbol &x, const FunctionSymbol &y, int high) { auto xname{x.get_name()}; auto yname{y.get_name()}; @@ -278,7 +278,7 @@ int main(int argc, char ** argv) return false; }; - return xy(x,y) or xy(y,x); + return xy(x,y, high) or xy(y,x, high); }}; array bvalues{1., 0.8, 0.8, -0.7, 0.6, -0.5, 0.4, -0.3, 0.2, -0.1}; @@ -350,7 +350,8 @@ int main(int argc, char ** argv) else if (out_mode == outmode::MSE) { auto ee = rcp_dynamic_cast( E(expand(pow(inn + v(k), 2_i ))) ); - mse_expanded = E.expand(ee); + // TODO: The expand should have a default value of 0! + mse_expanded = E.expand(ee, 0); seeds = states_vars(mse_expanded); } diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 09e088b..46c731a 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -5,3 +5,7 @@ add_test(NAME testall COMMAND testall) add_executable(testnode testnode.cpp) target_link_libraries(testnode PRIVATE stochastic) add_test(NAME testnode COMMAND testnode) + +add_executable(testmatrix testmatrix.cpp) +target_link_libraries(testmatrix PRIVATE stochastic) +add_test(NAME testmatrix COMMAND testmatrix) diff --git a/tests/testall.cpp b/tests/testall.cpp index 9e54229..c895609 100644 --- a/tests/testall.cpp +++ b/tests/testall.cpp @@ -83,8 +83,8 @@ void test_get_addtuple(){ void test_ExptectOperator_and_EquationSet(){ // u is independent of v but v is not independent of w. w and u are independent. - const ExpectedOperator E{[](const FunctionSymbol& x, const FunctionSymbol& y){ - auto xy = [](const FunctionSymbol &x, const FunctionSymbol &y) + const ExpectedOperator E{[](const FunctionSymbol& x, const FunctionSymbol& y, int high){ + auto xy = [](const FunctionSymbol &x, const FunctionSymbol &y, int high) { auto xname{x.get_name()}; auto yname{y.get_name()}; @@ -97,7 +97,7 @@ void test_ExptectOperator_and_EquationSet(){ return false; }; - return xy(x,y) or xy(y,x); + return xy(x,y, high) or xy(y,x, high); }}; const auto gamma{symbol("γ")}, sigma{symbol("σ")}, alpha{symbol("α")}; diff --git a/tests/testmatrix.cpp b/tests/testmatrix.cpp new file mode 100644 index 0000000..6997945 --- /dev/null +++ b/tests/testmatrix.cpp @@ -0,0 +1,15 @@ +#include +#include + + +using namespace std; +using nmatrix = NumMatrix::NumMatrix; + +int main() { + nmatrix a({{3,0,1},{2,2,2},{4,2,5}}); + nmatrix b({{1,0,0},{0,-2,0},{0,0,1.5}}); + + cout << NumMatrix::max_eigen_value(a) << endl; + cout << NumMatrix::max_eigen_value(b) << endl; + return 0; +}