From a9082f47c2c04ade424cacca685fbb2569fb50b7 Mon Sep 17 00:00:00 2001 From: Jean-Paul Pelteret Date: Sat, 30 Mar 2019 11:20:34 +0100 Subject: [PATCH] Add test for SD::Expression class --- tests/symengine/symengine_number_type_01.cc | 1142 +++++++++++++++++ .../symengine/symengine_number_type_01.output | 17 + 2 files changed, 1159 insertions(+) create mode 100644 tests/symengine/symengine_number_type_01.cc create mode 100644 tests/symengine/symengine_number_type_01.output diff --git a/tests/symengine/symengine_number_type_01.cc b/tests/symengine/symengine_number_type_01.cc new file mode 100644 index 0000000000..4b26361e1f --- /dev/null +++ b/tests/symengine/symengine_number_type_01.cc @@ -0,0 +1,1142 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2019 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Check that the wrapper for symengine expressions works as expected + +#include + +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include "../serialization/serialization.h" +#include "../tests.h" + +using namespace dealii; +namespace SD = Differentiation::SD; +namespace SE = ::SymEngine; + + +// Overload the comparison function in serialization.h +template <> +bool +compare(const SD::Expression &t1, const SD::Expression &t2) +{ + // Could also do the following if it ever becomes necessary: + // (t1.get_value().__cmp__(t2.get_value()) == 0) + return t1.get_expression() == t2.get_expression(); +} + + +int +main() +{ + initlog(); + + using SD_number_t = SD::Expression; + + deallog << "Constructors" << std::endl; + { + std::cout << "Constructor: Boolean" << std::endl; + { + const SD_number_t se_bool_true(true); + const SD_number_t se_bool_false(false); + Assert(SE::eq(se_bool_true.get_value(), *SE::boolTrue), + ExcMessage("Problem with constructor")); + Assert(SE::eq(se_bool_false.get_value(), *SE::boolFalse), + ExcMessage("Problem with constructor")); + } + + std::cout << "Constructor: Integer" << std::endl; + { + const int x = -1; + const SD_number_t se_number(x); + Assert(SE::eq(se_number.get_value(), *SE::integer(x)), + ExcMessage("Problem with constructor")); + } + + std::cout << "Constructor: Unsigned integer" << std::endl; + { + const unsigned int x = 2; + const SD_number_t se_number(x); + Assert(SE::eq(se_number.get_value(), *SE::integer(x)), + ExcMessage("Problem with constructor")); + } + + std::cout << "Constructor: Rational" << std::endl; + { + const unsigned int x = 2; + const unsigned int y = 3; + const SD_number_t se_number(x, y); + Assert(SE::eq(se_number.get_value(), *SE::rational(x, y)), + ExcMessage("Problem with constructor")); + } + + std::cout << "Constructor: Double" << std::endl; + { + const double x = 2.5; + const SD_number_t se_number(x); + Assert(SE::eq(se_number.get_value(), *SE::real_double(x)), + ExcMessage("Problem with constructor")); + } + + std::cout << "Constructor: Complex double" << std::endl; + { + const std::complex x(2.5, 1.5); + const SD_number_t se_number(x); + Assert(SE::eq(se_number.get_value(), *SE::complex_double(x)), + ExcMessage("Problem with constructor")); + } + + std::cout << "Constructor: Float" << std::endl; + { + const float x = 2.5; + const SD_number_t se_number(x); + Assert(SE::eq(se_number.get_value(), *SE::real_double(x)), + ExcMessage("Problem with constructor")); + } + + std::cout << "Constructor: Complex float" << std::endl; + { + const std::complex x(2.5, 1.5); + const SD_number_t se_number(x); + Assert(SE::eq(se_number.get_value(), *SE::complex_double(x)), + ExcMessage("Problem with constructor")); + } + + std::cout << "Constructor: Symbol" << std::endl; + { + const std::string f = "f"; + const bool parse_as_expression = false; + const SD_number_t se_number(f, parse_as_expression); + Assert(SE::eq(se_number.get_value(), *SE::symbol(f)), + ExcMessage("Problem with constructor")); + } + + std::cout << "Constructor: Symbol" << std::endl; + { + const std::string f = "func"; + const SD_number_t se_number(f.c_str()); + Assert(SE::eq(se_number.get_value(), *SE::symbol(f)), + ExcMessage("Problem with constructor")); + } + + std::cout << "Constructor: Symbol" << std::endl; + { + const std::string f = "func"; + const bool parse_as_expression = false; + const SD_number_t se_number(f, parse_as_expression); + Assert(SE::eq(se_number.get_value(), *SE::symbol(f)), + ExcMessage("Problem with constructor")); + } + + std::cout << "Constructor: Symbolic expression" << std::endl; + { + const std::string f = "x + y"; + const bool parse_as_expression = true; + const SD_number_t se_number(f, parse_as_expression); + Assert(SE::eq(se_number.get_value(), *SE::parse(f)), + ExcMessage("Problem with constructor")); + } + + std::cout << "Copy constructor" << std::endl; + { + const double x = 2.5; + const SD_number_t se_number_1(x); + const SD_number_t se_number_2(se_number_1); + Assert(SE::eq(se_number_2.get_value(), *SE::real_double(x)), + ExcMessage("Problem with constructor")); + } + } + + deallog << "Utility functions" << std::endl; + { + std::cout << "Parse" << std::endl; + { + const double x = 2.5; + SD_number_t se_number(x); + const std::string f = "x + y"; + se_number.parse(f); + Assert(SE::eq(se_number.get_value(), *SE::parse(f)), + ExcMessage("Problem with utility function: Parse")); + } + + std::cout << "Save/load" << std::endl; + { + // Output + const SD_number_t se_number_x_out("x"); + const SD_number_t se_number_y_out("y"); + const SD_number_t se_number_f_out("x + y + 2.5", + true /*parse_as_expression*/); + { + std::ofstream o_str("save_load.txt"); + Assert(o_str.is_open(), ExcInternalError()); + se_number_x_out.save(o_str); + se_number_y_out.save(o_str); + se_number_f_out.save(o_str); + o_str.close(); + } + + // Input + SD_number_t se_number_x_in; + SD_number_t se_number_y_in; + SD_number_t se_number_f_in; + { + std::ifstream i_str("save_load.txt"); + Assert(i_str.is_open(), ExcInternalError()); + se_number_x_in.load(i_str); + se_number_y_in.load(i_str); + se_number_f_in.load(i_str); + } + + Assert(static_cast(se_number_x_in == se_number_x_out), + ExcMessage("Problem with utility function: Save/load")); + Assert(static_cast(se_number_y_in == se_number_y_out), + ExcMessage("Problem with utility function: Save/load")); + Assert(static_cast(se_number_f_in == se_number_f_out), + ExcMessage("Problem with utility function: Save/load")); + } + + deallog.push("Serialization"); + std::cout << "Serialization" << std::endl; + { + const SD_number_t se_number_x_out("x"); + const SD_number_t se_number_y_out("y"); + const SD_number_t se_number_f_out("x + y + 2.5", + true /*parse_as_expression*/); + + // From serialization.h + SD_number_t se_number_x_in; + SD_number_t se_number_y_in; + SD_number_t se_number_f_in; + verify(se_number_x_out, se_number_x_in); + verify(se_number_y_out, se_number_y_in); + verify(se_number_f_out, se_number_f_in); + + Assert(static_cast(se_number_x_in == se_number_x_out), + ExcMessage("Problem with utility function: Serialization")); + Assert(static_cast(se_number_y_in == se_number_y_out), + ExcMessage("Problem with utility function: Serialization")); + Assert(static_cast(se_number_f_in == se_number_f_out), + ExcMessage("Problem with utility function: Serialization")); + } + deallog.pop(); + } + + deallog << "Comparison operators" << std::endl; + { + std::cout << "Equality (true)" << std::endl; + { + const int x = 2; + const int y = 2; + const SD_number_t se_number_1(x); + const SD_number_t se_number_2(y); + Assert(static_cast(se_number_1 == se_number_2) == true, + ExcMessage("Problem with comparison operator")); + } + + std::cout << "Equality (false)" << std::endl; + { + const int x = 2; + const int y = 3; + const SD_number_t se_number_1(x); + const SD_number_t se_number_2(y); + Assert(static_cast(se_number_1 == se_number_2) == false, + ExcMessage("Problem with comparison operator")); + } + + std::cout << "Non-equality (true)" << std::endl; + { + const int x = 2; + const int y = 3; + const SD_number_t se_number_1(x); + const SD_number_t se_number_2(y); + Assert(static_cast(se_number_1 != se_number_2) == true, + ExcMessage("Problem with comparison operator")); + } + + std::cout << "Non-equality (false)" << std::endl; + { + const int x = 2; + const int y = 2; + const SD_number_t se_number_1(x); + const SD_number_t se_number_2(y); + Assert(static_cast(se_number_1 != se_number_2) == false, + ExcMessage("Problem with comparison operator")); + } + + std::cout << "Less-than operator (true)" << std::endl; + { + const int x = 2; + const int y = 3; + const SD_number_t se_number_1(x); + const SD_number_t se_number_2(y); + Assert(static_cast(se_number_1 < se_number_2) == true, + ExcMessage("Problem with comparison operator")); + } + + std::cout << "Less-than operator (false)" << std::endl; + { + const int x = 3; + const int y = 2; + const int z = 3; + const SD_number_t se_number_1(x); + const SD_number_t se_number_2(y); + const SD_number_t se_number_3(z); + Assert(static_cast(se_number_1 < se_number_2) == false, + ExcMessage("Problem with comparison operator")); + Assert(static_cast(se_number_1 < se_number_3) == false, + ExcMessage("Problem with comparison operator")); + } + + std::cout << "Less-than-or-equals operator (true)" << std::endl; + { + const int x = 2; + const int y = 3; + const int z = 2; + const SD_number_t se_number_1(x); + const SD_number_t se_number_2(y); + const SD_number_t se_number_3(z); + Assert(static_cast(se_number_1 <= se_number_2) == true, + ExcMessage("Problem with comparison operator")); + Assert(static_cast(se_number_1 <= se_number_3) == true, + ExcMessage("Problem with comparison operator")); + } + + std::cout << "Less-than-or-equals operator (false)" << std::endl; + { + const int x = 3; + const int y = 2; + const SD_number_t se_number_1(x); + const SD_number_t se_number_2(y); + Assert(static_cast(se_number_1 <= se_number_2) == false, + ExcMessage("Problem with comparison operator")); + } + + std::cout << "Greater-than operator (true)" << std::endl; + { + const int x = 3; + const int y = 2; + const SD_number_t se_number_1(x); + const SD_number_t se_number_2(y); + Assert(static_cast(se_number_1 > se_number_2) == true, + ExcMessage("Problem with comparison operator")); + } + + std::cout << "Greater-than operator (false)" << std::endl; + { + const int x = 2; + const int y = 3; + const int z = 2; + const SD_number_t se_number_1(x); + const SD_number_t se_number_2(y); + const SD_number_t se_number_3(z); + Assert(static_cast(se_number_1 > se_number_2) == false, + ExcMessage("Problem with comparison operator")); + Assert(static_cast(se_number_1 > se_number_3) == false, + ExcMessage("Problem with comparison operator")); + } + + std::cout << "Greater-than-or-equals operator (true)" << std::endl; + { + const int x = 3; + const int y = 2; + const int z = 3; + const SD_number_t se_number_1(x); + const SD_number_t se_number_2(y); + const SD_number_t se_number_3(z); + Assert(static_cast(se_number_1 >= se_number_2) == true, + ExcMessage("Problem with comparison operator")); + Assert(static_cast(se_number_1 >= se_number_3) == true, + ExcMessage("Problem with comparison operator")); + } + + std::cout << "Greater-than-or-equals operator (false)" << std::endl; + { + const int x = 2; + const int y = 3; + const SD_number_t se_number_1(x); + const SD_number_t se_number_2(y); + Assert(static_cast(se_number_1 >= se_number_2) == false, + ExcMessage("Problem with comparison operator")); + } + } + + deallog << "Logical operators" << std::endl; + { + std::cout << "Logical not" << std::endl; + { + const bool x = true; + const SD_number_t se_number(x); + Assert(static_cast(!se_number) == (!x), + ExcMessage("Problem with logical operator")); + Assert(static_cast(!!se_number) == (!!x), + ExcMessage("Problem with logical operator")); + } + + std::cout << "Logical and" << std::endl; + { + auto test_logical_op = [](const bool x, const bool y) { + const SD_number_t se_number_1(x); + const SD_number_t se_number_2(y); + Assert(static_cast(se_number_1 & se_number_2) == (x & y), + ExcMessage("Problem with logical operator")); + }; + + test_logical_op(true, true); + test_logical_op(true, false); + test_logical_op(false, true); + test_logical_op(false, false); + } + + std::cout << "Logical inclusive or" << std::endl; + { + auto test_logical_op = [](const bool x, const bool y) { + const SD_number_t se_number_1(x); + const SD_number_t se_number_2(y); + Assert(static_cast(se_number_1 | se_number_2) == (x | y), + ExcMessage("Problem with logical operator")); + }; + + test_logical_op(true, true); + test_logical_op(true, false); + test_logical_op(false, true); + test_logical_op(false, false); + } + + std::cout << "Logical xor" << std::endl; + { + auto test_logical_op = [](const bool x, const bool y) { + const SD_number_t se_number_1(x); + const SD_number_t se_number_2(y); + Assert(static_cast(se_number_1 ^ se_number_2) == (x ^ y), + ExcMessage("Problem with logical operator")); + }; + + test_logical_op(true, true); + test_logical_op(true, false); + test_logical_op(false, true); + test_logical_op(false, false); + } + + std::cout << "And" << std::endl; + { + auto test_logical_op = [](const bool x, const bool y) { + const SD_number_t se_number_1(x); + const SD_number_t se_number_2(y); + Assert(static_cast(se_number_1 && se_number_2) == (x && y), + ExcMessage("Problem with logical operator")); + }; + + test_logical_op(true, true); + test_logical_op(true, false); + test_logical_op(false, true); + test_logical_op(false, false); + } + + std::cout << "Inclusive or" << std::endl; + { + auto test_logical_op = [](const bool x, const bool y) { + const SD_number_t se_number_1(x); + const SD_number_t se_number_2(y); + Assert(static_cast(se_number_1 || se_number_2) == (x || y), + ExcMessage("Problem with logical operator")); + }; + + test_logical_op(true, true); + test_logical_op(true, false); + test_logical_op(false, true); + test_logical_op(false, false); + } + } + + deallog << "Assignment operators" << std::endl; + { + std::cout << "Assignment operator (SD_number_t)" << std::endl; + { + const int x = 1; + const int y = 2; + SD_number_t se_number_1(x); + const SD_number_t se_number_2(y); + se_number_1 = se_number_2; + Assert(SE::eq(se_number_1.get_value(), *SE::integer(y)), + ExcMessage("Problem with assignment operator")); + } + + std::cout << "Assignment operator (arithmetic type)" << std::endl; + { + const int x = 1; + const double y = 2; + SD_number_t se_number(x); + se_number = y; + Assert(SE::eq(se_number.get_value(), *SE::real_double(y)), + ExcMessage("Problem with assignment operator")); + } + + std::cout << "Addition assignment (SD_number_t)" << std::endl; + { + const int x = 1; + const int y = 2; + SD_number_t se_number_1(x); + const SD_number_t se_number_2(y); + se_number_1 += se_number_2; + Assert(SE::eq(se_number_1.get_value(), *SE::integer(x + y)), + ExcMessage("Problem with assignment operator")); + } + + std::cout << "Addition assignment (arithmetic type; same type)" + << std::endl; + { + const int x = 1; + const int y = 2; + SD_number_t se_number(x); + se_number += y; + Assert(SE::eq(se_number.get_value(), *SE::integer(x + y)), + ExcMessage("Problem with assignment operator")); + } + + std::cout << "Addition assignment (arithmetic type; different type)" + << std::endl; + { + const int x = 1; + const double y = 2; + SD_number_t se_number(x); + se_number += y; + Assert(SE::eq(se_number.get_value(), *SE::real_double(x + y)), + ExcMessage("Problem with assignment operator")); + } + + std::cout << "Subtraction assignment (SD_number_t)" << std::endl; + { + const int x = 1; + const int y = 2; + SD_number_t se_number_1(x); + const SD_number_t se_number_2(y); + se_number_1 -= se_number_2; + Assert(SE::eq(se_number_1.get_value(), *SE::integer(x - y)), + ExcMessage("Problem with assignment operator")); + } + + std::cout << "Subtraction assignment (arithmetic type; same type)" + << std::endl; + { + const int x = 1; + const int y = 2; + SD_number_t se_number(x); + se_number -= y; + Assert(SE::eq(se_number.get_value(), *SE::integer(x - y)), + ExcMessage("Problem with assignment operator")); + } + + std::cout << "Subtraction assignment (arithmetic type; different type)" + << std::endl; + { + const int x = 1; + const double y = 2; + SD_number_t se_number(x); + se_number -= y; + Assert(SE::eq(se_number.get_value(), *SE::real_double(x - y)), + ExcMessage("Problem with assignment operator")); + } + + std::cout << "Multiplication assignment (SD_number_t)" << std::endl; + { + const int x = 2; + const int y = 3; + SD_number_t se_number_1(x); + const SD_number_t se_number_2(y); + se_number_1 *= se_number_2; + Assert(SE::eq(se_number_1.get_value(), *SE::integer(x * y)), + ExcMessage("Problem with assignment operator")); + } + + std::cout << "Multiplication assignment (arithmetic type; same type)" + << std::endl; + { + const int x = 2; + const int y = 3; + SD_number_t se_number(x); + se_number *= y; + Assert(SE::eq(se_number.get_value(), *SE::integer(x * y)), + ExcMessage("Problem with assignment operator")); + } + + std::cout << "Multiplication assignment (arithmetic type; different type)" + << std::endl; + { + const int x = 2; + const double y = 3; + SD_number_t se_number(x); + se_number *= y; + Assert(SE::eq(se_number.get_value(), *SE::real_double(x * y)), + ExcMessage("Problem with assignment operator")); + } + + std::cout << "Division assignment (SD_number_t)" << std::endl; + { + const int x = 8; + const int y = 2; + SD_number_t se_number_1(x); + const SD_number_t se_number_2(y); + se_number_1 /= se_number_2; + Assert(SE::eq(se_number_1.get_value(), *SE::integer(x / y)), + ExcMessage("Problem with assignment operator")); + } + + std::cout << "Division assignment (arithmetic type; same type)" + << std::endl; + { + const int x = 8; + const int y = 2; + SD_number_t se_number(x); + se_number /= y; + Assert(SE::eq(se_number.get_value(), *SE::integer(x / y)), + ExcMessage("Problem with assignment operator")); + } + + std::cout << "Division assignment (arithmetic type; different type)" + << std::endl; + { + const int x = 8; + const double y = 2; + SD_number_t se_number(x); + se_number /= y; + Assert(SE::eq(se_number.get_value(), *SE::real_double(x / y)), + ExcMessage("Problem with assignment operator")); + } + } + + deallog << "Math operators" << std::endl; + { + std::cout << "Addition" << std::endl; + { + const int x = 1; + const int y = 2; + const SD_number_t se_number_1(x); + const SD_number_t se_number_2(y); + const SD_number_t se_number_3 = se_number_1 + se_number_2; + Assert(SE::eq(se_number_3.get_value(), *SE::integer(x + y)), + ExcMessage("Problem with math operator")); + } + + std::cout << "Subtraction" << std::endl; + { + const int x = 1; + const int y = 2; + const SD_number_t se_number_1(x); + const SD_number_t se_number_2(y); + const SD_number_t se_number_3 = se_number_1 - se_number_2; + Assert(SE::eq(se_number_3.get_value(), *SE::integer(x - y)), + ExcMessage("Problem with math operator")); + } + + std::cout << "Multiplication" << std::endl; + { + const int x = 2; + const int y = 3; + const SD_number_t se_number_1(x); + const SD_number_t se_number_2(y); + const SD_number_t se_number_3 = se_number_1 * se_number_2; + Assert(SE::eq(se_number_3.get_value(), *SE::integer(x * y)), + ExcMessage("Problem with math operator")); + } + + std::cout << "Division" << std::endl; + { + const int x = 8; + const int y = 2; + const SD_number_t se_number_1(x); + const SD_number_t se_number_2(y); + const SD_number_t se_number_3 = se_number_1 / se_number_2; + Assert(SE::eq(se_number_3.get_value(), *SE::integer(x / y)), + ExcMessage("Problem with math operator")); + } + } + + deallog << "Differentiation" << std::endl; + { + std::cout << "Basic differentiation (SD_number_t)" << std::endl; + { + const SD_number_t x("x"); + SD_number_t f(1.0); + f *= x; + f *= x; // f = x^2 + const SD_number_t df_dx = f.differentiate(x); + Assert(static_cast(df_dx == (2.0 * x)), + ExcMessage("Problem with differentiaton")); + } + + std::cout << "Basic differentiation (Symbol)" << std::endl; + { + const auto symb_x = SE::symbol("x"); + const SD_number_t x(symb_x); + SD_number_t f(1.0); + f *= x; + f *= x; // f = x^2 + const SD_number_t df_dx = f.differentiate(symb_x); + Assert(static_cast(df_dx == SD_number_t(2.0 * x)), + ExcMessage("Problem with differentiaton")); + } + + std::cout << "Basic differentiation (Basic representing a symbolic type)" + << std::endl; + { + const auto symb_x_plus_y = SE::add(SE::symbol("x"), SE::symbol("y")); + const SD_number_t x_plus_y(symb_x_plus_y); + SD_number_t f(1.0); + f *= x_plus_y; + f *= x_plus_y; // f = (x+y)^2 + const SD_number_t df_dxpy = f.differentiate(symb_x_plus_y); + Assert(static_cast(df_dxpy == (2.0 * SD_number_t(symb_x_plus_y))), + ExcMessage("Problem with differentiaton")); + } + } + + deallog << "Substitution" << std::endl; + { + std::cout << "Basic substitution" << std::endl; + { + const SD_number_t x("x"); + const SD_number_t y("y"); + const SD_number_t x_plus_y = x + y; + const SD_number_t f = x_plus_y * x_plus_y; + Assert(static_cast(f == SD_number_t("(x + y)**2", true)), + ExcMessage("Problem with substitution")); + + const double x_val = 1.0; + const SD_number_t f_sub_x = f.substitute(x, SD_number_t(x_val)); + Assert(static_cast(f_sub_x == SD_number_t("(1.0 + y)**2", true)), + ExcMessage("Problem with substitution")); + + const double y_val = 2.0; + const SD_number_t f_sub_xy = f_sub_x.substitute(y, SD_number_t(y_val)); + Assert(static_cast(f_sub_xy) == + (x_val * x_val + 2.0 * x_val * y_val + y_val * y_val), + ExcMessage("Problem with substitution")); + } + + std::cout << "Templated-deduced substitution" << std::endl; + { + const SD_number_t x("x"); + const SD_number_t y("y"); + const SD_number_t x_plus_y = x + y; + const SD_number_t f = x_plus_y * x_plus_y; + Assert(static_cast(f == SD_number_t("(x + y)**2", true)), + ExcMessage("Problem with substitution")); + + const double x_val = 1.0; + const SD_number_t f_sub_x = f.substitute(x, x_val); + Assert(static_cast(f_sub_x == SD_number_t("(1.0 + y)**2", true)), + ExcMessage("Problem with substitution")); + + const double y_val = 2.0; + const SD_number_t f_sub_xy = f_sub_x.substitute(y, y_val); + Assert(static_cast(f_sub_xy) == + (x_val * x_val + 2.0 * x_val * y_val + y_val * y_val), + ExcMessage("Problem with substitution")); + } + + std::cout << "Map substitution" << std::endl; + { + const SD_number_t x("x"); + const SD_number_t y("y"); + const SD_number_t x_plus_y = x + y; + const SD_number_t f = x_plus_y * x_plus_y; + + const double x_val = 1.0; + const double y_val = 2.0; + SD::types::substitution_map sub_map; + sub_map[x] = SD_number_t(x_val); + sub_map[y] = SD_number_t(y_val); + const SD_number_t f_sub_xy = f.substitute(sub_map); + Assert(static_cast(f_sub_xy) == + (x_val * x_val + 2.0 * x_val * y_val + y_val * y_val), + ExcMessage("Problem with substitution")); + } + } + + deallog << "Evaluation" << std::endl; + { + std::cout << "Substitution with evaluation (integer)" << std::endl; + { + const SD_number_t x("x"); + const SD_number_t y("y"); + const SD_number_t x_plus_y = x + y; + const SD_number_t f = x_plus_y * x_plus_y; + + const int x_val = 1; + const int y_val = 2; + SD::types::substitution_map sub_map; + sub_map[x] = SD_number_t(x_val); + sub_map[y] = SD_number_t(y_val); + const int f_val = f.substitute_and_evaluate(sub_map); + Assert(f_val == (x_val * x_val + 2 * x_val * y_val + y_val * y_val), + ExcMessage("Problem with evaluation")); + } + + std::cout << "Substitution with evaluation (double)" << std::endl; + { + const SD_number_t x("x"); + const SD_number_t y("y"); + const SD_number_t x_plus_y = x + y; + const SD_number_t f = x_plus_y * x_plus_y; + + const double x_val = 1.0; + const double y_val = 2.0; + SD::types::substitution_map sub_map; + sub_map[x] = SD_number_t(x_val); + sub_map[y] = SD_number_t(y_val); + const double f_val = f.substitute_and_evaluate(sub_map); + Assert(f_val == (x_val * x_val + 2.0 * x_val * y_val + y_val * y_val), + ExcMessage("Problem with evaluation")); + } + + std::cout << "Substitution with evaluation (float)" << std::endl; + { + const SD_number_t x("x"); + const SD_number_t y("y"); + const SD_number_t x_plus_y = x + y; + const SD_number_t f = x_plus_y * x_plus_y; + + const float x_val = 1.0; + const float y_val = 2.0; + SD::types::substitution_map sub_map; + sub_map[x] = SD_number_t(x_val); + sub_map[y] = SD_number_t(y_val); + const float f_val = f.substitute_and_evaluate(sub_map); + Assert(f_val == (x_val * x_val + 2.0f * x_val * y_val + y_val * y_val), + ExcMessage("Problem with evaluation")); + } + + std::cout << "Substitution with evaluation (complex double)" << std::endl; + { + const SD_number_t x("x"); + const SD_number_t y("y"); + const SD_number_t x_plus_y = x + y; + const SD_number_t f = x_plus_y * x_plus_y; + + const std::complex x_val(1.0, 3.0); + const std::complex y_val(2.0, -2.0); + SD::types::substitution_map sub_map; + sub_map[x] = SD_number_t(x_val); + sub_map[y] = SD_number_t(y_val); + const std::complex f_val = + f.substitute_and_evaluate>(sub_map); + Assert(std::abs(f_val - (x_val * x_val + 2.0 * x_val * y_val + + y_val * y_val)) < 1e-12, + ExcMessage("Problem with evaluation")); + } + + std::cout << "Substitution with evaluation (complex float)" << std::endl; + { + const SD_number_t x("x"); + const SD_number_t y("y"); + const SD_number_t x_plus_y = x + y; + const SD_number_t f = x_plus_y * x_plus_y; + + const std::complex x_val(1.0, 3.0); + const std::complex y_val(2.0, -2.0); + SD::types::substitution_map sub_map; + sub_map[x] = SD_number_t(x_val); + sub_map[y] = SD_number_t(y_val); + const std::complex f_val = + f.substitute_and_evaluate>(sub_map); + Assert(std::abs(f_val - (x_val * x_val + 2.0f * x_val * y_val + + y_val * y_val)) < 1e-12, + ExcMessage("Problem with evaluation")); + } + + std::cout << "Substitution with evaluation (mixed arithmetic types)" + << std::endl; + { + const SD_number_t x("x"); + const SD_number_t y("y"); + const SD_number_t x_plus_y = x + y; + const SD_number_t f = x_plus_y * x_plus_y; + + const int x_val = 1; + const int y_val = 2; + SD::types::substitution_map sub_map; + sub_map[x] = SD_number_t(x_val); + sub_map[y] = SD_number_t(y_val); + const double f_val = f.substitute_and_evaluate(sub_map); + Assert(f_val == static_cast(x_val * x_val + 2.0f * x_val * y_val + + y_val * y_val), + ExcMessage("Problem with evaluation")); + } + + std::cout << "Substitution with evaluation (mixed complex types)" + << std::endl; + { + const SD_number_t x("x"); + const SD_number_t y("y"); + const SD_number_t x_plus_y = x + y; + const SD_number_t f = x_plus_y * x_plus_y; + + const std::complex x_val(1.0, 3.0); + const std::complex y_val(2.0, -2.0); + SD::types::substitution_map sub_map; + sub_map[x] = SD_number_t(x_val); + sub_map[y] = SD_number_t(y_val); + const std::complex f_val = + f.substitute_and_evaluate>(sub_map); + Assert(std::abs(f_val - static_cast>( + x_val * x_val + 2.0f * x_val * y_val + + y_val * y_val)) < 1e-12, + ExcMessage("Problem with evaluation")); + } + + std::cout << "Substitution with evaluation (single piecewise function)" + << std::endl; + { + const SD_number_t x("x"); + const SD_number_t y("y"); + const SD_number_t x_plus_y = x + y; + const SD_number_t x_minus_y = x - y; + const SD_number_t f_plus = x_plus_y * x_plus_y; + const SD_number_t f_minus = x_minus_y * x_minus_y; + const SD_number_t f((x > SD_number_t(0.0)), f_plus, f_minus); + + // Condition is met + { + const double x_val = 1.0; + const double y_val = 2.0; + SD::types::substitution_map sub_map; + sub_map[x] = SD_number_t(x_val); + sub_map[y] = SD_number_t(y_val); + const SD_number_t f_subs = f.substitute(sub_map); + const double f_val = static_cast(f_subs); + Assert(f_val == (x_val * x_val + 2.0 * x_val * y_val + y_val * y_val), + ExcMessage( + "Problem with evaluation of single piecewise function")); + } + + // Condition is not met + { + const double x_val = -1.0; + const double y_val = 2.0; + SD::types::substitution_map sub_map; + sub_map[x] = SD_number_t(x_val); + sub_map[y] = SD_number_t(y_val); + const double f_val = f.substitute_and_evaluate(sub_map); + Assert(f_val == (x_val * x_val - 2.0 * x_val * y_val + y_val * y_val), + ExcMessage( + "Problem with evaluation of single piecewise function")); + } + } + + std::cout << "Substitution with evaluation (double piecewise function)" + << std::endl; + { + const SD_number_t x("x"); + const SD_number_t y("y"); + const SD_number_t x_plus_y = x + y; + const SD_number_t x_minus_y = x - y; + const SD_number_t f_plus = x_plus_y * x_plus_y; + const SD_number_t f_plus_plus = 2.0 * x_plus_y * x_plus_y; + const SD_number_t f_minus = x_minus_y * x_minus_y; + const SD_number_t f((x > SD_number_t(0.0)), + SD_number_t(x > SD_number_t(2.0), + f_plus_plus, + f_plus), + f_minus); + + // Outer condition is met, inner condition is met + { + const double x_val = 3.0; + const double y_val = 2.0; + SD::types::substitution_map sub_map; + sub_map[x] = SD_number_t(x_val); + sub_map[y] = SD_number_t(y_val); + const SD_number_t f_subs = f.substitute(sub_map); + const double f_val = static_cast(f_subs); + Assert(f_val == (2.0 * x_val * x_val + 4.0 * x_val * y_val + + 2.0 * y_val * y_val), + ExcMessage( + "Problem with evaluation of double piecewise function")); + } + + // Outer condition is met, inner condition is not met + { + const double x_val = 1.0; + const double y_val = 2.0; + SD::types::substitution_map sub_map; + sub_map[x] = SD_number_t(x_val); + sub_map[y] = SD_number_t(y_val); + const SD_number_t f_subs = f.substitute(sub_map); + const double f_val = static_cast(f_subs); + Assert(f_val == (x_val * x_val + 2.0 * x_val * y_val + y_val * y_val), + ExcMessage( + "Problem with evaluation of double piecewise function")); + } + + // Condition is not met + { + const double x_val = -1.0; + const double y_val = 2.0; + SD::types::substitution_map sub_map; + sub_map[x] = SD_number_t(x_val); + sub_map[y] = SD_number_t(y_val); + const double f_val = f.substitute_and_evaluate(sub_map); + Assert(f_val == (x_val * x_val - 2.0 * x_val * y_val + y_val * y_val), + ExcMessage( + "Problem with evaluation of double piecewise function")); + } + } + + std::cout << "Substitution with evaluation (if-else_if-else type function)" + << std::endl; + { + const SD_number_t x("x"); + const SD_number_t f_1 = x; + const SD_number_t f_2 = 2 * x; + const SD_number_t f_3 = 3 * x; + const SD_number_t f_4 = 4 * x; + const SD_number_t f({{(x > SD_number_t(4.0)), f_1}, + {(x > SD_number_t(2.0)), f_2}, + {(x > SD_number_t(0.0)), f_3}}, + f_4); + + // Condition "if" is met + { + const double x_val = 5.0; + SD::types::substitution_map sub_map; + sub_map[x] = SD_number_t(x_val); + const SD_number_t f_subs = f.substitute(sub_map); + const double f_val = static_cast(f_subs); + Assert( + f_val == x_val, + ExcMessage( + "Problem with evaluation of 'if' branch of terminated piecewise function")); + } + + // Condition "else if (1)" is met + { + const double x_val = 3.0; + SD::types::substitution_map sub_map; + sub_map[x] = SD_number_t(x_val); + const SD_number_t f_subs = f.substitute(sub_map); + const double f_val = static_cast(f_subs); + Assert( + f_val == (2 * x_val), + ExcMessage( + "Problem with evaluation of 'else if (1)' branch of terminated piecewise function")); + } + + // Condition "else if (2)" is met + { + const double x_val = 1.0; + SD::types::substitution_map sub_map; + sub_map[x] = SD_number_t(x_val); + const SD_number_t f_subs = f.substitute(sub_map); + const double f_val = static_cast(f_subs); + Assert( + f_val == (3 * x_val), + ExcMessage( + "Problem with evaluation of 'else if (2)' branch of terminated piecewise function")); + } + + // Condition "else" is met + { + const double x_val = -1.0; + SD::types::substitution_map sub_map; + sub_map[x] = SD_number_t(x_val); + const SD_number_t f_subs = f.substitute(sub_map); + const double f_val = static_cast(f_subs); + Assert( + f_val == (4 * x_val), + ExcMessage( + "Problem with evaluation of 'else' branch of terminated piecewise function")); + } + } + + std::cout << "Substitution with evaluation (if-else_if type function)" + << std::endl; + { + const SD_number_t x("x"); + const SD_number_t f_1 = x; + const SD_number_t f_2 = 2 * x; + const SD_number_t f_3 = 3 * x; + const SD_number_t f({{(x > SD_number_t(4.0)), f_1}, + {(x > SD_number_t(2.0)), f_2}, + {(x > SD_number_t(0.0)), f_3}}); + + // Condition "if" is met + { + const double x_val = 5.0; + SD::types::substitution_map sub_map; + sub_map[x] = SD_number_t(x_val); + const SD_number_t f_subs = f.substitute(sub_map); + const double f_val = static_cast(f_subs); + Assert( + f_val == x_val, + ExcMessage( + "Problem with evaluation of 'if' branch of non-terminated piecewise function")); + } + + // Condition "else if (1)" is met + { + const double x_val = 3.0; + SD::types::substitution_map sub_map; + sub_map[x] = SD_number_t(x_val); + const SD_number_t f_subs = f.substitute(sub_map); + const double f_val = static_cast(f_subs); + Assert( + f_val == (2 * x_val), + ExcMessage( + "Problem with evaluation of 'else if (1)' branch of non-terminated piecewise function")); + } + + // Condition "else if (2)" is met + { + const double x_val = 1.0; + SD::types::substitution_map sub_map; + sub_map[x] = SD_number_t(x_val); + const SD_number_t f_subs = f.substitute(sub_map); + const double f_val = static_cast(f_subs); + Assert( + f_val == (3 * x_val), + ExcMessage( + "Problem with evaluation of 'else if (2)' branch of non-terminated piecewise function")); + } + + // No condition is met + { + const double x_val = -1.0; + SD::types::substitution_map sub_map; + sub_map[x] = SD_number_t(x_val); + const SD_number_t f_subs = f.substitute(sub_map); + const double f_val = static_cast(f_subs); + Assert( + numbers::is_nan(f_val) == true, + ExcMessage( + "Problem with evaluation of 'else' branch of non-terminated piecewise function")); + } + } + } + + deallog << "OK" << std::endl; +} diff --git a/tests/symengine/symengine_number_type_01.output b/tests/symengine/symengine_number_type_01.output new file mode 100644 index 0000000000..94923abb18 --- /dev/null +++ b/tests/symengine/symengine_number_type_01.output @@ -0,0 +1,17 @@ + +DEAL::Constructors +DEAL::Utility functions +DEAL:Serialization::0 0 1 x + +DEAL:Serialization::0 0 1 y + +DEAL:Serialization::0 0 11 2.5 + x + y + +DEAL::Comparison operators +DEAL::Logical operators +DEAL::Assignment operators +DEAL::Math operators +DEAL::Differentiation +DEAL::Substitution +DEAL::Evaluation +DEAL::OK -- 2.39.5